# Computation of inverses modulo p followup

In responding to

Fast computation of multiplicative inverse modulo q

I mentioned an algorithm for computing the inverse of $a \mod p$ different from the extended Euclidean algorithm, hoping that someone could tell me how its speed stacks up against other algorithms. Since no one did, I'm asking directly if someone can tell me.

To compute the inverse of a modulo p, you can run the Euclidean algorithm starting with $p^2$ and $ap+1$, comparing the size of each remainder with $p$. The first remainder less than $p$ that appears will be an inverse for $a \mod p$.

This will always take either the same number of steps to reach the inverse as it takes to reach $\gcd(a,p)=1$ using the Euclidean algorithm with $a$ and $p$, or else one additional step, depending on whether the least positive residue of an inverse of $a$ is greater than $p/2$ or less than $p/2$. Thus, it requires approximately the same number of computations as the first half of the extended Euclidean algorithm (albeit with bigger numbers initially), excepting an extra comparison with $p$ at each step.

Question: How does the speed of this algorithm compare to others?

Aside: Pedagogically, this is nice since the second half of the extended Euclidean algorithm is the one my students tend to mess up. However, assuming our ultimate goal is for students to understand why they're doing what they're doing, perhaps the extended Euclidean algorithm is preferable.

-
Quadratic complexity. Compared to standard quadratic algorithm: half the iterations, on twice-the-size numbers. Since the operations are more than linear in size, I would guess that the constant is larger than the constant of the standard algorithm. –  Dror Speiser Oct 12 '10 at 16:10
If you examine AVS's answer to the question you linked, you will see that your algorithm is substantially slower than Fast Euclidean. –  S. Carnahan Oct 12 '10 at 16:23
I didn't mention this in my answer to the linked question, but the half-gcd algorithm can be used to "jump" directly to any particular step of the extended Euclidean computation (not necessarily the last one) using a logarithmic number of recursive calls. One could perform a binary search to find the greatest remainder less than p using half-gcds, yielding a quasi-linear running time that is slower than the fast Euclidean algorithm by only a log factor, and even this might be avoided (at least on average) by jumping to the middle step and searching linearly from there. –  AVS Oct 12 '10 at 18:39
The algorithm sketched in my comment above was only meant to show that your approach could be implemented in quasi-linear time. But this comes at the cost of making the algorithm much more complicated (which may be a pedagogical flaw or feature, depending on one's point of view), and it will still be slower than just computing the extended gcd. –  AVS Oct 12 '10 at 18:45
Ah, thanks AVS! I wasn't expecting that performing the algorithm exactly as I wrote it would give a speed-up, but I was thinking that one of the other speed-ups might improve the algorithm I wrote to be comparable to others. One knows that the first remainder less than p will occur almost exactly half-way through the Euclidean algorithm with $p^2$ and $ap+1$, so jumping to the middle step will put you very close. –  Barry Oct 12 '10 at 18:45

Here is a tidy way to solve $ax+by=gcd(a,b)$. Start with the matrix $$\left(\begin{array}{ccc} 1&0&a\\ 0&1&b\\ \end{array}\right)$$ Suppose $a\ge b$. Then replace row 1 by row 1 minus $t$ times row 2, where $t=\lfloor a/b\rfloor$. Repeat this operation until the last entry in one of the two rows is zero. If the other row is $x,y,d$ then $ax+by=d$ and $d=gcd(a,b)$. This is very simple to program, and avoids the back-substitutions that students find confusing.

Most of the speedups of the Euclidean Algorithm described in Knuth's Art of Computer Programming v2 work here also.

-
I use this way to teach it as well. One has the invariant that a typical row is $r s ar+bs$. In case of an error one can zero in on the line where the error occurred by checking back and forth to see where the invariant failed. When working in $\mathbb{Z}[\sqrt{2}]$ or $\mathbb{Z}[i]$ that can be helpful. –  Aaron Meyerowitz Oct 13 '10 at 7:58
@Aaron: Strangely, I've never seen it in print. –  SJR Oct 13 '10 at 11:31

As you mention pedagogy: I solve $a x + b y = \gcd(a,b)$ by simple continued fractions, both by hand and in C++ code. For consecutive convergents $\frac{p_n}{q_n}$ and $\frac{p_{n+1}}{q_{n+1}},$ the standard relation is $p_n q_{n+1} - p_{n+1} q_n = (-1)^{n+1}$ where the starting point (from Khinchin's little book, Dover reprint) is $$\frac{p_{-2}}{q_{-2}} = \frac{0}{1}$$ and the fake but necessary $$\frac{p_{-1}}{q_{-1}} = \frac{1}{0}$$ For a rational number $\frac{a}{p}$ the final little 2 by 2 determinant may give the wrong sign (and does about half the time), in which case I negate everything.

So, this is the Euclidean algorithm but it all moves forward on the page, no back-substitutions or whatever is in the "second half" of extended gcd. And continued fractions are good for other things.

-