+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Calculators (and very old HP Computers) (/forum-3.html)
+--- Forum: General Forum (/forum-4.html)
+--- Thread: Calculating e^x-1 on classic HPs (/thread-5508.html)
Calculating e^x-1 on classic HPs - Dieter - 01-11-2016 10:20 PM
There is no doubt that the ln(1+x) and ex–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 ex–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)
u = Exp(x)
expm1 = (u - 1) + (x - Log(u)) * u
End FunctionOn a classic HP it could look like this:
Code:
ENTER
e^x
ENTER
ENTER
LN
CHS
R^
+
x
x<>y
1
-
+
RTNThis returns ex in Y and ex–1 in X.
What do you think?
Dieter
RE: Calculating e^x-1 on classic HPs - Thomas Klemm - 01-12-2016 01:41 AM
(01-11-2016 10:20 PM)Dieter Wrote: What do you think?
Brilliant!
Let's assume: \(u=e^x+\varepsilon\) where \(\varepsilon \ll 1\)
Then: \(e^x-1=u-1-\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^x-1\approx u-1+(x-\log(u))u\)
Cheers
Thomas
RE: Calculating e^x-1 on classic HPs - Claudio L. - 01-12-2016 05:13 AM
(01-11-2016 10:20 PM)Dieter Wrote: What do you think?
Dieter
I just checked at 2000 digits precision and the result was good.
I always thought e^x-1 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?
RE: Calculating e^x-1 on classic HPs - Dieter - 01-12-2016 02:02 PM
(01-12-2016 05:13 AM)Claudio L. Wrote: I just checked at 2000 digits precision and the result was good.
Fine. My results on 10- and 12-digit calculators did not show any results with larger errors than ±1 ULP, mostly because of different rounding. The vast majoriy was dead on.
(01-12-2016 05:13 AM)Claudio L. Wrote: I always thought e^x-1 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 ex–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.
(01-12-2016 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 15-16 digit accuracy a continued fraction method with 7 terms or a 16-term 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 hard-coded seven terms (k = 13, 11, 9, 7, 5, 3, 1):
Code:
Function expm1cf(x)
If (x < -0.1053605156578263) Or (x > 0.6931471805599453) Then
expm1cf = Exp(x) - 1
Else
s = 0
For k = 13 To 1 Step -2
s = x / (k - x / (s + 2))
Next k
expm1cf = s
End If
End FunctionDieter
RE: Calculating e^x-1 on classic HPs - ElectroDuende - 01-12-2016 03:34 PM
(01-12-2016 05:13 AM)Claudio L. Wrote: Question is: would a program computing terms of the power series in a loop be faster?
I suppose that the main reason not to use the program-loop approach is that not all the "classic" models are programmable...
RE: Calculating e^x-1 on classic HPs - Dieter - 01-12-2016 06:12 PM
(01-12-2016 01:41 AM)Thomas Klemm Wrote: Brilliant!
[...proof...]
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 E-6 and a 10-digit calculator.
Evaluating u = ex with 10 digit precision yields
u = 1,000003142
Accordingly, ex–1 is evaluated as
u-1 = 0,000003142
But that's not ex or ex-1. The correct values are
ex = 1,0000031415975888...
ex–1 = 0,0000031415975888...
The 10-digit result we got actually is not ex but eln u = e3,141995064E-6 (resp. this minus one). Simply because only eln u exactly equals u:
eln u = u = 1,000003142000000...
eln 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:
ex-1 ≈ e^(ln u) - 1 + (x - ln u) · d[ex–1]/dx
= (u - 1) + (x - ln u) · ex
...and since u agrees with ex to machine accuracy we finally get:
ex-1 ≈ (u - 1) + (x - ln u) · u
That's it.
Dieter
RE: Calculating e^x-1 on classic HPs - Dieter - 01-12-2016 06:30 PM
(01-12-2016 03:34 PM)ElectroDuende Wrote:(01-12-2016 05:13 AM)Claudio L. Wrote: Question is: would a program computing terms of the power series in a loop be faster?
I suppose that the main reason not to use the program-loop approach is that not all the "classic" models are programmable...
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 10-digit 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
RE: Calculating e^x-1 on classic HPs - Dieter - 01-12-2016 08:16 PM
(01-11-2016 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 ex and ex–1 in Y and X.
The following version returns ex–1 while it saves the Y-register. It also does not use R^ which is missing on some of the classic calculators.
Code:
e^x
ENTER
ENTER
LastX
x<>y
LN
-
x<>y
x
LastX
1
-
+
RTNAnd 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))*UDieter
RE: Calculating e^x-1 on classic HPs - Claudio L. - 01-12-2016 09:44 PM
(01-12-2016 06:30 PM)Dieter Wrote:(01-12-2016 03:34 PM)ElectroDuende Wrote: I suppose that the main reason not to use the program-loop approach is that not all the "classic" models are programmable...
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 10-digit 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
Great, now you have a very elegant solution, with hard proof that it's also speed and memory efficient. Good work.
RE: Calculating e^x-1 on classic HPs - Gerson W. Barbosa - 01-13-2016 03:46 PM
(01-12-2016 08:16 PM)Dieter Wrote: The following version returns ex–1 while it saves the Y-register. It also does not use R^ which is missing on some of the classic calculators.
Code:
e^x
ENTER
ENTER
LastX
x<>y
LN
-
x<>y
x
LastX
1
-
+
RTN
Very nice, thank you!
Just tested it on the HP-12C. I've tried an hyperbolic function approach on the HP-15C using the identity
\[e^{x}-1=\frac{2}{\coth \frac{x}{2}-1}\]
Code:
001-42.21.11 f LBL A
002- 2 2
003- 10 /
004-42.22.25 f HYP TAN
005- 15 1/x
006- 1 1
007- 30 -
008- 2 2
009- 34 x<>y
010- 10 /
011- 43 32 g RTNHowever 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.
RE: Calculating e^x-1 on classic HPs - Dieter - 01-13-2016 06:43 PM
(01-13-2016 03:46 PM)Gerson W. Barbosa Wrote: Just tested it on the HP-12C. I've tried an hyperbolic function approach on the HP-15C using the identity
\[e^{x}-1=\frac{2}{\coth \frac{x}{2}-1}\]
Ah, that's the inverse of the method I suggested for ln(1+x) in 2014, cf. this thread.
(01-13-2016 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 ex–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
RE: Calculating e^x-1 on classic HPs - Gerson W. Barbosa - 01-14-2016 02:47 AM
(01-13-2016 06:43 PM)Dieter Wrote: On the other hand a dedicated ex–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 W|A 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 10-digit calculators in the specified range (five terms for 12-digit 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 HP-49G/G+/50g version:
Code:
« 0 → x s
« 1 5
FOR i i DUP + DUP x SWAP ^ OVER IBERNOULLI * OVER ! / SWAP / 's' STO+
NEXT
s x 2 / + x LN + EXP
»
»2 LN EXM1 -> 1.
1.23E-8 EXM1 -> 1.23000000759E-8
-0.1 EXM1 RE -> -9.51625819644E-2
-0.1 EXPM -> -.095162581964 ; built-in EXPM function
1.23E-8 EXPM -> 1.23000000756E-8
Best regards,
Gerson.
Edited per Dieter's observation below (1.23E-8 EXPM -> 1.23000000756E-8, not 1.2300000076E-8)
RE: Calculating e^x-1 on classic HPs - Paul Dale - 01-14-2016 07:20 AM
A nice identity using hyperbolics is: ex-1 = tanh(x/2) (ex+1)
The 34S uses:
Code:
u = e^x
v = u - 1
if v = 0, return x
if v = -1, return -1
return v * x / ln(u)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
RE: Calculating e^x-1 on classic HPs - Dieter - 01-14-2016 01:44 PM
(01-14-2016 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 W|A missed it...
Bernoulli numbers? Waaaayyyy to complicated. ;-)
(01-14-2016 02:47 AM)Gerson W. Barbosa Wrote: Only four terms of the summand would suffice for 10-digit calculators in the specified range (five terms for 12-digit 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.
(01-14-2016 02:47 AM)Gerson W. Barbosa Wrote: 2 LN EXM1 -> 1.
1.23E-8 EXM1 -> 1.23000000759E-8
-0.1 EXM1 RE -> -9.51625819644E-2
-0.1 EXPM -> -.095162581964 ; built-in EXPM function
1.23E-8 EXPM -> 1.2300000076E-8
For the record: the exact 12-digit results are
1
1,23000000756 E-8
-9,51625819640 E-2
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 built-in EXPM function really return 1.2300000076 E-8, i.e. +4 ULP compared to the true result?
Dieter
RE: Calculating e^x-1 on classic HPs - Dieter - 01-14-2016 01:48 PM
(01-14-2016 07:20 AM)Paul Dale Wrote:Code:
u = e^x
v = u - 1
if v = 0, return x
if v = -1, return -1
return v * x / ln(u)
As William Kahan said of this: nearly full working relative accuracy despite cancellation no matter how tiny x may be.
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 vI 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
RE: Calculating e^x-1 on classic HPs - Gerson W. Barbosa - 01-14-2016 03:30 PM
(01-14-2016 01:44 PM)Dieter Wrote:(01-14-2016 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 W|A missed it...
Bernoulli numbers? Waaaayyyy to complicated. ;-)
You are quite right! Anyway, that was not supposed to be a valid method. I just liked the series, despite its uselessness :-)
(01-14-2016 01:44 PM)Dieter Wrote:(01-14-2016 02:47 AM)Gerson W. Barbosa Wrote: 2 LN EXM1 -> 1.
1.23E-8 EXM1 -> 1.23000000759E-8
-0.1 EXM1 RE -> -9.51625819644E-2
-0.1 EXPM -> -.095162581964 ; built-in EXPM function
1.23E-8 EXPM -> 1.2300000076E-8
For the record: the exact 12-digit results are
1
1,23000000756 E-8
-9,51625819640 E-2
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 built-in EXPM function really return 1.2300000076 E-8, i.e. +4 ULP compared to the true result?
I'd made a mistake. Fixed. Thanks!
Gerson.
RE: Calculating e^x-1 on classic HPs - emece67 - 01-14-2016 03:55 PM
(01-14-2016 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?
RE: Calculating e^x-1 on classic HPs - Dieter - 01-14-2016 06:42 PM
(01-14-2016 03:55 PM)emece67 Wrote:(01-14-2016 01:48 PM)Dieter Wrote: I did a few test runs on a 34s in SP mode with 100000 random numbers each. (...)
I'm really curious, how do you do such tests?
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. ;-)
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
ULP Kahan alternate ULP Kahan alternate
----------------------- -----------------------
≤-5 54 4 ≤-5 57 2
-4 58 86 -4 60 101
-3 316 80 -3 331 58
-2 4568 77 -2 4580 74
-1 22176 16796 -1 22325 16788
±0 45522 66011 ±0 45193 65769
+1 22215 16691 +1 22323 16935
+2 4617 71 +2 4667 69
+3 333 70 +3 326 78
+4 71 113 +4 64 126
≥+5 70 1 ≥+5 74 0Re. the meaning of the two seed values, cf. this thread. ;-)
Dieter
RE: Calculating e^x-1 on classic HPs - Gerson W. Barbosa - 01-14-2016 06:50 PM
(01-13-2016 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.
On the other hand a dedicated ex–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...
Here is my plan B for the hyperbolic function approach:
Code:
001- LBL B
002- e^x
003- LSTx
004- 2
005- /
006- HYP TAN
007- *
008- LSTx
009- +
010- RTNNo additional tests required. On the HP-15C I get the same results of the HP-12C 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
RE: Calculating e^x-1 on classic HPs - Dieter - 01-14-2016 07:38 PM
(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 results:
Code:
seed = 0,815 n = 100000
ULP Kahan alternate plan B
----------------------------------
-9 2 - -
-8 4 - -
-7 7 - -
-6 12 - -
-5 32 2 2
-4 60 101 229
-3 331 58 1852
-2 4580 74 8025
-1 22325 16788 20825
±0 45193 65769 38405
+1 22323 16935 20652
+2 4667 69 7973
+3 326 78 1792
+4 64 126 241
+5 34 - 4
+6 16 - -
+7 12 - -
+8 9 - -
+9 3 - -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