HP49-50G,VER17.5 Lambert Funktion Wk(x), k=k= 0, ±1, ±2, ±3, ±4..., x real or complex

+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Software Libraries (/forum-10.html)
+--- Forum: General Software Library (/forum-13.html)
+--- Thread: HP49-50G,VER17.5 Lambert Funktion Wk(x), k=k= 0, ±1, ±2, ±3, ±4..., x real or complex (/thread-21161.html)

HP49-50G,VER17.5 Lambert Funktion Wk(x), k=k= 0, ±1, ±2, ±3, ±4..., x real or complex - Gil - 01-14-2024 07:58 AM

My apologies, I deleted my whole previous post discussion relative to Lambert function W(x).

Here a partial summary of it
with the newest version 17.5 of W(x).

This version permits to calculate Wk(x), like WolframAlpha does, k being the branch number of W(x) and is equal to 0, ±1, ±2, ±3, ±4...

x may be a real or a complex number.

Input
a) k and x:
- k, in stack level 2, for the branch
- and argument x, in stack level 1
b) or simply x when by default k is supposed to be 0
(here, note the following, however: if there is already a number like 0, ±1, ±1., ±2, ±2., ±3, ±3., etc. in stack level 2, that number will be automatically taken as the branch level).

Output
a) Wk(x): w<
b) automatically 2 solution outputs
W0(x): w1
& W-1(x): w2, for x real -1/e<=x<=0

Observations
a) For initial guess of W-1 iteration process & -1/e<x<0, see Wikipedia;
b) when possible, we try to find W, with the built-in ROOT solver for 'W×EXP(W)-x';
c) for complex solution W, we try to use likewise the built-in MLSV solver for ['W×EXP(W)-x'];
d) for limit values, for instance by endless search near -1/e, we used instead one of the many methods proposed by Albert Chan in his own Lambert thread;
e) likewise, we used one of the many methods proposed by Albert Chan in his own Lambert thread when the output could be clearly much improved;
f) see also initial guesses according to Albert Chan, http://www.finetune.co.jp/~lyuka/technote/lambertw/clambertw.html, for W0 (and also for other integer values of k≠0);
g) general calculation of the different k branches, see https://www.johndcook.com/blog/2021/11/15/lambert-w-branch-number/
h) up to now, all test results, with k=0, ±1, ±2, ±3, ±4..., and x real or complex number, corresponded largely (with the exception, generally, of the last two digits) to WolframAlpha outputs;
i) however, contrarily to Wolfram Alpha, for inputs like {100 0}, ie k=100, argument=0, the program output is double, whereas Wolfram simplifies the result to -infinity:
for x—> -0, ie just negative, result here is a complex form:
"-infinity + i×200PI";
for x—> 0+, ie just positive, result will also present a complex form: "-infinity+i×199PI;
j) the same applies for instance for {50, infinity}, ie k=50 & x a "real" (9.99999999999E499) that tends to infinity, for which the program, here again, gives a complex output:
"infinity + i×100PI", when Wolfram simplifies the answer and just gives infinity;
k) extreme tiny input x ≠0, in particular ±1E-497, ±1E-498 ±1E-499, are accepted and calculated correctly;
l) entering the algebraic expression '-1/e' as input, you will get W0 & W-1 = exactly -1:
W0('-(1/e)'~-.367879441171): -1
& W-1('-(1/e)'~-.367879441171): -1;
but entering the approximate value of -1/e, ie the real number -.367879441171, will produce the double, different output W0(-.367879441171): -.999997710398
& W-1(-.367879441171): -1.00000141421.

RE: Lambert Funktion Wk(x), k=k= 0, ±1, ±2, ±3, ±4..., x real or complex - Gil - 01-19-2024 05:05 PM

New version, 8.02, with my thanks to Albert Chan.

The program accepts now all kind of argument : real(x) or imaginary (x) —>± infinity, for any integer branch k.

Input {k x}, x complex number or not.

