(28/48/50) Lambert W Function
04-01-2023, 05:36 PM
Post: #21
 John Keith Senior Member Posts: 1,039 Joined: Dec 2013
RE: (28/48/50) Lambert W Function
That sounds very promising. It will take me a while to wrap my head around it though.

As a bonus, we have an accurate Lambert W program for the HP-71B.

Meanwhile...
04-01-2023, 05:59 PM
Post: #22
 John Keith Senior Member Posts: 1,039 Joined: Dec 2013
RE: (28/48/50) Lambert W Function
Here is the current version of the Lambert W program for all RPL calculators. This version covers all representable numbers, real or complex, in any branch. In almost all cases the result will be accurate within 2 ULP's, except for a circle of radius .001 around -1/e.

This program replaces the first and last programs in post #1. As before, level 2 should have an integer representing the branch, and level 1 should have the real or complex number z.

Code:
 \<<   IF DUP ABS                                        @ z not 0?   THEN SWAP     IF DUP                                          @ k not 0?     THEN       IF DUP -1 SAME                                @ If branch -1       THEN         IF OVER RE DUP 0 < SWAP -2 > AND         3 PICK IM DUP 0 \>= SWAP .2 < AND AND       @ Near -1/e?         THEN DROP           IF DUP RE -.25 >                          @ Estimate for branch -1           THEN DUP NEG LN DUP NEG LN -              @ near -1/e           ELSE DUP -1. SWAP 1 EXP * 1 +             2 * \v/ -           END 6 0                                   @ Do 6 iterations         ELSE 1                                      @ Use complex code         END       ELSE 1                                        @ Use complex code       END     ELSE DROP DUP ABS .00000001       IF \>=                                        @ If |z| >= 1E-8       THEN DUP 1 EXP * 1 + \v/ DUP 1.14956131 *     @ Estimate for branch 0         1 + SWAP 1 + LN .4549574 * 1 + /         LN 2.036 * 1 - 3                            @ Do 3 iterations       ELSE DUP 1                                    @ DUP if |z| < 1E-8       END 0     END     IF NOT     THEN 1 SWAP                                     @ Iteration for branch 0       START DUP2 / LN 1 + SWAP DUP 1 + / *          @ and branch -1 near -1/e       NEXT     ELSE OVER SWAP (0, 6.28318530718) *             @ Estimate for       SWAP LN + DUP LN DUP2 DUP2 - ROT ROT SWAP / + @ complex branches       ROT ROT DUP 2 - * SWAP SQ 2 * / +       1 3       START DUP2 DUP EXP DUP ROT * ROT -            @ Halley iteration         3 PICK 1 + ROT * 3 PICK 2 + 3 PICK *         4 PICK 2 * 2 + / - / -       NEXT     END   ELSE                                              @ z = 0     IF OVER                                         @ If branch not 0, return     THEN INV NEG                                    @ negative infinity     END   END SWAP DROP \>>

Also, a shorter version for the HP 49 and 50 using UNROT, PICK3 and IFTE. The program can be used in approximate or exact mode, but the result will be approximate. Execution time on my 50g ranges from 110 to 490 ms depending on input.

Code:
 \<< \->NUM   IF DUP ABS   THEN SWAP IP     IF DUP     THEN DUP -1. SAME       { OVER RE DUP 0. < SWAP -2. > AND         PICK3 IM DUP 0. \>= SWAP .2 < AND AND         { DROP DUP RE -.25 >           { DUP NEG LN DUP NEG LN - }           { DUP -1. SWAP 1. EXP * 1. + 2. * \v/ -           } IFTE 6. 0.         } 1. IFTE       } 1. IFTE     ELSE DROP DUP ABS .00000001 \>=       { DUP 1. EXP * 1. + \v/ DUP 1.14956131 * 1. +         SWAP 1. + LN .4549574 * 1. + / LN 2.036 * 1. - 3.       }       { DUP 1. } IFTE 0.     END NOT     { 1. SWAP     START DUP2 / LN 1. + SWAP DUP 1. + / *     NEXT     }     { OVER SWAP (0.,6.28318530718) * SWAP LN +       DUP LN DUP2 DUP2 - UNROT SWAP / + UNROT       DUP 2. - * SWAP SQ 2. * / +       1. 3.       START DUP2 DUP EXP DUP ROT * ROT -         PICK3 1. + ROT * PICK3 2. + PICK3 *         4. PICK 2. * 2. + / - / -       NEXT      } IFTE   ELSE OVER { INV NEG } IFT   END SWAP DROP \>>
