Extended Euclidean algorithm

The extended Euclidean algorithm finds the greatest common divisor $g$ of two numbers $a$ and $b$ together with integers $x$ and $y$ such that a \cdot x + b \cdot y = g. With them you can solve equations such as 75x + 21y = 12 in integers and divide modulo any number. If you only need the GCD, the Euclidean algorithm is enough.

Explanation

A shop has a pan balance and many weights of 5 g and 3 g. Weights can go on either pan. To weigh 1 g of salt, put two 3 g weights on one pan and one 5 g weight on the other. The pans differ by 6 − 5 = 1 g, so 1 g of salt next to the 5 g weight makes them level.

Fig. 1. Two 3 g weights balance a 5 g weight and 1 g of salt.

Now the weights are 75 g and 21 g. Both are multiples of 3, so whatever we put on the pans, they differ by a whole number of 3 g: 3 is the greatest common divisor of 75 and 21. So 3 g is the least we can hope to weigh. Can we weigh exactly 3 g, and with how many weights of each kind?

The first idea is to try one 75 g weight, then two, and so on, and each time check whether some number of 21 g weights is exactly 3 g lighter or heavier. One 75 doesn't work: three 21s make 63, which is 12 g less, and four make 84, which is 9 g more. Two 75s make 150, and seven 21s make 147, exactly 3 g less. So two 75 g weights balance seven 21 g weights and 3 g of salt.

Fig. 2. Two 75 g weights balance seven 21 g weights and 3 g of salt: 150 − 147 = 3.

Two tries were enough here, but the number of tries grows with the weights. For weights of 1836311903 g and 1134903170 g the GCD is 1, and the first count that works is 433494437 of the first weight, against 701408733 of the second. Trying counts one by one takes about 4 \cdot 10^8 tries, several seconds by the rule of thumb of about 10^8 simple steps per second, while the Euclidean algorithm finds their GCD in 44 steps.

So the question is whether the steps that found 3 can also say how to weigh it. The Euclidean algorithm on 75 and 21 makes the remainders 12, 9 and 3, and each of them can be weighed. Since 75 = 3 · 21 + 12, 12 g is one 75 against three 21s. Next, 21 = 1 · 12 + 9, so 9 g is one 21 against one 12. There is no 12 g weight, so replace 12 by its weighing: 9 = 21 − 12 = 21 − (75 − 3 · 21) = 4 · 21 − 75. So 9 g is four 21s against one 75 (84 − 75 = 9).

Each new remainder is weighed the same way: take the weighing of the larger number of the pair and remove the weighing of the smaller number as many times as the quotient says. Write each weighing as two counts, the 75s and the 21s, with a minus for weights on the same pan as the load. The weight 75 itself is (1, 0) and 21 is (0, 1). Then 12 is (1, 0) − 3 · (0, 1) = (1, −3), and 9 is (0, 1) − (1, −3) = (−1, 4).

A table keeps each number with its counts. Each new row is the row two above it minus the quotient times the row just above it.

Fig. 3. The numbers of the Euclidean algorithm on 75 and 21, each with its counts of 75s and 21s.

This is the extended Euclidean algorithm: the Euclidean algorithm that keeps, next to each number, how many of each original number make it. When the remainder reaches 0, the GCD comes with its two counts. The fact that the GCD can always be made this way is called Bezout's identity, and the two counts are the Bezout coefficients.

The row for 3 is the answer: 3 = 2 · 75 − 7 · 21. On the way, the table also found a weighing for every remainder. The row for 0, seven 75s against twenty-five 21s, is a balanced set: 525 = 525. Adding it to a weighing doesn't change what it weighs, so 3 g can be weighed in infinitely many ways. For example, eighteen 21s against five 75s also give 378 − 375 = 3.

Equations a·x + b·y = c

To weigh 12 g, take the weighing for 3 g four times: eight 75s against twenty-eight 21s, and 600 − 588 = 12. To weigh 5 g is impossible, because every weighing is a whole number of 3 g. In general, the equation a x + b y = c has integer solutions exactly when g = \gcd(a, b) divides $c$. Then the Bezout coefficients multiplied by $c/g$ give one solution.

Adding or removing the balanced set gives the others. Seven 75s against twenty-five 21s is the smallest balanced set (525 is the least common multiple of 75 and 21), so the count of 75s changes by 7 = 21 / 3 from one solution to the next, and the count of 21s by 25 = 75 / 3. All solutions of 75x + 21y = 12 are x = 8 + 7k, y = -28 - 25k for integer $k$. In general, from one solution $x$, $y$ all of them are x + k \cdot \frac{b}{g}, y - k \cdot \frac{a}{g}.

Fig. 4. The solutions of 75x + 21y = 12 are 7 apart. Counting down from 8 by 7s reaches the smallest non-negative one, 1.

To get the smallest non-negative $x$, take the remainder of $x$ divided by $b/g$: here 8 divided by 7 leaves 1. Then y = -3, and one 75 against three 21s weighs 75 − 63 = 12 g.

