New Sum of Powers Log Function

[Albert Chan]

Even faster version, eliminated second LambertW, by taylor series approximation.

W(x+h) ≈ W(x) + W'(x) * h

w*e^w = x
(w+1)*e^w dw = dx
dw/dx = e^-w / (w+1) = w/(w*x+x)

Code:

def guess3(n,s):
    t = -log(n+.5)
    x = t/s
    w = lambertw(x,-1)
    p = w/t
    h = t/(s+0.5**p/p) - x
    return p-1 + p/(w*x+x) * h

To speed up search, solvex use log scale.

With secants method, and tolerance of 1e-5, we expected 5*1.6 = 8 "correct" decimals.
("correct" from solving root of formula point of view, not actual x)

Code:

RHS = lambda p,n: ((n+1)**(p+1)-1)/(p+1) - ((n+1)**p-1)/2 + p*((n+1)**(p-1)-1)/12
solvex = lambda n,s: findroot(lambda p: log(RHS(p,n)/s), log(s,n)-1, tol=1e-5)