04-02-2023, 11:12 PM
Post: #23
 Albert Chan Senior Member Posts: 2,703 Joined: Jul 2018
RE: (28/48/50) Lambert W Function
(03-26-2023 06:43 PM)Albert Chan Wrote:  Because solving for y = e^W(a) is easier, W code is not recommended.

W code is perfect for other branches! (expW code for W branch 0, -1)

Numbering the branches of the Lambert W function

Except for 2*k*pi*I term, it is basically same formula for W code!

lua> f = fn'x: x + I.log(x) - I.log(a) - 2*k*pi*I'
lua> df = fn'x: 1 + 1/x'

Let X = |x|, A = |a|

log(x) - log(a) = log(X/A) + I*(arg(x) - arg(a))
log(X/A) = 1 + log(X*(r+err)/A) ≈ 1 + log1p((X*err - (A+r) + r*(1+X)) / A)

This means we can evaluate f accurately around branch point!

f=0 guaranteed we get the correct branch, which is nice.
There is no need to use complicated W guesses.

lua> a, k = 3+4*I, 5
lua> x = 2*k*pi*I -- guess for Wk(a)
lua> repeat h=f(x)/df(x); x=x-h; print(x, I.abs(h)) until x == x+1e-6*h

(-1.815554248884793+30.714634540472343*I)       1.9462907525577031
(-1.817005891456528+30.71333413897562*I)        0.0019489253984594729
(-1.8170058918466274+30.713334137004896*I)      2.0089622485451688e-009
(-1.8170058918466274+30.713334137004896*I)      0

>>> from mpmath import *
>>> lambertw(3+4j, k=5)
mpc(real='-1.8170058918466274', imag='30.713334137004896')
04-03-2023, 07:24 PM
Post: #24
 John Keith Senior Member Posts: 1,039 Joined: Dec 2013
RE: (28/48/50) Lambert W Function
(04-02-2023 11:12 PM)Albert Chan Wrote:  f=0 guaranteed we get the correct branch, which is nice.
There is no need to use complicated W guesses.

lua> a, k = 3+4*I, 5
lua> x = 2*k*pi*I -- guess for Wk(a)
lua> repeat h=f(x)/df(x); x=x-h; print(x, I.abs(h)) until x == x+1e-6*h

(-1.815554248884793+30.714634540472343*I)       1.9462907525577031
(-1.817005891456528+30.71333413897562*I)        0.0019489253984594729
(-1.8170058918466274+30.713334137004896*I)      2.0089622485451688e-009
(-1.8170058918466274+30.713334137004896*I)      0

>>> from mpmath import *
>>> lambertw(3+4j, k=5)
mpc(real='-1.8170058918466274', imag='30.713334137004896')

Thanks for the John D. Cook link, I hadn't seen that one before.

My program gets exact or near-exact results for all branches, but it is slower for branches other than 0 and -1. I will try your code and see how it compares.
04-03-2023, 08:47 PM (This post was last modified: 04-03-2023 09:15 PM by Albert Chan.)
Post: #25
 Albert Chan Senior Member Posts: 2,703 Joined: Jul 2018
RE: (28/48/50) Lambert W Function
(04-03-2023 07:24 PM)John Keith Wrote:  My program gets exact or near-exact results for all branches, but it is slower for branches other than 0 and -1.

