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.]]>

]]>intaxbmodc(inta,intb,intc){returna ? (axbmodc(c % a, (a - b % a) % a, a) * c + b) / a : 0; }

for any prime p, and for any a not divisible by p:

a^{p-1} \equiv 1 (mod p)

therefore

a^{p-2} * a \equiv 1 (mod p)

and thus

a^{p-2} is the multiplicative inverse of a modulo p.]]>

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:

intinv(intx,intM) {inta = 1, b = x;while(b != 1) {intc = M / b; a *= c; a %= M; b *= c; b %= M;if(b > M/2) { a = M - a; b = M - b; } }returna; }

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.

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:

intbinom(intn,intp) { // binom(n, p) = C(n, p) [?]returnfact(n) / (fact(p) * fact(n-p)); // fact(k) = k! }

If all we're interested in is the value of

intbinom(intn,intp) { // binom(n, p) = C(n, p) [?]returnfact(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

intfact(intn) { // fact(k) = k! % 1000000007intans = 1;for(inti = 1; i <= n; i++) { ans *= i; ans %= 1000000007; }returnans; }

Finally, we have

intbinom(intn,intp) { // binom(n, p) = C(n, p) % 1000000007intnum = fact(n);intden = (fact(p) * fact(n-p)) % 1000000007;return(num * inv(den)) % 1000000007; }