Calculating e^x1 on classic HPs

01112016, 10:20 PM
(This post was last modified: 01112016 10:22 PM by Dieter.)
Post: #1




Calculating e^x1 on classic HPs
There is no doubt that the ln(1+x) and e^{x}–1 functions that appeared with the 41C are extremely useful as they allow more precise calculations for arguments close to zero, which otherwise would round to zero. The 15C Advanced Functions Handbook included a short routine for emulating ln(1+x) on calculators that no not feature this function, and another (IMHO slightly better) method has been posted in this forum.
On the other hand, the e^{x}–1 case seemed more challenging. Thomas Klemm posted a solution that included the hyperbolic sine for a short and elegant implementation of this function. However, not all HPs offer hyperbolic functions. Here is a method that should work on most classic HPs. I think it works fine but I did not do any extensive tests. So try it and see what you get. Coded in VBA, the algorithm is as follows: Code: Function expm1(x) On a classic HP it could look like this: Code: ENTER This returns e^{x} in Y and e^{x}–1 in X. What do you think? Dieter 

01122016, 01:41 AM
(This post was last modified: 01122016 01:49 AM by Thomas Klemm.)
Post: #2




RE: Calculating e^x1 on classic HPs
(01112016 10:20 PM)Dieter Wrote: What do you think? Brilliant! Let's assume: \(u=e^x+\varepsilon\) where \(\varepsilon \ll 1\) Then: \(e^x1=u1\varepsilon\) But since: \(u=e^x+\varepsilon=e^x(1+e^{x}\varepsilon)\) \(\log(u)=\log(e^x)+\log(1+e^{x}\varepsilon)\) Now we can use: \(log(1+x)\approx x\) for small \(x\). \(\log(u)\approx x+e^{x}\varepsilon\) Or: \(\varepsilon \approx (\log(u)x)e^x\) This leads to: \(e^x1\approx u1+(x\log(u))u\) Cheers Thomas 

01122016, 05:13 AM
(This post was last modified: 01122016 05:13 AM by Claudio L..)
Post: #3




RE: Calculating e^x1 on classic HPs
(01112016 10:20 PM)Dieter Wrote: What do you think? I just checked at 2000 digits precision and the result was good. I always thought e^x1 was the easiest of all to implement, since the power series expansion at 0 starts with a 1. All you have to do is skip the first term and with relatively few terms in the summation you get the answer. No need to worry about argument reduction because the whole point is to use it around zero. Question is: would a program computing terms of the power series in a loop be faster? 

01122016, 02:02 PM
(This post was last modified: 01122016 03:06 PM by Dieter.)
Post: #4




RE: Calculating e^x1 on classic HPs
(01122016 05:13 AM)Claudio L. Wrote: I just checked at 2000 digits precision and the result was good. Fine. My results on 10 and 12digit calculators did not show any results with larger errors than ±1 ULP, mostly because of different rounding. The vast majoriy was dead on. (01122016 05:13 AM)Claudio L. Wrote: I always thought e^x1 was the easiest of all to implement, since the power series expansion at 0 starts with a 1. All you have to do is skip the first term and with relatively few terms in the summation you get the answer. No need to worry about argument reduction because the whole point is to use it around zero. Right. This is true if you want to evaluate e^{x}–1 in terms of simple arithmetics, i.e. without the use of other transcendental functions. Which may be a good idea on a computer, but not on a calculator where speed and memory are much more limited. (01122016 05:13 AM)Claudio L. Wrote: Question is: would a program computing terms of the power series in a loop be faster? I'd say this depends on the implementation of the required exponential and logarithm functions. On most common processors this is done in hardware (the FPU provides these functions). On the required interval (ln 0,9 < x < ln 2) and 1516 digit accuracy a continued fraction method with 7 terms or a 16term series expansion will do. Now, how much faster are addition, multiplication and division compared to exp and log? The continued fraction method requires 14 divisions, 7 additions and 7 subtractions, plus a bit overhead for the range check and the loop counter. I suppose this is still faster than using exp and log. On a computer, that is. For the record, here is a VBA implementation using the continued fraction method and hardcoded seven terms (k = 13, 11, 9, 7, 5, 3, 1): Code: Function expm1cf(x) Dieter 

01122016, 03:34 PM
Post: #5




RE: Calculating e^x1 on classic HPs  
01122016, 06:12 PM
(This post was last modified: 01152016 02:12 PM by Dieter.)
Post: #6




