Post Reply 
Calculating e^x-1 on classic HPs
02-26-2021, 03:00 PM (This post was last modified: 02-28-2021 08:42 PM by Albert Chan.)
Post: #42
RE: Calculating e^x-1 on classic HPs
(01-14-2016 09:22 PM)Paul Dale Wrote:  Definitely Kahan's method not mine. I can't find a working link to it Sad
I've found where it used to be, but alas no more.

Kahan's article, How Java’s Floating-Point Hurts Everyone Everywhere mentioned expm1(x), although not directly.

JAVAhurt.pdf, page 29, estimate of A(X):=(X-1)/(e^(X-1)-1) Wrote:Real Function Â( Real X ) ;
Real Y, Z ;
Y := X – 1.0 ;
Z := EXP(Y) ;
If Z ≠ 1.0 then Z := LN(Z)/(Z – 1.0) ;
Return  := Z ;
End  .

This subtle program Â(X) is provably always accurate within a few ULPs unless Overflow occurs.

Let x = X-1. For Z = e^x ≠ 1, we have:

A(X) = x / expm1(x) ≈ ln(Z) / (Z-1)       → expm1(x) ≈ (Z-1)/ln(Z) * x

Update: I found another expm1(x) formula derivation, with explanation.
How Futile are Mindless Assessments of Roundoff in Floating-Point Computation ?
page 19, T(z) = expm1(z) / z ≈ (e^z-1)/ln(e^z)

Relative Error of numerator and denominator have same sign, almost cancelled out.

T(x+ε) ≈ (e^(x+ε)-1)/ln(e^(x+ε)) ≈ (e^x*(1+ε)-1) / (x+ε)

T(x+ε)/T(x) ≈ (1 + ε/(1-e^-x)) / (1 + ε/x) ≈ 1 + (ε/(1-e^-x) - ε/x)

If we use e^-x ≈ 1-x, relative error term totally cancelled.
Add 1 more term, using e^-x ≈ 1-x+x²/2, we have:

T(x+ε)/T(x) - 1 ≈ ε/(2-x) ≈ ε/2, for small x
Find all posts by this user
Quote this message in a reply
Post Reply 


Messages In This Thread
Calculating e^x-1 on classic HPs - Dieter - 01-11-2016, 10:20 PM
RE: Calculating e^x-1 on classic HPs - Albert Chan - 02-26-2021 03:00 PM



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