Choosing the best algorithm for approximating $e^{-x}$

Question

I am trying to learn numerical methods and understand why some algorithms are more efficient than others. My question os about the approximation of $\exp{-x}$. There are (at least) two possibilities. The first would be thinking of $e^{-x}$ as $$ \sum_{k=0}^\infty \frac{(-1)^k x^k}{k!} $$ and truncate this series to obtain the approximation.

The second option is think of $e^{-x}$ as

$$ \frac{1}{\sum_{k=0}^\infty \frac{x^k}{k!}} $$

and truncate the series.

I programmed both algorithms to compare them. I try with different values of $x$, seems that for small values ($x=0.5$) the first algorithm is better, while for bigger values ($x = 5$), the second is better than the first.

Could you help me to understand why is this?

Answer 1

There are probably multiple good answers for this question, but one reason is the relative computational error. If $x$ is large, then $\sum\limits_{k=0}^\infty \frac{(-1)^k x^k}{k!}$ adds and substracts large positive numbers at least in the beginning. This is bad as can be seen from the following toy example: If $a = 100$ and $b = 101$, then $b-a = 1$. A relative error of $1\%$ in $a$ or $b$ could give you up to $100\%$ error in the difference $b-a$. Addition (of positive numbers) does not have this problem; an error of $1\%$ in $b$ or $a$ will not give you more than $1\%$ error in the sum. Thus the alternating series should be avoided for large $x$. For small $x$ the terms get small very quickly so that you do not get this problem (as much).

Answer 2

Ok, so let's go down some rabbit holes that are actually facing:

Choosing the best algorithm for approximating $e^{−x}$

Spoiler: None of the real-world algorithms presented below avoids negative arguments or takes special measures depending on the sign. Reason is that the input domain is reduced to $|x|\leqslant(\ln2) / 2^{1+K}$ for some $K\in\Bbb N_0$, and one approach avoids Taylor altogether.

The first step in finding the "best" algorithm is of course to assess the requirements and what "good" means.

0. The Requirements

1. What is the maximal error you are willing to tolerate? Is it an absolute error or relative error? Relative errors are usually considered in the realm of floating-point algorithms / libraries like the IEEE-754 examples that will follow in the remainder.

2. Relates to the error is the required precision. It the precision known at design-time (like with the Newlib example below) or only known at run-time (like with the MPFR example below).

3. What are the constraints on run-time, code-memory and data-memory. Are they known at design time (which might be the case in deeply embedded systems)? How expensive are the basic arithmetic operations like addition, division, etc.?

In the remainder, I chose approaches from the floating-point + relative error + open source variety. One where the precision is known at design time, and one with arbitrary precision. This is just some random choice, but it shows the complexity of real world code.

1. Floating-point with known Precision (Newlib)

Take Newlib,a a C/C++ library, as an example. It implements$^1$ exp in newlib/libm/mathfp/s_exp.c.

In order to compute $\exp(x)$:

1. Compute $g = x\bmod (\ln 2)$ in the smallest remainder representation, i.e. $g = x-N\ln 2$ with $N\in\Bbb Z$ and $|g|\leqslant \frac12\ln2$.

2. Compute $R(g) = \exp(g) = 0.5 + P / (Q - P)$ where $P$ is an odd polynomial of degree 5 in $g$, and $Q$ is a polynomial of degree 3 in $g$.

3. Return $\exp(x) = 2^N \exp(g)$. As IEEE-754 is a binary floating-point format, the multiplication with $2^N$ is just an integer addition in the exponent; the mantissa is unaltered.

Bottom line:

• There is no "high level" computation involved: No power series or such, the essential step 2 is just evaluation of some polynomials with coefficients that appear to be arbitrary: The art lies in the determination of these coefficients$^2$; their theory and the effort to get them does not enter the very algorithm — except for the effectiveness of the outcome. These polynomial are special-purpose taylored to be optimal for the domain under consideration, $|g|\leqslant \frac12\ln2$ in the present case.

• In any case, whenever possible reduce the problem to a small, manageable interval where the function is well-behaved. The "reduce the problem" pattern is essential for an affactive approach, and you will see it it practical implementations for many other functions like: triginometric functions (exploit periodicity and symmetries), log and inverse triconometric functions (exploit functional equations like $\log (2^Nx) = 2^N + \log x$ or $\arctan(1/x)=\pi/2-\arctan x$, etc).

2. Floating-point with arbitrary Precision (MPFR)

MPFR is an arbitrary-precision$^3$ library written in C with bindings for many other languages. The entry in exp.c just deals with gory details and then dispatches to exp_2.c where it uses Brent's formula $$ \exp(x) = 2^n\cdot\left(1+r+\frac{r^2}{2!}+\frac{r^3}{3!}+\cdots\right)^{2^K} \quad\text{ where }\quad x = n\ln 2 + 2^Kr\tag 1$$

with $n\in\Bbb Z$, $K\in\Bbb N_0$, together with the Paterson-Stockmeyer ${\cal O}(t^{1/2})$ algorithm for the evaluation of power series. The series is then computed for the reduced value of

$$r = (x-n\ln 2)\cdot2^{-K}\tag 2$$

The additional parameter $K$ allows to dynamically adjust the domain for which the power series has to be evaluated, at the cost of $K$ additional squarings and increased internal precision.

In a different, unlikely route the main module uses a different, binary splitting of the argument$^4$ approach implemented in exp3.c.

Bottom line:

• Again, the major theme is to reduce the domain over which the actual "main" computation takes place. A considerable part of the code deals with working out what domain is optimal, which internal precision is needed, what the benefits and costs are, etc.

• An actual, real-world implementation of $\exp$ (or any other function for that matter) is likely to be much more complex and sophisticated than one might expect.

---

Notes

$^1$This is just the default, fallback implementation. Depending on the target hardware, other approaches might be better. If the hardware has an IEEE-754 compliant FPU that supports exp, then it's usually much preferred to use respective instructions.

$^2$The Newlib source refers to "Software Manual for the Elementary Functions" by William J. Cody, Jr. and William Waite, Prentice Hall, 1980. I don't have access to that source, and it wold be interesting to know theory behind the magic polynomials. Usually, one would use best (in terms of $\max$-norm) polynomials or rational functions for the domain in question, biu in the present case there is more about it.

$^3$"Arbitrary precision" means that the user of the library can specify the precision at run-time and change it as needed. Like with IEEE754, the semantics is that the result must be the same as if calculated with infinte precision, and then rounded according to the specified rounding mode (to nearest, to zero, to +∞, to −∞ etc.) for the precision active at that moment.

$^4$With that approach I am not familiar. Perhaps someone can explain the comments in that module or the general idea.