Calculating e^x-1 on classic HPs

[Gerson W. Barbosa]

(01-14-2016 07:38 PM)Dieter Wrote:  

(01-14-2016 06:50 PM)Gerson W. Barbosa Wrote:  Here is my plan B for the hyperbolic function approach:
(...)
No additional tests required.

Sorry, Gerson, but I simply could not resist. ;-)
So I ran the 34s test program once again with this method. Here are the resul

Hi, Dieter,

Yes, I know that has really sounded ambiguous, but what I meant were tests for x=0 and tests for the bounds between which the results had acceptable accuracy, as those that would be required in my previous 11-step program. Anyway, thanks for running the complete tests.

Gerson.

[emece67]

(01-14-2016 06:42 PM)Dieter Wrote:  Simple. Write a short test program on the 34s emulator (the "real thing" is quite fast, but too slow for tasks like this). Evaluate e^x–1 for a random argument within the desired domain, using the method you want to test. Then determine the difference to the exact result provided by the internal function. The result is an integer (the number of ULPs the approximation is off). Increment one of, say, 11 registers that count the number of occurences for ≤–5, –4, –3, ... , +4, ≥+5 ULP. Do this 100.000 times. ;-)

Thanks for the clarification. For no reason I assumed the results from the wp34s where compared against those of another "external" tool, such as Mathematica or similar. That's the reason for my surprise.

[Paul Dale]

(01-14-2016 01:48 PM)Dieter Wrote:  Where is this documented? Or is this "Dale's method" ?-)

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.


- Pauli

[Dieter]

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

An hour ago I happend to find a link to a Google books site where Kahan's method for calculating (e^x–1)/x was quoted in a book on the SIAM 100-digit challenge, with the original source listed among the references. Indeed it seems to be the method you posted.

Dieter

[Dieter]

(01-14-2016 09:14 PM)emece67 Wrote:  Thanks for the clarification. For no reason I assumed the results from the wp34s where compared against those of another "external" tool, such as Mathematica or similar. That's the reason for my surprise.

Well, the 34s internally works with 39 significant digits, and even XROM code uses 34 digit precision. So the result of the internal e^x–1 function can be expected to have 30+ digit accuracy. Which should be sufficient for an evaluation in 16-digit SP mode. ;-)

Edit: seems to be chapter 1.14.1 of this book.

Dieter

[Paul Dale]

I've got a copy of that book at home

- Pauli

[Paul Dale]

(01-14-2016 09:55 PM)Dieter Wrote:  Well, the 34s internally works with 39 significant digits, and even XROM code uses 34 digit precision. So the result of the internal e^x–1 function can be expected to have 30+ digit accuracy. Which should be sufficient for an evaluation in 16-digit SP mode. ;-)

e^x–1 is calculated using 39 digits.

I should switch to your method. It will possibly be slightly faster and it is more accurate.


Pauli

[matthiaspaul]

(01-14-2016 09:55 PM)Dieter Wrote:  Edit: seems to be chapter 1.14.1 of this book.

See also:
http://ftp.demec.ufpr.br/CFD/bibliografi...rithms.pdf
http://lib.mdp.ac.id/ebook/Karya%20Umum/...dition.pdf

Greetings,

Matthias

[Dieter]

(01-15-2016 01:19 AM)Paul Dale Wrote:  e^x–1 is calculated using 39 digits.

I should switch to your method. It will possibly be slightly faster and it is more accurate.

You should stay with the current code. Who cares about the last digit in a 39-digit calculation when only 34 digits are returned to the user? And who knows the accuracy of the exp and log functions? Finally, the current method is Kahan-proof. ;-)

But you can speed up the e^x–1 function by adding the test whether e^x>2 resp. v>1. This avoids, for larger arguments, a call of the log function which, as far as I remember, is particulary slow on the 34s.

Dieter

[Gerson W. Barbosa]

(01-15-2016 10:10 AM)matthiaspaul Wrote:  

(01-14-2016 09:55 PM)Dieter Wrote:  Edit: seems to be chapter 1.14.1 of this book.

See also:
http://ftp.demec.ufpr.br/CFD/bibliografi...rithms.pdf

Isn't this book Copyright 2002? (Oh, let's forget about that...)

HP-12C, second algorithm:

Code:


01- x=0
02- GTO 00
03- ENTER
04- e^x
05- 1
06- CHS 
07- x<>y 
08- + 
09- LSTx
10- LN
11- /
12- *
13- GTO 00

I'm currently running this on an ordinary 12C and Dieter's program on a 12C+, so I cannot compare speeds, but Dieter's should be slightly faster because it requires no tests.

Gerson.

P.S.: Thanks for the link! UFPR is just a walking distance from my place.

[Dieter]

(01-15-2016 06:07 PM)Gerson W. Barbosa Wrote:  I'm currently running this on an ordinary 12C and Dieter's program on a 12C+, so I cannot compare speeds, but Dieter's should be slightly faster because it requires no tests.

