lambertw, all branches
04-07-2023, 01:24 PM
Post: #1
 Albert Chan Senior Member Posts: 2,708 Joined: Jul 2018
lambertw, all branches
I am starting a new thread, to discuss issues of solving W(a, k=0), extended to complex numbers.

The thread is based on this post, which supposed guaranteed correct branch, if f=0
(04-02-2023 11:12 PM)Albert Chan Wrote:
(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 Wk(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')
04-07-2023, 02:47 PM
Post: #2
 Albert Chan Senior Member Posts: 2,708 Joined: Jul 2018
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.

(04-02-2023 11:12 PM)Albert Chan Wrote:  W code is perfect for other branches! (expW code for W branch 0, -1)

Again, for now on, we assumed k >= 0
f=0 may still have problem converging if k = 0 or 1, especially for small x
It does not matter if we have a really good guess.

It is why expW code is needed. It handle this issue nicely. (*)
expW code also handled accuracy issue elegantly, when a is close to branch point, -1/e

Problem is expW code lost phase information, y = e^(x + 2*pi*k*I) = e^x
expW can only get W(a, ±1) with the smaller imaginery part.

expW code can get this:

>>> lambertw(-1-.1j, 1)
(-0.373678737258224 + 1.41169809383182j)

But not this:

>>> lambertw(-1+.1j, 1)
(-2.04424331708213 + 7.48779784363809j)

Note that negative a imag part corresponded to smaller W(a, k=1) imag part.
It can be deduced from f

f = x + ln(x) - ln(a) - 2*k*pi*I = 0
x.imag + phase(x) = 2*k*pi + phase(a)

z.imag = abs(z) * sin(phase(z))

a.imag < 0      ⇒ phase(a) < 0      ⇒ RHS smaller      ⇒ x.imag smaller

(*) Based from tests, solve for y*ln(y)=a is better than x*e^x=a, even with help of Halley's method.
mpmath.lambertw make it work only because it throw in extra precisions to compensate.
And, it required a lot of code to get a good starting guess.
04-07-2023, 06:04 PM
Post: #3
 pier4r Senior Member Posts: 2,248 Joined: Nov 2014
RE: lambertw, all branches
Side observation. Wouldn't be better to format such posts a bit more and post them in the Article forum? I think they can be more valuable over time if they can get a bit more structured.

Wikis are great, Contribute :)
04-07-2023, 07:54 PM
Post: #4
 Albert Chan Senior Member Posts: 2,708 Joined: Jul 2018
RE: lambertw, all branches
Hi, pier4r

I am still learning, and math is not up to par into article section, where no feedback allowed.
Hopefully, someone reading this may suggest a better way.

(04-07-2023 02:47 PM)Albert Chan Wrote:  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)

Example, as I am writing code, I realized k "flipping" is unnecessary.
With equivalent starting guess, convergence is the same.

It is just here, "equivalent" may mean a difference of +0*I vs -0*I

Below is work in progress code, without flipping k to non-negative.
It simply solve f=0 by Newton's method, with a supplied guess x

lua> a, k = -0.1, -1
lua> x = I.log(-a)
lua> I.W(a, k, true, x)
(-2.3025850929940455+0*I)
(-3.7769077194622533-11.106812831380271*I)
(-3.949948702089913-0.9503184433887455*I)
(-3.558355491208183-0.04104686782028688*I)
(-3.577077533853747+8.297706605468053e-005*I)
(-3.5770735614781484-8.721290183674583*I)
...

lua> I.W(a, k, true, I.conj(x)) -- +0*I to -0*I fixed it!
(-2.3025850929940455-0*I)
(-3.7769077194622533-0*I)
(-3.579124176402325-0*I)
(-3.5771522746958055-0*I)
(-3.5771520639573-0*I)
(-3.577152063957297-0*I)

(04-07-2023 02:47 PM)Albert Chan Wrote:  f = x + ln(x) - ln(a) - 2*k*pi*I = 0
x.imag + phase(x) = 2*k*pi + phase(a)

Guess x with +0*I is a lousy one, with phase(x) = +pi
LHS = 0 + (+pi) = +pi
RHS = 2*(-1)*pi + (pi) = -pi

guess of conj(x) with -0*I is much better, LHS = 0 + (-pi) = -pi = RHS

The fix is simple, argument for log, instead of wrong implicit I.new(-a), use explicit -I.new(a)

lua> I.log(-I.new(a)) -- good W-1 guess
(-2.3025850929940455-0*I)
04-07-2023, 08:40 PM (This post was last modified: 04-14-2023 04:40 PM by Albert Chan.)
Post: #5
 Albert Chan Senior Member Posts: 2,708 Joined: Jul 2018
RE: lambertw, all branches
I am keeping code simple. Edge cases, example W0(0) = 0, are not implemented.

