Why is my C program for calculating Euler's constant of poor quality?

I was trying to post this code to a Wikipedia article, but it was soon removed. I then asked about the code in the page's talk section, and some other contributors said that is was "very poor" and that "there are clearly going to be underflows all over the place." What does that mean? What is an "underflow"?

Could you please tell me more about my program's flaws and how it can be improved? Please give me some constructive criticism. Here is where I originally posted the code.

#include 

long factorial(int n) {
    long result = 1;
    for (int i = 1; i <= n; ++i)
    result *= i;
    return result;
}

int main ()
{
    double n=0;
    int i;
    for (i=0; i<=32; i++) {
        n=(1.0/factorial(i))+n;
    }
    printf("%.32f\n", n);
}

Here is the result of the program 2.71828182845904553488480814849027

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 
#include 

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, 12, 123,... which is O(n^2), but computes 123... 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.

StackExchange Code Review Q#33015, answer score: 47