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

This is a discussion on How to compute the natural logarithm of a factorial: large numbers within the C Programming forums, part of the General Programming Boards category; Hello all, Is there a function in C to compute the natural logarithm of a factorial, say, log(factorial(646))? In most ...

  1. #1
    Registered User
    Join Date
    Aug 2010
    Posts
    20

    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
    Last edited by TiagoPereira; 08-14-2010 at 03:18 PM. Reason: To make the text clearer

  2. #2
    Registered User
    Join Date
    Aug 2010
    Posts
    20
    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;
          }
      }
    Last edited by TiagoPereira; 08-14-2010 at 03:45 PM.

  3. #3
    Registered User Char*Pntr's Avatar
    Join Date
    Sep 2007
    Location
    Lathrop, CA
    Posts
    198

    Cool "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 View Post
    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;
          }
      }
    Last edited by Char*Pntr; 08-14-2010 at 10:32 PM.

  4. #4
    Registered User
    Join Date
    Aug 2010
    Posts
    20
    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

  5. #5
    Algorithm Dissector iMalc's Avatar
    Join Date
    Dec 2005
    Location
    New Zealand
    Posts
    6,308
    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.
    My homepage
    Advice: Take only as directed - If symptoms persist, please see your debugger

    Linus Torvalds: "But it clearly is the only right way. The fact that everybody else does it some other way only means that they are wrong"

  6. #6
    Registered User Char*Pntr's Avatar
    Join Date
    Sep 2007
    Location
    Lathrop, CA
    Posts
    198

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

    You posted that for 845123 your result was: 10688479.004670

    My result was: 10688479.0046699200
    Last edited by Char*Pntr; 08-15-2010 at 12:14 AM. Reason: Update for clarification

  7. #7
    Registered User
    Join Date
    Aug 2010
    Posts
    20
    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.

  8. #8
    Registered User
    Join Date
    Jun 2010
    Posts
    19

    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

    Try that with your number.

  9. #9
    Registered User
    Join Date
    Aug 2010
    Posts
    20
    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.

Popular pages Recent additions subscribe to a feed

Similar Threads

  1. Large Numbers
    By rkjd2 in forum C++ Programming
    Replies: 2
    Last Post: 12-09-2002, 05:49 PM
  2. the definition of a mathematical "average" or "mean"
    By DavidP in forum A Brief History of Cprogramming.com
    Replies: 7
    Last Post: 12-03-2002, 10:15 AM
  3. Very large floating point numbers
    By bulsquare in forum C Programming
    Replies: 0
    Last Post: 04-08-2002, 04:10 PM
  4. printing large decimal numbers
    By Isometric in forum C++ Programming
    Replies: 5
    Last Post: 03-20-2002, 08:55 PM
  5. commas in large numbers
    By HKR in forum C Programming
    Replies: 7
    Last Post: 03-06-2002, 06:08 PM

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21