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

Since you are looking for nCr for multiple sequential values of n you can make use of the following:

    (n+1)Cr = (n+1)! / ((r!)*(n+1-r)!)
    (n+1)Cr = n!*(n+1) / ((r!)*(n-r)!*(n+1-r))
    (n+1)Cr = n! / ((r!)*(n-r)!) * (n+1)/(n+1-r)
    (n+1)Cr = nCr * (n+1)/(n+1-r)

This saves you from explicitly calling the factorial function for each i.

Furthermore, to save that first call to nCr you can use:

    nC(n-1) = n //where n in your case is 2^(k-1).

EDIT:
As Aki Suihkonen pointed out, (a/b) % m != a%m / b%m. So the method above so the method above won't work right out of the box. There are two different solutions to this:

1. 1000000009 is prime, this means that a/b % m == a*c % m where c is the inverse of b modulo m. You can find an explanation of how to calculate it here and follow the link to the Extended Euclidean Algorithm for more on how to calculate it.

2. The other option which might be easier is to recognize that since nCr * (n+1)/(n+1-r) must give an integer, it must be possible to write n+1-r == a*b where a | nCr and b | n+1 (the | here means divides, you can rewrite that as nCr % a == 0 if you like). Without loss of generality, let a = gcd(n+1-r,nCr) and then let b = (n+1-r) / a. This gives (n+1)Cr == (nCr / a) * ((n+1) / b) % MOD. Now your divisions are guaranteed to be exact, so you just calculate them and then proceed with the multiplication as before. EDIT As per the comments, I don't believe this method will work.

Another thing I might try is in your llu mulmod(llu a,llu b)

    llu mulmod(llu a,llu b)
    {
        llu q = a * b;
        if(q < a || q < b) // Overflow!
        {
            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);
        }
        else
        {
            return q % MOD;
        }
    }

That could also save some precious time.

Answer 2

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.