# Implementation of exp function in C using Taylor Series expansion

I am trying to write code to calcultate ex 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.

• Welcome to Code Review! Does your code do what you would like? I'm a little confused by your sample inputs. Mar 26 '16 at 21:48
• @SirPython In output of my program sum= is what my code generates, exp= is what math.h function exp() returns. My code works fine for values of x less than 21. After that difference between these two functions increases as x increases Mar 26 '16 at 21:51
• @SirPython In input: first number is x, the exponent to which e (2.718) is raised, second number is precision and is used for do-while condition, to end calculations. Mar 26 '16 at 21:59

### 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;
}

• I have tested it now on ideone and get the same answers as before. If I change type 'double' to 'long double' I get different answer than before, but just as wrong. Mar 27 '16 at 8:17
• @BojanPavlović Must be rounding errors. On my x86 based machine, I get the correct answer.
– JS1
Mar 27 '16 at 8:40

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.

• Right now I don't have time to look this with more attention, but I am going to look this when I finish with exams. Mar 28 '16 at 13:36