I was getting stupid. Of course you get near-exact result for other branches.

For branches besides (0, -1), a ≈ -1/e is not that important, since W(-1/e) ≠ -1
With huge complex parts in the way, there is no catastrophic cancellation issues.      (*)

The singularity is at only at 0, not -1/e

W-1 is the "hard" branch, with singularities at both places, 0 and -1/e.

Plot[{Re[LambertW[-1,x]], Im[LambertW[-1,x]]}, {x,-1,1}]

(*) Except for W1: lambertw(z, k) == conj(lambertw(conj(z), -k))

>>> lambertw(-.36787944+1e-99j, k=-1)
mpc(real='-1.0000798057615958', imag='-3.4063941215654215e-95')

>>> conj(lambertw(-.36787944-1e-99j, k=1))
mpc(real='-1.0000798057615958', imag='-3.4063941215654215e-95')
04-03-2023, 10:47 PM (This post was last modified: 04-04-2023 12:19 AM by Albert Chan.)
Post: #26
 Albert Chan Senior Member Posts: 2,703 Joined: Jul 2018
RE: (28/48/50) Lambert W Function
(04-02-2023 11:12 PM)Albert Chan Wrote:  lua> a, k = 3+4*I, 5
lua> x = 2*k*pi*I -- guess for Wk(a)

With Newton's rapid convergence, it might be overkill for better W guess ...

x = r cis(θ)
e^x = e^(r*cos(θ)) cis(r*sin(θ))
x*e^x = r*e^(r*cos(θ)) cis(θ+r*sin(θ)) = a

r*e^(r*cos(θ)) = |a|
θ+r*sin(θ) = arg(a) + 2*k*pi = T           // for W branch k

cos(θ) = ln(|a|/r)/r
sin(θ) = (T-θ)/r

We assume k positive --> both r and x imag part positive too.

lua> a, k = 3+4*I, 5
lua> A, T = I.abs(a), 2*k*pi + I.arg(a)
lua> A, T
5      32.34322175389954

lua> c = 0 -- = cos(θ)
lua> r = (T-acos(c)) / sqrt(1-c*c) -- = T - pi/2
lua> c = log(A/r)/r
lua> x = r * I.new(c, sqrt(1-c*c))
lua> x -- Wk(a) guess
(-1.81718109820366+30.71872424960138*I)

Rinse and repeat ...

lua> r = (T-acos(c)) / sqrt(1-c*c)
lua> c = log(A/r)/r
lua> x2 = r * I.new(c, sqrt(1-c*c))
lua> x2
(-1.8170057678766207+30.7133303234809*I)

lua> r = (T-acos(c)) / sqrt(1-c*c)
lua> c = log(A/r)/r
lua> x3 = r * I.new(c, sqrt(1-c*c))
lua> x3
(-1.8170058919343504+30.713334139703402*I)

lua> aitken(x, x2, x3)
(-1.8170058918466334+30.71333413700532*I)

We don't even need machine to support complex numbers.
We could just iterate for (r, c), improve with Aitken's delta-squared process
Assuming K is positve, after (r, c) converged, Wk(a) = complex(r*c, r*sqrt(1-c*c))
04-04-2023, 01:03 AM (This post was last modified: 04-04-2023 02:26 PM by Albert Chan.)
Post: #27
 Albert Chan Senior Member Posts: 2,703 Joined: Jul 2018
RE: (28/48/50) Lambert W Function
(04-03-2023 10:47 PM)Albert Chan Wrote:  cos(θ) = ln(|a|/r)/r
sin(θ) = (T-θ)/r

We assume k positive --> both r and x imag part positive too.

We can setup solver to get r, then Wk(a)
Let c = cos(θ), s = sin(θ)

s = (T-θ)/r
θ = T-r*s = T - r*√(1-c²)

