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

Question

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 <stdio.h>

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

Answer 1

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:

1. 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)).

2. 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.

Answer 2

Your algorithm uses intermediate numbers that are very large and very small, and they are therefore difficult for computers to work with. For example, 32! ≈ 2.6313 × 10³⁵, which is well beyond LONG_MAX (which may be as small as 2³¹ - 1 ≈ 2 × 10⁹).

Basically, once n is large enough (maybe 13, maybe 21), factorial(n) returns a "random" number. Since your result is long instead of unsigned long, most of the overflowed results will be very large positive or very large negative numbers, so their reciprocals will generally average to zero — but there's no guarantee of that.

So, what would be a good algorithm? You want something that avoids using very small or very large intermediate values. This question was posed as Stack Overflow question 3028312; this answer links to a paper describing many algorithms for computing e. Section 9 of the paper has an efficient algorithm (though it is code-golfed / mysterious):

This is a tiny C program from Xavier Gourdon to compute 9000 decimal digits of e on your computer. A program of the same kind exists for π and for some other constants defined by mean of hypergeometric series.

#include <stdio.h>
#define DIGITS 9000 /* decimal places (not including the '2') */
int main() {
    int N = DIGITS+9, a[DIGITS+9], x = 0;
    a[0] = 0;
    a[1] = 2;
    for (int n = 2; n < N; ++n) {
        a[n] = 1;
    }
    for ( ; N > 9; --N) {
        for (int n = N - 1; n > 0; --n) {
            a[n] = x % n;
            x = 10 * a[n-1] + x/n;
        }
        printf("%d", x);
    }
    return 0;
}

This program [when code-golfed] has 117 characters. It can be changed to compute more digits (change the value of DIGITS) and to be faster (change the constant 10 to another power of 10 and adjust the printf command). A not so obvious question is to find the algorithm used.

Answer 3

Technically the problem is not one of underflow, but of lost significance. Underflow is when the number is too small to be stored in the current format and it is replaced by 0.0. Lose of significance occurs when a (relatively) large number is added to a (relatively) small number. Recall that when fractional numbers are added the "decimal" (or binary) point must be lined up. For example in decimal we might add 2.5x10^0 and 2.5x10^-5, which expanded to line up the decimals looks like this...

 2.500000
+0.000025
----------
 2.500025

If we perform this calculation with only 4 places of accuracy the answer is simply 2.500 rather than 2.500025, which is what we would expect. The same thing happens with your program but of course you have more digits of accuracy available and of course it occurs in binary, not decimal, but the problem is the same.