Code:
r, err = 0.36787944117144233, -1.2428753672788363E-17 require'complex' function complex.log1p(x) local y = x-(-1); return I.log(y) - ((y-1-x)/y+0) end function complex.W01(a, K, verbal, y)   -- K = abs(W branch) = 0 or 1     K = K or 0     local h, s = a+r+err, function(y) return I.log1p((y-r-err)/r) end     if not y then         y = I.sqrt(2*r*h)*(1-2*abs(K))  -- assume close to branch point         y = r + y * I.sqrt(1+y/(3*r))   -- assume W imag part tiny     end     for i = 1, 10 do         if verbal then print(a/y, 'W01') end         h = y - (y+a) / s(y)         y = y - h -- Newton's correction         if y == y+h*0x1p-13 or not I.finite(y) then return a/y end     end     return NaN end function complex.W(a, k, verbal, x)     a, k = I.new(a), k or 0     local K, h, s = abs(k), a+r+err     if I.abs(h)<.25 and (K==0 or K==1 and signbit(k*I.imag(a))) then         return complex.W01(a, K, verbal), 'W01'     end     local A, T = I.abs(a), 2*k*pi+I.arg(a)     s = function(x) return x + I.new(log(I.abs(x)/A), I.arg(x)-T) end     if not x then     -- W rough guess         if K > 1 then    x = 2*k*pi*I         elseif K==1 then x = A==1 and 1 or I.log(-a)         else --[[K==0]]  x = A==1 and 1 or I.log1p(a)         end     end     for i = 1, 10 do         if verbal then print(x) end         h = x/(x-(-1)) * s(x)         x = x - h -- Newton's correction         if x == x+h*0x1p-13 or not I.finite(x) then return x end     end     return complex.W01(a, K, verbal, a/I.real(x)), 'W01' end

W01 is solving W branch 0 or an ambiguous ±1 branch (with the small imag part)
The code is solving for y*ln(y) = a --> y = e^W(a) --> W(a) = a/y

We can "pull out" conjugate: W01(conj(a),k) = conj(W01(a,k))

lua> I.W(3+4*I, 5, true)
(0+31.41592653589793*I)
(-1.815554248884793+30.714634540472343*I)
(-1.8170058914565277+30.71333413897562*I)
(-1.8170058918466272+30.713334137004896*I)
(-1.8170058918466274+30.713334137004896*I)

lua> I.W(3+4*I, -5, true)
(-0-31.41592653589793*I)
(-1.7565805546696363-28.861921295458465*I)
(-1.7547636407605829-28.85710103800315*I)
(-1.7547636279545493-28.85710104741841*I)
(-1.7547636279545493-28.85710104741841*I)

lua> I.W(-0.3,-1,true)
(-1.7824743623491122+0*I)      W01
(-1.7813374873395937+0*I)      W01
(-1.7813370234217047+0*I)      W01
(-1.781337023421628+0*I)        W01

lua> I.W(-0.03,-1,true)
(-3.506557897319982-0*I)
(-5.261733759227882-0*I)
(-5.144793893666586-0*I)
(-5.1444827237861634-0*I)
(-5.144482721515682-0*I)
(-5.144482721515682-0*I)

lua> I.W(-.1+1e-100*I, -1, true)
(-2.3025850929940455-1e-099*I)
(-3.7769077194622533+1.2592576024726166e-099*I)
(-3.579124176402325-8.545840728943592*I)
(-3.7689826579869887-0.8337217972782689*I)
(-3.5535464061313635-0.025256342802877585*I)
(-3.5771468383315437+6.535012002211599e-005*I)
(-3.5770902410326704-8.721224618630954*I)
(-3.7811926340275566-0.8444302692850991*I)
(-3.5536175209072614-0.02647513850865879*I)
(-3.5771431199160215+6.824715919089008e-005*I)
(-3.5770875001635236+0*I)      W01
(-3.5771520641833883-1.3880349425219542e-099*I)      W01
(-3.5771520639572962-1.3880252213229777e-099*I)      W01
(-3.5771520639572976-1.388025221322979e-099*I)        W01

W code f=0 have trouble converging. Switched to W01 help a lot.
W01 use W last guess, but dropped the unreliable imaginary part.
04-08-2023, 03:21 PM (This post was last modified: 04-15-2023 01:56 PM by Albert Chan.)
Post: #6
 Albert Chan Senior Member Posts: 2,708 Joined: Jul 2018
RE: lambertw, all branches
(04-07-2023 07:54 PM)Albert Chan Wrote:  It is just here, "equivalent" may mean a difference of +0*I vs -0*I

Just to show how sensitive this ±0*I business is ...

lua> I.W(-1e-9,0,true)
(-1.0000000005000001e-009+0*I)
(-1.000000001e-009+0*I)
(-1.0000000010000002e-009+0*I)

lua> I.W(I.conj(-1e-9),0,true)
(-1.0000000005000001e-009+0*I)
(-1.000000001e-009+6.283185316604365e-009*I)
(3.1450880362038825e-008-4.729032266472273e-010*I)
(-7.848729544901978e-008-9.717506161329017e-008*I)
...

2nd should be complex conjugate of first. W(-1E-9 ± 0*I) ≈ -1.000000001E-9 ± 0*I
Imaginary part sign of (a, W0(a)) should match.

At first, I thought it is also a case of bad W guess.
To keep x imag part intact, x+1 should be written as x-(-1)
I added +0 to correction too, just to be sure. Technically, it is not needed.

Signed zero: (±0) - 0 = (±0)      ,       (±0) + 0 = 0

Code:
< function complex.log1p(x) local y = x+1; return I.log(y) - (y-1-x)/y end > function complex.log1p(x) local y = x-(-1); return I.log(y) - ((y-1-x)/y+0) end

This get W0 guess phase correct, however iterations still diverges.
It is interesting Newton's 1st iteration essentially undo log1p patches.

lua> I.W(I.conj(-1e-9),0,true)
(-1.0000000005000001e-009-0*I)
(-1.000000001e-009+0*I)
(-1.0000000010000002e-009+6.283185319745956e-009*I)
(3.145088037733907e-008-4.729032324603788e-010*I)
...

The problem is Newton correction: 1/slope = x/(x+1). We fix it the same way:

Code:
< h = x/(x+1) * s(x) > h = x/(x-(-1)) * s(x)

lua> I.W(I.conj(-1e-9),0,true)
(-1.0000000005000001e-009-0*I)
(-1.000000001e-009-0*I)
(-1.0000000010000002e-009-0*I)

