Extended Euclidean Algorithm

Extended Euclidean Algorithm is an extension of standard Euclidean Algorithm for finding the GCD of two integers a and b. It also calculates the values of two more integers x and y such that: ax + by = gcd(a,b); where typically either x or y is negative. This algorithm is generally used to find multiplicative inverse in a finite field, because, if ax + by = gcd(a,b) = 1, i.e. a and be are co-primes, then x is the modular multiplicative inverse of a modulo b, and similarly, y is the modular multiplicative inverse of b modulo a.

This method computes expressions of the form r_i = ax_i + by_i for the remainder in each step i of the Euclidean algorithm. Each modulus can be written in terms of the previous two remainders and their whole quotient as follows:

By substitution, this gives:

The first two values are the initial arguments to the algorithm:

r₁ = a = a(1) + b(0)

r₂ = b = a(0) + b(1)

The expression for the last non-zero remainder gives the desired results since this method computes every remainder in terms of a and b, as desired.

Example: Compute the GCD of 120 and 23. Or, more formally, compute: x, y, g for 120x + 23y = g; where x, y are two integers and g is the gcd of 120 and 23...

The computation proceeds as follows:

Initial values:
Step 1: Reminder = 120;
Combine terms: 120 = 120 x 1 + 23 x 0

Step 2: Reminder = 23;
Combine terms: 23 = 120 x 0 + 23 x 1

Iterative steps:
Step 3: Quotient = 120 / 23 = 5; Reminder = 120 % 23 = 5;
5 = 120 - 23 x 5
=> 5 = (120 x 1 + 23 x 0) - (120 x 0 + 23 x 1) x 5 ;[from Step 1 and 2]
=> 5 = 120 x 1 + 23 x -5

Step 4: Quotient = 23 / 5 = 4; Reminder = 23 % 5 = 3;
3 = 23 - 5 x 4
=> 3 = (120 x 0 + 23 x 1) - (120 x 1 + 23 x -5) x 4 ;[from Step 2 and 3]
=> 3 = 120 x -4 + 23 x 21

Step 5: Quotient = 5 / 3 = 1; Reminder = 5 % 3 = 2;
2 = 5 - 3 x 1
=> 2 = (120 x 1 + 23 x -5) - (120 x -4 + 23 x 21) x 1 ;[from Step 3 and 4]
=> 2 = 120 x 5 + 23 x -26

Step 6: Quotient = 3 / 2 = 1; Reminder = 3 % 2 = 1;
1 = 3 - 2 x 1
=> 1 = (120 x -4 + 23 x 21) - (120 x 5 + 23 x -26) x 1 ;[from Step 4 and 5]
=> 1 = 120 x -9 + 23 x 47
End of Algorithm.

The last line reads 1 = −9×120 + 47×23, which is the required solution: x = −9 and y = 47, and obviously g = gcd(120,23) = 1

This also means that −9 is the multiplicative inverse of 120 modulo 23, and that 47 is the multiplicative inverse of 23 modulo 120.

−9 × 120 ≡ 1 mod 23 and also 47 × 23 ≡ 1 mod 120.

Algorithm:
By routine algebra of expanding and grouping like terms (refer to the previous example), the following algorithm for iterative method is obtained:

1. Apply Euclidean algorithm, and let q_n(n starts from 1) be a finite list of quotients in the division.
2. Initialize x₀, x₁ as 1, 0, and y₀, y₁ as 0,1 respectively.
3. Then for each i so long as q_i is defined,
• Compute x_i+1= x_i-1- q_ix_i
• Compute y_i+1= y_i-1- q_iy_i
• Repeat the above after incrementing i by 1.
4. The answers are the second-to-last of x_n and y_n.

Sample Program:
Here is a sample program written in C++ which implements the Extended Euclidean Algorithm:

#include <stdio.h>
/*
Takes a, b as input.
Returns gcd(a, b).
Updates x, y via pointer reference.
*/
int Extended_Euclid(int A, int B, int *X, int *Y)
{
    int x, y, u, v, m, n, a, b, q, r;
 
    /* B = A(0) + B(1) */
    x = 0; y = 1;
 
    /* A = A(1) + B(0) */
    u = 1; v = 0;
 
    for (a = A, b = B; 0 != a; b = a, a = r, x = u, y = v, u = m, v = n)
    {
        /* b = aq + r and 0 <= r < a */
        q = b / a;
        r = b % a;
 
        /* r = Ax + By - aq = Ax + By - (Au + Bv)q = A(x - uq) + B(y - vq) */
        m = x - (u * q);
        n = y - (v * q); 
    }   
 
    /* Ax + By = gcd(A, B) */
    *X = x; *Y = y;
 
    return b;
}

int main()
{
    int a, b, x, y, g;
    scanf("%d %d", &a, &b);
    g = Extended_Euclid(a, b, &x, &y);
    printf("X = %d; Y = %d; G = %d\n", x, y, g);
    return 0;
}

