lambertw, all branches

[Albert Chan]

(01-20-2024 10:52 PM)Gil Wrote:  For example, xo=ln(-a) failed to produce any answer for
W0
(
(1.E-100,-.0000000001)
)
: a bad guess.

I had to use instead the first formulae to get an answer for
W0
(
(1.E-100,-.0000000001)
)
: (9.99905166733E-21,-1.00000000005E-10).

W guess was explained in post 15.

lua> I.W(I(1e-100,-1e-10), 0, true)
(9.999999999990001e-21-1e-10*I)
(9.999999999990001e-21-1e-10*I)
(1.0000000000000001e-20-1e-10*I)

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

Quote:And what do you mean by I.log(-a)?

complex functions are stored in lua table complex
I = complex.I, and shared the same table.

lua> I.log, complex.log
function: 006acc20      function: 006acc20

[Gil]

A question, but sorry I do not understand how to interpret your code:

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

I take x=abs (a+1),
A)"if that result is less than 0.6"
B) A)"if that result is less than 0.6 AND that result is also less than a×EXP(a)
C) About T/2: it is an interval from [pi/4 pi/2]×pi, yes?
My question, how can I say that abs(a) is less than T/2 = an imaginary number with no real part?
D) "if the result abs(a) is not less than 0.6, then xo=I.log1p(2*a)/2

Sorry for these questions, but I am not a math guy, as you know, that however wishes to partially understand the criteria used in order to have his program work in any situation.

Thanks again.

[Albert Chan]

Hi, Gil

This is more lua language question than math question.
see http://lua-users.org/wiki/ExpressionsTutorial

x = A and B or C      → x = ((A and B) or C)      // parenthesis added to show precedence order

Code:

if A then
    local T = B
    if T then x=T else x=C end
else
    x = C  
end

Test expression, false or nil are false, everything else (e.g. 0.0, "", {}) are true.

Quote:< if k==0 then x = I.abs(a+1)<.6 and T/2 or I.log1p(a) end

Code:

if k==0 then
    if I.abs(a+1)<.6 then
        x = T/2
    else
        x = I.log1p(a)
    end
end

(01-21-2024 01:53 PM)Gil Wrote:  C) About T/2: it is an interval from [pi/4 pi/2]×pi, yes?

For k=0, I.W code ensured imaginery part is positive (x>0, y≥+0):

W₀(-x-y*I) = conj(W_-0(conj(-x-y*I))) = conj(W₀(-x+y*I))

T = (arg(-x+y*I) + 2*0*pi*I) = arg(-x+y*I)*I = (pi - atan(y/x))*I = (pi/2 .. pi)*I
T/2 = (pi/4 .. pi/2)*I

Quote:> if k==0 then x = I.log1p(2*a)/2 end

This is my 2nd iteration of W0 guess:
log1p(-1) = -∞ occurs when a = -0.5, inside lyuka branch |a+1/e|<.25, thus are safe.

if |a|→0, then x = log1p(2a)/2 ≈ a-a*a, matching W0(a) taylor series
if |a|→∞, then x = log1p(2a)/2 ≈ ln(2a)/2 ≈ ln(a)/2

Not quite match asymptote, but it will, for next iteration

f = x + ln(x) - ln(a) = ln(a)/2 - ln(ln(a)/2) - ln(a) ≈ -ln(a)/2
f' = 1 + 1/x ≈ 1

x - f/f' ≈ ln(a)/2 + ln(a)/2 = ln(a)

[Gil]

Thanks

And log1p(a)= (a*1) × EXP (a*1)?

Or rather, simply : ln(1*a).

My question is relative to the guess of xo,
when z=5+i×6, abs (z) >= 0.6

A simple guess for xo might be xo= ln(5+6i),
and not z×EXP(z)?

[Albert Chan]

(01-21-2024 04:35 PM)Gil Wrote:  And log1p(a)= (a*1) × EXP (a*1)?

I don't follow where EXP comes from
log1p(a) = ln(1+a), but evaluated more accurately when re(a) is tiny.

Quote:My question is relative to the guess of xo,
when z=5+i×6, abs (z) >= 0.6

A simple guess for xo might be xo= ln(5+6i),
and not z×EXP(z)?

Sure. W₀(z, |z|→∞) = ln(z)

[Albert Chan]

(01-21-2024 06:03 PM)Albert Chan Wrote:  Sure. W₀(z, |z|→∞) = ln(z)

For completeness, W_-1(z, |z|→0, arg(z)≥0) = ln(-z)
Note: z is complex number. -z flipped *both* complex parts sign.

ln(-z) = ln(z) - pi*I = (ln(z) - 2*pi*I) + pi*I = ln_-1(z) - sign(-1)*pi*I

This matched W_-1(0) asymptote.

