Get Time
Search | Watch Thread  |  My Post History  |  My Watches  |  User Settings
View: Flat (newest first)  | 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.)
Re: Modular division (response to post by ctgPi) | Reply
I don't think I really get the idea of your algorithm. Can you explain why it terminates? And is it better than using the Extended Euclidean algorithm?
Re: Modular division (response to post by dskloet) | Reply
... or better than using Fermat's little theorem?

for any prime p, and for any a not divisible by p:
a^{p-1} \equiv 1 (mod p)
a^{p-2} * a \equiv 1 (mod p)
and thus
a^{p-2} is the multiplicative inverse of a modulo p.
Re: Modular division (response to post by misof) | Reply
That'd get trickier if p isn't prime (would you have to calculate phi(p)?)
Anyway, nothing wrong with one line to solve ax == b (mod c):
int axbmodc(int a,int b,int c){ return a ? (axbmodc(c % a, (a - b % a) % a, a) * c + b) / a : 0; }
Re: Modular division (response to post by dskloet) | Reply
Since M % x < x, (M - x * M/x) = M % x < x, so a decreases every iteration, and therefore the algorithm always terminates. It's probably slightly worse than an extended GCD, but I can type it from memory.

I'll usually pull out my code library and copy the EGCD if I have to solve a*x == b (mod m) with a, m possibly not coprime, though.