I had applied above 2 patches to previous post.
04-08-2023, 05:54 PM (This post was last modified: 04-09-2023 08:31 PM by Albert Chan.)
Post: #7
 Albert Chan Senior Member Posts: 2,708 Joined: Jul 2018
RE: lambertw, all branches
(04-08-2023 03:21 PM)Albert Chan Wrote:  The problem is Newton correction: 1/slope = x/(x+1). We fix it the same way:

Code:
< h = x/(x+1) * s(x) > h = x/(x-(-1)) * s(x)

Patch is wrong. It get the correct result by accident.
If h imag part = -0, x = x - h may lose x sign information.

We could do x = x - (h + 0), but that is ugly.
I am sick of this ±0*I, and simply flip k to avoid it.

I.W(I.conj(-1e-9), 0) --> I.conj(I.W(-1e-9, -0))

We flip k to ensure x imag part always non-negative.
x = x - h will not generate x with imag part of -0.

Only update is complex.W code, the rest remains unaffected.
I am leaving 1/slope = x/(x-(-1)), but x/(x+1) will work as well.

Code:
function complex.W(a, k, verbal, x)     a, k = I.new(a), k or 0     local a0, h, s = a, a+r+err     local flip = k<0 or k==0 and signbit(I.imag(a))     flip = flip and I.conj or function(x) return x end     a, k = flip(a), abs(k) --> I.imag(x) non-negative         if I.abs(h)<.25 and (k==0 or k==1 and signbit(I.imag(a))) then         return complex.W01(a0, k, verbal), 'W01'     end     local A, T = I.abs(a), 2*k*pi+I.arg(a)     s = function(x) return x + I.new(log(I.abs(x)/A), I.arg(x)-T) end     if not x then     -- W rough guess         if k > 1 then    x = 2*k*pi*I         elseif k==1 then x = A==1 and 1 or I.log(-a)         else --[[k==0]]  x = A==1 and 1 or I.log1p(a)         end     end     for i = 1, 10 do         if verbal then print(flip(x)) end         h = x/(x-(-1)) * s(x)         x = x - h -- Newton's correction         if x == x+h*0x1p-13 or not I.finite(x) then return flip(x) end     end     return complex.W01(a0, k, verbal, a0/I.real(x)), 'W01' end
04-09-2023, 03:59 AM (This post was last modified: 04-22-2023 08:29 PM by Albert Chan.)
Post: #8
 Albert Chan Senior Member Posts: 2,708 Joined: Jul 2018
RE: lambertw, all branches
complex.W code version 2, with many improvements.
• W01 function removed, code placed inside W, for branch point cases only.
This allow for brute force confirmation that it always converge to correct branch.
Result guaranteed finite, NaN check (to get out of infinite loops) not needed.
• replaced hexfloat 0x1p-13, 0x1p-26 to more readable 1e-4, 1e-8
• W guess improved, for all branches.
• f formula rearranged, to cause catastrophic cancellation, on purpose!
Imag part (if small) are not very good, but at least it converge quickly.
We got very good estimated W real part, and a reasonable sized imag part.
• Added convergence confirmation, based from loop counts.
If we reached solution before max loops, we have convergence.
If not, result may still be good, but Wk(a) likely purely imaginary.

• if |Im(x)| < pi, we are on W branch 0, or branch ±1 with small imag part.
To be safe, I lower the limit, from pi to 3.

Within this region, we can improve result by y = e^W Newton step.
With good W real part, one correction is enough, to fix imag part.

y = (a+y) / (ln(y) + 1)
a/x = (a+a/x) / (ln(a/x) + 1)
x = x/(x+1) * (ln(a/x) + 1)

• Note that all iterations are using equivalent f = x - (ln(a) - ln(x)) - 2*k*pi*i
W01 code use accurate denominator around branch point, ln(y)+1 = log1p((y-1/e)*e)
W code use the general form of f, for all k's, but may cost of some accuracy, if |k|≤1
W final correction (with restriction, see next post) use accurate denominator around 0, ln(y)+1

Code:
r, err = 0.36787944117144233, -1.2428753672788363E-17 require'complex' function complex.log1p(x) local y = x-(-1); return I.log(y) - ((y-1-x)/y+0) end function complex.W(a, k, verbal)     a, k = I.new(a), k or 0     local flip = function(x) return x end   -- identity function     if k<0 or k==0 and signbit(I.imag(a)) then flip=I.conj; a=flip(a); k=-k end     local h, small_imag = a+r+err, (k==1 and signbit(I.imag(a)))     if I.abs(h)<.25 and (k==0 or small_imag) then         local y = I.sqrt(2*r*h) * (1-2*k)   -- close to branch point         y = r + y*I.sqrt(1+y/(3*r))         -- with small imag part         repeat             if verbal then print(flip(a/y)) end             h = y - (y+a) / I.log1p((y-r-err)/r)             y = y - h -- Newton's correction         until y == y+h*1e-4         return flip(a/y), true     end     local A, T = I.abs(a), (2*k*pi+I.arg(a))*I     local x = T -- W rough guess     if small_imag and A<.6 then x = I.log(-a) end     if k==0 then x = I.abs(a+1)<.6 and T/2 or I.log1p(a) end     k = 40      -- max loops     repeat         if verbal then print(flip(x)) end         h = x/(x+1) * (x + I.new(log(I.abs(x)/A),I.arg(x)) - T)         x, k = x-h, k-1     until k==0 or x == x+h*1e-8 or not I.finite(x)     if I.imag(x) >= 3 then return flip(x), k>0 end     if verbal then print(flip(x)) end     return flip(x/(x+1) * (I.log(a/x)+1)), k>0 end