(01-20-2024 04:48 PM)Albert Chan Wrote:  --> W_k≠0(0*cis(θ)) = ln_k(0*cis(θ)) - sign(k)*pi*I

[Gil]

Sorry for my ignorance.

I discover now the utility of this function that enables to get a result different of zero for cases ln(1+tiny), when normally ln(1+tiny) =ln(1) = zero.

The same useful function exists also in the HP50G.

[Gil]

One more question, if I may:

You wrote
"For k=0, I.W code ensured imaginery part is positive (x>0, y≥+0):

W0(-x-y*I) = conj(W-0(conj(-x-y*I))) = conj(W0(-x+y*I))

T = (arg(-x+y*I) + 2*0*pi*I) = arg(-x+y*I)*I = (pi - atan(y/x))*I = (pi/2 .. pi)*I
T/2 = (pi/4 .. pi/2)*I"

a) Suppose that k=0, s=(+1.E-100-i.0000000001).
—> abs(a) < 0.6
Then initial xo = T/2.
Could you give the different step to find the value of T/2?
I found, for xo =T/2, pi/4× PI, but that value leads to endless loops on my program, whereas xo =ln(1+a) gives the expected result.

But here I was lucky.
As I was not sure of the value of T/2,
b) for the other case k=0, a=(-.1,+i×1.E-100),
though abs(a) <0.6,
I tried as the above mentioned case,
xo=ln(1+a),
but then I got endless loops.
What would be the correct step to find here also T/2.

Meanwhile,Takayuki HOSODA's formulae for initial xo gives good, approximate results first digits exact) in both cases a) & b)

Thanks again for your comprehension and possible help.

[Gil]

In one of your posts, you wrote:
">>> 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)"

Shouldn't the correct answer rather be 0+i×PI/2?

Taking the first digits of PI, ie not the 'exact full form',
W0(-1.57079632679489661923132169163975144209858469968755291048) -2.78783334301039861107336040858316822982734956396028 × 10^-63 +1.5707963267948966192313216916397514420985846996875529104 i

With the exact form, with "full possible digits", real part —> 0.

[Albert Chan]

(01-21-2024 08:39 PM)Gil Wrote:  Suppose that k=0, s=(+1.E-100-i.0000000001).
—> abs(a) < 0.6

You misread W0 guess formula:

< if k==0 then x = I.abs(a+1)<.6 and T/2 or I.log1p(a) end

BTW, why did you pick this? The other version so much simpler!

> if k==0 then x = I.log1p(2*a)/2 end"

Quote:Then initial xo = T/2.
Could you give the different step to find the value of T/2?
I found, for xo =T/2, pi/4× PI, but that value leads to endless loops on my program, whereas xo =ln(1+a) gives the expected result.

For k=0, I.W code assumed a imag part positive (conj to positive if not true)

Python code does not do flips, perhaps you should use that as guide.
Python version (no flip arguments, no signed zero):

T = mpc(0,2*k*pi+arg(a)) = mpc(0, arg(a)) ≈ -pi/2*I
T/2 ≈ -pi/4*I

Of course, this is not applicable for x guess.

Quote:for the other case k=0, a=(-.1,+i×1.E-100),
though abs(a) <0.6,
I tried as the above mentioned case,
xo=ln(1+a),
but then I got endless loops.
What would be the correct step to find here also T/2.

Same issue, abs(a+1), not abs(a)

T = mpc(0,2*k*pi+arg(a)) = mpc(0, arg(a)) ≈ pi*I
T/2 ≈ pi/2*I

Again, this is not applicable for x guess.

To experiment, I patched I.W to allow user supplied guess. I don't see problem with log1p(a)

lua> a = I(-.1, 1e-100)

lua> I.W(a, 0, I.log1p(a))
(-0.1053605156578263+1.1111111111111111e-100*I)
(-0.11161906744205385+1.1848854610650873e-100*I)
(-0.11183232962808998+1.187433773120638e-100*I)
(-0.11183255915869776+1.1874365171431567e-100*I)
(-0.11183255915896298+1.1874365171463274e-100*I)
(-0.11183255915896298+1.2591382437197292e-100*I)

As expected, I.log1p(2a)/2 ≈ a-a*a is better, matching W0(a) taylor series.

lua> I.W(a, 0, I.log1p(2*a)/2)
(-0.11157177565710488+1.25e-100*I)
(-0.1118322166456466+1.2532842987580516e-100*I)
(-0.11183255915837244+1.2532886205550325e-100*I)
(-0.11183255915896304+1.2532886205624847e-100*I)
(-0.11183255915896297+1.259138243719729e-100*I)

[Gil]

For W0,
1) I can use now only xo=ln1p(2a)/2 for real numbers and the corresponding "ln(1+2a)/2" for complex numbers?
Right?

2) And, always for W0, completely ignore Takayuki HOSODA 's formula for initial guess xo?