s ≤ 1      → min(r) = T-θ > T-pi

lua> g = fn'r,c: c=log(A/r)/r; acos(c) - (T-r*sqrt(1-c*c))'

lua> S = require'solver'
lua> r = S.secant(g, T-pi, T+pi, nil, true)
29.20162910030975
35.48481440748934
30.766823278237727
30.76703434746101
30.767034374835603
30.767034374835603
lua> c = ln(A/r)/r
lua> r * I.new(c, sqrt(1-c*c)) -- = W5(3+4*I)
(-1.8170058918466274+30.713334137004896*I)

I setup function comparing angles. For |k| ≥ 2, r only have 1 real root.
04-04-2023, 07:09 PM
Post: #28
 John Keith Senior Member Posts: 1,039 Joined: Dec 2013
RE: (28/48/50) Lambert W Function
Your HP-71 program from post #19 is very good. It outperforms my program in both branches for -1/e < x <= -.2

I frankly don't feel that FNL from post #20 is worth the extra memory requirement since the post 19 program is always within 3 ULP's.

Here is my RPL translation. It is optimized for size rather than readability (lots of "stackrobatics") so I have added comments corresponding to the BASIC statements. As with the other programs, the level 2 argument is k, and the level 1 argument is x. This program is not intended for x > -.1, nor for k other than 0 or -1.

I'm not sure whether I will integrate it into the main program, which already takes over 600 bytes. For now I'll just leave it as-is for anyone who needs higher accuracy near the branch point.

Code:
 \<< -1. EXP DUP ROT +     4.42321595524E-13 + SWAP /    @ H=(A+R+R2)/R     DUP 2. * \v/ ROT DUP NOT + *  @ X=SQRT(2*H)*K      DUP 3. / 1. + \v/ *          @ X=X*SQRT(1+X/3)     DO DUP LNP1                   @ Z=LOGP1(X)       DUP2 - 4. PICK + SWAP /       OVER SWAP -                 @ Z=X-(X-Z+H)/Z       SWAP OVER -                 @ X=X-Z       SWAP .000001 * OVER +       @ X=X+Z*.000001     UNTIL OVER SAME     END SWAP DROP LNP1 1. -       @ LOGP1(X)-1 \>>
01-29-2024, 11:04 AM
Post: #29
 Gil Senior Member Posts: 637 Joined: Oct 2019
RE: (28/48/50) Lambert W Function
Hi, Albert.

You write 70 Y=(K=1)+K*A/4
80 REPEAT @ H=Y-(Y+A)/(LOG(Y)+1) @ Y=Y-H @ UNTIL Y=Y+H*.0001

What is the meaning of "(K=1)" in line 70 "
a) Y=k+k*abs(a)/4 in any case
b) or Y=k*abs(a)/4; Y=Y+(1, if k=1; 0, else)
c) or Y=?
01-29-2024, 02:47 PM (This post was last modified: 01-29-2024 02:58 PM by Gil.)
Post: #30
 Gil Senior Member Posts: 637 Joined: Oct 2019
RE: (28/48/50) Lambert W Function
I tried to translate your code in post 17, Albert.

{ "« x -.1 <=
IF
THEN .367879441171 'R' STO 4.42321595524E-13 'R2' STO 'sqrt(2.*R*(x+R+R2))*(k+k+1.)' —>NUM 'Y' STO 'R+Y*sqrt (1.+Y/(3.*R))' —) NUM 'Y' STO
DO Y 'Y+x' '(Y-R-R2)/R' —>NUM DUP TYPE 1. ==
IF
THEN 1. + LN
ELSE LNP1
END / - —>NUM 'H' STO 'Y-H' —>NUM 'Y' STO Y 'Y+H*.0001' —>NUM ==
UNTIL
END
ELSE 'k+k+1.' '(k+k+1.)*x/4.' + —>NUM 'Y' STO
DO 'Y-(Y+x)/(LN(Y)+1.)' —>NUM 'H' STO 'Y-H' —>NUM 'Y' STO
UNTIL Y 'Y+H*.0001' —>NUM ==
END
END x Y /

