Lambert W function / ProductLog without overflow

productlog.py

# https://en.wikipedia.org/wiki/Lambert_W_function
# https://code.activestate.com/recipes/577729-lambert-w-function/
import sys, math;

def w0_v1(x, epsilon=0.00000001): # Lambert W function using Newton's method
	w = float(x);
	while True:
		ew = math.exp(w);
		w_new = w - (w * ew - x) / (w * ew + ew);

		# Done condition
		if (abs(w - w_new) <= epsilon):
			return w_new;

		w = w_new;

def w0_v2(x, epsilon=0.00000001):
	w = math.log(float(x));
	max_exp = 30;
	max_exp_value = math.exp(max_exp);

	while True:
		# Expansion of "w_new = w - (w * math.exp(w) - x) / (w * math.exp(w) + math.exp(w))" that won't overflow
		x_div_w1 = x / (w + 1);
		w_exp = w;
		while (w_exp > max_exp):
			x_div_w1 /= max_exp_value;
			w_exp -= max_exp;
		x_div_w1 /= math.exp(w_exp);

		w_new = w - w / (w + 1) + x_div_w1;

		# Done condition
		if (abs(w - w_new) <= epsilon):
			return w_new;

		w = w_new;


# Examples
x = 20;
sys.stdout.write("w0_v1({0:d}) = {1:.16f}\n".format(x, w0_v1(x)));
sys.stdout.write("w0_v2({0:d}) = {1:.16f}\n\n".format(x, w0_v2(x)));

x = 60;
k = 0.000000005;

try:
	v = w0_v1(x / k);
	s = "{0:.16f}".format(v);
except OverflowError:
	s = "overflow";

sys.stdout.write("w0_v1({0:d} / {1:f}) = {2:s}\n".format(x, k, s));

try:
	v = w0_v2(x / k, 0.01);
	s = "{0:.16f}".format(v);
except OverflowError:
	s = "overflow";

sys.stdout.write("w0_v2({0:d} / {1:f}) = {2:s}\n".format(x, k, s));