3) And only use Takayuki HOSODA 's formula for k=1, ±2, ±3,± 4...
and for "k=-1 and (x<-1/e or x>0)".
Right?

Thanks again for confirming.

[Gil]

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


You write further
Wk(a → 0) = a if k==0 else lnk(a) - sign(k)*pi*I
Shouldn't it be rather

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

[Albert Chan]

Hi, Gil

If you want to implement W to HP50g, ignore lua I.W code, and use only Python code (post 16)
My guess is HP50g does not support signed zero anyway.

Python code, including blank lines, is only 34 lines long.
It is all in the code, but this is the gist of it ...

If |k| > 1, just use imaginery part as guess. There is no Newton's over-shoot issue.
If |k| = 1, and k and im(a) has same sign, again, use T as guess

A, T = abs(a), mpc(0,2*k*pi+arg(a))
x = T            # W rough guess

Now the hard part, we need good guess to avoid Newton over-shoot crossed discontinuity.

If |a+1/e| < 0.25, use lyuka's formula and y = e^W iteration formula, return a/y = x
This take cared of singularity at x = -1/e, W₀(-1/e) = W_-1(-1/e) = -1

If |k| = 1, guess x = LN(-a) + (0 if A<.5 else T/2)

We need customized LN() because we lost signed zero information.
If a ≈ -1, LN(-a) is very tiny, possibly 0, thus the need for adding T/2

If k = 0, guess x = log1p(a*s) / s, where s = (A+100)/(A+50) = (2 if A=0, 1 if A=∞)
We have good W0 asymptote estimate when A=0, and A=∞

repeat until loop does Newton's iteration. x = x - f/f'
f is designed to cause catastrophic cancellation, making x possibly not accurate.

Final touch with g = x*e^x - a, g' = (x+1)*e^x
--> h = g/g' = (x - a*e^-x) / (x+1)

Because of the exponential, h may not be finite.
If h is finite, we update for better x; if not, last x is the best we've got.

Quote:You write
Wk(+0) = 0 if k==0 else -∞ + (2*k-sign(k))*pi*I
OK

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

Wk(+0) is referring Python code, which does not support signed zero.

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

ln_k(0) = ln(|0|) + (arg(0) + 2*k*pi)*I = -∞ + 2*k*pi*I
ln_k(0) - sign(k)*pi*I = -∞ + (2*k-sign(k))*pi*I

With system that support signed zero, 0 may mean (±0 + ±0*I)
This is why ln_k(a) is left unevaluated. arg(0) may have multiple values.

lua> z = 0
lua> I(z,z):arg(), I(z,-z):arg()
0                                  -0
lua> I(-z,z):arg(), I(-z,-z):arg()
3.141592653589793      -3.141592653589793

[Gil]

You write
"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):“

From the above, I understand that
small_imag=1 if (k=1 and im(a) <0 or k=-1 and im(a) >≠0)
Else small_imag =0.

Then
For k=1, a=-1 —> im(a) = 0, ie im(a) is not < 0.
Consequently, small_imag=0.

But then you write
if small_imag: x = LN(-a) + (0 if A<.5 else T/2)
But a=-1 small_imag = 0
Then my question: for that case, what is the value for initial xo?

[Albert Chan]

(01-22-2024 06:09 PM)Gil Wrote:  Then
For k=1, a=-1 —> im(a) = 0, ie im(a) is not < 0.
Consequently, small_imag=0.

But then you write
if small_imag: x = LN(-a) + (0 if A<.5 else T/2)
But a=-1 small_imag = 0
Then my question: for that case, what is the value for initial xo?

We are left with easy case, x = T = mpc(0,2*k*pi+arg(a))

p2> W(-1, 1, True)
(0.0 + 9.42477796076938j)
(-2.05355675220527 + 7.63609249528568j)
(-2.06229155975165 + 7.58864445610365j)
(-2.06227772959937 + 7.58863117846966j)
(-2.06227772959828 + 7.58863117847251j)
(-2.06227772959828 + 7.58863117847251j)

[Gil]

1) You write
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.

What do you mean by more precision?
Can we get better result, at least regarding the sign of the real part?[/b]

2)You give the examples for
"W0(-pi /2), that should be (0+i×pi/2)
"W-1(-pi /2), that should be (0 -i×pi/2).

Is there a way to get the real part almost = to zero (±1E-100)?

The answer is have an initial guess xo with real part = 0 —> not to have ln(-a)+t/2, but only t/2, as ln(-a) = ln(pi/2) = real.

Could we say, instead of that code line
xo= ln(-a) + (0 of A<. 5, else t/2)

IF a>=0
THEN xo= ln(-a) + (0 of A<. 5, else t/2)
ELSE xo= always t/2
END

