Calculating e^x-1 on classic HPs

[Dieter]

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)
   u = Exp(x)
   expm1 = (u - 1) + (x - Log(u)) * u
End Function

On a classic HP it could look like this:

Code:

ENTER
e^x
ENTER
ENTER
LN
CHS
R^
+
x
x<>y
1
-
+
RTN

This returns e^x in Y and e^x–1 in X.
What do you think?

Dieter

[Thomas Klemm]

(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

[Claudio L.]

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

[Dieter]

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

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

Dieter

[ElectroDuende]

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

[Dieter]

(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 = e^x with 10 digit precision yields
  u = 1,000003142

Accordingly, e^x–1 is evaluated as
u-1 = 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 10-digit result we got actually is not e^x but e^ln u = e^{3,141995064E-6} (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

[Dieter]

(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

[Dieter]

(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 e^x and e^x–1 in Y and X.

The following version returns e^x–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

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

[Claudio L.]

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

[Gerson W. Barbosa]

(01-12-2016 08:16 PM)Dieter Wrote:  The following version returns e^x–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 RTN

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.

[Dieter]

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

[Gerson W. Barbosa]

(01-13-2016 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 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)

[Paul Dale]

A nice identity using hyperbolics is: e^x-1 = tanh(x/2) (e^x+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

[Dieter]

(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

[Dieter]

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

This may lead to an error in the last place for arguments where e^x = 10...11, 100...101, 1000...1001 etc. Here the last digit is lost. But this may also occur with the unchanged method:

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

[Gerson W. Barbosa]

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

[emece67]

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

[Dieter]

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

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           0

Re. the meaning of the two seed values, cf. this thread. ;-)

Dieter

[Gerson W. Barbosa]

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

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- RTN

No 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

[Dieter]

(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