lambertw, all branches
|
01-23-2024, 06:17 PM
Post: #41
|
|||
|
|||
RE: lambertw, all branches
(01-23-2024 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 Wk(-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 = lnk(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 |
|||
01-23-2024, 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) |
|||
01-23-2024, 11:00 PM
(This post was last modified: 01-24-2024 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 (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. |
|||
01-24-2024, 03:18 PM
(This post was last modified: 01-24-2024 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 |
|||
01-24-2024, 08:53 PM
(This post was last modified: 01-25-2024 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 Wk((2k)*pi*I) = (2k)*pi*I W1(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 Wk(-(2k+1/2)*pi) = (2k+1/2)*pi*I W1(-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 W1(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, W1(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 Wk(-odd*pi*I) = odd*pi*I W1(-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 Wk((odd+1/2)*pi) = (odd+1/2)*pi*I W-2(-5/2*pi) = -5/2*pi*I W1(3/2*pi) = 3/2*pi*I |
|||
01-25-2024, 12:37 AM
Post: #46
|
|||
|
|||
RE: lambertw, all branches
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? |
|||
01-25-2024, 01:10 AM
(This post was last modified: 01-25-2024 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 non-negative. But W-1(-1E-330)= -697.32277629546016099540752740546566360568199201428307824605649289390557552642032068023149199134240178698554915132 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. |
|||
01-25-2024, 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). |
|||
01-25-2024, 07:02 AM
Post: #49
|
|||
|
|||
RE: lambertw, all branches
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 = Wk(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)) --> Wk(a) = conj(W-k(conj(a))) |
|||
01-25-2024, 10:09 AM
(This post was last modified: 01-25-2024 10:13 AM by Gil.)
Post: #50
|
|||
|
|||
RE: lambertw, all branches
Thanks, Albert.
And how do you explain these different results? W-1(-1E-330)= -697.32277629546016099540752740546566360568199201428307824605649289390557552642032068023149199134240178698554915132 And W1(-1E-330)= -697.32281706295494653348757314785952050772713371449039168 +6.2922084418351405256893689216267218827742283825238692107 i Are they correct? |
|||
01-25-2024, 04:13 PM
Post: #51
|
|||
|
|||
RE: lambertw, all branches
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: Wk(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, 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 ... |
|||
01-25-2024, 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? |
|||
01-25-2024, 05:57 PM
Post: #53
|
|||
|
|||
RE: lambertw, all branches
(01-25-2024 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(-1E-330) == W+1(-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) |
|||
01-25-2024, 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. |
|||
01-28-2024, 11:18 PM
(This post was last modified: 01-29-2024 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(-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 ... |
|||
02-01-2024, 02:17 AM
Post: #56
|
|||
|
|||
RE: lambertw, all branches
(01-28-2024 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 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. |
|||
02-01-2024, 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 |
|||
02-02-2024, 11:49 AM
Post: #58
|
|||
|
|||
RE: lambertw, all branches
(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. |
|||
« Next Oldest | Next Newest »
|
User(s) browsing this thread: 5 Guest(s)