Calculating e^x-1 on classic HPs

[Dieter]

(02-04-2019 07:28 PM)Albert Chan Wrote:  A really short prove of above formulas: ...

Take a look at post #2 and #6 in this thread. ;-)
Yes, both formulas are based on a (one term) Taylor series.

Dieter

[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

[Albert Chan]

(02-26-2021 03:00 PM)Albert Chan Wrote:  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

Kahan's expm1 formula is good, but we can make it better.
Let h = x - ln(Z), rewrite it as corrections to Z-1:

expm1(x) ≈ (Z-1) + (Z-1)/(x-h) * h ≈ (Z-1) + √Z * h, when x is tiny

√Z version avoided testing divide-by-zero, when (Z-1)/ln(Z) = 0/0 = ?

Compare this against Dieter's version: expm1(x) ≈ (Z-1) + Z * h

Statistically, both versions seems to be equally good. This may be why.
For small x, Z ≈ 1+x, √Z ≈ 1+x/2:

Difference between two version is x*(h/2), correction too small to affect result.