lua> I.W(3+4*I, 5, true)
(0+32.34322175389954*I)
(-1.8166634488635778+30.71625713741827*I)
(-1.8170058961963056+30.713334138441088*I)
(-1.8170058918466274+30.713334137004896*I)      true

lua> I.W(3+4*I, -5, true)
(0-30.48863131789632*I)
(-1.754507767775116-28.860288696238875*I)
(-1.7547636338945982-28.857101048916533*I)
(-1.7547636279545493-28.85710104741841*I)      true

lua> I.W(-0.3678, 0, true)
(-0.9793607152032429+0*I)
(-0.9793607149578304+0*I)
(-0.9793607149578306+0*I)      true

lua> I.W(-.1+1e-100*I,-1,true)
(-2.3025850929940455-1e-099*I)
(-3.7769077194622533-5.084465616380436e-100*I)
(-3.579124176402325-5.180347740277737e-100*I)
(-3.5771522746958055-5.181454351722726e-100*I)
(-3.5771520639573-5.181454470168089e-100*I)
(-3.577152063957297-5.18145447016809e-100*I)
(-3.577152063957297-1.3880252213229779e-099*I)      true

lua> I.W(-.1+1e-100*I,0,true)
(-0.10536051565782631+1.1111111111111111e-100*I)
(-0.11161906744205384+1.1848854610650868e-100*I)
(-0.11183232962808998+1.1874337731206376e-100*I)
(-0.11183255915869776+1.1874365171431562e-100*I)
(-0.11183255915896297+1.1874365171463266e-100*I)
(-0.11183255915896298+1.259138243719729e-100*I)      true

lua> I.W(-.1+1e-100*I,1,true)
(0+9.42477796076938*I)
(-4.330508133087425+7.394500450320804*I)
(-4.448949088733658+7.307107936743056*I)
(-4.449098178536345+7.307060789152774*I)
(-4.44909817870089+7.3070607892176085*I)      true

lua> x = 5*I
lua> a = x * I.exp(x)
lua> x + x:log() - a:log()
(0+6.283185307179586*I)
lua> I.W(a, 1)
(3.1321400622130957e-068+5*I)      false
04-09-2023, 04:36 PM (This post was last modified: 04-10-2023 10:31 PM by Albert Chan.)
Post: #9
 Albert Chan Senior Member Posts: 2,708 Joined: Jul 2018
RE: lambertw, all branches
(04-09-2023 03:59 AM)Albert Chan Wrote:  [*] if |Im(x)| < pi, we are on W branch 0, or branch ±1 with small imag part.

(04-07-2023 02:47 PM)Albert Chan Wrote:  f = x + ln(x) - ln(a) - 2*k*pi*I = 0
x.imag + phase(x) = 2*k*pi + phase(a)

W code forced k non-negative, to guarantee x.imag ≥ 0
Since x.imag and phase(x) have same sign, we have:

k=0: (x.imag < phase(a) < pi)              // for k=0, we forced (0 ≤ phase(a) ≤ pi)
k=1: (x.imag < pi) ⇒ (LHS < pi + pi) ⇒ (phase(a) < 0), thus W branch ±1 with small imag part.

Quote:To be safe, I lower the limit, from pi to 3.

Experiments suggested |Im(x)| < pi is a hard limit. (why?)

lua> a = -100-92*I
lua> x = I.W(a, 1)
lua> x -- abs imag part < pi
(3.382739150014567+3.1375382604161515*I)
lua> x/(x+1) * (I.log(a/x) + 1) -- y = e^x Newton step
(3.382739150014567+3.1375382604161515*I)

If |Im(x)| > pi, y = e^W Newton step will not improve, but overshoot toward W branch 0.

lua> a = -100-93*I
lua> x = I.W(a, 1)
lua> x -- abs imag part > pi
(3.386390672646167+3.1426530733361844*I)
lua> x/(x+1) * (I.log(a/x) + 1) -- y = e^x Newton step
(4.064553974086096-2.1939787611722497*I)

Hard pi limit puzzle solved!

ln(a/x) ≈ ln(e^x) = Re(x) + i * smod(Im(x), 2*pi)         // signed mod to limit arg(z) within ±pi

If |Im(x)| < pi, ln(a/x) ≈ ln(e^x) = x
→ x/(x+1) * (ln(a/x) + 1) ≈ x/(x+1) * (x+1) = x            // Newton's step work as expected

If pi < |Im(x)| < 3*pi, ln(a/x) ≈ x ± 2*pi*i                       // sign opposite of Im(x)
→ x/(x+1) * (ln(a/x) + 1) ≈ x + x/(x+1)*(±2*pi*i)     // Newton's step overshoot toward W branch 0
04-10-2023, 04:44 PM
Post: #10
 Albert Chan Senior Member Posts: 2,708 Joined: Jul 2018
RE: lambertw, all branches
(04-09-2023 03:59 AM)Albert Chan Wrote:  [*] if |Im(x)| < pi, we are on W branch 0, or branch ±1 with small imag part.

Another way to look at this, if above is true, we can simplify f:

< f = x - (ln(a) - ln(x)) - 2*k*pi*i
> f = x - ln(a/x)

f general form has issues with catastrophic cancellation, simplifed f is much better.

The problem is we don't know |Im(x)| < pi, until we solve it.
Also, f simplifed form have multiple roots. A good guess is needed to get correct branch.
For f general form, with 1 root, guess good enough for convergence is all that is needed.

We use f general form for *good* guess, correct with using simplified f (if applicable)

