lambertw, all branches, Cas code

02062024, 01:44 PM
(This post was last modified: 02062024 10:36 PM by Albert Chan.)
Post: #3




RE: lambertw, all branches, Cas code
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, lolo%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((y1x)/y))  simulate HP Prime float operation. 214619445125684 Cas> m48(x) := trunc(x * 2^(48frexp(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)(y1x)/y) 214619445125684 Strange, lnp1(x) is worse than ln(1+x), even without correction. Update: I tried this on XCas 1.9.031 (mingw win32) I was expecting same result, since both use 48bits truncated float. XCas> log1p(.1)  ln(1+.1) → 4.88498130835e15 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(1e12) → [0., 39] 

« Next Oldest  Next Newest »

Messages In This Thread 
lambertw, all branches, Cas code  Albert Chan  02042024, 07:34 PM
RE: lambertw, all branches, Cas code  Albert Chan  02052024, 12:46 AM
RE: lambertw, all branches, Cas code  Albert Chan  02062024 01:44 PM

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