Post Reply 
lambertw, all branches
04-23-2023, 04:40 PM (This post was last modified: 02-01-2024 03:35 PM by Albert Chan.)
Post: #16
RE: lambertw, all branches
complex.W version 3, adapted for system without signed zero.
This time I made a customized LN(x) to honor signed zero.

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, im(a)<0
    small_imag = k==1 and s or k==-1 and not s
    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-9: return a/y
    s = 1 - 2 * (k<0 or k==0 and s)
    def LN(x, bad=-s*pi):
        x = ln(x)
        return conj(x) if 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 small_imag: 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
    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

Trivia: W(a < -1/e) is complex, W0(a) = conj(W-1(a))

x + ln(x) = lnk(a)

ln0(a) = ln(|a|) - pi*j + 2*pi*j = ln1(conj(a))

W0(a) = W1(conj(a)) = conj(W-1(a))

>>> a = -pi/2
>>> W(a, 0, True)
(0.375191803678016 + 1.59508875325509j)
(0.0128176527641888 + 1.58759064221342j)
(-5.58573510823744e-5 + 1.57084664181356j)
(6.00809277877961e-10 + 1.57079632603836j)
(5.72587285769658e-17 + 1.5707963267949j)
(1.30059124689582e-17 + 1.5707963267949j)
>>> W(a, -1, True)
(0.451582705289455 - 1.5707963267949j)
(0.0200654830289762 - 1.59320451813763j)
(-9.57300550001525e-5 - 1.57091376376939j)
(1.39642133562001e-9 - 1.5707963231279j)
(-3.14444301571107e-17 - 1.5707963267949j)
(-5.01152163329281e-17 - 1.5707963267949j)

Confirm x * exp(x) = a, both at the same time.

(±pi/2*j) * exp(±pi/2*j) = (±pi/2*j) * (±j) = pi/2 * j² = -pi/2

Udpate 1/19/24: added |a| = ∞ or 0 edge cases

p2> W(inf), W(-inf)
((+inf + 0.0j), (+inf + 3.14159265358979j))

p2> W(0, -1)
(-inf - 3.14159265358979j)

Default setup is 0 == +0, that is why we get W-1(+0) = -∞ - pi*I
To get W-1(0) = -∞, we need special case for k=-1, 0 == −0 + 0*I

Real part -0, Imag part +0 ... it is just too confusing.
I am keeping 0 == +0, for consistency. ln(0) = ln(+0) = -∞

--> Wk(+0) = 0 if k==0 else -∞ + (2*k-sign(k))*pi*I

Of course, we are free to patch it to make it pick what we wanted.
But then, what is Wk(0) in general ? Is 0 still meant +0 + 0*I ?
This is a weakness of not having signed zero.

Update 2/01/24