lua> a = 1e-100-1e-10*I
lua> I.W(a, 0, true)
(1e-100-1e-010*I)
(1.000000082740371e-020-1e-010*I)
(1.000000082740371e-020-1e-010*I)        -- general f=0 solution
(1.0000000000000002e-020-1e-010*I)      true

Here, initial guess log1p(a) ≈ a is not good (I.log code need work ...)
But, it doesn't matter. Guess is good enough for general f to converge.
f simplified form, with 1 Newton step, we get W(a) ≈ a - a*a
04-10-2023, 06:47 PM
Post: #11
 Albert Chan Senior Member Posts: 2,708 Joined: Jul 2018
RE: lambertw, all branches
(04-10-2023 04:44 PM)Albert Chan Wrote:  Here, initial guess log1p(a) ≈ a is not good (I.log code need work ...)
But, it doesn't matter. Guess is good enough for general f to converge.

As expected, better complex.log make little difference.
I considered Kahan's clog(z), but Python's cmath.log code is simpler.

Code:
function complex.log(z)     local x, y = I.real(z), I.imag(z)     local h = hypot(x, y)     if 0.71 <= h and h <= 1.73 then         if abs(x) < abs(y) then x,y = y,x end         h = log1p((x+1)*(x-1) + y*y)/2     else         h = log(h)     end     return I.new(h, I.arg(z)) end

lua> a = 1e-100-1e-10*I
lua> I.W(a, 0, true)
(5.0000000000000005e-021-1e-010*I)      -- log1p(a) ≈ a - a^2/2
(1.0000000413701856e-020-1e-010*I)
(1.0000000413701856e-020-1e-010*I)
(1.0000000000000001e-020-1e-010*I)      true
04-13-2023, 03:03 PM (This post was last modified: 04-17-2023 10:16 PM by Albert Chan.)
Post: #12
 Albert Chan Senior Member Posts: 2,708 Joined: Jul 2018
RE: lambertw, all branches
Here is equivalent complex.W version 2, adapted to machine without signed zero.

Without signed zero, this trick is out: W(a, k) = conj(W(conj(a), -k)
Instead, let W(a, k) = x = re(z) + s*im(z)*I , such that s = ±1, arg(z) ≥ 0

s = sign(arg(x)) = ±1            |x| = |z|            arg(x) = s * arg(z)

lua> I.W(-1e-100, -1, true)
(-230.25850929940458-0*I)
(-235.72143712473638-0*I)
(-235.72115887568603-0*I)
(-235.72115887568532-0*I)
(-235.72115887568535-0*I)      true

Above example, s=-1, all x iterations, arg(x) = s * pi
Almost identical, for efficiency: arg(x) ≠ bad = s * -pi

In other words, iteration overshoot allowed, where arg(z) ≥ 0 does not hold.
But, not to the extreme opposite direction, i.e. arg(z) ≠ -pi

We need to add arg(x) edge case, to make W (without signed zero) work the same as I.W
Code:
        h = arg(x)         if h == bad: h = -h

I use Python mpmath, which does not have 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 y == y+h*1e-4: return a/y, True     s = 1 - 2 * (k<0 or k==0 and s)     A, T, bad = abs(a), (2*k*pi+arg(a))*j, -s*pi     x = T   # W rough guess     if small_imag and A<.6: x = log(-a)     if k == 0: x = T/2 if abs(a+1)<.6 else log1p(a)     if s*im(x) < 0: x = conj(x)     # s = sign(im(x))     k = 40  # max loops     while 1:         if verbal: print x         h = arg(x)         if h == bad: h = -h         h = x/(x+1) * (x + mpc(log(abs(x)/A),h) - T)         x, k = x-h, k-1         if k==0 or x==x+h*1e-8 or not isfinite(x): break     if s*im(x) >= 3: return x, k>0     if verbal: print x     return x/(x+1) * (log(a/x)+1), k>0

>>> W(-1e-100, -1, True)
-230.258509299405
(-235.721437124736 + 0.0j)
(-235.721158875686 + 0.0j)
(-235.721158875685 + 0.0j)
(mpc(real='-235.72115887568535', imag='0.0'), True)
04-13-2023, 04:14 PM
Post: #13
 floppy Senior Member Posts: 389 Joined: Apr 2021
RE: lambertw, all branches
Nice. Since I see the lambert as usefull for infinite iteration, do you have any idea how to use the lambert (with "jumping" into branches?) for following type of problem
a) z = x + i *y
b) f(z) = sqrt(x*y) + i* (x+y)/2
c) where finally f(f(f..(z)))) = a + i*a (a=AGM(x,y). AGM arithemic geometric mean)

HP71B 4TH/ASM/Multimod, HP41CV/X/Y & Nov64d, PILBOX, HP-IL 821.62A & 64A & 66A, Deb11 64b-PC & PI2 3 4 w/ ILPER, VIDEO80, V41 & EMU71, DM41X
04-19-2023, 01:30 AM (This post was last modified: 02-01-2024 03:33 PM by Albert Chan.)
Post: #14
 Albert Chan Senior Member Posts: 2,708 Joined: Jul 2018
RE: lambertw, all branches
(04-07-2023 02:47 PM)Albert Chan Wrote:  (*) Based from tests, solve for y*ln(y)=a is better than x*e^x=a, even with help of Halley's method.

Not when W guess is extremely close to true root.
We can use Newton's method, f = x*e^x - a, for finishing touches

Since we are now using W loops for guess, we might as well make it slightly worse.
Instead of log(abs(x)/A), make it log(abs(x)) - log(A), literally, f = x + ln(x) - lnk(a)
Code is easier to read, and cheaper to compute. (last term is a constant)

