(28/48/50) Lambert W Function

[Albert Chan]

Hi, Gil

Post 18 HP71B code is short, I'll just copy/pasted in full

Code:

10 INPUT "a, k? ";A,K @ K=K+K+1
20 IF A<=-.1 THEN
30 R=.367879441171 @ R2=4.42321595524E-13
40 Y=SQRT(2*R*(A+R+R2))*K @ Y=R+Y*SQRT(1+Y/(3*R))
50 REPEAT @ H=Y-(Y+A)/LOGP1((Y-R-R2)/R) @ Y=Y-H @ UNTIL Y=Y+H*.0001
60 ELSE
70 Y=(K=1)+K*A/4
80 REPEAT @ H=Y-(Y+A)/(LOG(Y)+1) @ Y=Y-H @ UNTIL Y=Y+H*.0001
90 END IF
100 PRINT "e^W, W =";Y;A/Y

Line 10: K=K+K+1 map W branch (0,-1) into sign ±1
Line 40: this simplified lyuka's formula 2 branches ±√ term

But lyuka's formula is not suitable if A is too far from -1/e
Example, if -√ term keep growing, eventually we would get e^W guess = 0.
But, this does not match exactly at A = 0, which make W_-1 guess very bad.

For that, I just use simple guess (note: K = ±1)
Line 70: Y=(K=1)+K*A/4

K=-1 --> Y = -A/4            // --> X = ln(-A/4),     roughly ln(-A) if A → -0
K=+1 --> Y = 1 + A/4      // --> X = log1p(A/4), roughly ln(A) if A → INF

---

There was another way to get W, based from how lyuka's formula work?
Then, solve the same thing using accurate log(1+x)-x

Turns out, this super accuracy is not needed.
I made a version with FNL(x) = accurate log(1+x)-x, and it make little difference. (post #20)

(03-31-2023 10:07 PM)Albert Chan Wrote:  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

Gil Wrote:No solution is found: X is never equal to 'X+Z*.000001'.

Why this difference of "behaviour"?

Math assumed real numbers to begin with.
It may work with complex numbers, but you still need adjustments.