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

[Gil]

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/technot...bertw.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/1...ch-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.

[Albert Chan]

What algorithm were used for W_k(x) ?
Is the result guaranteed for k-th branch?

Numbering the branches of the Lambert W function



Perhaps of interest: thread lambertw, all branches

[Gil]

Yes, it's the algorithm that I found with the given reference.

The results found match well with Wolfram Alpha ones, at least for the numerous examples I tested.

Just a mistake: W('k≠0',x=0),
gives error instead of -infinity; will correct it soon.

Another special case W(k=-1,0) = -infinity, and not zero!

But, for the rest I would not guarantee anything... the Mr Maths being you!

But tell me otherwise how to implement a feasible solution for the different k branches.

Wolfram Alpha shows definitions with complex integrals.

Regards and thanks in advance for your insights, Albert.

Gil

[Albert Chan]

(01-14-2024 02:50 PM)Gil Wrote:  Yes, it's the algorithm that I found with the given reference.

If you meant given reference here, it is guaranteed to return correct branch.

x = W(a, k), we rewrite definition of branch k as function.

f(x) = x + log(x) - log(a) - 2*k*pi*I
f'(x) = 1 + 1/x

Newton: x ← x - f/f'

Solve for f(x)=0 thus guaranteed correct branch.

Without signed zero, code assumed real a == a+0*I
For illustration, here is lambertw code that support signed zeros.

lua> require'expW'
lua> I.W(-0.1, -1)
(-3.577152063957297-0*I)
lua> I.W(I.conj(-0.1), -1)
(-4.44909817870089-7.3070607892176085*I)

Technically, W_-1(-0.1) does not exist. (we have 2 different limits!)
By convention, we returned W_-1(-0.1+0*I), which does exist.

For comparison, W₀(-0.1) does exist.

lua> I.W(-0.1, 0)
(-0.11183255915896298+0*I)
lua> I.W(I.conj(-0.1), 0)
(-0.11183255915896298-0*I)

Note that 2nd result is conjugate of 1st.
In general: W(a,k) = conj(W(conj(a),-k))

[Gil]

Instead of Newton's method, I solved f(x) = 0, with the built-in HP50G MSLV command-solver, that is very fast,above all with EMU48 on the phone.

As initial value xo (complex value or not), I used the formula indicated in Takayuki HOSODA's paper.

Is it correct to use the same xo also, as I did, in the formulae (for the k branch) given in https://www.johndcook.com/blog/2021/11/1...ch-number/ ?

It seems that I get always the correct (Wolfram Alpha) output for the given k branch.

Thanks for you insight.

[Gil]

(01-14-2024 05:03 PM)Albert Chan Wrote:  

(01-14-2024 02:50 PM)Gil Wrote:  Yes, it's the algorithm that I found with the given reference.

If you meant given reference here, it is guaranteed to return correct branch.

x = W(a, k), we rewrite definition of branch k as function.

f(x) = x + log(x) - log(a) - 2*k*pi*I
f'(x) = 1 + 1/x

Newton: x ← x - f/f'

Solve for f(x)=0 thus guaranteed correct branch.

Without signed zero, code assumed real a == a+0*I
For illustration, here is lambertw code that support signed zeros.

lua> require'expW'
lua> I.W(-0.1, -1)
(-3.577152063957297-0*I)
lua> I.W(I.conj(-0.1), -1)
(-4.44909817870089-7.3070607892176085*I)

Technically, W_-1(-0.1) does not exist. (we have 2 different limits!)
By convention, we returned W_-1(-0.1+0*I), which does exist.

For comparison, W₀(-0.1) does exist.

lua> I.W(-0.1, 0)
(-0.11183255915896298+0*I)
lua> I.W(I.conj(-0.1), 0)
(-0.11183255915896298-0*I)

Note that 2nd result is conjugate of 1st.
In general: W(a,k) = conj(W(conj(a),-k))

[Gil]

“Technically, W-1(-0.1) does not exist. (we have 2 different limits!)
By convention, we returned W-1(-0.1+0*I), which does exist."

Sorry, but for the interval - 1/e<x<0, doesn't W(k=-1,x) represent —by convention ? — the lower part of the function?

For me, the solution z=(-4.4490981787,-7.30706078926)
is for W{k=-2,-.1) = (-4.4490981787,-7.30706078926).

[Albert Chan]

