Calculating e^x-1 on classic HPs

[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.