»-x =-.36787944118
W(k=->; (-.999999999985,-6.87304618656E-6)

Result ± ok for imaginary part.

Your translated code for more accuracy, post 19, same input values
« x -.1 
IF
THEN .367879441171 'R' STO 4.42321595524E-13 'R2' STO '(x+R+R2)/R' —>NUM 'H' STO 'sqrt(2.*H)*(2.*k+1.)' —>NUM 'X' STO 'X*sqrt(1.+X/3.)' —>NUM 'X' STO 'R+R*X' —>NUM 'Y' STO
DO X DUP TYPE 1. ==
IF
THEN 1. + LN
ELSE LNP1
END 'Z' STO 'X-(X-Z+H)/Z', —>NUM 'Z' STO X Z - 'X' STO X 'X+Z*.000001' —>NUM ==
UNTIL
END X DUP TYPE 1. ==
IF
THEN 1. + LN
ELSE LNP1
END 1. -
ELSE 'k+k+1.' '(k+k+1.)*x/4.' + —>NUM 'Y' STO
DO 'Y-(Y+x)/(LN(Y)+1.)' —>NUM 'H' STO 'Y-H' —>NUM 'Y' STO
UNTIL Y 'Y+H*.0001' —>NUM ==
END x Y /
END
»

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

Why this difference of "behaviour"?

And why, in code post 29, did you define Y=R+R*X, that seems never to be used again?
01-29-2024, 03:51 PM (This post was last modified: 01-29-2024 05:20 PM by Albert Chan.)
Post: #31
 Albert Chan Senior Member Posts: 2,703 Joined: Jul 2018
RE: (28/48/50) Lambert W Function
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.
01-29-2024, 06:46 PM
Post: #32
 Gil Senior Member Posts: 637 Joined: Oct 2019
RE: (28/48/50) Lambert W Function
Thanks.

And why, in code post 29, did you define Y=R+R*X, that seems never to be used again?
01-29-2024, 07:30 PM
Post: #33
 Albert Chan Senior Member Posts: 2,703 Joined: Jul 2018
RE: (28/48/50) Lambert W Function
(01-29-2024 06:46 PM)Gil Wrote:  And why, in code post 19, did you define Y=R+R*X, that seems never to be used again?

(R+R*X) is just another method to calculate   Y = e^W(A)
It is only for compare against Newton method Y = (Y+A)/LOGP1((Y-R-R2)/R)

I was trying to decide which way does the job better.

(12-30-2022 11:27 PM)Albert Chan Wrote:  When I try to understand lyuka's eW formula (previous post), I found a better one.

e^W(a) = y = (y+a) / (ln(y)+1)

Let r = 1/e
If a = -r + h (tiny h), we have y = r + z (tiny z)

e^W(a) = (r+z) = (h+z) / ln(1+e*z)

Let H = e*h (known), Z = e*z (unknown), we have:

(1+Z) = (H+Z) / ln(1+Z)

H = (1+Z) * ln(1+Z) - Z = Z^2/2 - Z^3/6 + Z^4/12 - Z^5/20 + ...
...

I am skipping details of how to estimate Z from H.
Once we have Z, e^W(a ≈ -1/e) problem is solved.

e^W(a) = y = r + z = r + Z/e = r + r*Z      // = R+R*X of HP71B code
01-29-2024, 09:50 PM
Post: #34
 Gil Senior Member Posts: 637 Joined: Oct 2019
RE: (28/48/50) Lambert W Function
What hawk eyes, Albert! Not only your brain works marvels, Albert, but your eyes too.
Sorry for having mixed up the post number!

More seriously: for k=-1, and a=-.36787944118 (>-1/e, about -.367879441171), how do I know that formulae post 19, more accurate, will fail (at least with the code I put in my calculator) and that formulae post 18, a wee bit less accurate, will work?

In other words: should I say, for "a" near -1/e, but just smaller — as in the example above —, I should always play safety and only use formulae in post 18?

Second question
k mappin — new k = ± k.
Does it mean that the derived formulae in posts 18 and 19 only work for k, initially = 0 or 1?

If yes, are there other pitfalls to look at for initial k≠{0, +-1} when "a" is near 0 or -1/e?

Regards
01-30-2024, 12:18 AM
Post: #35
 Albert Chan Senior Member Posts: 2,703 Joined: Jul 2018
RE: (28/48/50) Lambert W Function
(01-29-2024 09:50 PM)Gil Wrote:  for k=-1, and a=-.36787944118 (>-1/e, about -.367879441171), how do I know that formulae post 19, more accurate, will fail (at least with the code I put in my calculator) and that formulae post 18, a wee bit less accurate, will work?

They are equally accurate, if coded correctly.
They are really the same e^W formula, except we let Y = R+R*X
It is just scale and offset, convergence rate are exactly the same.

My preference is Y = (Y+A)/LOGP1((Y-R-R2)/R), because convergence is easily tested.
It is also easy to explain. Denominator is simply more accurate (LN(Y)+1)

The other version, X = (X-L+H)/L, where L = LOGP1(X), is harder to know when to stop.
Because of inaccurate X-L, X itself never really converged, only (1+X) = Y/R does.

OTTH, X formula is simpler than Y's (also, (R, R2) not used at all!)
With limited domain (A ≈ -1/e), fixed loops is simple to implement.

Quote:If yes, are there other pitfalls to look at for initial k≠{0, +-1} when "a" is near 0 or -1/e?

Yes! Only used this for k=0, or ±1 with small_imag part (k and im(a) of opposite sign)
This is because y = e^x lost branch information. Good guess is crucial to get right branch!

e^(2*k*pi*I) = cis(2*k*pi) = 1      // k is gone

What happen if k = ±1 and, k and im(a) have same sign?

x + ln(x) = lnk(a)

im(x) + arg(x) = (arg(a) + 2*k*pi) = [±2*pi .. ±3*pi]
im(x) = [±2*pi .. ±3*pi] - [0 .. ± pi] = [± pi .. ±3*pi]

Away from real line, no singularity around -1/e, easy to solve.
01-30-2024, 12:33 AM
Post: #36
 Gil Senior Member Posts: 637 Joined: Oct 2019
RE: (28/48/50) Lambert W Function
Thanks.

Strange anyway, a = -.36787944118, that in my calculator the until equality is never reached with method exposed in post 11.—> endless loops.
01-30-2024, 01:09 AM
Post: #37
 Albert Chan Senior Member Posts: 2,703 Joined: Jul 2018
RE: (28/48/50) Lambert W Function
(01-30-2024 12:18 AM)Albert Chan Wrote:  With limited domain (A ≈ -1/e), fixed loops is simple to implement.

For my Python code (~ 15 digits precision), it never get over 5 loops.
Perhaps you can use fixed Newton loops instead of testing convergence.

Also, you may use this to implement complex log1p(x)

y=1+x --> log1p(x) ≈ ln(y) - (y-1-x)/y

01-30-2024, 12:04 PM (This post was last modified: 01-30-2024 12:57 PM by Gil.)
Post: #38
 Gil Senior Member Posts: 637 Joined: Oct 2019
RE: (28/48/50) Lambert W Function
Hi, Albert.

With your function definition of post 20,

1) could I modify as follows?

