## arithmetic/src/algebra/misc/GCD.mg

```GENERIC MODULE GCD(R);
```
Arithmetic for Modula-3, see doc for details
```
IMPORT Arithmetic AS Arith;

<* UNUSED *>
CONST
Module = "GCD.";

PROCEDURE GCD (u, v: T; ): T =
(* returns the greatest common denominator for u and v. *)
(* use Euclid's algorithm *)
(* Arith.ErrorDivisionByZero cannot occur *)
<* FATAL Arith.Error *>
BEGIN
WHILE NOT R.IsZero(v) DO
WITH w = R.Mod(u, v) DO u := v; v := w; END;
END;
(*
WHILE v#0 DO
w:=u MOD v;
u:=v;
v:=w;
END;
*)
RETURN u;
END GCD;

PROCEDURE LCM (u, v: T; ): T =
BEGIN
TRY
RETURN R.Mul(R.Div(u, GCD(u, v)), v);
EXCEPT
Arith.Error (err) =>
<* ASSERT NOT ISTYPE(err, Arith.ErrorIndivisible) *>
RETURN R.Zero;
END;
END LCM;

PROCEDURE BezoutGCD
(u, v: T; VAR (*OUT*) c: ARRAY [0 .. 1], [0 .. 1] OF T; ): T =
(*
/ u \ - / q0 1 \ / q1 1 \ ... / gcd(u,v) \
\ v / - \ 1  0 / \ 1  0 /     \    0     /

... / 0  1  \ / 0  1  \ / u \ - / gcd(u,v) \
\ 1 -q1 / \ 1 -q0 / \ v / - \    0     /

/ a 1 \ / 0  1 \ - / 1 0 \
\ 1 0 / \ 1 -a / - \ 0 1 /
*)

(* Arith.ErrorDivisionByZero cannot occur *)
<* FATAL Arith.Error *>
VAR
next: ARRAY [0 .. 1] OF T;
BEGIN
c[0, 0] := R.One;
c[0, 1] := R.Zero;
c[1, 0] := R.Zero;
c[1, 1] := R.One;
WHILE NOT R.IsZero(v) DO
WITH qr = R.DivMod(u, v) DO
u := v;
v := qr.rem;
next[0] := R.Sub(c[0, 0], R.Mul(c[0, 0], qr.quot));
next[1] := R.Sub(c[0, 1], R.Mul(c[0, 1], qr.quot));
c[0] := c[1];
c[1] := next;
END;
END;
RETURN u;
END BezoutGCD;

PROCEDURE Bezout (u, v, w: T; VAR (*OUT*) c: ARRAY [0 .. 1], [0 .. 1] OF T)
RAISES {Arith.Error} =
(*
The former routine gives us c00, c01, c10, c11 for each u,v with
c00*u + c01*v = gcd(u,v)
c10*u + c11*v = 0
Now let
k = w / gcd(u,v)  (if this is indivisible the Bezout equation is not solvable)
then a solution is
k*c00*u + k*c01*v  = w
but the coefficients k*c00 and k*c01 are to large and can be reduced:
(k*c00-f*c10)*u + (k*c01-f*c11)*v = w
By dividing k*c00 by c10 with remainder one obtains the smallest possible coefficient for u and the coefficent for v cannot not to large since it is limitted by w and u and its coefficient.
*)
VAR
gcd            := BezoutGCD(u, v, c);
k              := R.Div(w, gcd);
fr : R.QuotRem;
BEGIN
c[0, 0] := R.Mul(c[0, 0], k);
c[0, 1] := R.Mul(c[0, 1], k);
(* reduce c[0,0] and c[0,1]*)
fr := R.DivMod(c[0, 0], c[1, 0]);
c[0, 0] := fr.rem;
c[0, 1] := R.Sub(c[0, 1], R.Mul(fr.quot, c[1, 1]));
END Bezout;

PROCEDURE MACDecompose (u, v: T; VAR (*OUT*) mac: MAC; ): T =
(* Arith.ErrorDivisionByZero cannot occur *)
<* FATAL Arith.Error *>
BEGIN
mac := NIL;
WHILE NOT R.IsZero(v) DO
WITH qr     = R.DivMod(u, v),
newmac = NEW(MAC, factor := qr.quot, next := mac) DO
mac := newmac;
u := v;
v := qr.rem;
END;
END;
RETURN u;
END MACDecompose;

BEGIN
END GCD.
```
```

```