(28/48/50) Lambert W Function

[Albert Chan]

(03-27-2023 11:30 PM)Albert Chan Wrote:  Let H=e*h, L = log(1+x) - x, accurately computed

f = (1+x)*log(1+x) - x - H
f' = log(1+x) + (1+x)/(1+x) - 1 = log(1+x)

x - f/f' = (H-L) / (x+L)

Note that numerator is really an addition, with both H, -L positive.
...
lua> r+r*x, log1p(x)-1 -- = e^W(a), W(a)

When we are close to branch point, x is tiny. Even rough x estimate suffice.
When we are far away, log1p(x)-x does not suffer catastrophic cancellation.

We may not need accurate log1p(x)-x after all.

Here, I added code to roughly solve above f=0 for x, then convert to W's
Old code with iterations of y remained (2nd e^W, W), for comparison.

>40 H=(A+R+R2)/R @ X=SQRT(2*H)*K @ X=X*SQRT(1+X/3) @ Y=R+R*X
>42 REPEAT @ Z=LOGP1(X) @ Z=X-(X-Z+H)/Z @ X=X-Z @ UNTIL X=X+Z*.000001
>44 PRINT "e^W, W =";R+R*X;LOGP1(X)-1

>run
a, k? -.367879441171,0
e^W, W = .367880011645 -.999998449291 ! true W = -.999998449288
e^W, W = .367880011646 -.999998449291
>run
a, k? -.36787944117,0
e^W, W = .367880471317 -.999997199776 ! true W = -.999997199775
e^W, W = .367880471317 -.999997199778
>run
a, k? -.3678794411,0
e^W, W = .367886691321 -.999980292244 ! true W = -.999980292245
e^W, W = .36788669132   -.999980292247
>run
a, k? -.367879441,0
e^W, W = .367890672444 -.999969470701 ! true W = -.999969470701
e^W, W = .367890672444 -.999969470702

...

>run
a, k? -.3678,0
e^W, W = .375551106332 -.979360714956 ! true W = -.979360714958
e^W, W = .37555110633   -.979360714962