Mathematicians has shown us that $e^x$ can be expressed as $\beta^{x/ln(\beta)}$, where $\beta$ is the number base of the computer, i.e. 2. This gives us that $e^x = 2^{x/ln(2)}$
It has also been shown that $e^x = \beta^{n_{1}}\beta^{1/2}\beta^{n_{2}}$ where $n_{1}$ is an integer, $-1/2 \le n_{2} < 1/2$, and $n_{1} + 1/2 + n_{2} = x/ln(\beta)$. So, $e^x = 2^{x/ln(2)} = 2^{n_{1}} \times 2^{1/2} \times 2^{n_{2}} = 2^{n_{1} + 1/2 + n_{2}}$
We will come back to this formula later on, but for now lets take a closer look at the expression $x/ln(2)$. Since $1/ln(2) = log_{2}(e)$, we can rewrite $x/ln(2)$ to be $x \times log_{2}(e)$. But $log_{2}(e)$ is a constant, $log_{2}(e) = 1.44269\ldots$, which means that we can create a new variable, let's call it $x_{0}$, which is equivalent to $x \times log_{2}(e)$. Thus, $x_{0} = x \times log_{2}(e)$, and $e^x = 2^{x_{0}}$.
In the code for exp() the variable $x_{0}$ is not created explicly, we simply multiply the variable x, which is the argument to the exp() function, with the constant log2e, and store the result of the multiplication back in x. The code for the assignment is simple:
x *= log2e;
From the exponent laws we know that $a^{x+y} = a^x \times a^y$. Using this law it is possible to form expressions such as $2^{6.2} = 2^6 \times 2^{0.2} = 73.51$. This means that we can split $2^{x_{0}}$ into two parts, one for the integer part of $x_{0}$ and one for the decimal part of $x_{0}$, let's name the parts $x_{1}$ and $x_{2}$. If we do the split we get $2^{x_{0}} = 2^{(x_{1} + x_{2})} = 2^{x_{1}} \times 2^{x_{2}}$.
To get the integral part of $x_{0}$ we use the floor() function, that is $x_{1}$ = floor($x_{0}$), and $x_{2} = x_{0} - x_{1}$. The code in exp() for this is shown below, remember that x in the code is equivalent to $x_{0}$.
x1 = floor(x);
At this moment we could compute the value of the integer part of
$2^{x_{1}} \times 2^{x_{2}}$, using the ldexp() function.
If we mix C-code and mathematics for a moment, we would get the following
expression
$e^{x} = 2^{xlog_{2}(e)} = 2^{x_{0}} = 2^{x_{1}} \times 2^{x_{2}} =$
ldexp(1, $x_{1}$) $\times 2^{x_{2}}$ =
ldexp(1, floor($x_{0}$)) $\times 2^{x_{2}}$
See ldexp.c for a detailed description of ldexp().
But, we will not do that, instead we will turn our attention to the decimal part, i.e. $2^{x_{2}}$. In the beginning we stated that $e^x = 2^{x/ln(2)} = 2^{xlog_2(e)} = 2^{n_{1}} \times 2^{1/2} \times 2^{n_{2}} = 2^{n_{1} + 1/2 + n_{2}}$.
If we replace $xlog_2(e)$ with $x_{0}$, we get $e^x = 2^{x/ln(2)} = 2^{x_{0}} = 2^{x_{1} + x_{2}} = 2^{n_{1}} \times 2^{1/2} \times 2^{n_{2}} = 2^{n_{1} + 1/2 + n_{2}}$.
Where $n_{1}$ is the integer part of $x_{0}$, and ${1/2} + {n_{2}}$ is the decimal part of $x_{0}$. If we remove the integer part from the $x_{0}$, then the decimal part can be expressed as $2^{x_{0} - x_{1}} = 2^{x_{2}} = 2^{1/2} \times 2^{n_{2}}$
This gives us that $n_{2} = x_{2} - 1/2$.
Let's rename $n_{2}$ in the last expression and give it the name $n$,
thus $n = x_{2} - 1/2$, or $2^{n} = {2^{x_{2}} \over 2^{1/2}}$.
The code in exp() to get the values for $x_{2}$ and $n$ is:
x2 = x - x1; n = x2 - 0.5;
Now that we have a value for $n$ we need to compute $2^{n}$, remember that $n$ is a decimal number, which means we cannot compute the value: for $2^{n}$ as we normally would if $n$ was an integer. So we use an arithmetic approximation for $2^n$ such that $2^{n} \approx {Q(n^2) + nP(n^2) \over Q(n^2) - nP(n^2)}$, where $0 \le n \le 1/2$.
Note that when we subtracted $1/2$ from $x_{2}$ to get the value of $n$, we also guaranteed that $0 \le n \le 1/2$. The next thing to do is to construct the two polynoms $P$ and $Q$.
$P = P_{2}(n^2)^2 + P_{1}(n^2)^1 + P_{0}(n^2)^0$
$Q = Q_{3}(n^2)^3 + Q_{2}(n^2)^2 + Q_{1}(n^2)^1 + Q_{0}(n^2)^0$
We get the coefficients for $P$ and $Q$ from table #1069 in Computer Approximations, Hart et al., 1968. Table #1069 also give us the order of the polynoms. Now that we have the polynoms we can compute the value for $2^n$, which is $2^{n} \approx {Q + n \times P \over Q - n \times P }$
In the code we use the macro Poly() to construct the polynoms:
nsquare = n * n; polyP = Poly(ORDERP, P, nsquare) * n; polyQ = Poly(ORDERQ, Q, nsquare);
See poly.c for a detailed description of Poly().
The resulting polynom that Poly() constructs looks like this:
polyP = (P[0] + nsquare * (P[1] + nsquare * P[2])) * n; polyQ = (Q[0] + nsquare * (Q[1] + nsquare * (Q[2] + nsquare * Q[3])));
Note that in the code we multiply polyP with n, i.e. $n \times P$, at the same time as we construct the polynom. Now that the polynoms are available to us in the code we can compute the value of $2^n$ simply by calculating (polyQ + polyP) / (polyQ - polyP). But, before we do that, remember what we have done so far. First we let $x_{0}$ denote the expression $xlog_2{e}$. Then we split $x_{0}$ into two parts, an integer part $x_{1}$, and a decimal part $x_{2}$. After that we saw that $2^{x_{1} + x_{2}} = 2^{n_{1}} \times 2^{n_{2}} \times 2^{1/2}$. We then removed the integer part which left us with $2^{x_2} = 2^{n_{2}} \times 2^{1/2}$. Then we calculated $n_{2}$ to be $x_{2} - 1/2$, and then we renamed $n_{2}$ to $n$. After that, we constructed the polynoms for the numerical approximation for $2^n$. So, now we have all the components to compute the value for the decimal part as follows: $2^{x_{2}} = 2^n + 2^{1/2} = 2^n \sqrt 2$.
If we add the integer part into the mix we get $e^x = 2^{xlog_2(e)} = 2^{x_{1} + x_{2}} = 2^{n_{1}} \times 2^{n_{2}} \times \sqrt 2 = 2^{n_{1}} \times 2^{n} \times \sqrt 2$, where $n_{1} = x_{0} - x_{2}$.
The expression $2^{n_{1}} \times 2^n \sqrt 2$ can be calculated with the ldexp() function. If we mix C-code and mathematics once again, we get the following expression: $e^x = 2^{n_{1}} \times 2^n \times \sqrt 2 = $ ldexp($2^n \times \sqrt 2$, $n_{1}$), but x1 $=$ floor(x) $ = n_{1}$ and $2^{n} = {Q(n^2) + nP(n^2) \over Q(n^2) - nP(n^2)}$ so, $e^x = $ ldexp(${Q(n^2) + nP(n^2) \over Q(n^2) - nP(n^2)} \times \sqrt 2$, x1)
The last line of exp() does the above calculation, and returns the result.
return ldexp((polyQ + polyP) / (polyQ - polyP) * sqrt2, x1);
/*
* exp.c
*
* Computes the exponential function of 'x'.
*/
#include "ieee754.h"
#include "math.h"
#define maxf 10000
#define log2e 1.4426950408889634073599247
#define sqrt2 1.4142135623730950488016887
static double P[] = {
0.2080384346694663001443843411e7,
0.3028697169744036299076048876e5,
0.6061485330061080841615584556e2
};
#define ORDERP ((int)(sizeof(P) / sizeof(P[0])) - 1)
static double Q[] = {
0.6002720360238832528230907598e7,
0.3277251518082914423057964422e6,
0.1749287689093076403844945335e4,
0.1e1
};
#define ORDERQ ((int)(sizeof(Q) / sizeof(Q[0])) - 1)
/*
* The exp() function computes the exponential function of 'x'.
* That is, e (2.71828...) raised to the 'x' power.
*
* Returns:
*
* The exp() function returns the exponential value.
*
* Internals:
*
* The coefficients for P and Q are from table #1069 in the book
* Computer Approximations, Hart et al. 1968.
*/
double exp(double x)
{
int x1; /* Integer part of x. */
double x2; /* Decimal part of x. */
double n, nsquare, polyP, polyQ;
if(x == 0){ return 1; }
if(x < -maxf){ return 0; }
if(x > maxf){ return inf(1); }
x *= log2e; /* Compute x0 */
x1 = floor(x);
x2 = x - x1;
n = x2 - 0.5;
nsquare = n * n;
polyP = Poly(ORDERP, P, nsquare) * n;
polyQ = Poly(ORDERQ, Q, nsquare);
return ldexp((polyQ + polyP) / (polyQ - polyP) * sqrt2, x1);
}