Let's try some large values of x (large for the real or for the imaginary part)

1) x=1E450, k=0

Input in HP50G
1E450

or
{0, 1E450}

Then press
LAMBERT

Output
:W0(1.E450): 1029.2267288

Does it go to +inf?
Try with Wolfram
Write LambertW (0,1E5000000)
Output: 1.15E7

Or write LambertW (0, infinity)
Output: infinity

And with HP50G?
Write sign for infinity
LS+0
Then press LAMBERT
Output: 'infinity'

2) x=1E450, k=15

Input in HP50G
{15 1E450}

Then press
LAMBERT

Output

:W15(1.E450): (1029.22256567,94.1565503694)

Does that result go to +inf (with no imaginary part, here = 94,16)
Try with Wolfram
Write LambertW (15,1E5000000)
Output: (1.15E7+94.24777)
And you see that 94.24777=30pi

But, if you try and write LambertW (0, infinity)
Output: infinity.
Well, let's admit it is correct.
Next example will put in question that kind of answer.

And with HP50G?
Write {15 infinity}
(LS+0 for infinty)
Then press LAMBERT
Output: "infinity + i(30pi)"


3) x=(-1E450), k=0,

Input in HP50G
(0,-1E450)

or
{0, - 1E450}

Then press
LAMBERT

Output
{ "W0(-inf):" inf-ipi" }

Does it go to +inf?
Try with Wolfram
Write LambertW(-9999E5000000)
Output: 1.15E7+i*3.14159238 (almost pi)

Or write LambertW (-i*infinity)
Output: infinity?!

And with HP50G?
Write sign for infinity
LS+0
Then press key ± (letter M)
Then press LAMBERT
Output: "infinity+ i(pi)"


4) x=(0, 1E450), k=0,
ie a large imaginary part

Input in HP50G
(0,1E450)

or
{0, (0,1E450)}

Then press
LAMBERT

Output
::W0((0.,1.E450)): (1029.22672764,1.56927161862)

Does it go to +inf?
Try with Wolfram
Write LambertW(0+i*9999E5000000)
Output: 1.15E7+i*1.570796

Or write LambertW (0+i*infinity)
Output: infinity?!

And with HP50G?
Write sign for infinity
LS+0
Then press LAMBERT
Output: "infinity+ i(pi/2)"

5) case x=(inf-i*inf), k= 30
Write 30
ENTER
key oo, ie LS+0[/code]
—>NUM
DUP
key ± (letter M)
R—>C (to get together inf, real part, and -inf, imaginary part, into a complex number)
2 —> list
You should get the following list
{30, (9.999999999E499, 9.999999999E499)}
Then press LAMBERT on the HP50G.

Output on HP50G
"oo +i(60pi).

Approximation with Wolfram
LambertW(30, 1E10000+i*E10000)
Result: 23015.807+i×188.487, about i×(60 pi)[code]

RE: HP49-50G, VER 12 Lambert Funktion Wk(x), k=k= 0, ±1, ±2, ±3, ±4..., x real or complex - Gil - 01-27-2024 04:44 PM

New version 16.2 of W(x).

For x real between -1/e and 0, two results: one for the upper branch, the default branch W0(x); and the other, for the lower branch branch, W-1.

For W-1, x—>0 is a special case: if x—>0- (x "slightly negative") —> W-1 = -infinity; and if x—>0+ (x "slightly positive"), W-1—> "- infinity - iPI“, but simplified in Wolfram to -infinity.

To get the answer for W (k=100, z=i×infinity), procede as follows:
1) 100 ENTER
2) 0 (for the real part)
3) LS +" number zero", ie white arrow + "key right to On-key" (to get the infinity character or corresponding the infinity number 9.99999999999E499)
4) —>NUM (to convert infinity character to the real number 9.99999999999E499)
5) command R—>C (to convert the two numbers in the stack, 0 and 9.99999999999E499, into a complex (0, 9.99999999999E499))
6) 2 —>LIST (to get a list with two objects: 100, =k the branch number, as first object of the list ; and (0, 9.99999999999E499), the complex z, as second object of the list)
7) Having the appropriate list {100 (0, 9.99999999999E499)} in stack level 1, launch the program.

