Lambert W function (for HP Prime)

[Albert Chan]

(11-02-2020 09:53 AM)lyuka Wrote:  The convergence test by Urabe's theorem is used to determine the convergence.
It is a method to judge that the calculation limit has been reached and terminate
the iteration when the correction value does not change or increases.

We are assuming guess is good enough that we never get a false positive.
With your guess function, I think we are safe, but hard to be sure.

Example, if we try iterating for W(1e99), guess = 200, we get this:

Code:

n  O(4) correction   W(1e99)
1  2.98522167315105  202.985221673151  
2  2.98543581371237  205.970657486863  
3  2.98564322111156  208.956300707975  
4  2.98583215049592  211.942132858471  
5  2.98576178491479  214.927894643386  
6  2.98060044381505  217.908495087201  
7  2.87730531320572  220.785800400406  
8  1.69794037402345  222.483740774430  
9  0.06702817368893  222.550768948119  
10 0.00000000763134  222.550768955750 <-- converged

Supplied guess is 10% low. "Chaos" did not stop until 5th iteration.
Trivia: 4-th order correction limited to ±3, similar to 3-rd order Halley ±2

It is better if we apply convergence test after we reached half-precision.
Actually, half-precision test alone is enough here, since next iteration does nothing.

Even better, flatten the curve (see my previous post)
Newton's method beat Householder's method of order 4, converged in 3 iterations.

>>> nest(w_newton2, 200, 1e99)
0 200
1 222.544882427724
2 222.550768955402
3 222.55076895575
4 222.55076895575
5 222.55076895575

---

A few suggestions to LW0H(x) code

1) 0 := t may be wrong, should be t := 0

Better yet, t := 9, with a more robust test:

< UNTIL p == y OR 0 == t OR (t >= r AND 0 <> r);
> UNTIL p == y OR (r < 1 AND t >= r);

Since r < 1, we are in the safe zone.
Convergence test will not generate false positives, even with bad guess.

I removed 0 == t, since p == y already covered it.

2) if we detected chaos, updated iteration likely worse than before.

< return y;
> return p;

3) correction numerator can be made more accurate

< v := y + 2;
< t := s + s;
< t := 3 * (u * (u * v + t) - x * (x * v + t));

We expected u = y*exp(y) ≈ x, thus u-x is exact

> t := 3 * (u-x) * ((u+x)*(y+2) + 2*s)

Update:

4) correction denominator can be made more accurate, as polynomial of (u-x)

< v := (y + 3) * ( u * u + x * (x + 4 * u)) + 6 * s * (u + x + x + s)
> v := 6*((y+1)*u*x + (s+x)*(s+x+u)) + (u-x) * (3*(u+x) + y*(u-x))

The constant term is written this way because around singularity x≈-1/e, y≈-1, s≈-x