JOIN
Get Time
forums   
Search | Watch Thread  |  My Post History  |  My Watches  |  User Settings
View: Flat  | Threaded  | Tree
Previous Thread  |  Next Thread
Modular division | Reply
This article lacked an overview of modular division -- in practice, that's the most useful number-theoretic trick.

For those who have never heard about it: given a modulus M, the inverse of a number x is a number a for which a*x == 1 (mod M) (a is unique (mod M), as long as x and M are coprime). One can obtain it with an algorithm like this:

int inv(int x, int M) {
  int a = 1, b = x;
  while (b != 1) {
    int c = M / b;
    a *= c; a %= M;
    b *= c; b %= M;
    if (b > M/2) { a = M - a; b = M - b; }
  }
  return a;
}


Basically, what this does is start with the equation x == x (mod M), and keep multiplying by suitable factors so the right-hand side becomes 1.

Example: Take x == 3 (mod 5). In the first step, c = 1 and the if triggers, so the equation becomes 4x == 2 (mod 5). In the second step, c = 2 and the if triggers, so the equation becomes 2x == 1 (mod 5), and we conclude that the inverse of 3 (mod 5) is 2.

(The code above can fail for composite M; fixing it wuch that this is not the case anymore is left as an exercise for the reader.)

The cool part about modular inverses is that they're a drop-in replacement for division. A naive implementation of a function that returns binomial numbers can read like this:

int binom(int n, int p) { // binom(n, p) = C(n, p) [?]
  return fact(n) / (fact(p) * fact(n-p)); // fact(k) = k!
}


If all we're interested in is the value of binom() modulo, say, 1000000007, we could write binom(n, p) % 1000000007, but that will fail due to overflows for n as small as 21, even using 64-bit integers. With modular inverses, we can replace the division by fact(p) * fact(n-p) by a multiplication by inv(fact(p) * fact(n-p)):

int binom(int n, int p) { // binom(n, p) = C(n, p) [?]
  return fact(n) * inv(fact(p) * fact(n-p)); // fact(k) = k!
}


Of course, the problem of overflow still remains. But now, we're doing only additions, subtractions and multiplications, so we need not operate with the full values of fact(), but merely with their reductions (mod 1000000007):

int fact(int n) { // fact(k) = k! % 1000000007
  int ans = 1;
  for (int i = 1; i <= n; i++) {
    ans *= i; ans %= 1000000007;
  }
  return ans;
}


Finally, we have binom() returning correct values for all values of n and p:

int binom(int n, int p) { // binom(n, p) = C(n, p) % 1000000007
  int num = fact(n);
  int den = (fact(p) * fact(n-p)) % 1000000007;
  return (num * inv(den)) % 1000000007;
}


(Of course, many of the ints here should actually be long longs, and one should precalculate fact(), inv(), etc., but I digress.)
Subject Author Date
Modular division ctgPi Mar 5, 2007 at 11:54 PM EST
Re: Modular division dskloet Mar 6, 2007 at 2:59 AM EST
Re: Modular division misof Mar 6, 2007 at 4:09 AM EST
Re: Modular division sql_lall Mar 6, 2007 at 4:28 AM EST
Re: Modular division ctgPi Mar 6, 2007 at 10:04 AM EST
RSS