lambertw, all branches

01232024, 06:17 PM
Post: #41




RE: lambertw, all branches
(01232024 04:57 PM)Gil Wrote: We accept... 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) (04092023 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. (04192023 01:30 AM)Albert Chan Wrote: We can use Newton's method, f = x*e^x  a, for finishing touches 

01232024, 06:44 PM
Post: #42




RE: lambertw, all branches
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) 

01232024, 11:00 PM
(This post was last modified: 01242024 03:23 PM by Gil.)
Post: #43




RE: lambertw, all branches
Again, W(k=1,a=+PI/2 )
seems impossible to give — at least with my HP50G program —the nice answer (0i×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: :W1000('i*3998*(pi/2)'=(0.,6280.04371455)): (3.6616915941E12,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×42)) = 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) < 1E10} THEN result = i × imaginary (result) END And more generally With W0(pi/2)=0+i×Pi/2 With W1(pi/2)=0i×Pi/2 IF {imaginary (result) = integer≠0 × PI/2 and abs (real (result) < 1E10} 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 (M0M1)<0.0001 and abs (M2) < 1E11} 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 =(4k1)×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 (M0M1)<0.0001 and abs (M2) < 1e11 are fulfilled. 

01242024, 03:18 PM
(This post was last modified: 01242024 03:36 PM by Gil.)
Post: #44




RE: lambertw, all branches
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 

01242024, 08:53 PM
(This post was last modified: 01252024 12:27 PM by Albert Chan.)
Post: #45




RE: lambertw, all branches
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_{1}(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_{1}(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_{1}(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_{1}(3/2*pi) = 3/2*pi*I too, because it matched last formula form. odd = 2n  sign(n) // n ≠ 0 x = odd*pi*I = (2nsign(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_{1}(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 = (2nsign(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 = ((2n1/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_{1}(3/2*pi) = 3/2*pi*I 

01252024, 12:37 AM
Post: #46




RE: lambertw, all branches
I try with Wolfram
W0 (1.2E323) and get the correct answer, ie about 1.2E323. 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.2E324), I don't get now the expected approximate value 1.2E324, but directly zero. Why this distinction at 1E323/E324? 

01252024, 01:10 AM
(This post was last modified: 01252024 01:21 AM by Gil.)
Post: #47




RE: lambertw, all branches
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 nonnegative. But W1(1E330)= 697.32277629546016099540752740546566360568199201428307824605649289390557552642032068023149199134240178698554915132 And W1(1E330)= 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. 

01252024, 03:04 AM
Post: #48




RE: lambertw, all branches
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). 

01252024, 07:02 AM
Post: #49




RE: lambertw, all branches
Hi, Gil
Complex Conjugate formula is always true, if we have signed zero. With signedzero, 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))) 

01252024, 10:09 AM
(This post was last modified: 01252024 10:13 AM by Gil.)
Post: #50




RE: lambertw, all branches
Thanks, Albert.
And how do you explain these different results? W1(1E330)= 697.32277629546016099540752740546566360568199201428307824605649289390557552642032068023149199134240178698554915132 And W1(1E330)= 697.32281706295494653348757314785952050772713371449039168 +6.2922084418351405256893689216267218827742283825238692107 i Are they correct? 

01252024, 04:13 PM
Post: #51




RE: lambertw, all branches
Quote:With signedzero, arg(z) + arg(conj(z)) = 0 Without signedzero, 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("1E330", "1E9999") p2> W(z, 1) (766.494908744112  1.00130634441663e9669j) p2> conj(W(conj(z), 1)) (766.494908744112  1.00130634441663e9669j) 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, Still, I think better solution is to properly add signedzero to mpmath. fredrikjohansson 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 ... 

01252024, 05:14 PM
Post: #52




RE: lambertw, all branches
Thanks, but what about the results given by Wolfram?
I ask, because I never used your flip trick for k a negative integer. Fortunately? 

01252024, 05:57 PM
Post: #53




RE: lambertw, all branches
(01252024 05:14 PM)Gil Wrote: Thanks, but what about the results given by Wolfram? Wolfram's answer is correct. Why would you expect W_{1}(1E330) == W_{+1}(1E330) ? 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) 

01252024, 06:19 PM
Post: #54




RE: lambertw, all branches
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. 

01282024, 11:18 PM
(This post was last modified: 01292024 01:25 AM by Albert Chan.)
Post: #55




RE: lambertw, all branches
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): My Python code used customized LN, to track signed zero, to give correct imag parts. Example, to simulate ln(10*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.790747928e7j) (3.577152064 + 1.841873128e14j) (3.577152064  4.31950186e29j) (3.577152064 + 0.0j) (3.5771520639573 + 0.0j) We can use better solver, without worrying about overshooting! 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.173902656e6j) (3.577152064  1.483011857e16j) (3.577152064  2.350988702e38j) (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 ... 

02012024, 02:17 AM
Post: #56




RE: lambertw, all branches
(01282024 11:18 PM)Albert Chan Wrote: But let's try a *bad* guess, x0=B+j, crossed discontinuity! 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): For debugging purpose, lets add a way to input usersupplied 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.48877473305296e13j) (3.5771520639573 + 5.31405873378874e16j) (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. 

02012024, 04:16 PM
Post: #57




RE: lambertw, all branches
Python W code (post 16)
Quote:s = (im(a) < 0) It is very convoluted. Let's make a table Code: old s updated s 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 

02022024, 11:49 AM
Post: #58




RE: lambertw, all branches
(02012024 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 signedzero 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. 

« Next Oldest  Next Newest »

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