Is there a mathematical technique for "n choose r modulus m"?

Question

How to calculate nCr % m (i.e. "(n choose r) modulus m") in programming when n, r, m are very large numbers?

Answer 1

There are techniques to reduce the complexity of computing binomial coefficients modulo an integer m, provided that none of m's prime factors is too large resp. divides m with a too high power.

The first step is factorising m,

m = ∏ (p_i ^ e_i)

then one computes the binomial coefficient modulo each of these prime powers, and combines the result with the Chinese remainder theorem.

The computation of a binomial coefficient modulo a prime (e_i == 1) can be computed a little easier than the general case, cf. e.g. this answer, but it can also be subsumed there.

For a prime p and an n >= 0, let us define

F(p, n) =   ∏ k   = p^(n/p) * (n/p)!
         1<=k<=n
           p | k

and

G(p, n) =    ∏ k
          1<=k<=n
         gcd(k,n)=1

Then we have

n! = F(p, n) * G(p, n)

and iteratively, using the same splitting for the (n/p)! appearing in F(p, n),

           m
n! = p^K * ∏ G(p, n/(p^j))
          j=0

where p^m <= n < p^(m+1). All factors in G(p, x) are coprime to p^e, so the corresponding factors in the denominator of the binomial coefficient can be inverted modulo p^e, and if we find an efficient way to compute a G(p, x) modulo p^e, we have an efficient way to compute the binomial coefficient modulo p^e.

For the binomial coefficient, we then have

n! / (r! * (n-r)!) = p^M * (∏ G(p, (n/p^j)) * [ ∏ G(p, r/(p^j)) * ∏ G(p, (n-r)/(p^j)) ]^(-1)

Let H(p, e, n) = G(p, n) % (p^e). The crucial point is that the product of all numbers coprime to p^e not exceeding p^e is quite simple. It is congruent to -1 modulo p^e, unless p = 2 and e > 2, in which case it is congruent to 1.

So

H(p, e, n) ≡ (-1)^(n/(p^e)) * H(p, e, n % (p^e)) (mod p^e)

(unless p = 2 and e > 2, in which case the first factor is 1), and we only need to compute H(p, e, k) for 0 <= k < p^e, then we can look up the result.

Code:

// invert k modulo p, k and p are supposed coprime
unsigned long long invertMod(unsigned long long k, unsigned long long p) {
    unsigned long long q, pn = 1, po = 0, r = p, s = k;
    unsigned odd = 1;
    do {
        q = r/s;
        q = pn*q + po;
        po = pn;
        pn = q;
        q = r%s;
        r = s;
        s = q;
        odd ^= 1;
    }while(pn < p);
    return odd ? p-po : po;
}

// Calculate the binomial coefficient (n choose k) modulo (prime^exponent)
// requires prime to be a prime, exponent > 0, and 0 <= k <= n,
// furthermore supposes prime^exponent < 2^32, otherwise intermediate
// computations could have mathematical results out of range.
// If k or (n-k) is small, a direct computation would be more efficient.
unsigned long long binmod(unsigned long long prime, unsigned exponent,
                          unsigned long long n, unsigned long long k) {
    // The modulus, prime^exponent
    unsigned long long ppow = 1;
    // We suppose exponent is small, so that exponentiation by repeated
    // squaring wouldn't gain much.
    for(unsigned i = 0; i < exponent; ++i) {
        ppow *= prime;
    }
    // array of remainders of products
    unsigned long long *remainders = malloc(ppow * sizeof *remainders);
    if (!remainders) {
        fprintf(stderr, "Allocation failure\n");
        exit(EXIT_FAILURE);
    }
    for(unsigned long long i = 1; i < ppow; ++i) {
        remainders[i] = i;
    }
    for(unsigned long long i = 0; i < ppow; i += prime) {
        remainders[i] = 1;
    }
    for(unsigned long long i = 2; i < ppow; ++i) {
        remainders[i] *= remainders[i-1];
        remainders[i] %= ppow;
    }

    // Now to business.
    unsigned long long pmult = 0, ntemp = n, ktemp = k, mtemp = n-k,
                       numer = 1, denom = 1, q, r, f;
    if (prime == 2 && exponent > 2) {
        f = 0;
    } else {
        f = 1;
    }
    while(ntemp) {
        r = ntemp % ppow;
        q = ntemp / ppow;
        numer *= remainders[r];
        numer %= ppow;
        if (q & f) {
            numer = ppow - numer;
        }
        ntemp /= prime;
        pmult += ntemp;
    }
    while(ktemp) {
        r = ktemp % ppow;
        q = ktemp / ppow;
        denom *= remainders[r];
        denom %= ppow;
        if (q & f) {
            denom = ppow - denom;
        }
        ktemp /= prime;
        pmult -= ktemp;
    }
    while(mtemp) {
        r = mtemp % ppow;
        q = mtemp / ppow;
        denom *= remainders[r];
        denom %= ppow;
        if (q & f) {
            denom = ppow - denom;
        }
        mtemp /= prime;
        pmult -= mtemp;
    }
    // free memory before returning, we don't use it anymore
    free(remainders);
    if (pmult >= exponent) {
        return 0;
    }
    while(pmult > 0) {
        numer = (numer * prime) % ppow;
        --pmult;
    }
    return (numer * invertMod(denom, ppow)) % ppow;
}

which computes n choose k modulo p^e in O(p^e + log n) steps.

Answer 2

I assume you want to calculate nCr kind of stuff for large numbers and searching for some library functions which may be optimized. In that case you can look at GNU Scientific Library.

Answer 3

First of all, it depends on what you mean by "very large". If your problem is simply that the intermediate values become too large for standard 64-bit integers, then you could use something like gmp or BigInteger.

It's possible that the numbers involved will get so big that you will run out of memory or patience, and you won't be able to calculate the intermediate values to complete precision. In this case, your best bet is to first determine the prime factorization of each factorial, use these factorizations to determine the factorization of the binomial, then multiply this factorization out using the modulus at each intermediate step.

You will need a list of all primes up to n.

The following is pseudocode. I'm using int, but you should replace it with whatever large number library you're using.

int factorial_prime_power(int f, int p) {
    // When p is prime, returns k such that p^k divides f!
    int k = 0;
    while (f > p) {
        f = f / p;
        k = k + f;
    }
    return k;
}
int binomial_prime_power(int n, int r, int p) {
    // when p is prime, returns k such that p^k divides nCr
    return factorial_prime_power(n,p) - factorial_prime_power(r,p) - factorial_prime_power(n-r,p);
}
int powmod(int p, int k, int m) {
    // quickly calculates p^k mod m
    int res = 1;
    int q = p;
    int j = k;
    while (j > 0) {
        // invariant:  p^k is congruent to res * q^j
        if (j is odd) {
            res = (res * q) % m;
            j = (j-1)/2;
        } else {
            j = j / 2;
        }
        q = (q * q) % m;
    }
    return res;
}
int big_binomial(int n, int r, int m) {
    if (n < r or r < 0) {
        return 0;
    }
    int res = 1;
    for(p in all primes from 2 to n) {
        k = binomial_prime_power(n,r,p);
        res = (res * powmod(p,k,m)) % m;
    }
    return res;
}

Answer 4

Do it in chunks.

For example: n!/(r!(n-r)!) = 1*2*...*n/(1*2*...r*1*2*...*(n-r))=1/1*1 * 2/2*2 * 3/3*3 * ... *... .

Each chunk is easily computable thus you should avoid overflows in your computations. Compute chunk and multiply your current result by it.

Also, it's worth doing modulo on n and r before that.