Python implementation:

def Extended_Euclid(a, b):
    x, last_x = 0, 1
    y, last_y = 1, 0
 
    while b:
        quotient = a / b
        a, b = b, a % b
        x, last_x = last_x - quotient*x, x
        y, last_y = last_y - quotient*y, y
 
    return (last_x, last_y, a)

For complete reading, click here. It is also useful to have a look at the Chinese Remainder Theorem.

Hope this will help...

Largest impossible score

Problem:
If two possible points in a game is p and q, what is the largest impossible score of the game which can't be achieved? Or it can't be determined? In other words, for integers p, q ≥ 0, what is the maximum number m which can't be written as a linear combination of p, q, i.e. px + qy, where, x, y ≥ 0.

Solution:
The score can only be determined when p and q are co-primes, i.e. g = GCD(a,b) = 1. Otherwise, we can only make the multiples of g and the remaining numbers can't be made which means the maximum number can't be determined, or infinite! This also covers the case if exactly one of p or q is 0.

One other required condition is p,q > 1. If anyone of p or q is 1, there is no such score which can't be achieved, i.e. answer will be 0. If p and q both are 0, we can't even make any score...

Now, we can arrange the numbers from 0 to p*q-1 in a rectangular array as follows:

   0         1        .......       q-1
   q        q+1       .......       2q-1
  ...       ...       .......       ...
  ...       ...       .......       ...
(p-1)q   (p-1)(q-1)   .......       pq-1

For example, if p = 3 and q = 7, we have the following table:

(0)    1     2   (3)   4     5   (6)
 7     8    (9)   10   11   (12)  13
 14   (15)   16   17  (18)   19   20

We have circled the multiples of p in the above table and claim these observations:

• There will be q circled numbers. As we have exactly p numbers in each of q columns, so each column will have exactly one circled number, i.e. multiple of p, which implies, there will be a total of q circled number in this table.


• No two circled numbers are in same same column and there will be only one circled number in any column. For the same reason stated above, exactly one circled number in each column.


• Any circled number is possible. This is trivial, because, the circled numbers are multiple of one of the two given numbers, i.e. in the form p*x + q*0 or p*0 + q*x.


• Any number below the circled one in the same column is possible. We have seen any circled number is possible. A circled number is the first multiple of p in the respective column, i.e. n*p, where n ≥1, so the following numbers in the column will be: n*p + q, n*p + 2*q. n*p + 3*q, ... ... ... which are always possible.


• Any number above the circled one in the same column is impossible. This is also explained in the previous step. As any circled number in a column is the first one which is a multiple of one element, the number above this one will be neither the multiple of any of the two numbers. So it is impossible.


• Any number ≥ p*q is possible. As we complete p rows of the table with q columns, we already have one circled number in each column, so, any number > p*q-1 (i.e. the last element of the table) will be surely following any of the q circled numbers, so any number starting from the p+1^th row, starting from p*q will be possible.

Now considering all these observations, we can see that, largest impossible number is the number just above the last circled number. The last circled number is p*q-p, i.e. 21-3 = 18 and number above it is (p*q-p)-q, i.e. 18-7 = 11. So in general, largest impossible number is: p*q - p - q.

Conclusion:
If you are given p, q :: two positive numbers, the largest number, m, which can't be expressed as a linear combination of p, q i.e. in the form px + qy, where x and y are two non-negative numbers, is p * q - p - q, when p, q > 1 and GCD(p, q) is 1 i.e. p, q are co-primes.

Binary GCD algorithm

Here is a sample code for calculating GCD bitwise:

typedef unsigned int ui;

ui gcd(ui u, ui v)
{
    ui shift, diff;
    
    if (u == 0 || v == 0)
        return u | v;
    
    /* Let shift := lg K,
    where K is the greatest power of 2 dividing both u and v. */

    for (shift = 0; ((u | v) & 1) == 0; ++shift)
    {
        u >>= 1;
        v >>= 1;
    }
    
    while ((u & 1) == 0)
        u >>= 1;
        
    /* From here on, u is always odd. */

    do {
        while ((v & 1) == 0)
            v >>= 1;
            
        /* Now u and v are both odd, so diff(u, v) is even.
        Let u = min(u, v), v = diff(u, v)/2. */

        if (u < v)
            v -= u;
        else
        {
            diff = u - v;
            u = v;
            v = diff;
        }
        v >>= 1;
    
    } while (v != 0);
    
    return u << shift;
}

Actually, this algorithm outperforms the basic euclidean algorithm, so their asymptotic performance are same. According to Wikipedia, it is sometimes 60% faster than normal euclidean algorithm [proven].

But who needs it when he can do the following:

int gcd(int a, int b)
{
    while(b) b ^= a ^= b ^= a %= b;
    return a;
}

Stunned? oh that's normal...