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

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;
      }
  }

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!

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;
      }
  }

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

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.

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
Please confirm.

You posted that for 845123 your result was: 10688479.004670

My result was: 10688479.0046699200

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.

Any math handbook or google can tell you the following formula for large N

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

Try that with your number.

No, no! Please, don't use that approximation for statistical tests!

If you want to implement a routine for other purposes, where you precision is not that important you can use it. If you need precision, however, (i.e. to compute exact P-values for some statistical tests), the best compromise between efficiency and accuracy is to use the gamma approximation:

Code:

#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);
      }
  }

The code above calculates the natural logarithm of n! for very large numbers of n using the gamma approximation. For me, (32-bits, linux) the maximum n can be 2147483646, since 2147483647 is the limit for a 32-bits system.