lambertw, all branches

[Albert Chan]

(01-23-2024 04:57 PM)Gil Wrote:  We accept...
You decided to accept that the real part in question is indeed zero...

That's not what I mean.

I specifically say real part will *not* be zero, even with more precision.
This is because of numerical limits, a = 3/2*pi is only an approximation.

I accepted limited precision limits, and confirm symbolic result (if pi is true pi, ε is true 0 ...)

see https://en.wikipedia.org/wiki/Experimental_mathematics

---

Test for W_k(-pi/2), k=0 or -1, is to make sure we don't get stuck in loops.

That's the reason for f = x + ln(x) - T, where T = ln_k(a) = (ln(a) + 2*k*pi*I)

(04-09-2023 03:59 AM)Albert Chan Wrote:  * f formula rearranged, to cause catastrophic cancellation, on purpose!

With cancellation, final x might not be as accurate, but loops will terminate.
This is then used as an extremely good guess, with another f.

(04-19-2023 01:30 AM)Albert Chan Wrote:  We can use Newton's method, f = x*e^x - a, for finishing touches

[Gil]

Any real tiny value will always do the job:
(i×b+ real_tiny) × EXP ((i×b+ real_tiny)
=(i×b+ real_tiny) × EXP (i×b) × EXP (real_tiny)
=(i×b+ real_tiny) × EXP (i×b) × 1
=(i×b+ real_tiny) × EXP (i×b)
=(i×b) × EXP (i×b)
=B×EXP(B)

W(B×EXP(B) =B=(i×b)
==(i×b+ real_tiny)

[Gil]

Again, W(k=-1,a=+PI/2 )
seems impossible to give — at least with my HP50G program —the nice answer (0-i×PI),
returning "tiny real - i×PI"

And also W(k=-1, a=-4×PI/2 )
HP output is "tiny - i×2PI"

And also W(k=-2,a=i×6PI/2)
HP output is "tiny - i×3PI"

LambertW (k=-3,0+10*i*pi/2)
HP output is "tiny - i×5PI"

LambertW (k=-1000,i*3998 *pi/2)
Output Wolfram Alpha
i ×
6280.04371452599668368682412317572626551014162935083653612891424=

Output HP:
:W-1000
(
'i*3998*(pi/2)'=(0.,6280.04371455)
)
: (3.6616915941E-12,-6280.04371452)

And LambertW (-5, -20*i*pi/2)
=-i*10*pi/2 =-i*20*pi/2 = a


Pattern
W(k=n, n<0, a=i×(-n×4-2)) = a —> complex with no real part
W(k=n, n>0, a=i×(n×4)) = a —> complex with no real part

All the results are quite OK, though unsatisfactory, as already said.

Rule of thumb:
IF {imaginary (result) = integer≠0 × PI
and abs (real (result) < 1E-10}
THEN result = i × imaginary (result)
END

And more generally
With W0(-pi/2)=0+i×Pi/2
With W-1(-pi/2)=0-i×Pi/2
IF {imaginary (result) = integer≠0 × PI/2
and abs (real (result) < 1E-10}
THEN result = i × imaginary (result)
END

Is that deduction or assumption correct?

Of course, because of the roundings and the fact that we never have the real value of pi, we might never — or rather will never — have a real integer multiple of pi/2.

Final program steps/test
M0=abs(Imaginary part of the final result) /(pi/2)
M1= round (M0, 0) —> to get the nearest integer to M0
M2= real part of the final result
IF {abs (M0-M1)<0.0001 and abs (M2) < 1E-11}
THEN
M0 is considered as a multiple of (pi/2)
& real part of the result, found to be "tiny",
is to be set to exactly = zero:
M2=0
final result = M2+i×(Imaginary part of the final result)
(ELSE
The found result is OK and is not to be "tampered":
final result = final result)
END

Check a=real, multiple of pi/2, with no imaginary part
Tests with Wolfram
Conclusion
W(k, k<0 or k>0, a= real =(-4k-1)×pi/2=-i×a (and zero for the real part)

Overall conclusion
Too many tests about the input on k & a=±i × multiple of pi/2.
Easier to test the final result and correct the real part of that final result (tiny into exactly zero), as described above in Final program steps/test, if necessary or wished, ie if & only if when both conditions abs (M0-M1)<0.0001 and abs (M2) < 1e-11 are fulfilled.

[Gil]

Always regarding the Albert Chan's
result = "tiny real part" ± i × n(PI/2), n≠0

here is a possible "patch" for HP50G to fix up such an issue,
where next, final result will be then: 0 ± i × n(PI/2).

Code:

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

[Albert Chan]

Gil's Trivia (previous posts) summary

x = (2n)*pi*I
a = x*e^x = x*cis(0) = x
k = (im(x) + arg(x) - arg(a)) / (2*pi) = (2n)*pi / (2*pi) = n

W_k((2k)*pi*I) = (2k)*pi*I

W₁(2*pi*I) = 2*pi*I
W_-1(-2*pi*I) = -2*pi*I      // conjugate symmetry

---

x = (2n+1/2)*pi*I
a = x*e^x = x*cis(pi/2) = x*I

if n≥0 then k = ((2n+1/2)*pi + pi/2 - pi) / (2*pi) = n
if n<0 then k = ((2n+1/2)*pi − pi/2 − 0) / (2*pi) = n

W_k(-(2k+1/2)*pi) = (2k+1/2)*pi*I

W₁(-5/2*pi) = 5/2*pi*I
W_-1(3/2*pi) = -3/2*pi*I

We can use conjugate symmetry for version with RHS = -x

W_-k(-(2k+1/2)*pi - 0*I) = -(2k+1/2)*pi*I

W_-1(-5/2*pi - 0*I) = -5/2*pi*I
W₁(3/2*pi - 0*I) = 3/2*pi*I

But this required signed zero. I had another version without using signed zero (last formula)
Example, W₁(3/2*pi) = 3/2*pi*I too, because it matched last formula form.

---

odd = 2n - sign(n)      // n ≠ 0

x = odd*pi*I = (2n-sign(n))*pi*I
a = x*e^x = x*cis(odd*pi) = −x

if n>0 then k = ((2n−1)*pi + pi/2 + pi/2) / (2*pi) = n
if n<0 then k = ((2n+1)*pi − pi/2 − pi/2) / (2*pi) = n

a = -x should have -0 real part, but here, it does not matter.

k = (odd+sign(odd))/2
W_k(-odd*pi*I) = odd*pi*I

W₁(-pi*I) = pi*I
W_-1(pi*I) = -pi*I      // conjugate symmetry

---

For im(x) mod (2*pi) = 3/2*pi --> im(x)/pi = odd + 1/2

x = (odd+1/2)*pi*I = (2n-sign(n)+1/2)*pi*I      // n ≠ 0
a = x*e^x = x*cis(3/2*pi) = x/I                        // im(a) = +0

if n>0 then k = ((2n-1/2)*pi + pi/2 - 0) / (2*pi) = n
if n<0 then k = ((2n+3/2)*pi - pi/2 - pi) / (2*pi) = n

k = (odd+sign(odd))/2
W_k((odd+1/2)*pi) = (odd+1/2)*pi*I

W_-2(-5/2*pi) = -5/2*pi*I
W₁(3/2*pi) = 3/2*pi*I

[Gil]

I try with Wolfram
W0 (-1.2E-323)
and get the correct answer, ie about -1.2E-323.
Curiously, it offers the possibility to get more digits, two times, but does not show them at the second request.

More strange:
when I ask W0 (-1.2E-324),
I don't get now the expected approximate value -1.2E-324, but
directly zero.

Why this distinction at -1E-323/E-324?

[Gil]

In your second post, Albert, you wrote

RE: lambertw, all branches
f = x + ln(x) - ln(a) - 2*k*pi*I

Solve for f=0 is easy if x is big.
For |k| > 1, we can solve directly, since x imaginery part ≈ 2*k*pi (see OP example)

For k=-1, solve for f=0 is especially difficult.
Fortunately, we can flip it, to solve for k=1 instead.
W(z, k) == conj(W(conj(z), -k)
In fact, I flip all negative k's. So now, k is non-negative.

But W-1(-1E-330)=
-697.32277629546016099540752740546566360568199201428307824605649289390557552642032068023​149199134240178698554915132

And W1(-1E-330)=
-697.32281706295494653348757314785952050772713371449039168 +6.2922084418351405256893689216267218827742283825238692107 i

It seems I missed a point when I may use
W(z, k) == conj(W(conj(z), -k).

Thanks for your lights, Albert.

[Gil]

I think that
W(z, k) == conj(W(conj(z), -k)
is true when z is not a real x,
x/z belonging to [-1/e, 0).

[Albert Chan]

Hi, Gil

Complex Conjugate formula is always true, if we have signed zero.
With signed-zero, arg(z) + arg(conj(z)) = 0

Let z = r * cis(θ)      // polar form, θ = arg(z)

ln(conj(z)) = ln(r*cis(-θ)) = ln(r) - θ*I = conj(ln(r) + θ*I) = conj(ln(z))

Read from right to left, conj can go inside ln

x = W_k(a)
x + ln(x) = ln(a) + 2*k*pi*I
conj(x + ln(x)) = conj(ln(a) + 2*k*pi*I)
conj(x) + ln(conj(x)) = ln(conj(a)) + 2*(-k)*pi*I
conj(x) = W_-k(conj(a))

--> W_k(a) = conj(W_-k(conj(a)))

[Gil]

Thanks, Albert.

And how do you explain these different results?

W-1(-1E-330)=
-697.32277629546016099540752740546566360568199201428307824605649289390557552642032068023​149199134240178698554915132

And W1(-1E-330)=
-697.32281706295494653348757314785952050772713371449039168 +6.2922084418351405256893689216267218827742283825238692107 i

Are they correct?

[Albert Chan]

Quote:With signed-zero, arg(z) + arg(conj(z)) = 0

Without signed-zero, above might not hold.
If a is negative real, conj(a) returns itself, we get sum = pi + pi ≠0

In this case, this does not hold either: W_k(a) = conj(W_-k(conj(a)))

One way to get sum=0 is to to use ±eps for ±0 (★)

p2> z = mpc("-1E-330", "1E-9999")

p2> W(z, -1)
(-766.494908744112 - 1.00130634441663e-9669j)
p2> conj(W(conj(z), 1))
(-766.494908744112 - 1.00130634441663e-9669j)

p2> W(z, 1)
(-766.494942472642 + 6.2913931263013j)
p2> conj(W(conj(z), -1))
(-766.494942472642 + 6.2913931263013j)


(★) with eps hack, we may get into another ugly situation.

What if eps underflowed to zero? Then, the sign is lost!

Luckily, mpmath (by design) will not overflow/underflow.

https://mpmath.org/doc/current/technical.html Wrote:Mpmath uses arbitrary precision integers for both the mantissa and the exponent,
so numbers can be as large in magnitude as permitted by the computer’s memory.

Still, I think better solution is to properly add signed-zero to mpmath.

fredrik-johansson Wrote:Should negative zero (-0) be supported? It is virtually supported already; the object 'fnzero' is defined in libmpf.py (though not used anywhere). It just needs to be supported in standard operations, and elsewhere "x ==fzero" needs to be replaced with something else so that -0 is not forwarded to places where it shouldn't go ...

[Gil]

Thanks, but what about the results given by Wolfram?

I ask, because I never used your flip trick for k a negative integer.

Fortunately?

[Albert Chan]

(01-25-2024 05:14 PM)Gil Wrote:  Thanks, but what about the results given by Wolfram?

I ask, because I never used your flip trick for k a negative integer.

Wolfram's answer is correct.
Why would you expect W_-1(-1E-330) == W₊₁(-1E-330) ?

x + ln(x) = ln(a) + 2*k*pi*I

With k = ±1, same a, RHS have difference of 4*pi*I ! Of course solved x is different.

I already explained why flip trick don't work in this case (see my previous post)

[Gil]

I meant the same real part.

As I said, though reading your posts, I stuck to the original branch sign for calculation — and did well.

[Albert Chan]

Note that W curve is only in 1st and 3rd quadrant.
If x is in ℝ (implied a in ℝ too), they have same sign, then ln(x/a) = ln(|x|) - ln(|a|)

This give me an idea!

Code:

def LN(x, bad=-s*pi):
    x = ln(x)
    return conj(x) if im(x)==bad else x

My Python code used customized LN, to track signed zero, to give correct imag parts.
Example, to simulate ln(-1-0*I), set s=-1 --> LN(-1) = conj(pi*j) = -pi*j

But, what if re(a) were never negative to begin with?
Then, LN is not needed! Unlike ±pi, ±0 are numerically the same thing!

Let s = -1 if re(a)<0 else 1

x * e^x = a
(s*x) * e^x = s*a
ln(s*x) + x = ln(s*a) = B

p2> f = lambda x: show(x) + ln(s*x) - B
p2> df = lambda x: 1 + 1/x
p2> def show(x): print nstr(x, n=10); return x
...
p2> a = mpc(-0.1)
p2> s = -1 if re(a)<0 else 1
p2> B = ln(s*a)

p2> W(a, -1)
(-3.5771520639573 + 0.0j)
p2> findroot(f, x0=B, df=df, solver='newton')
(-2.302585093 + 0.0j)
(-3.776907719 + 0.0j)
(-3.579124176 + 0.0j)
(-3.577152275 + 0.0j)
(-3.577152064 + 0.0j)
(-3.577152064 + 0.0j)
(-3.5771520639573 + 0.0j)

This is nothing special, iterations are exactly the same as W code.
But let's try a *bad* guess, x0=B+j, crossed discontinuity!
(this guess for for W(a,-1) would get stuck in endless loops.)

p2> findroot(f, x0=B+j, df=df, solver='newton')
(-2.302585093 + 1.0j)
(-3.44870968 - 0.2167163462j)
(-3.575061443 + 0.002993043143j)
(-3.577151814 - 6.790747928e-7j)
(-3.577152064 + 1.841873128e-14j)
(-3.577152064 - 4.31950186e-29j)
(-3.577152064 + 0.0j)
(-3.5771520639573 + 0.0j)

We can use better solver, without worrying about over-shooting!

p2> findroot(f, x0=B, df=df, solver='halley')
(-2.302585093 + 0.0j)
(-3.48604165 + 0.0j)
(-3.577141461 + 0.0j)
(-3.577152064 + 0.0j)
(-3.577152064 + 0.0j)
(-3.5771520639573 + 0.0j)

p2> findroot(f, x0=B+j, df=df, solver='halley')
(-2.302585093 + 1.0j)
(-3.730387658 - 0.002761341752j)
(-3.577193828 - 2.173902656e-6j)
(-3.577152064 - 1.483011857e-16j)
(-3.577152064 - 2.350988702e-38j)
(-3.577152064 + 0.0j)
(-3.5771520639573 + 0.0j)

Comment:

What we really wanted is s = -1 if re(x)<0 else 1
This keep |arg(s*x)| closer to 0 then pi.

It happens that if final x is real, sign(x) = sign(a)
Unfortunately, for complex a, we don't know what s is ...

[Albert Chan]

(01-28-2024 11:18 PM)Albert Chan Wrote:  But let's try a *bad* guess, x0=B+j, crossed discontinuity!
(this guess for for W(a,-1) would get stuck in endless loops.)

This is a safe way to fix stuck in loops problem for bad guess.
We know ahead of time final sign of imag part is s or 0 (if x is real)
If any iterations does not match the sign, LN will nudge it back.

Code:

def LN(x, bad=-s):
    x = ln(x)
    return conj(x) if sign(im(x))==bad else x

For debugging purpose, lets add a way to input user-supplied guess.
Put this right before Newton while 1 loop

Code:

if verbal and verbal!=True: x = verbal

Repeat the same test, with bad guess

p2> a = mpc(-0.1)
p2> B = ln(-a)
p2> W(a, -1, B+j)

p2> W.W(a, -1, B+j)
(-2.30258509299405 + 1.0j)
(-3.14484755109194 - 1.43195611385943j)
(-3.45450302775194 - 0.017653789718918j)
(-3.57799843751793 + 0.000256987554174169j)
(-3.57715209366168 - 0.000199346248602327j)
(-3.57715206180199 - 8.48877473305296e-13j)
(-3.5771520639573 + 5.31405873378874e-16j)
(-3.5771520639573 + 0.0j)

It is hard to nudge back to real line, but it work!

I think guess is good enough without this trick, so it is not updated to version 3.
Perhaps I would put it to next version.

[Albert Chan]

Python W code (post 16)

Quote:s = (im(a) < 0)
small_imag = k==1 and s or k==-1 and not s
...
s = 1 - 2 * (k<0 or k==0 and s)

It is very convoluted. Let's make a table

Code:

old s       updated s
im(a)<0  k>0   k<0   k=0                   
True     +1    -1    -1
False    +1    -1    +1

It is simpler to have old s = - sign(im(a) or 1) = ±1
Then, small_imag = (k==s), too simple to even need a variable for it.

Previous post customized LN and ability to enter W guess also included.

Python W, version 4

Code:

from mpmath import *    # no signed zero support

r, err = 0.36787944117144233, -1.2428753672788363E-17

def W(a, k=0, verbal=False):
    h, s = a+r+err, -sign(im(a) or 1)
    if abs(h)<.25 and (k==0 or k==s):
        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
            x = y
            y = (y+a)/log1p((y-r-err)/r)
            if abs((y-x)/y) < 1e-9: return a/y
    def LN(x, bad = -sign(k) or s):
        x = ln(x)
        return conj(x) if sign(im(x))==bad else x
    A, T = abs(a), mpc(0,2*k*pi+arg(a))
    x = T          # W rough guess
    if k==0: x = (A+100)/(A+50); x = log1p(a*x)/x
    if k==s: x = LN(-a) + (0 if A<.5 else T/2)
    T = T + log(A) # = ln_k(a)
    if not isfinite(A): return T
    if A==0: return a if k==0 else T - sign(k)*pi*j
    if verbal and verbal!=True: x = verbal
    while 1:
        if verbal: print x
        h = x/(x+1) * (x + LN(x) - T)
        x = x - h
        if not (abs(h/x) >= 1e-6): break
    if verbal: print x
    h = (x-a*exp(-x))/(x+1) # final touch
    return x-h if isfinite(h) else x

[Albert Chan]

(02-01-2024 04:16 PM)Albert Chan Wrote:  Previous post customized LN and ability to enter W guess also included.

Python W nudge imag part from inside LN, because it is free.
Customized LN need to fix ±0j effect to phase angle anyway.

Lua I.W (post 14) work differently, by forcing positive x imag part, with "flips"
The reason is signed-zero arithmetic preferred +0. It is much easier to maintain +0 than -0.

If I were to implement equivalent safe feature, I would just force this for each iteration.
This is a safe patch, since it is what is assumed. Just place it before Newton correction.

Code:

x = I(I.real(x), abs(I.imag(x)))

For the ability to user to enter W guess, add this before Newton loops.

Code:

if verbal and verbal~=true then x = flip(verbal) end

With these patches to post 14, this is Lua I.W version 4

The patch does not mean we can use any guess.
The problem is again ±0*I, with negative real part.

We get phase angle of ±pi, difference of ±2*pi, exactly matched branch ±1
This is why branch 0, ±1 with small imag parts *must* have good guess.