3

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?

talonmies
  • 70,661
  • 34
  • 192
  • 269
Zagorax
  • 11,440
  • 8
  • 44
  • 56
  • Are you aware of this paper on calculating Integer Reciprocals (http://ipa.ece.illinois.edu/mif/pubs/web-only/Frank-RawMemo12-1999.html)? I believe that it may be useful to you. – RBarryYoung Dec 30 '12 at 18:45

3 Answers3

12

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.

Community
  • 1
  • 1
Daniel Fischer
  • 181,706
  • 17
  • 308
  • 431
8

The normal way to compute the modular inverse is by the extended Euclidean algorithm:

function inverse(x, m)
    a, b, u := 0, m, 1
    while x > 0
        q, r := divide(b, x)
        x, a, b, u := b % x, u, x, a - q * u
    if b == 1 return a % m
    error "must be coprime"
user448810
  • 17,381
  • 4
  • 34
  • 59
  • That's normal? Says who? Normal means nothing in science. It's a valid alternative, definitely, but you can't say this is the "normal way". – IVlad Dec 30 '12 at 18:59
  • I highly doubt those were his exact words, but either way, it doesn't mean it's right. It's not any more or less normal than using Fermat's theorem. – IVlad Dec 30 '12 at 19:07
  • @IVlad The extended euclidean algorithm is the normal way of computing the inverse modulo a prime. Where for "normal" I(and probably user448810 too) mean what is actually used in any serious work. Just go read some papers that talk about prime inverses and you'll see that everybody uses the above algorithm, since it is much better then the modular-exponentiation. – Bakuriu Dec 30 '12 at 21:24
  • @Bakuriu - uh, no. Nobody in a serious paper talks about which algorithm is normal and which isn't. They choose one over another for its advantages in that context. Saying one is normal means nothing. It's not "much better then[sic]", each one has its uses, advantages and disadvantages. – IVlad Dec 31 '12 at 00:42
  • @IVlad What I meant is this: Let us redefine "normal" in a clearer way: we say that an algorithm A is the normal approach to solve problem P if more than 51% of the papers and programs(ever written in the history of the humankind) that deal with problem P use algorithm A to solve it. With this definition of "normal" the extended euclidean algorithm is the normal algorithm to solve the problem of finding the inverse modulo a prime. – Bakuriu Dec 31 '12 at 08:52
1

There is no reason this isn't a good algorithm for computers, you just need to be careful about the implementation, which isn't exactly trivial I guess, but it's not difficult either.

Just use exponentiation by squaring, then it most likely won't matter how big p is.

a^n = a^(n / 2) * a^(n / 2) for n even
    = a*a^(n - 1) for n odd
IVlad
  • 43,099
  • 13
  • 111
  • 179
  • 2
    I'm pretty sure you wrote that incorrectly...`a^n = a^(n / 2) for n even` doesn't make any sense, and it's not what that link you posted says either – nbrooks Dec 30 '12 at 20:19