New Sum of Powers Log Function

[Albert Chan]

Revised guess3(n,s), with less code.

Code:

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

Note that guess3 had a weekness, when t/s < -1/e, LambertW returned complex values.
Even if w is real, if n ≫ s > 1, we have p=x+1 ≈ 0. Two terms taylor series estimate is still bad.

Code:

from mpmath import *    
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)
exactx = lambda n,s: findroot(lambda x: log(fsum(k**x for k in range(1,n+1))/s), log(s,n)-1)

>>> n, s = 100, 13       # t/s ≈ -0.35463, slightly above -1/e ≈ -0.36788
>>> print exactx(n,s), solvex(n,s), guess3(n,s), log(s,n)-1
-0.622429308602986    -0.622506385172213    -0.544275714709962    -0.443028323846582

---

Thus, if we want fully converged x, we should start with solvex.
We can use y = log(fsum(k**x for k in range(1,n+1))/s), to estimate next x.
Note: x correction is opposite the sign of y

>>> n, s = 1000, 5000 # previous post example
>>> x = solvex(n,s)
>>> x1, y1 = x, log(fsum(k**x for k in range(1,n+1))/s)
>>> x1, y1
(mpf('0.26718762849802312'), mpf('7.9143542807225195e-8'))

>>> x = x1 - abs(x1)*y1/2
>>> x2, y2 = x, log(fsum(k**x for k in range(1,n+1))/s)
>>> x2, y2
(mpf('0.26718761792493534'), mpf('1.4440531777112248e-8'))

>>> x = x2 - y2 * (x2-x1)/(y2-y1)
>>> x
mpf('0.26718761556521503')