RE: Calculating e^x1 on classic HPs
(01122016 01:41 AM)Thomas Klemm Wrote: Brilliant! Thank you very much. Indeed the idea behind this is rather trivial. I'd like to explain this method with an example. Assume x = 3,141592654 E6 and a 10digit calculator. Evaluating u = e^{x} with 10 digit precision yields u = 1,000003142 Accordingly, e^{x}–1 is evaluated as u1 = 0,000003142 But that's not e^{x} or e^{x}1. The correct values are e^{x} = 1,0000031415975888... e^{x}–1 = 0,0000031415975888... The 10digit result we got actually is not e^{x} but e^{ln u} = e^{3,141995064E6} (resp. this minus one). Simply because only e^{ln u} exactly equals u: e^{ln u} = u = 1,000003142000000... e^{ln u } – 1 = 0,000003142000000... So we got an exact result for an argument that is slightly off. Precisely, we want the result for an argument that it is off by dx = x – ln u Now the rest is easy. With the first terms of a Taylor series we get f(x+dx) ≈ f(x) + dx · f'(x) and thus: e^{x}1 ≈ e^(ln u)  1 + (x  ln u) · d[e^{x}–1]/dx = (u  1) + (x  ln u) · e^{x} ...and since u agrees with e^{x} to machine accuracy we finally get: e^{x}1 ≈ (u  1) + (x  ln u) · u That's it. Dieter 

01122016, 06:30 PM
(This post was last modified: 01122016 07:48 PM by Dieter.)
Post: #7




RE: Calculating e^x1 on classic HPs
(01122016 03:34 PM)ElectroDuende Wrote:(01122016 05:13 AM)Claudio L. Wrote: Question is: would a program computing terms of the power series in a loop be faster? An iterative approach on such a classic calculator would be much slower and require more memory than the suggested method. The two transcendental functions are evaluated faster (and each in one single line of code) than, say, 12 loops of a power series (for 10digit accuracy), plus the required overhead. Addendum: I just tried the CF approach in a program for the 41 series. Ten digits require just four loops, but nevertheless the program needs > 3 seconds for the calculation. Compared to this, the method suggested in my initial post returns its result in less than a second – and occupies only 1/3 of the memory. Dieter 

01122016, 08:16 PM
(This post was last modified: 01122016 08:43 PM by Dieter.)
Post: #8




RE: Calculating e^x1 on classic HPs
(01112016 10:20 PM)Dieter Wrote: On a classic HP it could look like this: The code shown in the initial post uses the complete stack, OTOH it returns both e^{x} and e^{x}–1 in Y and X. The following version returns e^{x}–1 while it saves the Yregister. It also does not use R^ which is missing on some of the classic calculators. Code: e^x And using the 35s in equation mode (as well as one data register) it can even be done in a single line: Code: EXP(REGX)►U–1+(REGX–LN(U))*U Dieter 

01122016, 09:44 PM
Post: #9




RE: Calculating e^x1 on classic HPs
(01122016 06:30 PM)Dieter Wrote:(01122016 03:34 PM)ElectroDuende Wrote: I suppose that the main reason not to use the programloop approach is that not all the "classic" models are programmable... Great, now you have a very elegant solution, with hard proof that it's also speed and memory efficient. Good work. 

01132016, 03:46 PM
Post: #10




RE: Calculating e^x1 on classic HPs
(01122016 08:16 PM)Dieter Wrote: The following version returns e^{x}–1 while it saves the Yregister. It also does not use R^ which is missing on some of the classic calculators. Very nice, thank you! Just tested it on the HP12C. I've tried an hyperbolic function approach on the HP15C using the identity \[e^{x}1=\frac{2}{\coth \frac{x}{2}1}\] Code:
However this will require two additional steps to handle x=0. Also, for x > 3 the results start to lose accuracy as tanh(x/2) approaches 1. Congratulations for your brilliant method! Gerson. 

01132016, 06:43 PM
(This post was last modified: 01132016 06:57 PM by Dieter.)
Post: #11