(01-14-2024 07:23 PM)Gil Wrote:  “Technically, W-1(-0.1) does not exist. (we have 2 different limits!)
By convention, we returned W-1(-0.1+0*I), which does exist."

Sorry, but for the interval - 1/e<x<0, doesn't W(k=-1,x) represent —by convention ? — the lower part of the function?

Yes, if you limited x to real, -1 branch is the lower part.
But, if we expand x to complex numbers, sign of 0 matters.

lua> x, k = -0.1, -1
lua> w = I.W(x, k)
lua> f = w + I.log(w) - I.log(x) - 2*k*pi*I
lua> f, w
(-4.440892098500626e-16+0*I)      (-3.577152063957297-0*I)

Note that w has -0*I, with arg(w) = -pi
Sign matters! If it were +0*I, we would get f ≈ 2*pi*I ≠ 0

Unfortunately, most calculators does not support signed zeros.
Without signed zero, my code compensate for this with customized LN(x).

We pick a result, based on some convention. W(-0.1, -1) --> W(-0.1+0*I, -1)
But it is all arbitrary! We could have pick -0*I, with +1 branch:

W(-0.1-0*I, +1) = (-3.577152063957297+0*I)

This is not that far-fetched! Example, ASIN/ACOS picked -0*I branch. It is a mess.

lua> I.acos(2)
(0-1.3169578969248166*I)
lua> I.acos(I.conj(2))
(0+1.3169578969248166*I)

Free42: 2 ACOS      → 0+1.316957896924816708625046347307968i

Branch Cuts for Complex Elementary Functions or Much Ado About Nothing' s Sign

[Gil]

Well, well.
Thanks for your explanations!, Albert!

[Albert Chan]

(01-14-2024 08:51 PM)Albert Chan Wrote:  Without signed zero, my code compensate for this with customized LN(x).

Customized LN(x) is needed to account for lost of signed zero.
This is because of the discontinuity of W_-1(-0.1 ± 0*I)

p2> from mpmath import * # functions without signed zeros
p2> mp.pretty = 1

p2> a, k = -0.1, -1
p2> x = ln(-a)             # guess = -2.302585...
p2> f = x + ln(x) - ln(a) - 2*k*pi*1j
p2> f
(0.834032445247956 + 6.28318530717959j)

f has complex parts, (x - f/f') also have complex parts, game over!
Newton's step will get across the discontinuity, unable to converge.

Customized LN(x) accounted for this, by knowing final x imag part sign.

p2> s = (im(a)<0)      # im(a)=0 treated as +0
p2> s = 1 - 2 * (k<0 or k==0 and s)
p2> s                         # final x imag part sign
-1
p2> def LN(x, bad=-s*pi): x=ln(x); return conj(x) if im(x)==bad else x
...
p2> f = x + LN(x) - ln(a) - 2*k*pi*1j
p2> f
(0.834032445247956 + 0.0j)

f still on real line, so does (x - f/f'). Iterations stay on real line, not crossing discontinuity.

---

To illustrate what happens without customized LN, we can use HP Prime

Cas> arg(-1.), arg(conj(-1.))      → [pi, pi], no signed zero support

Code:

#cas
nestlist(f,x,n)
BEGIN
  x := [x];
  WHILE((n-=1) >= 0) DO
    x := append(x, f(x[0]));
  END;
  RETURN x;
END;
#end

Cas> f := x + ln(x) - ln(a) - 2*k*pi*i
Cas> a, k := -0.1, -1
Cas> nestlist(unapply(x-f/f',x), ln(-a), 5)

[−2.30258509299,
−3.77690771946−11.1068128314*i,
−3.94994870209−0.950318443389*i,
−3.55835549121−4.10468678203e−2*i,
−3.57707753385+8.29770660431e−5*i,
−3.57707356148−8.72129018367*i]

For comparison, this is Python code with customized LN(x)

p2> x = W(-0.1, -1, True)
-2.30258509299405
(-3.77690771946225 + 0.0j)
(-3.57912417640232 + 0.0j)
(-3.57715227469581 + 0.0j)
(-3.5771520639573 + 0.0j)

[Gil]

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]

[Gil]

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.

[Gil]

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.

[Gil]

New version 16.2

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

(previously it was a list)

[Gil]

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 W
k
(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':
   . W
0
 (x): 
w1
 
   . W-
1
(x): 
w2

     
Formulae initial guess
.Albert Chan
.WIKIP: W-1 & -1/e\<=x\<=0

For branch 
k
 of W
k
(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
  \>>
\>>