Instantly share code, notes, and snippets.

# nutbread/productlog.pySecret Last active Aug 29, 2015

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));