HP49-50G,VER16.2 Lambert Funktion Wk(x), k=k= 0, ±1, ±2, ±3, ±4..., x real or complex
01-14-2024, 01:58 PM (This post was last modified: 02-11-2024 09:13 PM by Gil.)
Post: #1
 Gil Senior Member Posts: 590 Joined: Oct 2019
HP49-50G,VER16.2 Lambert Funktion Wk(x), k=k= 0, ±1, ±2, ±3, ±4..., x real or complex
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/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.

Attached File(s)
Lambert.16.2.hp (Size: 5.8 KB / Downloads: 0)
01-14-2024, 02:22 PM
Post: #2
 Albert Chan Senior Member Posts: 2,555 Joined: Jul 2018
RE: Lambert Funktion Wk(x), k=k= 0, ±1, ±2, ±3, ±4..., x real or complex
What algorithm were used for Wk(x) ?
Is the result guaranteed for k-th branch?

Numbering the branches of the Lambert W function

Perhaps of interest: thread lambertw, all branches
01-14-2024, 02:50 PM (This post was last modified: 01-14-2024 09:39 PM by Gil.)
Post: #3
 Gil Senior Member Posts: 590 Joined: Oct 2019
RE: Lambert Funktion Wk(x), k=k= 0, ±1, ±2, ±3, ±4..., x real or complex
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
01-14-2024, 05:03 PM
Post: #4
 Albert Chan Senior Member Posts: 2,555 Joined: Jul 2018
RE: Lambert Funktion Wk(x), k=k= 0, ±1, ±2, ±3, ±4..., x real or complex
(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))
01-14-2024, 05:45 PM (This post was last modified: 01-14-2024 06:24 PM by Gil.)
Post: #5
 Gil Senior Member Posts: 590 Joined: Oct 2019
RE: Lambert Funktion Wk(x), k=k= 0, ±1, ±2, ±3, ±4..., x real or complex
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.
01-14-2024, 07:13 PM
Post: #6
 Gil Senior Member Posts: 590 Joined: Oct 2019
RE: Lambert Funktion Wk(x), k=k= 0, ±1, ±2, ±3, ±4..., x real or complex
(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))
01-14-2024, 07:23 PM (This post was last modified: 01-14-2024 07:31 PM by Gil.)
Post: #7
 Gil Senior Member Posts: 590 Joined: Oct 2019
RE: Lambert Funktion Wk(x), k=k= 0, ±1, ±2, ±3, ±4..., x real or complex
“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).
01-14-2024, 08:51 PM
Post: #8
 Albert Chan Senior Member Posts: 2,555 Joined: Jul 2018
RE: Lambert Funktion Wk(x), k=k= 0, ±1, ±2, ±3, ±4..., x real or complex
(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
01-14-2024, 09:52 PM
Post: #9
 Gil Senior Member Posts: 590 Joined: Oct 2019
RE: Lambert Funktion Wk(x), k=k= 0, ±1, ±2, ±3, ±4..., x real or complex
Well, well.
Thanks for your explanations!, Albert!
01-15-2024, 01:18 PM
Post: #10
 Albert Chan Senior Member Posts: 2,555 Joined: Jul 2018
RE: Lambert Funktion Wk(x), k=k= 0, ±1, ±2, ±3, ±4..., x real or complex
(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)
01-19-2024, 11:05 PM (This post was last modified: 02-11-2024 09:08 PM by Gil.)
Post: #11
 Gil Senior Member Posts: 590 Joined: Oct 2019
RE: Lambert Funktion Wk(x), k=k= 0, ±1, ±2, ±3, ±4..., x real or complex
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-ipi" }

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]

Attached File(s)
Lambert.16.2.hp (Size: 5.8 KB / Downloads: 0)
01-27-2024, 10:44 PM (This post was last modified: 02-11-2024 09:11 PM by Gil.)
Post: #12
 Gil Senior Member Posts: 590 Joined: Oct 2019
RE: HP49-50G, VER 12 Lambert Funktion Wk(x), k=k= 0, ±1, ±2, ±3, ±4..., x real or complex
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.

Attached File(s)
Lambert.16.2.hp (Size: 5.8 KB / Downloads: 0)
02-09-2024, 12:43 AM (This post was last modified: 02-11-2024 09:12 PM by Gil.)
Post: #13
 Gil Senior Member Posts: 590 Joined: Oct 2019
RE: HP49-50G, VER 15 Lambert Funktion Wk(x), k=k= 0, ±1, ±2, ±3, ±4..., x real or complex
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.
02-11-2024, 09:17 PM
Post: #14
 Gil Senior Member Posts: 590 Joined: Oct 2019
RE: HP49-50G,VER16.2 Lambert Funktion Wk(x), k=k= 0, ±1, ±2, ±3, ±4..., x real or complex
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 Wk(x):     . w    . or '\oo \177 in\pi/2'      for ext.big comp.      (\=/Wolfram\-> \oo)    . or '-\oo \177 in\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   \<< 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>)" + "i1\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   \>> \>
>

Attached File(s)
Lambert.16.2.hp (Size: 5.8 KB / Downloads: 0)
 « Next Oldest | Next Newest »

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