You should get, as a result, "infinity + i × 401×pi/2“.

Approximate check with Wolfram
LambertW (100, i*1E500000)
gives1.1512785901131729534744838666712253701249035476598766913241 × 10^6 +629.88877992373306983756183247568673861816336497433660193071 i

And 629.8887799/(pi/2)
gives 400.99965169 —> 401

Note that, if you enter in Wolfram
LambertW (100, i*infinity),
you get infinity, and not "infinity + i × 401×pi/2“.

For more details about the LambertW function and its subtilities, look at the recent special thread in this forum on LambertW function by Albert Chan.

RE: HP49-50G, VER 15 Lambert Funktion Wk(x), k=k= 0, ±1, ±2, ±3, ±4..., x real or complex - Gil - 02-08-2024 06:43 PM

New version 16.2

Initially, only ROOT or MLSV built-in commands of the calculator were used —> and the answer then might take a while to appear — or even never appear at all.

Because of many "pitfalls", all nicely solved by Albert Chan in his thread LambertW, the program had to be modified, hopefully more or less correctly, thanks to and according to Albert Chan's implementations.

The resulting program, with its possible apparent lack of continuity in the structure and its redundancies, is not supposed to be elegant or always most efficient; it just works fine, above all with Android phone application EMU48, and reflects, with more or less success, its successive "improvements" according to the different situations.

Don't hesitate to compare the exactness of the results with Wolfram Alpha or with the very fast and compact program by John Keith after Albert Chan's suggestions.

Regarding the accuracy of the results
Entering the algebraic expression '-1/e' as input, you will get W0 & W-1 = exactly -1:
W0('-(1/e)'~-.367879441171): -1
& W-1('-(1/e)'~-.367879441171): -1.
But entering the approximate value of -1/e, ie the real number -.367879441171, will produce the double output ≠-1 W0(-.367879441171): -.999997710398
& W-1(-.367879441171): -1.00000141421.

RE: HP49-50G,VER16.2 Lambert Funktion Wk(x), k=k= 0, ±1, ±2, ±3, ±4..., x real or complex - Gil - 02-11-2024 03:17 PM

New version 16.2

Main change:
- Branch k in stack level 2
- Argument x in stack level 1

(previously it was a list)

RE: HP49-50G,VER16.2 Lambert Funktion Wk(x), k=k= 0, ±1, ±2, ±3, ±4..., x real or complex - Gil - 12-17-2024 09:15 PM

New version 17.5

Corrections depending on the given arguments.

Code:
\<< "Lambert (ver 17.5)

 2 inputs 
 k & x with branch k
   =0, \1771, \1772, \1773, \1774

 Or 1 input x possible
  for k = 0 by default
 when stack level 2
  .\=/0, \1771, \1772, \1773, \1774
  .or empty

 .With x=W(x)*e^W(x),
  W(x) unknown
 .x may be a complex #
   (a,b), even ext.big
   complex, i.e a or b
   =\1779.99999999999E499

 Output(s) 
 -1 output Wk(x): 
   . w
   . or '\oo \177 in\pi/2'
     for ext.big comp.
     (\=/Wolfram\-> \oo)
   . or '-\oo \177 in\pi'
     for x\->-0 or x\->0+
     (\=/Wolfram\-> -\oo)
 -Or 2 outputs,
  if -1/e \<=x \<=0
      & 'k=0 or k=-1':
   . W0 (x): w1 
   . W-1(x): w2
     
Formulae initial guess
.Albert Chan
.WIKIP: W-1 & -1/e\<=x\<=0

