(28/48/50) Lambert W Function

[Albert Chan]

(03-21-2023 01:53 PM)John Keith Wrote:  I imagine that the best compromise would be to use extended precision reals, but SysRPL is beyond my skill set.

Just to be clear, we don't need extended precison numbers.
My Lua code only use the same 64-bits float throughout.

All we need is more precise constant 1/e, so that (y-1/e) does not suffer cancellation errors.

Quote:The estimation code I used for branch 0 is amazingly accurate over the entire numerical range representable by Saturn-based HP's. Also, the recurrence that I used converges faster than Newton's method while not being much larger in code size.

I am unfamilar with RPL, can you post actual formula for W(a) guess, and the iteration formula?

Convergence can be affected greatly by what the iteration formula look like.
Newton's method is not necessarily much slower than higher order method.

To give an example, lets do W(a = 1e10), with guess a0 = log1p(a) ≈ 23.0259
Note: nestlist is my defined function, not part of XCas

XCas> a := 1e10; a0 := log1p(a)
XCas> f := x*e^x - a
XCas> nestlist(unapply(x-f/f',x), a0, 3)                      → [23.0259, 22.1091, 21.2606, 20.568]
XCas> nestlist(unapply(x-f/(f'-f''/2*f/f'),x), a0, 3)      → [23.0259, 21.2714, 20.1788, 20.029]

This is a bad setup, with huge derivatives, f' = (x+1)*e^x, f'' = (x+2)*e^x
If a is bigger , bad setup is going to take longer to converge.

Let's try a better setup.

XCas> f := x - ln(a/x)
XCas> nestlist(unapply(x-f/f',x), a0, 3)                      → [23.0259, 20.0198, 20.0287, 20.0287]
XCas> nestlist(unapply(x-f/(f'-f''/2*f/f'),x), a0, 3)      → [23.0259, 20.0279, 20.0287, 20.0287]

Here, there is not much need for complicated halley's method. Newton is almost as good.
If a is bigger , better setup is going to converge faster! Perhaps in 1 iteration!

Quote:Of more interest to me is computing W(z) for branches other than 0 and -1. I have not been able to find any methods for doing so. Any insight in this area would be greatly appreciated.

Python mpmath lambertw(z, k=0) support all branches.

https://mpmath.org/doc/0.19/functions/po...w-function
https://github.com/mpmath/mpmath/blob/ma...ns.py#L464