200 DEF FNL(X) !
« —> X
«210 IF ABS(X)>=.4
THEN X=X/(SQRT(1+X)+1) @ L=FNL(X)*2-X*X
ELSE
220 X2=X/(X+2) @ X4=X2*X2
230 X4=X4*(5005-X4*(5082-X4*969))/(15015-X4*(24255-X4*(11025-X4*1225)))
240 L=X2*(X4+X4-X)
250 END
»
»

for my HP50G calculator program.

2) Could I, consequently, simplify "your" line 42 as follows?

42 DO FNL(X) @Z=X-(H-L)/(X+L) @ X=X-Z @ UNTIL X=X+Z*.0001 END
44 PRINT "e^W, W =";R+R*X;LOGP1(X)-1

3) And I should start the program with line 40 defined as before
>40 H=(A+R+R2)/R @ X=SQRT(2*H)*K @ X=X*SQRT(1+X/3) @ Y=R+R*X

Correct?

In fact, it seems that L cannot be equal to @ L=FNL(X)*2-X*X

Or the defined function should be rather
200 DEF FNL(X) !
« —> X
«210 IF ABS(X)>=.4
THEN X=X/(SQRT(1+X)+1) @ L=X*2-X*X
ELSE
220 X2=X/(X+2) @ X4=X2*X2
230 X4=X4*(5005-X4*(5082-X4*969))/(15015-X4*(24255-X4*(11025-X4*1225)))
240 L=X2*(X4+X4-X)
250 END