Hm, seems you forgot an essential test here: if x is so small that e^x rounds to 1 (here –5 E–11 < x < 5 E–10), the term (u–1)/ln u becomes 0/0... #-) In such a case the program should return x. This additional test would also handle the x=0 case that's now tested in the first line.

BTW, W. Kahan's method is supposed to calculate (e^x–1)/x. Even if this quotient is evaluated exactly, i.e. for small x very close to 1, there is always a potential error of half an ULP or 5 E–10 on a ten-digit device. This translates to a possible error of up to 5 ULP in the final result, i.e. after the multiplication with x.

Example: x=9 E–9. The exact quotient is 1,000000045... which on a ten-digit calculator is rounded to 1,000000005. This indeed is the correctly rounded ten-digit value. But the final result is returned as 9,000000045 E–9 instead of 9,0000000405000...E–9, i.e. there is an error of +4,5 ULP.

Dieter

[Gerson W. Barbosa]

(01-15-2016 06:46 PM)Dieter Wrote:  

(01-15-2016 06:07 PM)Gerson W. Barbosa Wrote:  I'm currently running this on an ordinary 12C and Dieter's program on a 12C+, so I cannot compare speeds, but Dieter's should be slightly faster because it requires no tests.

Hm, seems you forgot an essential test here: if x is so small that e^x rounds to 1 (here –5 E–11 < x < 5 E–10), the term (u–1)/ln u becomes 0/0... #-) In such a case the program should return x. This additional test would also handle the x=0 case that's now tested in the first line.

BTW, W. Kahan's method is supposed to calculate (e^x–1)/x. Even if this quotient is evaluated exactly, i.e. for small x very close to 1, there is always a potential error of half an ULP or 5 E–10 on a ten-digit device. This translates to a possible error of up to 5 ULP in the final result, i.e. after the multiplication with x.

Example: x=9 E–9. The exact quotient is 1,000000045... which on a ten-digit calculator is rounded to 1,000000005. This indeed is the correctly rounded ten-digit value. But the final result is returned as 9,000000045 E–9 instead of 9,0000000405000...E–9, i.e. there is an error of +4,5 ULP.

So, it's better to stick as close as possible to the modified original algorithm:

Code:


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

Code:


