Open main menu

Wikibooks β

Algorithm Implementation/Mathematics/Extended Euclidean algorithm

< Algorithm Implementation‎ | Mathematics

Contents

CEdit

Iterative algorithmEdit

#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, r;
    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);
}

Recursive algorithmEdit

//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;		
}

JavaEdit

Iterative algorithmEdit

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 BigIntegerEdit

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() + "]");
    }
}

PythonEdit

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 algorithmEdit

def xgcd(b, n):
    x0, x1, y0, y1 = 1, 0, 0, 1
    while n != 0:
        q, b, n = b // n, n, b % n
        x0, x1 = x1, x0 - q * x1
        y0, y1 = y1, y0 - q * y1
    return  b, x0, y0

Recursive algorithmEdit

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

# return (g, x, y) a*x + b*y = gcd(x, y)
def egcd(a, b):
    if a == 0:
        return (b, 0, 1)
    else:
        g, x, y = egcd(b % a, a)
        return (g, y - (b // a) * x, x)

Modular inverseEdit

An application of extended GCD algorithm to finding modular inverses:

# x = mulinv(b) mod n, (x * b) % n == 1
def mulinv(b, n):
    g, x, _ = egcd(b, n)
    if g == 1:
        return x % n

HaskellEdit

UnextendedEdit

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

ExtendedEdit

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)