Calculating e^x-1 on classic HPs
01-14-2016, 07:50 PM (This post was last modified: 01-14-2016 07:51 PM by Gerson W. Barbosa.)
Post: #21
 Gerson W. Barbosa Senior Member Posts: 1,343 Joined: Dec 2013
RE: Calculating e^x-1 on classic HPs
(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:
(...)

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.
01-14-2016, 09:14 PM
Post: #22
 emece67 Senior Member Posts: 363 Joined: Feb 2015
RE: Calculating e^x-1 on classic HPs
(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 ex–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.

César - Information must flow.
01-14-2016, 09:22 PM (This post was last modified: 01-14-2016 09:35 PM by Paul Dale.)
Post: #23
 Paul Dale Senior Member Posts: 1,653 Joined: Dec 2013
RE: Calculating e^x-1 on classic HPs
(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
01-14-2016, 09:49 PM
Post: #24
 Dieter Senior Member Posts: 2,397 Joined: Dec 2013
RE: Calculating e^x-1 on classic HPs
(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 (ex–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
01-14-2016, 09:55 PM (This post was last modified: 01-14-2016 10:11 PM by Dieter.)
Post: #25
 Dieter Senior Member Posts: 2,397 Joined: Dec 2013
RE: Calculating e^x-1 on classic HPs
(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 ex–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
01-15-2016, 01:15 AM
Post: #26
 Paul Dale Senior Member Posts: 1,653 Joined: Dec 2013
RE: Calculating e^x-1 on classic HPs
I've got a copy of that book at home

- Pauli
01-15-2016, 01:19 AM
Post: #27
 Paul Dale Senior Member Posts: 1,653 Joined: Dec 2013
RE: Calculating e^x-1 on classic HPs
(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 ex–1 function can be expected to have 30+ digit accuracy. Which should be sufficient for an evaluation in 16-digit SP mode. ;-)

ex–1 is calculated using 39 digits.

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

Pauli
01-15-2016, 10:10 AM (This post was last modified: 01-15-2016 08:05 PM by matthiaspaul.)
Post: #28
 matthiaspaul Senior Member Posts: 385 Joined: Jan 2015
RE: Calculating e^x-1 on classic HPs
(01-14-2016 09:55 PM)Dieter Wrote:  Edit: seems to be chapter 1.14.1 of this book.
http://ftp.demec.ufpr.br/CFD/bibliografi...rithms.pdf
http://lib.mdp.ac.id/ebook/Karya%20Umum/...dition.pdf

Greetings,

Matthias

--
"Programs are poems for computers."
01-15-2016, 02:00 PM (This post was last modified: 01-15-2016 02:03 PM by Dieter.)
Post: #29
 Dieter Senior Member Posts: 2,397 Joined: Dec 2013
RE: Calculating e^x-1 on classic HPs
(01-15-2016 01:19 AM)Paul Dale Wrote:  ex–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 ex–1 function by adding the test whether ex>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
01-15-2016, 06:07 PM (This post was last modified: 01-15-2016 06:09 PM by Gerson W. Barbosa.)
Post: #30
 Gerson W. Barbosa Senior Member Posts: 1,343 Joined: Dec 2013
RE: Calculating e^x-1 on classic HPs
(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.
http://ftp.demec.ufpr.br/CFD/bibliografi...rithms.pdf

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.
01-15-2016, 06:46 PM (This post was last modified: 01-15-2016 07:13 PM by Dieter.)
Post: #31
 Dieter Senior Member Posts: 2,397 Joined: Dec 2013
RE: Calculating e^x-1 on classic HPs
(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 ex 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 (ex–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
01-15-2016, 10:45 PM
Post: #32
 Gerson W. Barbosa Senior Member Posts: 1,343 Joined: Dec 2013
RE: Calculating e^x-1 on classic HPs
(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 ex 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 (ex–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.
01-16-2016, 08:19 AM (This post was last modified: 01-16-2016 09:10 AM by Dieter.)
Post: #33
 Dieter Senior Member Posts: 2,397 Joined: Dec 2013
RE: Calculating e^x-1 on classic HPs
(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 (ex–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
01-16-2016, 12:40 PM (This post was last modified: 01-16-2016 12:42 PM by Gerson W. Barbosa.)
Post: #34
 Gerson W. Barbosa Senior Member Posts: 1,343 Joined: Dec 2013
RE: Calculating e^x-1 on classic HPs
(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.
01-16-2016, 01:32 PM (This post was last modified: 01-16-2016 02:02 PM by Dieter.)
Post: #35
 Dieter Senior Member Posts: 2,397 Joined: Dec 2013
RE: Calculating e^x-1 on classic HPs
(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 ex – 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
02-02-2019, 07:22 PM
Post: #36
 Albert Chan Senior Member Posts: 1,159 Joined: Jul 2018
RE: Calculating e^x-1 on classic HPs
(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.

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.
02-02-2019, 07:32 PM
Post: #37
 Albert Chan Senior Member Posts: 1,159 Joined: Jul 2018
RE: Calculating e^x-1 on classic HPs
(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.
02-03-2019, 06:46 PM
Post: #38
 Dieter Senior Member Posts: 2,397 Joined: Dec 2013
RE: Calculating e^x-1 on classic HPs
(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
02-03-2019, 08:28 PM (This post was last modified: 02-03-2019 08:31 PM by Albert Chan.)
Post: #39
 Albert Chan Senior Member Posts: 1,159 Joined: Jul 2018
RE: Calculating e^x-1 on classic HPs
Hi, Dieter

I like your expm1 code very much ! Thanks

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
02-04-2019, 07:28 PM
Post: #40
 Albert Chan Senior Member Posts: 1,159 Joined: Jul 2018
RE: Calculating e^x-1 on classic HPs
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.
 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)