Thanks again, Albert, for your possible and always precious help.

Gil
01-30-2024, 01:47 PM
Post: #39
 Albert Chan Senior Member Posts: 2,703 Joined: Jul 2018
RE: (28/48/50) Lambert W Function
Hi, Gil

FNL(x) = accurate ln(1+x)-x, was designed for real number.
It may need work to extend for complex numbers.

Also, it is defined recursively, until x is small, not a IF THEN ELSE structure.
see thread: Accurate x - log(1+x) (for consisteny, I switched order, and do log1p_sub(x))

Anyway, accurate FNL(x) is not needed here.
Here is a direct replacement snippet for my Python W code

I am using formula to get log1p, instead of calling built-in log1p
Quote:y=1+x --> log1p(x) ≈ ln(y) - (y-1-x)/y

Because of inaccurate x-log1p(x), x does not converge, but (1+x) does.
(again, you may forget about checking convergence, and just loop a few times)

Code:
    if abs(h)<.25 and (k==0 or small_imag):         h = h/r         x = sqrt(2*h) * (1-2*k*k)   # close to branch point         x *= sqrt(1 + x/3)         while 1:             if verbal: print a/(r+r*x)             y = 1+x             z = ln(y) - (y-1-x)/y   # = log1p(x)             x = (x-z+h)/z             z = 1+x             if abs((z-y)/z) < 1e-12: return a/(r+r*x)

With identical termination condition, results compared well to W original version.
Code:
    if abs(h)<.25 and (k==0 or small_imag):         y = sqrt(2*r*h) * (1-2*k*k) # close to branch point         y = r + y*sqrt(1+y/(3*r))   # with small imag part         while 1:             if verbal: print a/y             h = y - (y+a) / log1p((y-r-err)/r)             y = y - h             if abs(h/y) < 1e-12: return a/y
01-30-2024, 02:52 PM
Post: #40
 Gil Senior Member Posts: 637 Joined: Oct 2019
RE: (28/48/50) Lambert W Function
I tried and changed the code fo my function FNL relative to ln(1+x)-x

FNL=
«
WHILE x .4 >=
REPEAT 'x/(sqrt(1+x)+1)' ->NUM 'x' STO x FNL 2 * x SQ - 'fnl' STO
END 'x/(x+2)' ->NUM 'X2' STO X2 DUP * 'X4' STO 'X4*(5005-X4*(5082-X4*969))/(15015-X4*(24255-X4*(11025-X4*1225)))' ->NUM 'X4' STO 'X2*(X4+X4-x)' ->NUM 'fnl' STO
»

But, for initial x=4, I get a final fnl = -4.77474540488E-2 that is incorrect, of course.
 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)