How to compute the natural logarithm of a factorial: large numbers

• 08-14-2010
TiagoPereira
How to compute the natural logarithm of a factorial: large numbers
Hello all,

Is there a function in C to compute the natural logarithm of a factorial, say, log(factorial(646))?

In most statistical packages, such as R, the factorial can be computed up to 177 (e.g. factorial(177) or so. But for statistical purposes, all formulae can be worked in the logarithm scale. As a result, these packages compute very large lnfactorials, such as lnfactorial(854511124).

I will be really grateful for any tips.

[it is important to note that I do not need approximations, since the lnfactorials will be used to compute statistical estimates that need to be extremely precise up to the 7th decimal place).

Cheers!

Tiago
• 08-14-2010
TiagoPereira
I guess I have solved this question with a simple function by noting the additive effects on the logarithm scale:

Code:

``` double lnfactorial( int a)   {           int y; double z;       if (a == 1)           return 0;       else       {               z = 0;                       for (y = 2; y<=a; y++ )                                 z = log(y)+z;                                           return z;       }   }```
• 08-14-2010
Char*Pntr
"int y" should be declared "double y" ?
Hello,

I tried the code below to verify the solution... I got the following error from my compiler:

"more than one instance of overloaded function "log" matches the argument list"

In the function definition, I noticed the "int y" which is used in " z = log( y ) + z ; " and when I changed the variable "y" to a type double, the code below compiled with no errors.

I'm using Microsoft VS 2010 compiler, and I have my doubts whether it's ANSI compliant in the first place. Especially when it comes to mathematical work. So if this is the case, please accept my apology in advance.

I'm very interested in seeing how this program works, is the reason why I tried it in the first place. Cheers!

Quote:

Originally Posted by TiagoPereira
I guess I have solved this question with a simple function by noting the additive effects on the logarithm scale:

Code:

``` double lnfactorial( int a)   {           int y; double z;       if (a == 1)           return 0;       else       {               z = 0;                       for (y = 2; y<=a; y++ )                                 z = log(y)+z;                                           return z;       }   }```

• 08-14-2010
TiagoPereira
Nice! Did you include the math.h library? I am under linux, so the commands to compile were:

gcc -Wall -I/usr/local/include -c lnfactorial.c

and then:

gcc -L/usr/local/lib lnfactorial.o -lm -o lnfactorial

I have tried the following code and it worked like a charm:

Code:

``` #include <stdio.h> #include <math.h> int x; double f; double lnfactorial( int a); int main(void) {                         printf("Enter an integer :");                 scanf("%d" , &x);                                         f = lnfactorial(x);                 printf("\nThe ln factorial of %d is %f\n", x, f);                 return 0;         }  double lnfactorial( int a)   {                     int y;           double z;       if (a == 1)           return 1;       else       {               z = 0;                       for (y = 2; y<=a; y++ )                                 z = log(y)+z;                                           return (z);       }   }```
On terminal, typing ./lnfactorial gives me, for example:

./lnfactorial
Enter an integer number:845123

The ln factorial of 845123 is 10688479.004670
• 08-14-2010
iMalc
Well the formatting is a little off, but that looks like it would work (on some compilers at least). Given log is usually overloaded for different types (float, double) and you're giving it an int, standards conformant compilers should complain about the ambiguity of passing log an int I believe. Making y a double as noted, would be better.
You can also write z += log(y) for simplicity.
Well done on solving your own question btw.
• 08-14-2010
Char*Pntr
Confirm Data
Yes, I included the <math.h> library file. In your initial post, you used a very large number: 854511124.

This is the number that I used to test your code. It took my computer aprox. 33 seconds to compute this with a 25% cpu load. (Quad core AMD 9850-2500 black edition CPU with high speed ram) & Windows 7.

The result was:16719398910.9735430000

You posted that for 845123 your result was: 10688479.004670

My result was: 10688479.0046699200
• 08-15-2010
TiagoPereira
Yes, your results look fine, as long as you don't need an extreme precision. I have used this simple function to compute the exact P-value for a couple of statistical tests, and they match exactly those reported in papers.

In R, I got:

> x<- lfactorial(854511124)
> print(x, digits=20)
[1] 16719398910.97366

Taking the code above, and changing everthing to the long double format, I got:

./lnfactorial
Enter an integer :854511124

The ln factorial of 854511124 is 16719398910.973660741933

For me, it takes around 42 seconds (old dual core, on linux). Hence, it stands to logic that that function is not efficient for large numbers. For example, the R function takes 0.05 seconds. I believe, however, that it uses a mathematical approximation, whereas we are computing the result based on a "brute force" way.

For statistical purposes, I guess your implementation is OK.

I am eager to see more efficient alternatives.
• 08-15-2010
mr.wu
asymptotic expression if you don't need high precision
Any math handbook or google can tell you the following formula for large N

ln(N!) \approx [ N ln(N) ] -N

``` #include <stdio.h> #include <math.h> int x; long double f; long double lnfactorial( int a); int main(void) {                         printf("Enter an integer :");                 scanf("%d" , &x);                                         f = lnfactorial(x);                 printf("\nThe ln factorial of %d is %9.12Lf\n", x, f);                 return 0; }  long double lnfactorial( int a)   {                   long double z;       if (a == 1)           return 0;       else       {               z = lgamma(a+1);                             return (z);       }   }```