lyuka branch relative limit reduced, from 1e-12 to 1e-9 (see post #14)
Find all posts by this user
Quote this message in a reply
Post Reply 


Messages In This Thread
lambertw, all branches - Albert Chan - 04-07-2023, 01:24 PM
RE: lambertw, all branches - Albert Chan - 04-07-2023, 02:47 PM
RE: lambertw, all branches - Albert Chan - 04-19-2023, 01:30 AM
RE: lambertw, all branches - pier4r - 04-07-2023, 06:04 PM
RE: lambertw, all branches - Albert Chan - 04-07-2023, 07:54 PM
RE: lambertw, all branches - Albert Chan - 04-08-2023, 03:21 PM
RE: lambertw, all branches - Albert Chan - 04-08-2023, 05:54 PM
RE: lambertw, all branches - Albert Chan - 04-07-2023, 08:40 PM
RE: lambertw, all branches - Albert Chan - 04-09-2023, 03:59 AM
RE: lambertw, all branches - Albert Chan - 04-09-2023, 04:36 PM
RE: lambertw, all branches - Albert Chan - 04-10-2023, 04:44 PM
RE: lambertw, all branches - Albert Chan - 04-10-2023, 06:47 PM
RE: lambertw, all branches - Albert Chan - 04-13-2023, 03:03 PM
RE: lambertw, all branches - floppy - 04-13-2023, 04:14 PM
RE: lambertw, all branches - Albert Chan - 04-23-2023, 02:49 PM
RE: lambertw, all branches - Albert Chan - 04-23-2023 04:40 PM
RE: lambertw, all branches - Albert Chan - 01-19-2024, 04:14 PM
RE: lambertw, all branches - Albert Chan - 01-20-2024, 04:48 PM
RE: lambertw, all branches - Gil - 01-20-2024, 10:52 PM
RE: lambertw, all branches - Albert Chan - 01-21-2024, 01:14 AM
RE: lambertw, all branches - Albert Chan - 01-21-2024, 01:54 AM
RE: lambertw, all branches - Gil - 01-21-2024, 01:53 PM
RE: lambertw, all branches - Albert Chan - 01-21-2024, 04:19 PM
RE: lambertw, all branches - Gil - 01-21-2024, 04:35 PM
RE: lambertw, all branches - Albert Chan - 01-21-2024, 06:03 PM
RE: lambertw, all branches - Albert Chan - 01-21-2024, 07:01 PM
RE: lambertw, all branches - Gil - 01-21-2024, 07:30 PM
RE: lambertw, all branches - Gil - 01-21-2024, 08:39 PM
RE: lambertw, all branches - Albert Chan - 01-21-2024, 10:06 PM
RE: lambertw, all branches - Gil - 01-21-2024, 09:51 PM
RE: lambertw, all branches - Gil - 01-21-2024, 10:56 PM
RE: lambertw, all branches - Albert Chan - 01-22-2024, 01:34 AM
RE: lambertw, all branches - Gil - 01-21-2024, 11:15 PM
RE: lambertw, all branches - Gil - 01-22-2024, 06:09 PM
RE: lambertw, all branches - Albert Chan - 01-22-2024, 07:29 PM
RE: lambertw, all branches - Gil - 01-22-2024, 11:33 PM
RE: lambertw, all branches - Albert Chan - 01-23-2024, 02:32 AM
RE: lambertw, all branches - Gil - 01-23-2024, 02:35 PM
RE: lambertw, all branches - Albert Chan - 01-23-2024, 03:54 PM
RE: lambertw, all branches - Gil - 01-23-2024, 04:57 PM
RE: lambertw, all branches - Albert Chan - 01-23-2024, 06:17 PM
RE: lambertw, all branches - Gil - 01-23-2024, 06:44 PM
RE: lambertw, all branches - Gil - 01-23-2024, 11:00 PM
RE: lambertw, all branches - Gil - 01-24-2024, 03:18 PM
RE: lambertw, all branches - Albert Chan - 01-24-2024, 08:53 PM
RE: lambertw, all branches - Gil - 01-25-2024, 12:37 AM
RE: lambertw, all branches - Gil - 01-25-2024, 01:10 AM
RE: lambertw, all branches - Gil - 01-25-2024, 03:04 AM
RE: lambertw, all branches - Albert Chan - 01-25-2024, 07:02 AM
RE: lambertw, all branches - Gil - 01-25-2024, 10:09 AM
RE: lambertw, all branches - Albert Chan - 01-25-2024, 04:13 PM
RE: lambertw, all branches - Gil - 01-25-2024, 05:14 PM
RE: lambertw, all branches - Albert Chan - 01-25-2024, 05:57 PM
RE: lambertw, all branches - Gil - 01-25-2024, 06:19 PM
RE: lambertw, all branches - Albert Chan - 01-28-2024, 11:18 PM
RE: lambertw, all branches - Albert Chan - 02-01-2024, 02:17 AM
RE: lambertw, all branches - Albert Chan - 02-01-2024, 04:16 PM
RE: lambertw, all branches - Albert Chan - 02-02-2024, 11:49 AM



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