Fast N Choose K Mod P for Large N

Fast n choose k mod p for large n?

So, here is how you can solve your problem.

Of course you know the formula:

comb(n,k) = n!/(k!*(n-k)!) = (n*(n-1)*...(n-k+1))/k! 

(See http://en.wikipedia.org/wiki/Binomial_coefficient#Computing_the_value_of_binomial_coefficients)

You know how to compute the numerator:

long long res = 1;
for (long long i = n; i > n- k; --i) {
  res = (res * i) % p;
}

Now, as p is prime the reciprocal of each integer that is coprime with p is well defined i.e. a^-1 can be found. And this can be done using Fermat's theorem a^p-1=1(mod p) => a*a^p-2=1(mod p) and so a^-1=a^p-2.
Now all you need to do is to implement fast exponentiation(for example using the binary method):

long long degree(long long a, long long k, long long p) {
  long long res = 1;
  long long cur = a;

  while (k) {
    if (k % 2) {
      res = (res * cur) % p;
    }
    k /= 2;
    cur = (cur * cur) % p;
  }
  return res;
}

And now you can add the denominator to our result:

long long res = 1;
for (long long i = 1; i <= k; ++i) {
  res = (res * degree(i, p- 2)) % p;
}

Please note I am using long long everywhere to avoid type overflow. Of course you don't need to do k exponentiations - you can compute k!(mod p) and then divide only once:

long long denom = 1;
for (long long i = 1; i <= k; ++i) {
  denom = (denom * i) % p;
}
res = (res * degree(denom, p- 2)) % p;

EDIT: as per @dbaupp's comment if k >= p the k! will be equal to 0 modulo p and (k!)^-1 will not be defined. To avoid that first compute the degree with which p is in n*(n-1)...(n-k+1) and in k! and compare them:

int get_degree(long long n, long long p) { // returns the degree with which p is in n!
  int degree_num = 0;
  long long u = p;
  long long temp = n;

  while (u <= temp) {
    degree_num += temp / u;
    u *= p;
  }
  return degree_num;
}

long long combinations(int n, int k, long long p) {
  int num_degree = get_degree(n, p) - get_degree(n - k, p);
  int den_degree = get_degree(k, p);

  if (num_degree > den_degree) {
    return 0;
  }
  long long res = 1;
  for (long long i = n; i > n - k; --i) {
    long long ti = i;
    while(ti % p == 0) {
      ti /= p;
    }
    res = (res * ti) % p;
  }
  for (long long i = 1; i <= k; ++i) {
    long long ti = i;
    while(ti % p == 0) {
      ti /= p;
    }
    res = (res * degree(ti, p-2, p)) % p;
  }
  return res;
}

EDIT: There is one more optimization that can be added to the solution above - instead of computing the inverse number of each multiple in k!, we can compute k!(mod p) and then compute the inverse of that number. Thus we have to pay the logarithm for the exponentiation only once. Of course again we have to discard the p divisors of each multiple. We only have to change the last loop with this:

long long denom = 1;
for (long long i = 1; i <= k; ++i) {
  long long ti = i;
  while(ti % p == 0) {
    ti /= p;
  }
  denom = (denom * ti) % p;
}
res = (res * degree(denom, p-2, p)) % p;

Finding binomial coefficient for large n and k modulo m

Just use the fact that

(n, k) = n! / k! / (n - k)! = n*(n-1)*...*(n-k+1)/[k*(k-1)*...*1]

so you actually have just 2*k=2*10^5 factors. For the inverse of a number you can use suggestion of kfx since your m is prime.

How can we compute N choose K modulus a prime number without overflow?

To compute (n choose k) % M, you can separately compute the nominator (n!) modulus M and the denominator (k!*(n - k)!) modulus M and then multiply the nominator by the denominator's modular multiplicative inverse (in M). Since M is prime, you can use Fermat's Little Theorem to calculate the multiplicative inverse.

There is a nice explanation, with sample code, on the following link (problem SuperSum):

http://www.topcoder.com/wiki/display/tc/SRM+467

nCk modulo p when n % p or k % p == 0

You can use Lucas' theorem to reduce the problem to ceil(log_P(N)) subproblems with k, n < p: Write n = n_m * p^m + ... + n_0 and k = k_m * p^m + ... + k_0 in base p (n_i, k_i < p are the digits), then we have

C(n,k) = PROD(i = 0 to m, C(n_i, k_i))  (mod p)

The subproblems are easy to solve, because every factor of k! has an inverse modulo p. You get an algorithm with runtime complexity O(p log(n)), which is better than that of Ivaylo's code in case of p << n, if I understand it correctly.

int powmod(int x, int e, int p) {
    if (e == 0) return 1;
    if (e & 1) return (long long)x * powmod(x, e - 1, p) % p;
    long long rt = powmod(x, e / 2, p);
    return rt * rt % p;
}

int binom_coeff_mod_prime(int n, int k, int p) {
    long long res = 1;
    while (n || k) {
        int N = n % p, K = k % p;
        for (int i = N - K + 1; i <= N; ++i)
            res = res * i % p;
        for (int i = 1; i <= K; ++i)
            res = res * powmod(i, p - 2, p) % p;
        n /= p;
        k /= p;
    }
    return res;
}

Finding nCk that can have huge values of n and k

Because 2^(p-1)==1 mod p, you can do all the calculations of the exponents modulo p-1.

How do I find mod of large C(n , r)

I know algorithm with complexity O(r*log_n)
Firstly look at algorithm to calc C(n,r) without mod k:

int res = 1;
for(int i=1; i<=r; i++){
  res*=(n+1-i);
  res/=i;
}

In your case, you can't divide, because you use modular arithmetics. But you can multiply on the modular multiplicative inverse element, information about it you can find here https://en.wikipedia.org/wiki/Modular_multiplicative_inverse.
You code will be like this:

int res = 1;
for(int i=1; i<=r; i++){
  res*=(n+1-i);
  res%=k;
  res*=inverse(i,k);
  res%=k;
}

---

Related Topics