New Sum of Powers Log Function

[Albert Chan]

(03-31-2021 01:27 PM)Namir Wrote:  Here is the SopLog implementation in Matlab ...

Code:

    f0 = sumFx(N, S, x);
    fp = sumFx(N, S, x + h);
    fm = sumFx(N, S, x - h);
    diff = 2 * h * f0 / (fp - fm);
    x = x - diff;
    ...

Without function for the derivative, Newton's method is very inefficient.
We should consider secant's method (or, improved secant)

Also, curve is very close to LambertW, we should use log scale.

Example: soplog(1000,5000), start with 2 guesses.

>>> n, s = 1000, 5000
>>> x = log(s,n)-1
>>> x1, y1 = x, log(fsum(k**x for k in range(1,n+1))/s)
>>> x1, y1
(mpf('0.23299000144533966'), mpf('-0.20890713759105706')

>>> x = guess3(n,s)
>>> x2, y2 = x, log(fsum(k**x for k in range(1,n+1))/s)
>>> x2, y2
(mpf('0.26718816370743487'), mpf('3.3544082779438775e-6'))

>>> x = x2 - y2 * (x2-x1)/(y2-y1)
>>> x3, y3 = x, log(fsum(k**x for k in range(1,n+1))/s)
>>> x3, y3
(mpf('0.26718761459859175'), mpf('-5.9153427886634265e-9'))

>>> x = x3 - y3 * (x3-x2)/(y3-y2)
>>> x4, y4 = x, log(fsum(k**x for k in range(1,n+1))/s)
>>> x4, y4
(mpf('0.26718761556521503'), mpf('-2.2204460492503136e-16'))

---

Turns out mpmath implemented Euler–Maclaurin Asymptotic expansion of sums
Summation will continue, as long as additional term improve the sum, by at least 1 decimal.

>>> x = solvex(n, s) # my version, without higher derivatives
>>> x, fsum(k**x for k in range(1,n+1))
(mpf('0.26718762849802313'), mpf('5000.0003957177305'))

>>> x = findroot(lambda x: log(sumem(lambda k: k**x, [1,n])/s), log(s,n)-1)
>>> x, fsum(k**x for k in range(1,n+1))
(mpf('0.26718761309749891'), mpf('4999.9999244928886'))