(28/48/50) Lambert W Function

[John Keith]

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...

[John Keith]

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
\>>

[Albert Chan]

(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 W_k(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')

[John Keith]

(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 W_k(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.

[Albert Chan]

(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 W₁: 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')

[Albert Chan]

(04-02-2023 11:12 PM)Albert Chan Wrote:  lua> a, k = 3+4*I, 5
lua> x = 2*k*pi*I -- guess for W_k(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 -- W_k(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, W_k(a) = complex(r*c, r*sqrt(1-c*c))

[Albert Chan]

(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 W_k(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)) -- = W₅(3+4*I)
(-1.8170058918466274+30.713334137004896*I)

I setup function comparing angles. For |k| ≥ 2, r only have 1 real root.

[John Keith]

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
\>>

[Gil]

In your post 18,
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=?

[Gil]

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?

[Albert Chan]

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.

[Gil]

Thanks.

And why, in code post 29, did you define Y=R+R*X, that seems never to be used again?

[Albert Chan]

(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

[Gil]

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

[Albert Chan]

(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.

There is no wrong answer. Pick your preferred choice.

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) = ln_k(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.

[Gil]

Thanks.

Strange anyway, a = -.36787944118, that in my calculator the until equality is never reached with method exposed in post 11.—> endless loops.

[Albert Chan]

(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

https://www.hpmuseum.org/forum/thread-55...#pid111682

[Gil]

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

[Albert Chan]

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

[Gil]

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.