# 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])))
```

### 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)
```