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("%.0Lf\n",factorial(25)); // this result is correct
printf("%.0Lf\n",factorial(26));
//correct result is 403291461126605635584000000,but it return 403291461126605635592388608
return 0;
}
Back of the envelope calculation: 25! is slightly more than 1025; 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 = 310 · 56 · 73 · 112 · 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
log2(f) = 10·log23 + 6·log25 + 3·log27 + 2·log211 + log213 + log217 + log219 + log223 = 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.