For branch k of Wk(x):
https://www.johndcook.com/blog/2021/11/15/lambert-w-branch-number/
" DROP
  IFERR OVER
  THEN 0
  ELSE \-> k
    \<< k TYPE 0 == k TYPE 28 == OR
      IF
      THEN k FP 0 == { SWAP } 0 IFTE
      ELSE 0
      END
    \>>
  END SWAP DUP \->NUM 0 RCLF 0 0 0 0 .367879441171 4.42321595524E-13 0 0 0 0 0 \-> k x0 x f fg xR xi t x0\175 R R2 H X Y Z h
  \<< -3 CF x IM DUP 'xi' STO 0 ==
    IF
    THEN x RE 'x' STO
    END x RE 0 \<= x RE R NEG \>= AND k -1 == AND
    IF
    THEN 0 'k' STO
    END -103 SF x RE DUP 'xR' STO ABS \->NUM \oo \->NUM == xi ABS \->NUM \oo \->NUM == OR
    IF
    THEN "W" k R\->I + "((" +
      CASE xi \oo \->NUM NEG ==
        THEN xR ABS \oo \->NUM ==
          IF
          THEN xR 0 > "\oo" "-\oo" IFTE
          ELSE xR DUP FP 0 == { R\->I } IFT
          END + ",-\oo))" + 4 k * 1 - DUP 0 > "+i\183" "-i\183" IFTE SWAP ABS R\->I + "\pi/2" +
        END xi \oo \->NUM ==
        THEN xR ABS \oo \->NUM ==
          IF
          THEN xR 0 > "\oo" "-\oo" IFTE
          ELSE xR DUP FP 0 == { R\->I } IFT
          END + ",\oo))" + 4 k * 1 + DUP 0 > "+i\183" "-i\183" IFTE SWAP ABS R\->I + "\pi/2" +
        END xR \oo \->NUM NEG == xi 0 < AND
        THEN "-\oo," + xi DUP FP 0 == { R\->I } IFT + "))" + 2 k * 1 - DUP 0 > "+i\183" "-i\183" IFTE SWAP ABS R\->I + "\pi" +
        END xR \oo \->NUM NEG ==
        THEN "-\oo," + xi DUP FP 0 == { R\->I } IFT + "))" + 2 k * 1 + DUP 0 > "+i\183" "-i\183" IFTE SWAP ABS R\->I + "\pi" +
        END "\oo," + xi DUP FP 0 == { R\->I } IFT + "))" + 2 k * DUP 0 > "+i\183" "-i\183" IFTE SWAP ABS R\->I + "\pi" +
      END "i\1830" "" SREPL 1 ==
      IF
      THEN DROP '\oo'
      ELSE "\oo " SWAP + "\pi" "\183\pi" SREPL DROP "\1831\183" "\183" SREPL DROP " " " <U>" SREPL DROP "<U>" + "<U>" 19 CHR 3 CHR OVER + + SREPL DROP
      END SWAP \->TAG
    ELSE -105 SF RAD '2*k*\pi' x ARG + i * \->NUM 't' STO
      CASE xi 0 == xR R NEG < AND xR -.38 > AND
        THEN '\v/(2*R*(x+R+R2))*(k+k+1)' \->NUM 'Y' STO 'R+Y*\v/(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 x Y / "W" k R\->I + "(" + x0 TYPE 9 ==
          IF
          THEN x0 + "=" +
          END x xR FP 0 == xi 0 == AND xR ABS 1.E13 < AND { R\->I } IFT + ")" + \->TAG fg STOF KILL
        END k 1 == xi 0 < AND k -1 == xi 0 \>= AND OR
        THEN x NEG LN x ABS .5 <
          IF
          THEN 0
          ELSE t 2 /
          END +
        END k 0 ==
        THEN xi 0 == xR 0 > AND
          IF
          THEN x 2 * LNP1
          ELSE x 2 * 1 + DUP 0 ==
            IF
            THEN DROP 1.6 DUP R\->C
            ELSE LN
            END
          END 2 /
        END t
      END 1 \->ARRY xi 0 == xR 0 \<= AND xR e INV NEG \>= AND k 0 == k -1 == OR AND
      IF
      THEN 'W*e^W' x - DUP 'f' STO 'W' ROT 1 GET ROOT
        CASE x e INV NEG \->NUM == x0 TYPE 9 == AND
          THEN DROP -1
          END x -.3615 \<= x e NEG INV > AND
          THEN DROP '(x+R+R2)/R' \->NUM 'H' STO '\v/(2*H)*(2*k+1)' \->NUM 'X' STO 'X*\v/(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 -
          END
        END DUP FP 0 == { R\->I } IFT "W0(" x0 TYPE 9 ==
        IF
        THEN x0 + x R NEG == "~" "=" IFTE +
        END x xR FP 0 == xi 0 == AND { R\->I } IFT + ")" + \->TAG
        CASE x 0 < x e INV NEG \>= AND
          THEN x 0 < x ABS 1.E-495 < AND
            IF
            THEN x NEG LN 1 3
              START 'x0\175' STO '(1-(LN(-x0\175)-LN(-x)))/(1+1/x0\175)' \->NUM
              NEXT
            ELSE f 'W'
              IF x 4 INV NEG \<=
              THEN '-1-\v/(2*(1+e*x))'
              ELSE 'LN(-x)-LN(-LN(-x))'
              END \->NUM ROOT x e NEG INV \->NUM == x0 TYPE 9 == AND
              IF
              THEN DROP -1
              END
            END
          END x 0 ==
          THEN -105 CF "-\oo (-\oo<U>-i\pi<U>)" "<U>" 19 CHR 3 CHR OVER + + SREPL DROP
          END -1
        END x -.3636 \<= x e NEG INV > AND
        IF
        THEN -1 'k' STO DROP '(x+R+R2)/R' \->NUM 'H' STO '\v/(2*H)*(2*k+1)' \->NUM 'X' STO 'X*\v/(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 -
        END "W-1(" x0 TYPE 9 ==
        IF
        THEN x0 + x R NEG == "~" "=" IFTE +
        END x xR FP 0 == xi 0 == AND { R\->I } IFT xR 0 ==
        IF
        THEN DROP "-0 (0+)"
        END + ")" + \->TAG 'W' PURGE
      ELSE x 0 \=/
        IF
        THEN 'k*2*\pi*i+LN(x)' 'W+LN(W)' - 1 \->ARRY [ 'W' ] ROT MSLV 1 GET DUP IM 0 == { RE } IFT UNROT DROP2 DUP \-> r
          \<< TYPE 1 ==
            IF
            THEN r C\->R DUP ABS DUP .0001 + \pi 2 / \>= SWAP \pi 2 / / \->NUM DUP 0 RND - ABS .0001 \<= AND ROT ABS .00000000001 < AND
              IF
              THEN 0 SWAP R\->C
              ELSE DROP r
              END
            ELSE r
            END
          \>>
        ELSE DROP "-\oo<U>" 2 k * k SIGN - DUP 0 >
          IF
          THEN "+i"
          ELSE "-i"
          END \-> kk i
          \<< i + kk ABS k 0 > 1 -1 IFTE + R\->I + "\pi<U>" + " (-\oo<U>" + i + kk ABS + "\pi<U>)" + "i1\pi" "i\pi" SREPL DROP "<U>" 19 CHR 3 CHR OVER + + SREPL DROP
          \>>
        END "W" k R\->I + "(" + x0 TYPE 9 ==
        IF
        THEN x0 + "=" +
        END x xR 0 == xi 0 == AND
        IF
        THEN DROP "\->-0 (\->0+)"
        ELSE xR FP 0 == xi 0 == AND xR ABS 1.E13 < AND { R\->I } IFT
        END + ")" + \->TAG
      END
    END fg STOF
    IFERR 0 DOERR
    THEN
    END
  \>>
\>>