New Sum of Powers Log Function

[Namir]

Albert,

Your guess3() function provides good approximation to the SopLog function.

Here is the SopLog implementation in Matlab:

Code:


function x = soplog(N, S)
%SOPLOG implements the SopLog function.

  toler=1e-8; % default tolerance
  x = 1; % initial guess for x
  maxiter = 1000;
  for iter =1:maxiter
    h = 0.001 * (1 + abs(x));
    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;
    if abs(diff) < toler, break; end
  end
end

function y =sumFx(N,S,x)
  sumx = 1;
  for i=2:N
    sumx = sumx + i^x;
  end
  y = S - sumx;
end

Here is my copy of your Python guiess3() + some test calls:

Code:

from mpmath import *
import numpy as np

def guess3(n,s):
    t = -np.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


print(guess3(100,1500))
print(guess3(1000,5000))
print(guess3(1000,500))

Here is my Matlab version of your guess3() function which I call soplogw3:

Code:


function x = soplogw3(n,s)
%SOPLOGW Summary of this function goes here
  t = -log(n+.5);
  x = t/s;
  w = lambertw(-1,x);
  p = w/t;
  h = t/(s+0.5^p/p) - x;
  x = real(p-1 + p/(w*x+x) * h);
end

I do beleive we must use the solution path of -1 with the Lambert function. The W0(x) does not give the same answers.

Here are three test values:

Code:

soplog(100,1500)  -> 0.701661431056288
print(guess3(100,1500)) -> 0.70166598312237
soplogw3(200,2500)       -> 0.701665983122370

soplog(1000,5000)           -> 0.267187615565215
print(guess3(1000,5000)) -> 0.267188163707435
soplogw3(1000,5000)       -> 0.267188163707435

soplog(1000,500)           -> -0.118482712209941
soplogw3(1000,500)       -> -0.118485900562627
print(guess3(1000,500)) -> -0.118485900562627

The soplog() is my original Matlab implementation. The guess3() is you latest implementation, and the soplogw3() is the Matlab equivalent of your guess3(). The results show that using the Lambert W function yields results that match to 5 or 6 decimals. Pretty good!