Algorithm Implementation/Mathematics/Extended Euclidean algorithm
 | Wikipedia has related information at Extended Euclidean algorithm. |
C
[edit | edit source]Iterative algorithm
[edit | edit source]#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; for (;;) { q = a / b; a = a % b; aa[0] = aa[0] - q*aa[1]; bb[0] = bb[0] - q*bb[1]; if (!a) { 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
[edit | edit source]//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
[edit | edit source]Iterative algorithm
[edit | edit source]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
[edit | edit source]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
[edit | edit source]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
[edit | edit source]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
[edit | edit source]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
[edit | edit source]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
[edit | edit source]Extended
[edit | edit source](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
[edit | edit source]Unextended
[edit | edit source]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
[edit | edit source]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)