RE: Calculating e^x1 on classic HPs
(01132016 03:46 PM)Gerson W. Barbosa Wrote: Just tested it on the HP12C. I've tried an hyperbolic function approach on the HP15C using the identity Ah, that's the inverse of the method I suggested for ln(1+x) in 2014, cf. this thread. (01132016 03:46 PM)Gerson W. Barbosa Wrote: However this will require two additional steps to handle x=0. Also, for x > 3 the results start to lose accuracy as tanh(x/2) approaches 1. Accuracy problems have also been addressed in the mentioned 2014 thread. That's why back then I suggested another ln(1+x) method using sinh. On the other hand a dedicated e^{x}–1 function is only required for ln 0,9 < x < ln 2, and here your method works fine. One might add two additional tests to check whether the argument is within these bounds or not. Which again makes the method less elegant... Dieter 

01142016, 02:47 AM
(This post was last modified: 01142016 03:14 PM by Gerson W. Barbosa.)
Post: #12




RE: Calculating e^x1 on classic HPs
(01132016 06:43 PM)Dieter Wrote: On the other hand a dedicated e^{x}–1 function is only required for ln 0,9 < x < ln 2, and here your method works fine. One might add two additional tests to check whether the argument is within these bounds or not. Which again makes the method less elegant... Still not elegant, but I decided to try a series approach. Starting with an empirical approximation (I noticed e^x  1 ~ e^(ln(x) + x/2) for small x), I came up with a new Taylor series  Well, at least apparently WA missed it... \[e^{x}1=e^{\ln(x)+p(x)}\] \[\ln(e^{x}1)=\ln(x)+p(x)\] \[p(x)=\ln(e^{x}1)\ln(x)\] \[p(x)=\ln\left (\frac{e^{x}1}{x}\right)\] \[\cdots\] \[p(x)=\frac{x}{2}+\frac{x^{2}}{24}\frac{x^{4}}{2880}+\frac{x^{6}}{181440}\frac{x^{8}}{9676800}+\frac{x^{10}}{479001600}\frac {691\cdot x^{12}}{15692092416000}+\cdots\] \[p(x)=\frac{x}{2}+\sum_{k=1}^{\infty} \frac{B_{2k}\cdot x^{2k}}{2k\cdot (2k)!}\] \[\therefore e^{x}1=e^{\ln(x)+\frac{x}{2}+\sum_{k=1}^{\infty} \frac{B_{2k}\cdot x^{2k}}{2k\cdot (2k)!} }\] Only four terms of the summand would suffice for 10digit calculators in the specified range (five terms for 12digit calculators). This would take too many steps on the Voyagers, however. Anyway, the first two or three terms might be useful for certain applications if x is always positive and small enough. Just for the sake of illustration and testing, here is an HP49G/G+/50g version: Code:
2 LN EXM1 > 1. 1.23E8 EXM1 > 1.23000000759E8 0.1 EXM1 RE > 9.51625819644E2 0.1 EXPM > .095162581964 ; builtin EXPM function 1.23E8 EXPM > 1.23000000756E8 Best regards, Gerson. Edited per Dieter's observation below (1.23E8 EXPM > 1.23000000756E8, not 1.2300000076E8) 

01142016, 07:20 AM
Post: #13




RE: Calculating e^x1 on classic HPs
A nice identity using hyperbolics is: e^{x}1 = tanh(x/2) (e^{x}+1)
The 34S uses: Code:
As William Kahan said of this: nearly full working relative accuracy despite cancellation no matter how tiny x may be. Dieter's method avoids the conditionals which is a plus for some machines.  Pauli 

01142016, 01:44 PM
Post: #14




RE: Calculating e^x1 on classic HPs
(01142016 02:47 AM)Gerson W. Barbosa Wrote: Still not elegant, but I decided to try a series approach. Starting with an empirical approximation (I noticed e^x  1 ~ e^(ln(x) + x/2) for small x), I came up with a new Taylor series  Well, at least apparently WA missed it... Bernoulli numbers? Waaaayyyy to complicated. ;) (01142016 02:47 AM)Gerson W. Barbosa Wrote: Only four terms of the summand would suffice for 10digit calculators in the specified range (five terms for 12digit calculators). Maybe you can take a look at the CF method in post #4. No transcendetal functions (except when e^x–1 is calculated the trivial way outside of the critical domain) and four loops will do for 10 digits as well. (01142016 02:47 AM)Gerson W. Barbosa Wrote: 2 LN EXM1 > 1. For the record: the exact 12digit results are 1 1,23000000756 E8 9,51625819640 E2 The method I suggested happens to return exactly these results. I did some further tests with several runs of 100.000 random numbers in [ln 0,9, ln 2] on a 34s in SP mode (16 digits) which showed that about 2/3 of the results were dead on an 99,5% within 1 ULP. However, there are rare cases where the result may be off by several ULP. BTW, does the builtin EXPM function really return 1.2300000076 E8, i.e. +4 ULP compared to the true result? Dieter 

