Calculating e^x-1 on classic HPs

[Albert Chan]

(01-14-2016 09:22 PM)Paul Dale Wrote:  Definitely Kahan's method not mine. I can't find a working link to it
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