1

For fun, I was trying to evaluate the Gaussian integral from 0 to 1 using a series expansion. For this reason, I wrote a factorial function which works well up to 20!(I checked) and then I wrote this:

int main(){
    int n;
    long double result=0;
    for(n=0; n<=5; n++){
        if(n%2==0){
            result+=(((long double) 1/(long double)(factorial(n)*(2*n+1))));
        } else {
            result-=(((long double) 1/(long double)(factorial(n)*(2*n+1))));
        }
    }        
    printf("The Gaussian integral from 0 to 1 is %Lf\n", result);
}

This gives me a strange negative number which is obviously not even close. I suspect the problem is with the cast, but I don't know what it is. Any thoughts? This is not the first thing I tried. I tried converting anything in the expression and putting the explicit cast at the beginning, but it didn't work.

Gennaro Marco Devincenzis
  • 141
  • 1
  • 5

3 Answers3

3

You are using the MinGW compiler (port of gcc for Windows), which has issues with the long double type. This is due to conflicts between GCC's implementation of long double and Microsoft's C library. See also this question.

According to this question, defining __USE_MINGW_ANSI_STDIO may solve this. If not, using double instead will work.

Community
  • 1
  • 1
interjay
  • 107,303
  • 21
  • 270
  • 254
2

In (long double)(factorial(n)*(2*n+1), the multiplications are integer multiplications and the first one could overflow if the result of factorial is already close to the limit of the integer type used.

Write ((long double)(factorial(n))*(2*n+1) so that the first multiplication is a floating-point multiplication.

Pascal Cuoq
  • 79,187
  • 7
  • 161
  • 281
  • `n<=5` here, so it shouldn't overflow. – interjay Aug 28 '14 at 14:41
  • @interjay Yes. I have to assume that the code shown in the question is not exactly the one that overflows. However, assuming a factorial function that works for the values passed to it and returns some integer type (not specified in the question), the integer multiplication by `2*n+1` is the obvious red flag in the question, and converting immediately to floating-point is the simplest remedy. – Pascal Cuoq Aug 28 '14 at 14:44
  • 1
    OP claims "factorial function which works well up to 20!". I doubt that as it would overflow even a 64 bit integer. So I agree with Pascal - this is not an SSCCE. So +1 answer (even if mine is cuter). – Bathsheba Aug 28 '14 at 14:45
  • @Bathsheba 20! ~ 2.432902e+18 is close to the limit of 64-bit signed integers. But if `factorial` returns a 64-bit signed integer, why would `factorial(5)` cause any trouble? Hmmm... – Pascal Cuoq Aug 28 '14 at 14:47
  • Regarding the 20! thing: it's a naive algoritm, designed for not many digits of precision. However, I checked on Wolfram and 20! is correct. – Gennaro Marco Devincenzis Aug 28 '14 at 14:48
2

You're almost certainly overflowing your integer type. In C this is technically undefined behaviour.

For 32 bit unsigned integer, 13! will overflow. On 64 bit, 21! will overflow.

Your algorithm will survive a little longer if you use a floating point double type or an extension like __uint128 (gives you, I think, up to 34!) if your compiler supports it.

Another problem that you have is that you are progressively adding terms of decreasing size to your total. That's never a good idea when working with floating point types. If you run your for loop in the reverse order then the result will be more accurate.

Bathsheba
  • 231,907
  • 34
  • 361
  • 483