Try W(k=-1)
For a=-2, and put zero for the real part of the "found xo initial guess. The iterating process finds all the same the good solution.
For any number a <0, for instance -0.4, and then procede as above, taking xo = t/2,and the expected answer seems to be always found!
Of course, it's not a mathematical solution, but a possible "trick answer" for dummies. To be thought over?

[Albert Chan]

(01-22-2024 11:33 PM)Gil Wrote:  1) You write
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.

What do you mean by more precision?

Final x is basically the best we can do with IEEE double.
More Newton steps are not able to get any better, not even sign of real part.

p2> from mpmath import *
p2> a = 4.794621373315692+1.4183109273161312j
p2> x = x0 = -6.474631941441244e-017+5j
p2> for k in range(4): x -= (x-a*exp(-x))/(x+1); print complex(x)
...
(1.15986489997e-16+5j)
(-4.99631494016e-17+5j)
(1.30201076547e-16+5j)
(-3.62952777182e-17+5j)

Now, with more bits

p2> mp.prec = 113 # binary128
p2> x = x0
p2> for k in range(4): x -= (x-a*exp(-x))/(x+1); print complex(x)
...
(4.51973914236e-18+5j)
(4.51973914236e-18+5j)
(4.51973914236e-18+5j)
(4.51973914236e-18+5j)

Quote:2)You give the examples for
"W0(-pi /2), that should be (0+ipi/2)
"W-1(-pi /2), that should be (0 -pi/2).

Is there a way to get the real part almost = to zero (±1E-100)?

You can throw in more precisions, or confirm symbolically.

x = ±pi/2*I
x * e^x = ±pi/2*I * ±I = -pi/2 = a

k = im(x + ln(x) - ln(a)) / (2*pi)

if x = +pi/2*I, then k = (pi/2 + pi/2 - pi) / (2*pi) = 0
if x = −pi/2*I, then k = (-pi/2 - pi/2 - pi) / (2*pi) = -1

--> W₀(-pi/2) = +pi/2*I
--> W_-1(-pi/2) = -pi/2*I

Actually, we only need to prove one, the other will follow.

(04-23-2023 04:40 PM)Albert Chan Wrote:  Trivia: W(a < -1/e) is complex, W₀(a) = conj(W_-1(a))

x + ln(x) = ln_k(a)

ln₀(a) = ln(|a|) - pi*j + 2*pi*j = ln₁(conj(a))

W₀(a) = W₁(conj(a)) = conj(W_-1(a))

(01-22-2024 11:33 PM)Gil Wrote:  The answer is have an initial guess xo with real part = 0 —> not to have ln(-a)+t/2, but only t/2

Symbolically yes, but numerically no.

Final touch, x -= (x-a*exp(-x)) / (x+1)

pi is just an approximation of true pi = 3.14159...
exp(±pi/2*I) will not get back exactly ±I, but with tiny real part.
Numerically, we can't get cos(x) to 0, see smallest |cos(x)| ?

(01-22-2024 11:33 PM)Gil Wrote:  Try W(k=-1)
For a=-2, and put zero for the real part of the "found xo initial guess.
The iterating process finds all the same the good solution.

If final x has decent amount of complex parts, good guess is not as critical.
You may play around for better guess. Just make sure it *always* work!

[Gil]

For W (k=-1, a=+3*PI/2),
the result for the real part is tiny, but not 0 — zero as the exact result is according to Wolfram.

Could you check, Albert, what you find with your program(s)?

PS: as you mentioned, with roundings, cos 3pi/2 ≠ 0, though sin 3*pi/2 = exactly 1.

Thanks in advance for your answer.

Regards

[Albert Chan]

(01-23-2024 02:35 PM)Gil Wrote:  For W (k=-1, a=+3*PI/2),
the result for the real part is tiny, but not 0 — zero as the exact result is according to Wolfram.

Could you check, Albert, what you find with your program(s)?

p2> a = 3*pi/2
p2> a
mpf('4.7123889803846897')
p2> W(a, -1)
mpc(real='-1.4351953215276557e-16', imag='-4.7123889803846897')

Note that a is not *exactly* \( 3\pi/2 \), but only an approximation.
Even with more precisions, W(a, -1) will have tiny real part.

We accept this numerical limits, and confirm result instead.

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

x * e^x = -a*I * I = a

[Gil]

We accept...
You decided to accept that the real part in question is indeed zero.
But according to your other examples, for much more tinier results in the real part output you did not round — and did correctly — that real part to zero.
The rules should be set up before the "game" starts.
I understand that tiny is tiny, but sometimes we round it to zero, sometimes not, after checking the results. In a way, it's not very fair.

But I understand that, in that case, we cannot do much better.
The exception being for the previous case W0/W-1 (-pi/2), where we can "cheat" in advance in the program — if we want to, of course — fixing up the issue for a<=0, and get the the correct answer, as I developed in a previous post.