How to Efficiently Calculate N! Using a Fixed-Point Bignumber Library?

Problem:

Determine an efficient method to calculate N! (N factorial) with no precision loss using a fixed-point bignumber library. Specifically, find the formula for:

T2 = { (2N+1).(2N+3).(2N+5)...(4N-1) }.(2^N)/(N!)

Solution:

1. Compute T1:

   T1 = T2 * N!

   Copy after login

   where N! is determined using the following formula:

   N! = ((2N)!).((2N)!).{ (2N+1).(2N+3).(2N+5)...(4N-1) }

   Copy after login

   This allows for the computation of N! from T1.

2. Calculate the exponent e for prime i:

   e = sum(j=1,2,3,4,5,...) of (4N/(i^j))-(2N/(i^j)) )

   Copy after login

   where i is any prime <= 4N
   and j is any integer such that i^j <= 4N

3. Compute T2:

   T2 = multiplication(i=all primes<=4N) of [i^e]

   Copy after login

   where e is the exponent calculated in Step 2.

Code Snippet for New Equation:

// edit beg:
// Sorry, forget to copy sorted list of all primes up to max n here it is
// end of table is marked with 0
// Primes are in DWORDs so they only 4Byte per number
// so the table is very small compared with lookup table for the same max n!
// and also primes are needed for many other routines in bignum
// can compute n! for n <= max prime in table
DWORD _arithmetics_primes[]={2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,0};
// edit end.

longnum fact(const DWORD &x)
    {
    if (x<=4)
        {
        if (x==4) return 24;
        if (x==3) return  6;
        if (x==2) return  2;
        if (x==1) return  1;
        if (x==0) return  1;
        }
    int N4,N2,p,i,j,e; longnum c,pp;
    N4=(x>>2)<<2;
    N2=N4>>1;
    c=fact(N2); c*=c;                 // c=((2N)!)^2;
    for (i=0;;i++)                    // c*= T2
        {
        p=_arithmetics_primes[i];
        if (!p) break;
        if (p>N4) break;
        for (e=0,j=N4;j;e+=j&1,j/=p);
        if (e)                        // c*=p^e
            {
            if (p==2) c<<=e;
            else for (pp=p;;)
                {
                if (int(e&1)) c*=pp;
                e>>=1; if (!e) break;
                pp*=pp;
                }
            }
        }
    for (i=N4+1;i<=x;i++) { c*=i; } c.round();
    return c;
    }

Time Measurements:


N!
Time (ms)



1!
0.001


2!
0.000


3!
0.000


4!
0.006


5!
0.006


6!
0.007


7!
0.005


8!
0.006


9!
0.007


10!
0.008


11!
0.012


12!
0.013


13!
0.014


14!
0.016


15!
0.014


16!
0.015


17!
0.017


18!
0.019


19!
0.016


20!
0.017


21!
0.019


22!
0.021


23!
0.023


24!
0.025


25!
0.027


26!
0.029


27!
0.032


28!
0.034


29!
0.037


30!
0.039


31!
0.034


32!
0.037


33!
0.039


34!
0.041


35!
0.039


36!
0.041


37!
0.044


38!
0.046


39!
0.041


40!
0.044


41!
0.046


42!
0.049


43!
0.048


44!
0.050


45!
0.054


46!
0.056


47!
0.056


48!
0.060


49!
0.063


50!
0.066


51!
0.065


52!
0.069


53!
0.072


54!
0.076


55!
0.077


56!
0.162


57!
0.095


58!
0.093


59!
0.089


60!
0.093


61!
0.098


62!
0.096


63!
0.090


64!
0.100


65!
0.104


66!
0.111


67!
0.100


68!
0.121


69!
0.109


70!
0.119


71!
0.104


72!
0.124


73!
0.113


74!
0.118


75!
0.118


76!
0.123


77!
0.129


78!
0.133

79



The above is the detailed content of How to Efficiently Calculate N! Using a Fixed-Point Bignumber Library?. For more information, please follow other related articles on the PHP Chinese website!