01- ENTER
02- e^x
03- 1
04- CHS
05- x<>y 
06- +
07- x=0
08- GTO 14
09- LSTx
10- LN
11- /
12- *
13- x<>y
14- Rv
15- GTO 00

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. The program is now two stpes longer and slightly slower than the previous one. Anyway, a good exercise after two cups of cheap Argentinian wine (they're still better than our finest ones :-)

Quoting Professor Kahan, as found in the book, perhaps applied to your proposed method:

"There will always be a small but steady demand for error-analysts to ...
expose bad algorithms' big errors and, more important,
supplant bad algorithms with provably good ones.
- WILLIAM M. KAHAN, Interval Arithmetic Options in the
Proposed IEEE Floating Point Arithmetic Standard (1980)"

Gerson.

[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

[Gerson W. Barbosa]

(01-16-2016 08:19 AM)Dieter Wrote:  

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

Code:


01- ENTER
02- e^x
03- ENTER
04- 1/x
05- -
06- x=0
07- GTO 15
08- x<>y
09- CLx
10- LSTx
11- 1
12- +
13- /
14- x<>y
15- Rv
16- GTO 00

This appeared to work, but some cancellation occurs at a critical region:

+1E+0 -> +1.718281829E+0 (+1.718281828E+0) 0001 ULP
+1E-1 -> +1.051709180E-1 (+1.051709181E-1) 0001 ULP
+1E-2 -> +1.005016702E-2 (+1.005016708E-2) 0006 ULP
+1E-3 -> +1.000500000E-3 (+1.000500167E-3) 0167 ULP !
+1E-4 -> +1.000050000E-4 (+1.000050002E-4) 0002 ULP
+1E-5 -> +1.000000000E-5 (+1.000005000E-5) 5000 ULP !
+1E-6 -> +1.000000500E-6 (+1.000000050E-6) 0000 ULP

-1E+0 -> -6.321205588E-1 (-6.321205588E-1) 0000 ULP
-1E-1 -> -9.516258195E-2 (-9.516258196E-2) 0001 ULP
-1E-2 -> -9.950166234E-3 (-9.950166251E-3) 0017 ULP !
-1E-3 -> -9.995001001E-4 (-9.995001666E-4) 0665 ULP !
-1E-4 -> -9.999500001E-5 (-9.999500017E-5) 0016 ULP !
-1E-5 -> -9.999950000E-6 (-9.999950000E-6) 0000 ULP
-1E-6 -> -9.999995000E-7 (-9.999995000E-6) 0000 ULP

What's a few ULPs among friends anyway? :-)

At least it gets those right:

+9E-09 -> +9.000000041E-09 (+9.000000041E-09) 0000 ULP
+5E-11 -> +5.000000000E-11 (+5.000000000E-11) 0000 ULP

That should be expected by just looking at the identity I have used:

\[e^{x}-1=\frac{e^{x}-e^{-x}}{e^{-x}+1}\]

This resembles a bit the HP-35S trigonometry bug.

Gerson.

[Dieter]

(01-16-2016 12:40 PM)Gerson W. Barbosa Wrote:  This appeared to work, but some cancellation occurs at a critical region:
(...)
What's a few ULPs among friends anyway? :-)

Or a few hundred, or a few thousand... ;-)

(01-16-2016 12:40 PM)Gerson W. Barbosa Wrote:  (...)
That should be expected by just looking at the identity I have used:

\[e^{x}-1=\frac{e^{x}-e^{-x}}{e^{-x}+1}\]

It works better after replacing e^x – e^–x with 2 sinh(x) – provided an accurate sinh function is available. I tried this on the 34s and in 100.000 random samples there were no errors beyond ±7 ULP. About 75% of the results were within ±1 ULP.

Finally I did a test run with one million random values in [ln 0.9, ln 2] with the method I suggested, and the results are essentially the same as before:

Code:

seed = 47,11   n = 1000000

ULP      #cases
--------------------------
 -5          12 ─── 0,001%
 -4        1092 ─┐
 -3         674  ├─ 0,24%
 -2         676 ─┘
 -1      167703 ─┐
 ±0      659710  ├─ 99,5%
 +1      167784 ─┘
 +2         666 ─┐
 +3         677  ├─ 0,23%
 +4         992 ─┘
 +5          14 ─── 0,001%

Again, there are two peaks at +4 and –4 ULP that also showed up in earlier tests. I wonder where these come from. Any idea?

As already mentioned, for arguments in [ln10, ln11[ or [ln100, ln101[ etc. the last digit is lost, which leads to an interesting error distribution. Here are some results for x in [ln10, ln11]:

Code:

seed = 47,11    n = 100000

ULP      #cases
--------------------------
 -5       4990 ─── 5%
 -4      10025 ─┐
 -3       9952  ├─ 30%
 -2       9975 ─┘
 -1      10002 ─┐
 ±0      10024  ├─ 30%
 +1      10170 ─┘
 +2       9888 ─┐
 +3       9941  ├─ 30%
 +4      10042 ─┘
 +5       4991 ─── 5%

As expected, the errors are evenly distributed in this interval. Again, they do not exceed ±5 ULP.

Dieter

[Albert Chan]

(01-16-2016 12:40 PM)Gerson W. Barbosa Wrote:  \[e^{x}-1=\frac{e^{x}-e^{-x}}{e^{-x}+1}\]

(01-16-2016 01:32 PM)Dieter Wrote:  ... there are two peaks at +4 and –4 ULP that also showed up in earlier tests. I wonder where these come from. Any idea?

My guess is that error goes at opposite direction, thus make things worse.

exp(x) = 1/exp(-x), so if exp(x) error is positive, exp(-x) is negative.

This made numerator even bigger.
Worse, denominator shrink. Errors reinforcing each other, making bad peaks !

---
This is very different than Kahan's calculating (exp(x)-1)/x via (u-1)/log(u)

For u close to 1, both numerator and denominator error goes the same way.
A bonus is that the size of error also similar size, cancelling each other.

Although calculated values of u-1 and log(u) are not accurate, the ratio is very good.

[Albert Chan]

(01-13-2016 03:46 PM)Gerson W. Barbosa Wrote:  \[e^{x}-1=\frac{2}{\coth \frac{x}{2}-1}\]
... this will require two additional steps to handle x=0

Is there any reason to use a possible divide-by-0 coth(x/2) ?

If t = tanh(x/2), we have exp(x) - 1 = 2 t / (1 - t), avoided the singularity.

[Dieter]

(02-02-2019 07:32 PM)Albert Chan Wrote:  

(01-13-2016 03:46 PM)Gerson W. Barbosa Wrote:  \[e^{x}-1=\frac{2}{\coth \frac{x}{2}-1}\]
... this will require two additional steps to handle x=0

Is there any reason to use a possible divide-by-0 coth(x/2) ?

If t = tanh(x/2), we have exp(x) - 1 = 2 t / (1 - t), avoided the singularity.

Both methods suffer from digit cancellation for large x where coth and tanh both tend towards 1. That's why I prefer the original method in the first post. But maybe this can also cause numeric problems – did you find any?

Dieter

[Albert Chan]

Hi, Dieter

I like your expm1 code very much ! Thanks

It is amazing about its symmetry with your log1p code:

expm1(x) = (u-1) - (ln(u) - x) * u, where u = exp(x), rounded

log1p( x ) = ln(u) - ((u-1) - x) / u, where u = 1+x, rounded

[Albert Chan]

Quote:expm1(x) = (u-1) - (ln(u) - x) * u, where u = exp(x), rounded

log1p( x ) = ln(u) - ((u-1) - x) / u, where u = 1+x, rounded

A really short prove of above formulas: If z ≈ x, f(x) ≈ f(z) + (x-z) * f'(z)

expm1(x) = expm1(z) + (x-z) * exp(z)
log1p(x) = log1p(z) + (x-z) * 1/(1+z)

Let z = ln(exp(x)), rounded for expm1 eqn, and z = (1+x)-1, rounded for log1p, we get the quoted formulas.