Algorithm Implementation/Mathematics/Extended Euclidean algorithm

C

Iterative algorithm

#include <stdio.h>
#include <stdlib.h> // strtol() function

void xgcd(long *result, long a, long b){
    long aa[2]={1,0}, bb[2]={0,1}, q;
    while(1) {
        q = a / b; a = a % b;
        aa[0] = aa[0] - q*aa[1];  bb[0] = bb[0] - q*bb[1];
        if (a == 0) {
            result[0] = b; result[1] = aa[1]; result[2] = bb[1];
            return;
        };
        q = b / a; b = b % a;
        aa[1] = aa[1] - q*aa[0];  bb[1] = bb[1] - q*bb[0];
        if (b == 0) {
            result[0] = a; result[1] = aa[0]; result[2] = bb[0];
            return;
        };
    };
}
 
int main(int argc, char** argv){
    long a, b, c[3];
    char *end;
    if(argc < 3) {printf("Usage: %s a b  (compute xgcd(a,b))\n",argv[0]); exit(-1);};
    printf("long int is %lu bytes (or %lu bits)\n",sizeof(long),8*sizeof(long));
    a = strtol(argv[1], &end, 10);
    b = strtol(argv[2], &end, 10);
    printf("xgcd(%ld,%ld) = ",a,b);
    xgcd(c,a,b);
    printf("%ld = %ld*%ld + %ld*%ld\n",c[0],c[1],a,c[2],b);
return 0;
}

Recursive algorithm

//https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm
//clrs p.937

#include <stdio.h>
#include <stdlib.h>

int mod (int a, int b){
	return a %b; 
}	

int *extendedEuclid (int a, int b){
	int *dxy = (int *)malloc(sizeof(int) *3);

	if (b ==0){
		dxy[0] =a; dxy[1] =1; dxy[2] =0;
		
		return dxy;
	}
	else{
		int t, t2;
		dxy = extendedEuclid(b, mod(a, b));
		t =dxy[1];
		t2 =dxy[2];
		dxy[1] =dxy[2];
		dxy[2] = t - a/b *t2;

		return dxy;
	}
}

int main(void)
{
	int a =99, b =78;
	int *ptr;

	ptr =extendedEuclid (a, b);
	printf("%d = %d * %d + %d * %d \n",ptr[0], a, ptr[1], b, ptr[2] );

	return 0;		
}

Java

Iterative algorithm

public class xgcd {
    public static long[] xgcd(long a, long b){
		long[] retvals={0,0,0};
		long aa[]={1,0}, bb[]={0,1}, q=0;
		while(true) {
    		q = a / b; a = a % b;
	    	aa[0] = aa[0] - q*aa[1];  bb[0] = bb[0] - q*bb[1];
    		if (a == 0) {
      			retvals[0] = b; retvals[1] = aa[1]; retvals[2] = bb[1];
      			return retvals;
	    	};
    		q = b / a; b = b % a;
    		aa[1] = aa[1] - q*aa[0];  bb[1] = bb[1] - q*bb[0];
	    	if (b == 0) {
    	  		retvals[0] = a; retvals[1] = aa[0]; retvals[2] = bb[0];
      			return retvals;
	    	};
    	}
    }
    
    public static void main(String[] args){
	if(args.length < 2){
	    System.out.println("Usage: xgcd(a,b)"); System.exit(-1);};
	long a = Integer.parseInt(args[0]);
	long b = Integer.parseInt(args[1]);
	long[] retvals;
	retvals = xgcd(a,b);
	System.out.println("xgcd("+args[0]+", "+args[1]+") = ["
		+String.valueOf(retvals[0])+", "+
		String.valueOf(retvals[1])+", "+
		String.valueOf(retvals[2])+"]");
	System.out.println("   "+String.valueOf(retvals[0])+" = "
		+args[0]+"*("+String.valueOf(retvals[1])+") + "
		+args[1]+"*("+String.valueOf(retvals[2])+")");
    };
}

Iterative with BigInteger

import java.math.BigInteger;