Some problems need both x \ge 0 and y \ge 0. As $k$ grows, $x$ grows and $y$ falls, so the good values of $k$ form a range, and it is enough to check its ends. For 75x + 21y = 12 the range is empty: at the smallest non-negative $x$, y = -3 is already negative.

Modular inverse

The inverse of $a$ modulo $m$ is the number $x$ such that a \cdot x leaves remainder 1 when divided by $m$. Dividing by $a$ modulo $m$ means multiplying by its inverse. To find the inverse of 7 modulo 25, run the algorithm on 7 and 25. The remainders are 4, 3 and 1, and the counts of 1 give 7 \cdot (-7) + 25 \cdot 2 = 1. So 7 \cdot (-7) is 1 more than a multiple of 25, and in modular arithmetic that makes $-7$ an inverse of 7 modulo 25. The usual answer is between 0 and 24: -7 + 25 = 18, and indeed 7 \cdot 18 = 126 = 5 \cdot 25 + 1.

This works for any $m$, prime or not, but only when \gcd(a, m) = 1. For example, 21 has no inverse modulo 75: 21x - 75y is always a multiple of 3, so it is never 1.

Looking for $x$ such that a \cdot x leaves remainder $c$ modulo $m$ is the same equation a x + m y = c from the previous subsection.

Implementation

In code, the pair is a and b, as in the Euclidean algorithm. The counts of a are xa and ya, and the counts of b are xb and yb. On the input 75 21, a always equals 75 * xa + 21 * ya, and b equals 75 * xb + 21 * yb. At the start a is 75 with counts 1 and 0, and b is 21 with counts 0 and 1. The quotient is a / b, integer division that drops the remainder.

1. Take the remainder: q = a / b, r = a % b.

2. Update the counts: the counts of r are xa - q * xb and ya - q * yb.

3. Shift the pair: b and its counts move into a, and r and its counts move into b.

4. Stop at zero: when b == 0, a is the GCD, and xa, ya are its Bezout coefficients.

For positive a and b the Bezout coefficients are small: |x| <= b / g and |y| <= a / g. They are still kept in long long, because the code that uses them multiplies them by a or b.

The function ext_gcd(a, b, x, y) returns the GCD and writes the Bezout coefficients into x and y. They are passed by reference (&), so the function changes the caller's variables. The program reads two positive numbers and prints the GCD and the Bezout coefficients. The input for the weights is:

75 21

#include <iostream>
using namespace std;

// returns gcd(a, b) and writes x, y with a * x + b * y == gcd(a, b)
long long ext_gcd(long long a, long long b, long long& x, long long& y) {
    long long xa = 1, ya = 0;   // a = 1 · a + 0 · b
    long long xb = 0, yb = 1;   // b = 0 · a + 1 · b

    while (b != 0) {            // stop at zero
        // take the remainder
        long long q = a / b, r = a % b;
        // update the counts
        long long xr = xa - q * xb, yr = ya - q * yb;
        // shift the pair
        a = b;  xa = xb;  ya = yb;
        b = r;  xb = xr;  yb = yr;
    }

    x = xa;
    y = ya;
    return a;
}

int main() {
    long long a, b;
    cin >> a >> b;
    long long x, y;
    long long g = ext_gcd(a, b, x, y);
    cout << g << ' ' << x << ' ' << y << '\n';
}

If a or b can be negative, run the function on their absolute values and flip the sign of the coefficient whose number was negative.

The loop makes the same steps as the Euclidean algorithm, so it takes O(\log \min(a, b)) time and $O(1)$ memory.

Solving a·x + b·y = c

The function solve finds the solution with the smallest non-negative x, or returns false if there is none. It scales the Bezout coefficients by c / g and then takes x modulo b / g. The y that goes with that x is (c - a * x) / b, and the division is exact, because c - a * x equals b * y. solve expects positive a and b. If one of them is 0, the equation has one unknown and needs a separate check. With a, b and |c| up to about 10^9, every product fits in long long.

// finds x, y with a * x + b * y == c and the smallest x >= 0; false if there are none
bool solve(long long a, long long b, long long c, long long& x, long long& y) {
    long long g = ext_gcd(a, b, x, y);
    if (c % g != 0) return false;
    x *= c / g;                     // one solution
    long long d = b / g;
    x = (x % d + d) % d;            // the smallest x >= 0
    y = (c - a * x) / b;
    return true;
}

On 75 21 12 it gives x = 1, y = -3.

Modular inverse

The function inverse(a, m) returns the inverse of a modulo m between 0 and m - 1, or −1 if there is none. The algorithm can return a negative x, and in C++ -7 % 25 is -7, so (x % m + m) % m moves it into range.

// the inverse of a modulo m, or -1 if gcd(a, m) != 1
long long inverse(long long a, long long m) {
    long long x, y;
    if (ext_gcd(a, m, x, y) != 1) return -1;
    return (x % m + m) % m;
}

inverse(7, 25) returns 18. To divide c by a modulo m, compute c % m * inverse(a, m) % m. For m up to 2 \cdot 10^9 that product can reach 4 \cdot 10^{18}, so it needs long long.

The recursive form

