If you're to implement something like this, you should first learn about how these things are done. I hope this doesn't sound too harsh. To explain it better, here is my variant of your code, with comparison to what C computes as e using expl(1)
:
#include <stdio.h>
#include <math.h>
int main ()
{
long double n = 0, f = 1;
int i;
for (i = 28; i >= 1; i--) {
f *= i; // f = 28*27*...*i = 28! / (i-1)!
n += f; // n = 28 + 28*27 + ... + 28! / (i-1)!
} // n = 28! * (1/0! + 1/1! + ... + 1/28!), f = 28!
n /= f;
printf("%.64llf\n", n);
printf("%.64llf\n", expl(1));
printf("%llg\n", n - expl(1));
printf("%d\n", n == expl(1));
}
Output:
2.7182818284590452354281681079939403389289509505033493041992187500
2.7182818284590452354281681079939403389289509505033493041992187500
0
1
There are two important changes to your code:
This code doesn't compute 1, 1*2, 1*2*3,... which is O(n^2), but computes 1*2*3*... in one pass (which is O(n)).
It starts from smaller numbers. Let us, for a moment, assume that your factorials are correct (see below). When you compute
1/1 + 1/2 + 1/6 + ... + 1/20!
and try to add it 1/21!, you are adding
1/21! = 1/51090942171709440000 = 2E-20,
to 2.something, which has no effect on the result (double holds about 16 significant digits). This effect is called underflow.
However, had you started with these numbers, i.e., if you computed 1/32!+1/31!+... they would all have some impact.
Notice that 32! needs 118 bits, i.e., 15 bytes. Your long type doesn't hold as much (the standard says at least 32 bits; it's not likely it will hold more than 64). If we add printf("%ld\n", factorial(i));
to your main for
-loop, we see the factorials:
1
2
6
24
120
720
5040
40320
362880
3628800
39916800
479001600
6227020800
87178291200
1307674368000
20922789888000
355687428096000
6402373705728000
121645100408832000
2432902008176640000
-4249290049419214848
-1250660718674968576
8128291617894825984
-7835185981329244160
7034535277573963776
-1569523520172457984
-5483646897237262336
-5968160532966932480
-7055958792655077376
-8764578968847253504
4999213071378415616
-6045878379276664832
See the negatives? That's when your (long) integer grew too much and restarted from its lowest possible value. Read about overflows; it's important to understand them.
By the way, I'm getting the same result (on my computer) even if I replace long
with long long
(and apply the proper printf
format %lld
).
Computing this kind of stuff is complex, and takes quite a bit of knowledge about how the numbers are stored in the computer, as well as the numerical mathematics and the mathematical analysis. I don't consider myself an expert, so so there might be more to this problem than what I've wrote. However, my solution seems in accordance to what C computes with its expl
function, on my 64bit machine, compiled with gcc 4.7.2 20120921.
double
's precision is silly as well. \$\endgroup\$ – CodesInChaos Oct 22 '13 at 13:35