+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Software Libraries (/forum-10.html)
+--- Forum: General Software Library (/forum-13.html)
+--- Thread: HP49-50G,VER16.2 Lambert Funktion Wk(x), k=k= 0, ±1, ±2, ±3, ±4..., x real or complex (/thread-21161.html)
HP49-50G,VER16.2 Lambert Funktion Wk(x), k=k= 0, ±1, ±2, ±3, ±4..., x real or complex - Gil - 01-14-2024 01:58 PM
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 16.2 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 - Albert Chan - 01-14-2024 02:22 PM
What algorithm were used for Wk(x) ?
Is the result guaranteed for k-th branch?
Numbering the branches of the Lambert W function
![[Image: unwindw.svg]](https://www.johndcook.com/unwindw.svg)
Perhaps of interest: thread lambertw, all branches
RE: Lambert Funktion Wk(x), k=k= 0, ±1, ±2, ±3, ±4..., x real or complex - Gil - 01-14-2024 02:50 PM
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
RE: Lambert Funktion Wk(x), k=k= 0, ±1, ±2, ±3, ±4..., x real or complex - Albert Chan - 01-14-2024 05:03 PM
(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, W0(-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))
RE: Lambert Funktion Wk(x), k=k= 0, ±1, ±2, ±3, ±4..., x real or complex - Gil - 01-14-2024 05:45 PM
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/15/lambert-w-branch-number/ ?
It seems that I get always the correct (Wolfram Alpha) output for the given k branch.
Thanks for you insight.
RE: Lambert Funktion Wk(x), k=k= 0, ±1, ±2, ±3, ±4..., x real or complex - Gil - 01-14-2024 07:13 PM
(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, W0(-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))
RE: Lambert Funktion Wk(x), k=k= 0, ±1, ±2, ±3, ±4..., x real or complex - Gil - 01-14-2024 07:23 PM
“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).
RE: Lambert Funktion Wk(x), k=k= 0, ±1, ±2, ±3, ±4..., x real or complex - Albert Chan - 01-14-2024 08:51 PM
(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
RE: Lambert Funktion Wk(x), k=k= 0, ±1, ±2, ±3, ±4..., x real or complex - Gil - 01-14-2024 09:52 PM
Well, well.
Thanks for your explanations!, Albert!
RE: Lambert Funktion Wk(x), k=k= 0, ±1, ±2, ±3, ±4..., x real or complex - Albert Chan - 01-15-2024 01:18 PM
(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;
#endCas> 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)
RE: Lambert Funktion Wk(x), k=k= 0, ±1, ±2, ±3, ±4..., x real or complex - Gil - 01-19-2024 11: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 10: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-09-2024 12:43 AM
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 09: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)
Code:
\<< "Lambert (ver 16.2)
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
\<< x IM DUP 'xi' STO 0 ==
IF
THEN x RE 'x' STO
END x 0 \<= x 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 DUP FP 0 == { 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 DUP FP 0 == { 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 DUP FP 0 == { R\->I } IFT x 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 DUP 0 ==
IF
THEN DROP "\->-0 (\->0+)"
ELSE DUP FP 0 == { R\->I } IFT
END + ")" + \->TAG
END
END fg STOF
IFERR 0 DOERR
THEN
END
\>>
\>