Ytukyg

I try to write a factorial function to compute a large number(factorial(105)),its result have 168 digit, so use long double, but it seems to be a error, can't it use like this?

#include <stdio.h>



long double factorial(long double n,long double base = 1){

  if (n <= 1){

    return 1*base;

  }else{

    return factorial(n-1,n * base);

  }

}  

int main(){  

    printf("%.0Lfn",factorial(25));     // this result is correct



    printf("%.0Lfn",factorial(26));    

    //correct result is 403291461126605635584000000,but it return 403291461126605635592388608

    return 0;

}

Back of the envelope calculation: 25! is slightly more than 10²⁵; three orders of magnitude is roughly 10 bits, so you would need roughly 83 bits of mantissa even just to represent precisely the result.

Given that a long double, on platforms that support it, is generally 80 bits for the whole value (not just the mantissa!), apparently there's no way you have enough mantissa to perform that calculations of that order of magnitude with integer precision.

However: factorial here is a bit magic, as many of the factors contain powers of two, which just add binary zeroes to the right, which do not require mantissa (they end up in the exponent). In particular:

25! = 2   4   2   8   2    4    2    16    2    4     2    8    = 2²² · m

        3   5 3 7   9 5 11 3 13 7 15    17 9 19 5 21 11 23 3 25 

(m being the product of all non-2 factors, namely m = 3¹⁰ · 5⁶ · 7³ · 11² · 13 · 17 · 19 · 23, so effectively the data we have to store in the mantissa)

Hence, our initial estimate exceeds the actual requirements by 22 bits.

It turns out that

log₂(f) = 10·log₂3 + 6·log₂5 + 3·log₂7 + 2·log₂11 + log₂13 + log₂17 + log₂19 + log₂23 = 61.68

which is indeed just below the size of the mantissa of 80 bit long double (64 bit). But when you multiply it by 26 (so, excluding the factor 2, which ends up in the exponent, by 13), you are adding log2(13) = 3.7 bits. 61.7+3.7 is 65.4, so from 26! onwards you no longer have the precision to perform your calculation exactly.

Firstly, nobody knows what long double can or cannot represent. Specific properties of the format depend on the implementation.

Secondly, X86 extended precision floating-point format has 64-bit significand with explicit leading 1, which means that it can represent contiguous integers in ±2⁶⁴ range. Beyond that range representable integers are non-contiguous (i.e. they begin to "skip" with wider and wider gaps). Your factorials fall far outside that range, meaning that it is perfectly expected that they won't be represented precisely.

Since 105! is a huge number that does not fit in a single word (or two of them), you want arbitrary precision arithmetic, also known as bignums. Use the Stirling's approximation to get an idea of how big 105! is and read the wikipage on factorials.

Standard C (read n1570 to check) don't have bignums natively, but you'll find many libraries for that.

I recommend GMPlib. BTW, some of its code is hand-written assembly for performance reason (when coding bignum addition, you want to take advantage of add with carry machine instructions).

I recommend to avoid writing your own bignum operations. The existing libraries use very clever algorithms (and you'll need to make some PhD work to get something better). If you try coding your own bignum library, it probably will be much worse than competitors (unless you spend years of work).