The Euclid Code
| Table of contents |
Greatest common divisor(gcd), Least common multiple(lcm), Modular inverses, Extended Euclidean Algorithm
All of these functions from number theory can be obtained with the extended euclidean algorithm (EEA). This algorithm, given integers a and b, solves the following for x,y,and g.
One we have x,y,g then gcd(a,b) = g, and lcm(a,b) = a * b / g.
In the case g=1, notice ax = 1 (mod b) and x is the modular inverse of a. The code below makes x positive, when possible, for convenience.
Following is C++ code written by landen. This is extracted from a larger program which uses an arbitrary precision integer class Big in MIRACL from Shamus Software. You can typdef in whatever type you like.
void ExtendedEuclideanKB(const Big &a, const Big &b, Big &x, Big &y, Big &g)
{
// solves ax+by = g. for x,y,g. g= gcd(a,b). x is positive if possible
// if g is 1 then ax = 1 (mod b) so x is multiplicative inverse of a
// Big can me any signed integer type you like
// Uses the iterative method of Knuth with the speedup of Bradley. Knuth, Seminumerical Algorithms page 325
// of seminumerical algorithms
Big t1,t2,t3,u1,u2,u3,v1,v2,v3,q; //Knuth notation
const Big& u = a, v = b;
u1 = 1; u3 = u; v1 = 0; v3 = v; // initialize
while (v3 != 0) {
q = u3/v3;
t1 = u1 - v1*q; t3 = u3 - v3*q;
u1 = v1; u3 = v3; //
// u1.swap(v1); u3.swap(v3); // use this line instead of preceding if you have swap
v1 = t1; v3 = t3;
// v1.swap(t1); v3.swap(t3); // use this line instead of preceding if you have swap
}
if (b == 0) { //special case
g = u;
x = 1;
y = 0;
return;
}
else
u2 = (u3 - a*u1)/v;
// end of Knuth's algorithm. rest is just to fix up signs
if (u3 >= 0)
{ x = u1; y=u2; g=u3;
}
else {
x = -u1; y=-u2; g=-u3;
}// cosmetic to force g positive
if(x < 0) { // this is just cosmetic to get x > 0; y can sometimes be negative
if (g <= 1) {
if (b > 0) {
x = x+b;
y = y-a;
}
else {
x = x-b;
y = y+a;
}
}
else {
if (b > 0) {
x = x+b/g;
y = y-a/g;
}
else {
x = x-b/g;
y = y+a/g;
}
}
}
}
Maxima Code for GCD with Gaussian integers
gintp(z):= integerp(realpart(z))and integerp(imagpart(z)); /* checks to see if z is a gaussian integer */
grem(z1,z2):= /* finds remainder on division of two gaussian integers, |r|<|z2| */
block([q,qi,r],
if not(gintp(z1) and gintp(z2)) then return("arguments to grem must be gaussian integers"),
q:rectform(z1/z2),
qi:floor(realpart(q)+1/2)+%i*floor(imagpart(q)+1/2),
r:rectform(z1-qi*z2),
return(r)
);
canonical_ans(a) := /* puts gcd answer in a unique canonical form in the first quadrant */
block(
if not (gintp(a)) then return ("canonical_ans must be called with a gaussian integer argument"),
if realpart(a)=0 then a:%i*a,
if realpart(a)<0 then a:-a,
if imagpart(a)<0 then a:%i*a,
return(rectform(a))
);
ggcd(z1,z2):= /* finds gcd in gaussian integers, puts the answer in a unique canonical form */
block([ans,k],
if not(gintp(z1) and gintp(z2)) then return("arguments to ggcd must be gaussian integers"),
if z2=0 then return(canonical_ans(z1)),
k:0, while true do ( k:k+1, if k > 1000000 then return(ans:"probably looping"), if (z1:grem(z1,z2))=0 then (ans:z2,return(z2)), if (z2:grem(z2,z1))=0 then (ans:z1,return(z1)) ), if not(gintp(ans)) then return (ans), return(canonical_ans(ans)) );
Fun Code for the mIRC client
landen found this script very hard to write. Put the following lines in a file called euclid.ini and put it in your MIRC directory. Then /euclid n1 n2 will compute the EEA.
/euclid {
; Let's check they typed in 2 parameters exactly
/echo -a given: Euclid $1-
; start of checking the input parameters for problems
if ($1 == $null) goto needhelp
if ($2 == $null) goto needhelp
if ($3 != $null) goto needhelp
; But were both parameters integers? If they were then +0 will do nothing
%a = $1 + 0
%b = $2 + 0
if (%a != $1) goto needhelp
if (%b != $2) goto needhelp
%a = $abs(%a)
%b = $abs(%b)
%frac = %a % 1
if (%frac != 0) goto needhelp
%frac = %b % 1
if (%frac != 0) goto needhelp
; Did they type in a 0 parameter or two?
if (%a != 0) {
if (%b == 0) {
%re = %a
%xe = 1
%ye = 0
goto odd_success
}
}
if (%a == 0) {
if (%b != 0) {
%re = %b
%xe = 0
%ye = 1
goto odd_success
}
else goto needhelp
}
; End of stuff for 0 parameters
; End of input check
; even numbered row in a two row table to start
%k = 0
%ae = %a
%be = %b
%qe = 0
%re = %ae
%xe = 1
%ye = 0
;//echo row: %k %ae %qe %be %re %xe %ye
; odd numbered row in the two row table to start
%k = 1
%ao = %be
%bo = %re
%ro = %ao % %bo
%qo = %ao - %ro
%qo = %qo / %bo
%temp = %qo * %xe
%xo = -1 * %temp
%temp = %qo * %ye
%yo = 1 - %temp
;//echo row: %k %ao %qo %bo %ro %xo %yo
if (%ro == 0) goto odd_success
:loop
; compute a new even row
%k = %k + 1
%ae = %bo
%be = %ro
%re = %ae % %be
%qe = %ae - %re
%qe = %qe / %be
%temp = %qe * %xo
%xe = %xe - %temp
%temp = %qe * %yo
%ye = %ye - %temp
;//echo row: %k %ae %qe %be %re %xe %ye
if (%re == 0) goto even_success
; compute a new odd row
%k = %k + 1
%ao = %be
%bo = %re
%ro = %ao % %bo
%qo = %ao - %ro
%qo = %qo / %bo
%temp = %qo * %xe
%xo = %xo - %temp
%temp = %qo * %ye
%yo = %yo - %temp
;//echo row: %k %ao %qo %bo %ro %xo %yo
if (%ro == 0) goto odd_success
else goto loop
:even_success
; use the previous row in the table
%r = %ro
%x = %xo
%y = %yo
goto Tada
halt
:odd_success
; use the previous row in the table
%r = %re
%x = %xe
%y = %ye
goto Tada
halt
:Tada
if (%r < 0) {
%r = -1 * %r
%x = -1 * %x
%y = -1 * %y
}
if (%x < 0) {
%temp = %b / %r
%x = %x + %temp
%temp = %a / %r
%y = %y - %temp
}
//echo ans: %a * %x + %b * %y = %r
halt
:needhelp
/echo OOps! Use exactly two integer arguments. One must be nonzero. Eg: /euclid -75 35
halt
}
Another mIRC script from toad_
alias -l eeuclidr {
return $eeuclidc($1,$$2)
}
alias -l eeuclidc {
var %a = $1,%b = $$2
if (%b == 0) return 1 0
else {
var %temp,%x,%y
%temp = $eeuclidr(%b , $calc(%a % %b))
%x = $gettok(%temp,1,32)
%y = $gettok(%temp,-1,32)
return %y $calc( %x - %y * ((%a - (%a % %b))/ %b ))
}
}
alias eeuclid {
if ($int($1) == $1) && ($int($$2) == $2) {
var %temp = $eeuclidc($1,$2), %gcd
%gcd = $abs($calc($1 * $gettok(%temp,1,32) + $2 * $gettok(%temp,2,32))))
echo -a GCD: %gcd = $1 * $gettok(%temp,1,32) + $2 * $gettok(%temp,2,32)
}
else echo -a Usage: /eeuclid a b where a and b are integers.
}