Create a gist now

Instantly share code, notes, and snippets.

@nutbread /productlog.py Secret
Last active Aug 29, 2015

What would you like to do?
Lambert W function / ProductLog without overflow
# 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));
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment