Implementation of exp function in C using Taylor Series expansion

Question

I am trying to write code to calcultate e^x using:

$$e^x = \sum_{n=0}^\infty \frac{x^n}{n!} = 1 + x + \frac{x^2}{2!} + \frac{x^3}{3!} + \frac{x^4}{4!} + \cdots $$

This is the code I have, which works for small values of x. Please do not pay attention to the variable e, specifically to the fact that user enters it.

#include <stdio.h>
#include <math.h>

int main(void) {
    double x0,x1,e,x,sum=0.0;
    int n=0;
    scanf("%lf %lf",&x,&e);
    x1=1;
    do{
        sum+=x1;
        x0=x1;
        x1=(x/++n)*x0;
    }
    while(x1/x0>e);

    printf("sum=%lf, \nexp=%lf",sum,exp(x));
    return 0;
}

For input 50 0.00001 number of decimal places for e doesn't matter. it generates:

sum=5184705528587076239360.000000, 
exp=5184705528587072045056.000000

For input 21 0.00000001:

sum=1318815734.483214, 
exp=1318815734.483215

The difference only increases as x increases. I assume that condition is wrong but see no other way of doing it.

Answer 1

Termination condition not correct

To get the precision you request, you should stop the loop when x1 <= e, not when x1 / x0 <= e. If you do that, you can now get rid of x0 because it no longer serves a purpose. I adjusted your program to use the new termination condition and it gave me the results I believe you are looking for. I also renamed e to precision because it was confusing for a variable to be named e when it didn't actually hold the value of the constant \$e\$.

#include <stdio.h>
#include <math.h>

int main(void) {
    double x1,precision,x,sum=0.0;
    int n=0;
    scanf("%lf %lf",&x,&precision);

    // Sums: 1 + x + x^2/2! + x^3/3! + x^4/4! + x^5/5! + ...
    //
    // To get an approximation of e^x.
    //
    // Stops when the next term to add becomes smaller than precision.
    x1=1;
    do {
        sum += x1;
        x1  *= (x/++n);
    } while (x1 > precision);

    printf("sum=%lf,\nexp=%lf\n",sum,exp(x));
    return 0;
}

Answer 2

Short answer: Your computation is numerically unsound.

Long answer:

The termination condition looks funny indeed. From a pure mathematical standpoint x1/x0 is equal to x/n, so the condition actually means x/n > e, and one would expect the loop to run until n = x/e, in your case to 5000000. However, the unfriendly nature of floating point strikes again, and the loop breaks much earlier, at n = 539. Upon loop termination the debugger shows that x1 is literally 0. For the record, x0 at this point is 4.94E-324, and on the next iteration double ran out of mantissa bits.

It may look inconceivable that numbers as small as 4.94E-324 may significantly affect the result. In fact your problems start much earlier than that.

Look how terms of the sequence behave: they grow while n is less then x, and start decreasing after that. The maximum is obviously achieved at \$n \approx x\$, and is equal to \$\frac{x^x}{x!} \approx \frac{x^x}{\sqrt{2\pi x}}\$ (by Stirling approximation). That is, for \$x = 50\$ the 50'th term is just 1/17'th of the result, and the accumulated sum is already in order of E21 of magnitude. The terms meanwhile start to decrease, and you find yourself in a very unfortunate situation:

you add small values to large ones.

To add two floating point values, the summator must align them to the same exponent. This means that the smaller value inevitably loses some bits of mantissa. There goes precision, and what's even worse, the whole mantissa may disappear: that is from the certain n terms stop to contribute whatsoever. For x = 50 underflow happen at n = 119, and the first completely lost term is about 269909. Considering further lost terms, that it is quite in line with the 9437184 error you've got.

Solution:

The huge intermediate result defeats small terms one by one. Together they can defend. Start from the small numbers (that is large n) and work your way down to n = 0. The Horner schema was specifically designed to deal with issues you are facing.