01142016, 01:48 PM
(This post was last modified: 01142016 02:53 PM by Dieter.)
Post: #15




RE: Calculating e^x1 on classic HPs
(01142016 07:20 AM)Paul Dale Wrote: Ah, interesting. I have been looking for an equivalent to Kahan's ln(1+x) method, but I could not find anything. Where is this documented? Or is this "Dale's method" ?) Adding another test (before the other two) may speed up the code as it avoids the log in the last line: Code: if v >= 1, return v I did a few test runs on a 34s in SP mode with 100000 random numbers each. About 45% were exact, and another 45% was within ±1 ULP. Another 9% was ±2 ULP and the rest had larger errors. So accuracy varies a bit more than with the method I suggested (90% vs. 99,5% within ±1 ULP). Note: I originally posted that the results using your method were at most –1 ULP off. Sorry, this is wrong, there was an error in the test routine. #) Dieter 

01142016, 03:30 PM
Post: #16




RE: Calculating e^x1 on classic HPs
(01142016 01:44 PM)Dieter Wrote:(01142016 02:47 AM)Gerson W. Barbosa Wrote: Still not elegant, but I decided to try a series approach. Starting with an empirical approximation (I noticed e^x  1 ~ e^(ln(x) + x/2) for small x), I came up with a new Taylor series  Well, at least apparently WA missed it... You are quite right! Anyway, that was not supposed to be a valid method. I just liked the series, despite its uselessness :) (01142016 01:44 PM)Dieter Wrote:(01142016 02:47 AM)Gerson W. Barbosa Wrote: 2 LN EXM1 > 1. I'd made a mistake. Fixed. Thanks! Gerson. 

01142016, 03:55 PM
Post: #17




RE: Calculating e^x1 on classic HPs
(01142016 01:48 PM)Dieter Wrote: I did a few test runs on a 34s in SP mode with 100000 random numbers each. About 45% were exact, and another 45% was within ±1 ULP. Another 9% was ±2 ULP and the rest had larger errors. So accuracy varies a bit more than with the method I suggested (90% vs. 99,5% within ±1 ULP). I'm really curious, how do you do such tests? 

01142016, 06:42 PM
(This post was last modified: 01142016 07:10 PM by Dieter.)
Post: #18




RE: Calculating e^x1 on classic HPs
(01142016 03:55 PM)emece67 Wrote:(01142016 01:48 PM)Dieter Wrote: I did a few test runs on a 34s in SP mode with 100000 random numbers each. (...) 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. ;) Using the same seed for the random number generator, this is what I got for the Kahan method resp. the one suggested in this thread. The results of the two different runs are quite similar. Code: seed = 4711 n = 100000 seed = 0,815 n = 100000 Re. the meaning of the two seed values, cf. this thread. ;) Dieter 

01142016, 06:50 PM
(This post was last modified: 01142016 06:53 PM by Gerson W. Barbosa.)
Post: #19




RE: Calculating e^x1 on classic HPs
(01132016 06:43 PM)Dieter Wrote: Accuracy problems have also been addressed in the mentioned 2014 thread. That's why back then I suggested another ln(1+x) method using sinh. Here is my plan B for the hyperbolic function approach: Code:
No additional tests required. On the HP15C I get the same results of the HP12C running your program, except for occasional differences of one digit in the last significant digit. I won't post the formula here because it's very difficult to edit it in this tiny 4" screen, but it's apparent from the listing. Windows 10 refuses to start on my desktop computer. I guess the "CPU Over Voltage Error!" I've been receiving when attempting to boot has something to do with that... Gerson 

01142016, 07:38 PM
(This post was last modified: 01152016 09:42 PM by Dieter.)
Post: #20




RE: Calculating e^x1 on classic HPs
(01142016 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 results: Code: seed = 0,815 n = 100000 Kahan = method posted by Pauli, probably due to W. Kahan alternate = method I suggested in this thread Plan B = hyperbolic method suggested by Gerson All results for random arguments in [ln 0.9, ln 2]. Edit: Table extended to show errors up to ±9 ULP. Larger errors were not detected. Dieter 

« Next Oldest  Next Newest »

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