algorithm - Find the inverse of a number modulo a prime

Question

I was thinking about an algorithm to solve the congruence ax = 1 mod p with p prime. I was thinking about using Fermat theorem. Since I know that

a ^ (p-1) = 1 mod p

and that

a ^ (p-1) = a * (a ^ (p-2))

It means that a ^ (p-2) mod p is the solution. Unfortunately this solution, although mathematically correct, isn't good for computer since for big primes I have to do a ^ (p-2) which is usually not calculable.

Which algorithm is good for computer science?

Answer 1

since for big primes I have to do a ^ (p-2) which is usually not calculable.

You need modular exponentiation, so with the exponentiation by squaring mentioned by IVlad you only need Θ(log p) modular multiplications of numbers of size at most p-1. The intermediate results are bounded by p^2, so despite a^(p-2) not being calculable for large primes, (a ^ (p-2)) % p usually is. That method is simple to implement:

unsigned long long invert_mod(unsigned long long a, unsigned long long p) {
    unsigned long long ex = p-2, result = 1;
    while (ex > 0) {
        if (ex % 2 == 1) {
            result = (result*a) % p;
        }
        a = (a*a) % p;
        ex /= 2;
    }
    return result;
}

but has a few drawbacks. (p-1)^2 must be representable in the used type (not a problem [except for huge p] if arbitrary precision integers are used), or you get invalid results due to overflow, and it always uses at least log (p-2)/log 2 modular multiplications.

The extended Euclidean algorithm, as suggested by user448810, or equivalently the continued fraction method, never produces intermediate values larger than p, thus avoiding all overflow problems if p is representable, and usually needs fewer divisions. Additionally, it computes the modular inverse not only for primes, but for any two coprime numbers.

unsigned long long invert_mod(unsigned long long a, unsigned long long p) {
    unsigned long long new = 1, old = 0, q = p, r, h;
    int pos = 0;
    while (a > 0) {
        r = q%a;
        q = q/a;
        h = q*new + old;
        old = new;
        new = h;
        q = a;
        a = r;
        pos = !pos;
    }
    return pos ? old : (p - old);
}

The code is slightly longer, but an optimising compiler ought to compile it to a short loop using only one division per iteration.