lambertw, all branches

01222024, 01:34 AM
Post: #33




RE: lambertw, all branches
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 overshoot 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 overshoot 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_{0}(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) 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*ksign(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 

« Next Oldest  Next Newest »

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