Instead of comparing *both* complex parts converged, we test for relative errors.
This solved the problem if any complex parts is bouncing around zero.
For symmetry, I used relative error test for Y loop as well.

With relative error test, it should not get stuck in infinte loops → loop limits removed

complex.W code, version 3
Code:
r, err = 0.36787944117144233, -1.2428753672788363E-17 require'complex' function complex.W(a, k, verbal)     a, k = I.new(a), k or 0     local flip = function(x) return x end   -- identity function     if k<0 or k==0 and signbit(I.imag(a)) then flip=I.conj; a=flip(a); k=-k end     local h = I(a:real() + r + err, a:imag())     local small_imag = (k==1 and signbit(I.imag(a)))     if I.abs(h)<.25 and (k==0 or small_imag) then         local y = I.sqrt(2*r*h) * (1-2*k)   -- close to branch point         y = r + y*I.sqrt(1+y/(3*r))         -- with small imag part         repeat             if verbal then print(flip(a/y)) end             h = y - (y+a) / I.log1p((y-r-err)/r)             y = y - h         until I.abs(h/y) < 1e-9         return flip(a/y)     end     local A, T = I.abs(a), I.new(0,2*k*pi+I.arg(a))     local x = T     -- W rough guess     if k==0  then x = (A+100)/(A+50); x = I.log1p(a*x) / x end     if small_imag then x = I.log(-a) + (A<.5 and 0 or T/2) end     T = T + log(A)  -- = ln_k(a)     if not finite(A) then return flip(T) end     if A==0 then return flip(k==0 and a or T-pi*I) end     if verbal and verbal~=true then x = flip(verbal) end     repeat         if verbal then print(flip(x)) end         h = x/(x+1) * (x + I.log(x) - T)         x = x - h     until not (I.abs(h/x) >= 1e-6)     if verbal then print(flip(x)) end     h = (x-a*I.exp(-x))/(x+1)   -- final touch     return flip(I.finite(h) and x-h or x) end

(04-09-2023 03:59 AM)Albert Chan Wrote:  lua> x = 5*I
lua> a = x * I.exp(x)
lua> x + x:log() - a:log()
(0+6.283185307179586*I)
lua> I.W(a, 1)
(3.1321400622130957e-068+5*I)      false

complex.W version 2 was not able to converge after 40 loops.
Redo with version 3:

lua> a
(4.794621373315692+1.4183109273161312*I)
lua> I.W(a, 1, true)
(0+6.570796326794897*I)
(-0.033367069043767406+4.994921913968373*I)
(1.97521636071743e-005+5.000010560789135*I)
(3.753082249819371e-012+5.000000000009095*I)
(1.0141772755203484e-016+5*I)
(-6.474631941441244e-017+5*I)

We need more precisions to get real part right.
For reference, with hexfloat of a: W1(a) ≈ 4.51974e-018 + 5.00000 I

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

lua> I.W(inf), I.W(-inf)
(inf+0*I)      (inf+3.141592653589793*I)

lua> z = 0
lua> I.W(-z,-1), I.W(z,-1)
(-inf-0*I)      (-inf-3.141592653589793*I)

Update 1/29/24

< h = (a + r + err)
> h = I(a:real() + r + err, a:imag())      -- keep -0*I sign

Update 2/01/24

lyuka branch relative limit reduced, from 1e-12 to 1e-9
Although convergence near -1/e is hard, lyuka (improved) formula is very accurate there.
Further away, guess formula is not as good, but convergence O(rel_err^2) constant < 1

I also removed lua implementation of log1p(x). it had moved to C code
04-23-2023, 02:49 PM (This post was last modified: 04-25-2023 10:32 AM by Albert Chan.)
Post: #15
 Albert Chan Senior Member Posts: 2,708 Joined: Jul 2018
RE: lambertw, all branches
Previous post, complex.W version 3 old guess code (OK, but now replaced)
Note: k flipped to non-negative, to force Im(x) non-negative too, to avoid -0*I issue.
Code:
    local A, T = I.abs(a), (2*k*pi+I.arg(a))*I     local x = T     -- W rough guess     if k==0 then x = I.abs(a+1)<.6 and T/2 or I.log1p(a) end     if small_imag and A<.6 then x = I.log(-a) end

k=0 and Im(a)≥+0, W0(a) guess of log1p(a) is good, except around log1p singularity, a ≈ -1
If a is real, inside (-1,-1/e), iterations may get stuck on real line, or take a while to get out.
That's why I had added a "hole" around -1, with purely complex guess.

x + ln(x) = ln(|a|) + T
Im(x + ln(x)) = T/i = [pi/2, pi]

a ≈ -1      → guess x = T/2 = [pi/4, pi/2]*i

We could hide log1p singularity inside branch region, for simpler and better W0 guess
Code:
< if k==0 then x = I.abs(a+1)<.6 and T/2 or I.log1p(a) end > if k==0 then x = I.log1p(2*a)/2 end

New guess is better for small a.

W0(ε)          = ε - ε² + 3/2*ε³ - ...
log(1+2ε)/2 = ε - ε² + 4/3*ε³ - ...

New guess is worse for huge a. (factor of 2) But, it doesn't matter.
x ≫ ln(x), curve is very flat. Cost is at most 1 extra iteration.

Update: perhaps we can adjust constant 2, to have best of both worlds.
Constant is 2 when a is small, 1 when a is huge (getting back log1p(a))

Code:
< if k==0 then x = I.log1p(2*a)/2 end > if k==0 then x = (A+100)/(A+50); x = I.log1p(a*x) / x end

