38
\$\begingroup\$

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

\$\endgroup\$
  • 2
    \$\begingroup\$ I wouldn't say that the code is poor at all. I would probably mention the inconsistent braces, and you self asign result in your factorial method, and you don't in your main method, but I wouldn't say it was poor because of that. \$\endgroup\$ – Robert Snyder Oct 21 '13 at 19:34
  • 5
    \$\begingroup\$ @RobertSnyder The int-overflow is a critical flaw IMO. printing more digits than double's precision is silly as well. \$\endgroup\$ – CodesInChaos Oct 22 '13 at 13:35
46
\$\begingroup\$

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.

\$\endgroup\$
  • 1
    \$\begingroup\$ +1 Your putting the division by the common denominator 28! at the end is probably quite a performance improvement and allows for storing a the approximation as a fraction as well. \$\endgroup\$ – Tobias Kienzler Oct 22 '13 at 7:45
  • 3
    \$\begingroup\$ This algorithm is not mentioned in this similar SO question, maybe you want to post it there as well \$\endgroup\$ – Tobias Kienzler Oct 22 '13 at 7:56
  • \$\begingroup\$ Where did you get that expl code from? For comparison, here's GNU libc's 64-bit implementation (it's nuts!) and 32-bit implementation (relatively sane). There are also implementations for other architectures in the ieee754 directory. \$\endgroup\$ – 200_success Oct 22 '13 at 8:39
  • 2
    \$\begingroup\$ I didn't get the code. I merely called it, and compared the results (== should mean that all that bits of the two long double outputs are identical). My GNU libc is version 2.15-59.fc17.x86_64. This was intended as a simple test. True test would involve the estimate of the top bound of the error, but I felt that this would be beyond what was asked here. \$\endgroup\$ – Vedran Šego Oct 22 '13 at 13:13
  • \$\begingroup\$ Time complexity does not matter, please use factorials for better readibility. \$\endgroup\$ – Caridorc Aug 10 '15 at 14:00
8
\$\begingroup\$

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 × 1035, which is well beyond LONG_MAX (which may be as small as 231 - 1 ≈ 2 × 109).

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.

\$\endgroup\$
  • \$\begingroup\$ It looks like this is an HTML version containing the code in its original form: numbers.computation.free.fr/Constants/TinyPrograms/… \$\endgroup\$ – mwfearnley Oct 21 '13 at 20:59
  • 2
    \$\begingroup\$ 1+1/n doesn't seem to make any sense for integers, but it does: it's like n>1?1:2, but 2 characters shorter. \$\endgroup\$ – pts Oct 21 '13 at 21:07
  • 1
    \$\begingroup\$ If anyone can further de-golf the code or offer an explanation, please feel free to edit this answer! \$\endgroup\$ – 200_success Oct 22 '13 at 6:24
  • \$\begingroup\$ Those two whiles would probably read better as for loops, or at least the inner one: while (N-- > 9) { for (int n = N; n > 0; n--) { ... \$\endgroup\$ – Tobias Kienzler Oct 22 '13 at 7:37
  • 1
    \$\begingroup\$ I edited in your quote to stackoverflow.com/a/3028312/321973, since that link-only answer otherwise were poor following today's SE standards \$\endgroup\$ – Tobias Kienzler Oct 22 '13 at 8:09
4
\$\begingroup\$

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.

\$\endgroup\$
  • 1
    \$\begingroup\$ You're right that technically the problem is not underflow. However, lost significance isn't the issue either. As @VedranŠego points out, overflow of the long value leads to random wrong numbers being added. \$\endgroup\$ – 200_success Oct 22 '13 at 8:18
  • 1
    \$\begingroup\$ The question was "what does underflow mean"? \$\endgroup\$ – Dwayne Towell Oct 22 '13 at 14:57
  • \$\begingroup\$ @200_success Underflow (or loss of significance, whatever you want to call it) is a problem. Not as big as the overflow in the factorials, but it can contribute to the wrong solution. This is one of the reasons that I have reversed the order of the computation. See my point 2. \$\endgroup\$ – Vedran Šego Oct 22 '13 at 23:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.