How to calculate (n!)%1000000009

Question

I need to find n!%1000000009. n is of type 2^k for k in range 1 to 20. The function I'm using is:

#define llu unsigned long long
#define MOD 1000000009

llu mulmod(llu a,llu b) // This function calculates (a*b)%MOD caring about overflows
{
    llu x=0,y=a%MOD;
    while(b > 0)
    {
        if(b%2 == 1)
        {
            x = (x+y)%MOD;
        }
        y = (y*2)%MOD;
        b /= 2;
    }
    return (x%MOD);
}

llu fun(int n)   // This function returns answer to my query ie. n!%MOD
{
    llu ans=1;
    for(int j=1; j<=n; j++)
    {
        ans=mulmod(ans,j);
    }
    return ans;
}

My demand is such that I need to call the function 'fun', n/2 times. My code runs too slow for values of k around 15. Is there a way to go faster?

EDIT: In actual I'm calculating 2*[(i-1)C(2^(k-1)-1)]*[((2^(k-1))!)^2] for all i in range 2^(k-1) to 2^k. My program demands (nCr)%MOD caring about overflows.

EDIT: I need an efficient way to find nCr%MOD for large n.

Answer 1

The mulmod routine can be speeded up by a large factor K.

1) '%' is overkill, since (a + b) are both less than N.
- It's enough to evaluate c = a+b; if (c>=N) c-=N;
2) Multiple bits can be processed at once; see optimization to "Russian peasant's algorithm"
3) a * b is actually small enough to fit 64-bit unsigned long long without overflow

Since the actual problem is about nCr mod M, the high level optimization requires using the recurrence

(n+1)Cr mod M = (n+1)nCr / (n+1-r) mod M.

Because the left side of the formula ((nCr) mod M)*(n+1) is not divisible by (n+1-r), the division needs to be implemented as multiplication with the modular inverse: (n+r-1)^(-1). The modular inverse b^(-1) is b^(M-1), for M being prime. (Otherwise it's b^(phi(M)), where phi is Euler's Totient function.)

The modular exponentiation is most commonly implemented with repeated squaring, which requires in this case ~45 modular multiplications per divisor.

If you can use the recurrence

nC(r+1) mod M = nCr * (n-r) / (r+1) mod M

It's only necessary to calculate (r+1)^(M-1) mod M once.