# Algorithm Implementation/Mathematics/Extended Euclidean algorithm

## 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, a):
x0, x1, y0, y1 = 1, 0, 0, 1
while a != 0:
q, b, a = b // a, a, b % a
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
```

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