k=1 and Im(a)≤-0, W1(a) guess of log(-a) is good, except a ≈ -1 (log(1)=0, terrible W guess)
If a is real < -1/e, iterations may get stuck on real line, or take a while to get out.

My solution was only use log(-a) for a ≈ 0, and purely imaginary T as default W guess.

T/i = 2*pi + arg(a)

re(a)≤0:      T/i = 2*pi .+ (-pi, -pi/2] = (pi, 3/2*pi]
re(a)≥0:      T/i = 2*pi .+ [-pi/2, -0.] = [3/2*pi, 2*pi]

Another approach is to add reasonable complex part to log(-a).
if Re(h = a+1/e) < 0, W is complex.
Code:
< if small_imag and A<.6 then x = I.log(-a) end > if small_imag then x = I.log(-a) + (I.real(h)<0 and T/2 or 0) end

New guess, Re(a) < -1/e, which covered a ≈ -1:

arg(a) = (-pi, -pi/2]
im(log(-a) + T/2) = (arg(a)+pi) + (pi+arg(a)/2) = (0, pi/2] + (pi/2, 3/4*pi] = (pi/2, 5/4*pi]

I've patched new W guess to complex.W, version 3. Better W guess welcome ...

Update: it seems to fit better if we again treat a ≈ 0 as the hole.
Code:
< if small_imag then x = I.log(-a) + (I.real(h)<0 and T/2 or 0) end > if small_imag then x = I.log(-a) + (A<.5 and 0 or T/2) end

Outside the hole, with Re(a) > 0, we have same formula, but different domain.

arg(a) = [-pi/2, -0]
im(log(-a) + T/2) = (arg(a)+pi) + (pi+arg(a)/2) = [pi/2, pi] + [3/4*pi, pi] = [5/4*pi, 2*pi]

Added T/2 match asymptotes imaginary part better.

lua> I.W(1e300-I, 1, true)
(690.7755278982137+6.283185307179586*I) -- 2*pi
(684.2471217875482+6.274017030794511*I)
(684.2471666555265+6.274016338969102*I)
(684.2471666555265+6.274016338969102*I) -- ≈ 2*pi
lua> I.W(1-1e300*I, 1, true)
(690.7755278982137+3.9269908169872414*I) -- 5/4*pi
(684.2471409153251+4.705523410032926*I)
(684.2471850188219+4.705512170223622*I)
(684.2471850188219+4.705512170223621*I) -- ≈ 3/2*pi
04-23-2023, 04:40 PM (This post was last modified: 02-01-2024 03:35 PM by Albert Chan.)
Post: #16
 Albert Chan Senior Member Posts: 2,708 Joined: Jul 2018
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)
01-19-2024, 04:14 PM
Post: #17
 Albert Chan Senior Member Posts: 2,708 Joined: Jul 2018
RE: lambertw, all branches
I have added w = Wk(a), for edge cases, |a| = ∞ or 0 (updated to post 14 and 16)

For |a| → ∞, we have Wk(a) = lnk(a)

For |a| → 0, k=0

XCas> series(LambertW(x))      → x - x^2 + 3/2*x^3 - 8/3*x^4 + 125/24*x^5

For |a| → 0, k≠0, we use W branch k definition

w + ln(w) = lnk(a) = ln(|a|) + (arg(a) + 2*k*pi)*I

Let w = r + θ*I, matching complex parts

--> (r + re(ln(r)) ≈ r) = ln(|a|) = -∞
--> (θ + sign(θ)*pi) = (arg(a) + 2*k*pi)

--> sign(θ + sign(θ)*pi) = sign(arg(a) + 2*k*pi)
--> Since k≠0,   sign(θ) = sign(k) = ±1

Wk(a → 0) = a if k==0 else lnk(a) - sign(k)*pi*I

This depends on meaning of arg(a) = atan2(im(a), re(a)) = atan2(±0, ±0)
IEEE starndard assumed real part is (ε → 0), but imaginery part is truly 0.

IEEE: arg(±0 ± 0*I) → arg(±ε ± 0*I) → ±0 or ±pi

Technically, we don't really know arg(0), except its full range (-pi, pi]
WA avoid all these, ∞ + (finite number, real or complex) = ∞
WA also don't have signed zero, ±0 ±0*I collapsed to plain 0

WA: Wk(0) = 0 if k==0 else -∞
01-20-2024, 04:48 PM
Post: #18
 Albert Chan Senior Member Posts: 2,708 Joined: Jul 2018
RE: lambertw, all branches
Newton's method convergence to prove previous post.

x * e^x = a               // x = Wk(a)
x + ln(x) = lnk(a) = c

If c = ±∞ + θ*I, then (|x| → ∞) ⇒ (slope = 1+1/x ≈ 1)

Iteration formula is practically equivalent to Newton's method.

x = c - ln(x)

x1 = c                      // guess
x2 = c - ln(x1)
x3 = c - ln(x2)          // ... until x's converged

if a = ∞*cis(θ), r = re(c) = ∞

c = r + θ*I ≈ |r| * cis(θ/r)
ln(c) ≈ ln(|r|) + (θ/r)*I
x2 ≈ c - ln(c) ≈ (r-ln(|r|)) + (θ-(θ/r))*I ≈ (r + θ*I) = x1

--> Wk(∞*cis(θ)) = lnk(∞*cis(θ))

if a = 0*cis(θ) and k≠0, r = -∞, sign(θ) = sign(k)

c = r + θ*I ≈ |r| * cis(sign(k)*(pi - |θ/r|))
ln(c)   ≈ ln(|r|) + (sign(k)*(pi - |θ/r|))*I  ≈ ln(|r|) + sign(k)*pi*I
x2 = c - ln(c) ≈ (r-ln(|r|)) + (θ-sign(k)*pi)*I ≈ r + θ2*I

