Calculating e^x-1 on classic HPs

[Dieter]

(01-15-2016 10:45 PM)Gerson W. Barbosa Wrote:  So, it's better to stick as close as possible to the modified original algorithm:

The original algorithm calculated (e^x–1)/x, so here a multiplication by x is required which is missing in the first "if" branch. Correction:

Code:

% Algorithm 2.
y = e^x
if y = 1
   f = x
else
   f = (y - 1)/(log y)*x
end if

OTOH the 12C code seems to be OK.

Well, almost. ;-) Something is missing both in "algorithm 2" and the 12C code: the case where y underflows to zero, e.g. for x=–300. This would cause an attempt at calculating ln 0. In this case the result –1 has to be returned. So it's more like this:

Code:

% Algorithm 2a.
y = e^x
if y = 0 then
   f = -1
else
   if y = 1 then
      f = x
   else
      f = (y - 1)/(log y)*x
   end if
end if

(01-15-2016 10:45 PM)Gerson W. Barbosa Wrote:  I'd imagined the final multiplication by x might be an issue, as you've pointed out, but in my (few) tests I didn't find significant differences compared to your results.

You'll find them after a few hundred thousand samples. ;-)

In post #20 I provided a table with the error distribution in ULPs. I have updated this table to show differences up to 9 ULP which occured with my implementation of the Kahan method.

Dieter