lambertw, all branches, Cas code

[Albert Chan]

I experimented ln / lnp1 for real agrument.
First, some reference numbers, using lua

lua> p48 = fn'x: local hi,lo = bit.split(x); bit.join(hi, lo-lo%32)' -- Lua to Prime float
lua> m48 = fn'x: trunc(frexp(x) * 2^48)' -- HP Prime mantissa bits
lua> x = p48(.1)
lua> y = p48(1+x)
lua> m48(x)
225179981368524
lua> m48(log1p(x))
214619445125685
lua> m48(log(y))
214619445125674
lua> m48(p48(log(y)) - p48((y-1-x)/y)) -- simulate HP Prime float operation.
214619445125684

Cas> m48(x) := trunc(x * 2^(48-frexp(x)[2]))
Cas> x := .1
Cas> y := 1+x
Cas> m48(x)
225179981368524
Cas> m48(lnp1(x))
214619445125630
Cas> m48(ln(y))
214619445125674
Cas> m48(ln(y)-(y-1-x)/y)
214619445125684

Strange, lnp1(x) is worse than ln(1+x), even without correction.

Update:
I tried this on XCas 1.9.0-31 (mingw win32)
I was expecting same result, since both use 48-bits truncated float.

XCas> log1p(.1) - ln(1+.1)      → 4.88498130835e-15
XCas> ans() * 2^51                → 11.

XCas log1p(.1) gives exactly lua reference, 74 + 11 = 85.
lnp1(.1) bug is only for HP Prime, error = 85 - 30 = 55 ULP

---

I discovered 2 more HP Prime bugs, during experimentation.

1.) round(x) likely defined as floor(x+0.5), even if x is an integer.
I was expecting symmetry, round(x) = sign(x)*round(|x|), like C++ round
Not matching others is not a bug, but round an integer to something else is.

Cas> m := float(m48(x))
Cas>  round( m)       → 225179981368524
Cas> -round(-m)      → 225179981368523

2.) frexp(x) may lose its mantissa.
This is why m48 not defined as trunc(frexp(x)[1] * 2^48)

Cas> frexp(1e-12)      → [0., -39]