public class XGCD {

    // Assumes a and b are positive
    public static BigInteger[] xgcd(BigInteger a, BigInteger b) {
    	BigInteger x = a, y=b;
    	BigInteger[] qrem = new BigInteger[2];
    	BigInteger[] result = new BigInteger[3];
    	BigInteger x0 = BigInteger.ONE, x1 = BigInteger.ZERO;
    	BigInteger y0 = BigInteger.ZERO, y1 = BigInteger.ONE;
    	
    	while (true) {
    	    qrem = x.divideAndRemainder(y);
    	    x = qrem[1];
    	    x0 = x0.subtract(y0.multiply(qrem[0]));
    	    x1 = x1.subtract(y1.multiply(qrem[0]));
    	    if (x.equals(BigInteger.ZERO)) {
    	        result[0] = y;
    	        result[1] = y0;
    	        result[2] = y1;
    	        return result;
    	    }
    	    
    	    qrem = y.divideAndRemainder(x);
    	    y = qrem[1];
    	    y0 = y0.subtract(x0.multiply(qrem[0]));
    	    y1 = y1.subtract(x1.multiply(qrem[0]));
    	    if (y.equals(BigInteger.ZERO)) {
    	        result[0] = x;
    	        result[1] = x0;
    	        result[2] = x1;
    	        return result;
    	    }
    	}
    }

    public static void main(String[] args) {
    	BigInteger a = new BigInteger(args[0]);
    	BigInteger b = new BigInteger(args[1]);
    	BigInteger[] result = new BigInteger[3];
    	System.out.println("a = " + a.toString() + ";   b = " + b.toString());
    	result = xgcd(a,b);
    	System.out.println("xgcd = " + result[0].toString() + 
    			   "  [" + result[1].toString() + ", " + result[2].toString() + "]");
    }
}

Python

Both functions take positive integers a, b as input, and return a triple (g, x, y), such that ax + by = g = gcd(a, b).

Iterative algorithm

from typing import Tuple
def xgcd(a: int, b: int) -> Tuple[int, int, int]:
    """return (g, x, y) such that a*x + b*y = g = gcd(a, b)"""
    x0, x1, y0, y1 = 0, 1, 1, 0
    while a != 0:
        (q, a), b = divmod(b, a), a
        y0, y1 = y1, y0 - q * y1
        x0, x1 = x1, x0 - q * x1
    return b, x0, y0

Recursive algorithm

import sys
from typing import Tuple
sys.setrecursionlimit(1000000)  # long type, 32bit OS 4B, 64bit OS 8B (1bit for sign)

def egcd(a: int, b: int) -> Tuple[int, int, int]:
    """return (g, x, y) such that a*x + b*y = g = gcd(a, b)"""
    if a == 0:
        return (b, 0, 1)
    else:
        b_div_a, b_mod_a = divmod(b, a)
        g, x, y = egcd(b_mod_a, a)
        return (g, y - b_div_a * x, x)

Modular inverse

An application of extended GCD algorithm to finding modular inverses:

def modinv(a: int, b: int) -> int:
    """return x such that (x * a) % b == 1"""
    g, x, _ = egcd(a, b)
    if g != 1:
        raise Exception('gcd(a, b) != 1')
    return x % b

Clojure

Extended

(defn xgcd
  "Extended Euclidean Algorithm. Returns [gcd(a,b) x y] where ax + by = gcd(a,b)."
  [a b]
  (if (= a 0)
    [b 0 1]
    (let [[g x y] (xgcd (mod b a) a)]
      [g (- y (* (Math/floorDiv b a) x)) x])))

Haskell

Unextended

euclid :: Integral a => a -> a -> a
euclid 0 b = abs b
euclid a 0 = abs a
euclid a b = euclid b $ rem a b

Extended

eGCD :: Integer -> Integer -> (Integer,Integer,Integer)
eGCD 0 b = (b, 0, 1)
eGCD a b = let (g, s, t) = eGCD (b `mod` a) a
           in (g, t - (b `div` a) * s, s)