sign(θ2) = sign(θ) = sign(k)      // 1-sided convergence (★)
ln(x2) ≈ ln(|r|) + (sign(k)*(pi - |θ2/r|))*I ≈ ln(|r|) + sign(k)*pi*I

x3 = c - ln(x2) ≈ c - ln(c) = x2

--> Wk≠0(0*cis(θ)) = lnk(0*cis(θ)) - sign(k)*pi*I

(★) Without slope correction, convergence is 1-sided, not crossing discontinuity.

lua> c = I(-99999*log(10), -pi)      -- = ln-1(-10^-99999)
lua> c
(-230256.2067143116-3.141592653589793*I)
lua> c - I.log(_)
(-230268.55366222185-1.3643899976489848e-05*I)
lua> c - I.log(_)
(-230268.55371584298-5.925215873503475e-11*I)
lua> c - I.log(_)
(-230268.5537158432+0*I)          -- W-1(-10^-99999) conjugate

Note that final imaginery part should be -0
This is why I.W(a, k) flip to positive k, iterate, then flip back for negative k

lua> c = c:conj()
lua> c
(-230256.2067143116+3.141592653589793*I)
lua> c - I.log(_)
(-230268.55366222185+1.3643899976489848e-05*I)
lua> c - I.log(_)
(-230268.55371584298+5.925215873503475e-11*I)
lua> c - I.log(_)
(-230268.5537158432+0*I)
lua> _:conj()
(-230268.5537158432-0*I)          -- W-1(-10^-99999)
01-20-2024, 10:52 PM (This post was last modified: 01-20-2024 11:49 PM by Gil.)
Post: #19
 Gil Senior Member Posts: 637 Joined: Oct 2019
RE: lambertw, all branches
I looked your case k=-1 & k=1, x=-.1 + i ×1E-99.

With my initial guess=
'LN(((-.1,1.E-100)+1./e)*(e-\v/2.-1.)+\v/(((-.1,1.E-100)+1./e)*2./e)+1./e)'

= (-.11284787341,1.26802807724E-100),

the built-in multisolver from HP50 runs endless.

I looked in your programs, but could not understand which initial Xo formulae I should use for tiny imaginary part.

I tried, from what I guessed from your code, as a guess xo=ln(-a), and it worked.

When should I used that guess formula for my initial guess?

For example, xo=ln(-a) failed to produce any answer for

W0((1.E-100,-.0000000001)): (9.99905166733E-21,-1.00000000005E-10).

And what do you mean by I.log(-a)?
Is I.log(-x) a special function or is it i×ln(-x).

Could you help me on that matter?

01-21-2024, 01:14 AM (This post was last modified: 01-24-2024 07:27 PM by Albert Chan.)
Post: #20
 Albert Chan Senior Member Posts: 2,708 Joined: Jul 2018
RE: lambertw, all branches
(01-20-2024 10:52 PM)Gil Wrote:  I looked your case k=-1 & k=1, x=-.1 + i ×1E-99.

With my initial guess=
'LN(((-.1,1.E-100)+1./e)*(e-\v/2.-1.)+\v/(((-.1,1.E-100)+1./e)*2./e)+1./e)'

= (-.11284787341,1.26802807724E-100),

the built-in multisolver from HP50 runs endless.

x = Wk(a)

Newton's method required a good guess to start. (iterations cannot cross discontinuity!)

lyuka's e^W formula is only used when a ≈ -1/e, and (k=0 or k=±1 and x with small imag part)
Above example probably too far away from -1/e (my code use lyuka if |a+1/e| < 0.25)

Also, your formula is only good for k=0. For k=-1, you need to flip sqrt sign.

lua> a, e = I(-0.1, 1e-100), exp(1)
lua> I.log((a+1/e)*(e-sqrt(2)-1) + I.sqrt((2/e)*(a+1/e)) + 1/e) -- x guess for k=0
(-0.1128478734106899+1.268028077236754e-99*I)
lua> I.log((a+1/e)*(e-sqrt(2)-1) I.sqrt((2/e)*(a+1/e)) + 1/e) -- x guess for k=-1
(-5.225138594746179-9.75118016335884e-98*I)

If x imag part is big, it just use that.

-- local A, T = I.abs(a), I.new(0,2*k*pi+I.arg(a))
-- local x = T -- W rough guess

lua> I.W(a, 1, true)
(0+9.42477796076938*I)
(-4.330508133087424+7.394500450320804*I)
(-4.448949088733658+7.307107936743056*I)
(-4.449098178536345+7.307060789152774*I)
(-4.44909817870089+7.3070607892176085*I)
(-4.44909817870089+7.3070607892176085*I)

If k=-1 and x imag part is small, we wanted real guess when a is negative real.
Also, x→-∞ if a→0, log(-a) fit the bill.

But if a is too negative (left of -1/e), code added some imaginery part.
Imag part is needed because a=-1 --> log(-a) = log(1) = 0, very bad x guess!

-- if small_imag then x = I.log(-a) + (A<.5 and 0 or T/2) end

lua> I.W(a, -1, true)
(-2.3025850929940455-1e-99*I)
(-3.776907719462253-5.084465616380437e-100*I)
(-3.5791241764023245-5.180347740277738e-100*I)
(-3.577152274695805-5.181454351722727e-100*I)
(-3.5771520639572993-5.18145447016809e-100*I)
(-3.577152063957297-1.388025221322979e-99*I)
 « Next Oldest | Next Newest »

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