The algorithm is also often written as a recursive function, like the short gcd. It calls itself on the pair b, a % b, which gives x1 and y1 with b * x1 + (a % b) * y1 == g. The counts for a and b are built from them on the way back. When b == 0, the answer is a = 1 * a + 0 * b.

long long ext_gcd(long long a, long long b, long long& x, long long& y) {
    if (b == 0) { x = 1; y = 0; return a; }
    long long x1, y1;
    long long g = ext_gcd(b, a % b, x1, y1);
    x = y1;
    y = x1 - (a / b) * y1;
    return g;
}

The calls go down to (3, 0), where b = 0: it returns 3 with x = 1, y = 0, because 3 = 1 · 3 + 0 · 0. Under each call, the grey number is its quotient q = a / b and the green pair is its x, y.

(9, 3) gets x1 = 1, y1 = 0 from (3, 0). x = y1 = 0, y = x1 − q · y1 = 1 − 3 · 0 = 1, and 9 · 0 + 3 · 1 = 3.

(12, 9) gets x1 = 0, y1 = 1 from (9, 3). x = y1 = 1, y = x1 − q · y1 = 0 − 1 · 1 = −1, and 12 − 9 = 3.

(21, 12) gets x1 = 1, y1 = −1 from (12, 9). x = y1 = −1, y = x1 − q · y1 = 1 − 1 · (−1) = 2, and −21 + 24 = 3.

(75, 21) gets x1 = −1, y1 = 2 from (21, 12). x = y1 = 2, y = x1 − q · y1 = −1 − 3 · 2 = −7, and 150 − 147 = 3.

The first call (75, 21) returns the answer: 3 = 2 · 75 − 7 · 21, the same counts as in the table.

1/6

Fig. 5. Step through the recursive form with the buttons or the arrow keys. Each call takes x1, y1 from the call on its right, swaps them and subtracts the quotient times y1.

It can replace the loop version above and gives the same Bezout coefficients as the table: 2 and −7 on 75 21. The recursion is as deep as the number of steps, so it uses O(\log \min(a, b)) memory.

Correctness proof

The invariant: each number of the pair equals A \cdot x + B \cdot y for its counts $x$, $y$, where $A$ and $B$ are the input numbers. At the start it holds: A = A \cdot 1 + B \cdot 0 and B = A \cdot 0 + B \cdot 1. A step makes the new number r = a - q \cdot b and its counts by the same subtraction, x_a - q \cdot x_b and y_a - q \cdot y_b. If a = A x_a + B y_a and b = A x_b + B y_b, then subtracting $q$ times the second from the first gives r = A (x_a - q x_b) + B (y_a - q y_b). So the invariant holds for $r$ too. The pair goes through the same numbers as in the Euclidean algorithm, so it ends with the GCD, and the GCD carries its counts.

Fig. 6. The number and its counts are computed by the same subtraction, so the new row passes the check whenever the two rows above it do.

$g$ is also the smallest positive value of A x + B y. Every such value is a multiple of $g$, because $g$ divides $A$ and $B$, and the algorithm reaches $g$ itself.

Any two solutions of A x + B y = c differ by a balanced set, a solution of A x + B y = 0. If A x + B y = c and A x' + B y' = c, then A (x - x') = B (y' - y). Divide both sides by $g$: \frac{A}{g} (x - x') = \frac{B}{g} (y' - y). The GCD of \frac{A}{g} and \frac{B}{g} is 1, so by Bezout's identity \frac{A}{g} u + \frac{B}{g} v = 1 for some integers $u$, $v$. Multiply by x - x': x - x' = \frac{A}{g} (x - x') \cdot u + \frac{B}{g} (x - x') \cdot v, and both terms are multiples of \frac{B}{g}. So the values of $x$ are exactly \frac{B}{g} apart, and each $x$ fixes its $y$.

Common mistakes

• Using the Fermat formula for the inverse, a^{m-2} modulo $m$, when $m$ isn't prime: it holds only for a prime $m$. Modulo 25, 2^{23} leaves 8, and 2 \cdot 8 = 16, not 1. The real inverse is 13, since 2 \cdot 13 = 26.

• Shifting solutions by b instead of b / g. For 75x + 21y = 12 scaling gives x = 8, and 8 % 21 keeps 8, but the smallest non-negative $x$ is 1. Take x modulo b / g, here 7.

• Overwriting x before using it in the recursive form. With the same variables, x = y; y = x - (a / b) * y; uses the new x in the second line and returns 0 and 0 on 75 21. Keep x1 and y1 from the call, or use a temporary.

Similar algorithms

The Euclidean algorithm is this algorithm without the counts. Use it when only the GCD is needed, for example to check whether a x + b y = c has a solution at all.

For a prime modulus, the modular inverse is often computed as a^{m-2} with binary exponentiation, which is shorter to write. Both take O(\log m) steps. The modular inverse article also builds the inverses of all numbers from 1 to $n$ with a linear-time table.

The Chinese remainder theorem finds a number with given remainders modulo several numbers. Merging two such conditions means solving an equation a x + b y = c, and it uses this algorithm inside.

Practice