Calculating e^x1 on classic HPs

01142016, 07:50 PM
(This post was last modified: 01142016 07:51 PM by Gerson W. Barbosa.)
Post: #21




RE: Calculating e^x1 on classic HPs
(01142016 07:38 PM)Dieter Wrote:(01142016 06:50 PM)Gerson W. Barbosa Wrote: Here is my plan B for the hyperbolic function approach: 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 11step program. Anyway, thanks for running the complete tests. Gerson. 

01142016, 09:14 PM
Post: #22




RE: Calculating e^x1 on classic HPs
(01142016 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. 

01142016, 09:22 PM
(This post was last modified: 01142016 09:35 PM by Paul Dale.)
Post: #23




RE: Calculating e^x1 on classic HPs
(01142016 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 

01142016, 09:49 PM
Post: #24




RE: Calculating e^x1 on classic HPs
(01142016 09:22 PM)Paul Dale Wrote: Definitely Kahan's method not mine. I can't find a working link to it 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 100digit challenge, with the original source listed among the references. Indeed it seems to be the method you posted. Dieter 

01142016, 09:55 PM
(This post was last modified: 01142016 10:11 PM by Dieter.)
Post: #25




RE: Calculating e^x1 on classic HPs
(01142016 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 16digit SP mode. ;) Edit: seems to be chapter 1.14.1 of this book. Dieter 

01152016, 01:15 AM
Post: #26




RE: Calculating e^x1 on classic HPs
I've got a copy of that book at home
 Pauli 

01152016, 01:19 AM
Post: #27




RE: Calculating e^x1 on classic HPs
(01142016 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 16digit 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 

01152016, 10:10 AM
(This post was last modified: 01152016 08:05 PM by matthiaspaul.)
Post: #28




RE: Calculating e^x1 on classic HPs
(01142016 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  "Programs are poems for computers." 

01152016, 02:00 PM
(This post was last modified: 01152016 02:03 PM by Dieter.)
Post: #29




RE: Calculating e^x1 on classic HPs
(01152016 01:19 AM)Paul Dale Wrote: e^{x}–1 is calculated using 39 digits. You should stay with the current code. Who cares about the last digit in a 39digit 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 Kahanproof. ;) 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 

01152016, 06:07 PM
(This post was last modified: 01152016 06:09 PM by Gerson W. Barbosa.)
Post: #30




RE: Calculating e^x1 on classic HPs
(01152016 10:10 AM)matthiaspaul Wrote:(01142016 09:55 PM)Dieter Wrote: Edit: seems to be chapter 1.14.1 of this book.See also: Isn't this book Copyright 2002? (Oh, let's forget about that...) HP12C, second algorithm: Code:
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. 

01152016, 06:46 PM
(This post was last modified: 01152016 07:13 PM by Dieter.)
Post: #31




RE: Calculating e^x1 on classic HPs
(01152016 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 tendigit 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 tendigit calculator is rounded to 1,000000005. This indeed is the correctly rounded tendigit 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 

01152016, 10:45 PM
Post: #32




RE: Calculating e^x1 on classic HPs
(01152016 06:46 PM)Dieter Wrote:(01152016 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. So, it's better to stick as close as possible to the modified original algorithm: Code:
Code:
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 erroranalysts 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. 

01162016, 08:19 AM
(This post was last modified: 01162016 09:10 AM by Dieter.)
Post: #33




RE: Calculating e^x1 on classic HPs
(01152016 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. 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. (01152016 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 

01162016, 12:40 PM
(This post was last modified: 01162016 12:42 PM by Gerson W. Barbosa.)
Post: #34




RE: Calculating e^x1 on classic HPs
(01162016 08:19 AM)Dieter Wrote:(01152016 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. Code:
This appeared to work, but some cancellation occurs at a critical region: +1E+0 > +1.718281829E+0 (+1.718281828E+0) 0001 ULP +1E1 > +1.051709180E1 (+1.051709181E1) 0001 ULP +1E2 > +1.005016702E2 (+1.005016708E2) 0006 ULP +1E3 > +1.000500000E3 (+1.000500167E3) 0167 ULP ! +1E4 > +1.000050000E4 (+1.000050002E4) 0002 ULP +1E5 > +1.000000000E5 (+1.000005000E5) 5000 ULP ! +1E6 > +1.000000500E6 (+1.000000050E6) 0000 ULP 1E+0 > 6.321205588E1 (6.321205588E1) 0000 ULP 1E1 > 9.516258195E2 (9.516258196E2) 0001 ULP 1E2 > 9.950166234E3 (9.950166251E3) 0017 ULP ! 1E3 > 9.995001001E4 (9.995001666E4) 0665 ULP ! 1E4 > 9.999500001E5 (9.999500017E5) 0016 ULP ! 1E5 > 9.999950000E6 (9.999950000E6) 0000 ULP 1E6 > 9.999995000E7 (9.999995000E6) 0000 ULP What's a few ULPs among friends anyway? :) At least it gets those right: +9E09 > +9.000000041E09 (+9.000000041E09) 0000 ULP +5E11 > +5.000000000E11 (+5.000000000E11) 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 HP35S trigonometry bug. Gerson. 

01162016, 01:32 PM
(This post was last modified: 01162016 02:02 PM by Dieter.)
Post: #35




RE: Calculating e^x1 on classic HPs
(01162016 12:40 PM)Gerson W. Barbosa Wrote: This appeared to work, but some cancellation occurs at a critical region: Or a few hundred, or a few thousand... ;) (01162016 12:40 PM)Gerson W. Barbosa Wrote: (...) 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 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 As expected, the errors are evenly distributed in this interval. Again, they do not exceed ±5 ULP. Dieter 

02022019, 07:22 PM
Post: #36




RE: Calculating e^x1 on classic HPs
(01162016 12:40 PM)Gerson W. Barbosa Wrote: \[e^{x}1=\frac{e^{x}e^{x}}{e^{x}+1}\] (01162016 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 (u1)/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 u1 and log(u) are not accurate, the ratio is very good. 

02022019, 07:32 PM
Post: #37




RE: Calculating e^x1 on classic HPs  
02032019, 06:46 PM
Post: #38




RE: Calculating e^x1 on classic HPs
(02022019 07:32 PM)Albert Chan Wrote:(01132016 03:46 PM)Gerson W. Barbosa Wrote: \[e^{x}1=\frac{2}{\coth \frac{x}{2}1}\] 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 

02032019, 08:28 PM
(This post was last modified: 02032019 08:31 PM by Albert Chan.)
Post: #39




RE: Calculating e^x1 on classic HPs
Hi, Dieter
I like your expm1 code very much ! Thanks It is amazing about its symmetry with your log1p code: expm1(x) = (u1)  (ln(u)  x) * u, where u = exp(x), rounded log1p( x ) = ln(u)  ((u1)  x) / u, where u = 1+x, rounded 

02042019, 07:28 PM
Post: #40




RE: Calculating e^x1 on classic HPs
Quote:expm1(x) = (u1)  (ln(u)  x) * u, where u = exp(x), rounded A really short prove of above formulas: If z ≈ x, f(x) ≈ f(z) + (xz) * f'(z) expm1(x) = expm1(z) + (xz) * exp(z) log1p(x) = log1p(z) + (xz) * 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)