e^x without math.h

Question

I'm trying to find e^x without using math.h. My code gives wrong anwsers when x is bigger or lower than ~±20. I tried to change all double types to long double types, but it gave some trash on input.

My code is:

#include <stdio.h>

double fabs1(double x) {
    if(x >= 0){
        return x;
    } else {
        return x*(-1);
    }
}

double powerex(double x) {
    double a = 1.0, e = a;
    for (int n = 1; fabs1(a) > 0.001; ++n) {
        a = a * x / n;
        e += a;
    }
    return e;
}

int main(){
    freopen("input.txt", "r", stdin);
    freopen("output.txt", "w", stdout);
    int n;
    scanf("%d", &n);
    for(int i = 0; i<n; i++) {
        double number;
        scanf("%lf", &number);
        double e = powerex(number);
        printf("%0.15g\n", e);
    }
    return 0;
}

Input:

8
0.0
1.0
-1.0
2.0
-2.0
100.0
-100.0
0.189376476361643

My output:

1
2.71825396825397
0.367857142857143
7.38899470899471
0.135379188712522
2.68811714181613e+043
-2.91375564689153e+025
1.20849374134639

Right output:

1
2.71828182845905
0.367879441171442
7.38905609893065
0.135335283236613
2.68811714181614e+43
3.72007597602084e-44
1.20849583696666

You can see that my answer for e⁻¹⁰⁰ is absolutely incorrect. Why does my code output this? What can I do to improve this algorithm?

Answer 1

When x is negative, the sign of each term alternates. This means each successive sum switches widely in value rather than increasing more gradually when a positive power is used. This means that the loss in precision with successive terms has a large effect on the result.

To handle this, check the sign of x at the start. If it is negative, switch the sign of x to perform the calculation, then when you reach the end of the loop invert the result.

Also, you can reduce the number of iterations by using the following counterintuitive condtion:

e != e + a

On its face, it appears that this should always be true. However, the condition becomes false when the value of a is outside of the precision of the value of e, in which case adding a to e doesn't change the value of e.

double powerex(double x) {
    double a = 1.0, e = a;
    int invert = x<0;
    x = fabs1(x);
    for (int n = 1; e != e + a ; ++n) {
        a = a * x / n;
        e += a;
    }
    return invert ? 1/e : e;
}

We can optimize a bit more to remove one loop iteration by initializing e with 0 instead of a, and calculating the next term at the bottom of the loop instead of the top:

double powerex(double x) {
    double a = 1.0, e = 0;
    int invert = x<0;
    x = fabs1(x);
    for (int n = 1; e != e + a ; ++n) {
        e += a;
        a = a * x / n;
    }
    return invert ? 1/e : e;
}

Answer 2

For values of x above one or so, you may consider to handle the integer part separately and compute powers of e by squarings. (E.g. e^9 = ((e²)²)².e takes 4 multiplies)

Indeed, the general term of the Taylor development, x^n/n! only starts to decrease after n>x (you multiply each time by x/k), so the summation takes at least x terms. On another hand, e^n can be computed in at most 2lg(n) multiplies, which is more efficient and more accurate.

So I would advise

• to take the fractional part of x and use Taylor,
• when the integer part is positive, multiply by e raised to that power,
• when the integer part is zero, you are done,
• when the integer part is negative, divide by e raised to that power.

---

You can even spare more by considering quarters: in the worst case (x=1), Taylor requires 18 terms before the last one becomes negligible. If you consider subtracting from x the immediately inferior multiple of 1/4 (and compensate multiplying by precomputed powers of e), the number of terms drops to 12.

E.g. e^0.8 = e^(3/4+0.05) = 2.1170000166126747 . e^0.05