module Data.Nat.Divisible.GCD where
Greatest common divisorsπ
The greatest common divisor of a pair of natural numbers is the largest number which divides them both. The definition we use is slightly unorthodox, since it requires only that the GCD be the greatest in the divisibility ordering, but the divisibility ordering is finer than the usual ordering, so this will turn out to imply they are greatest in magnitude, as well.
We start by defining what it means for a number to be a GCD β that is, we start by defining the graph of the function . Only after will we show that such a function exists.
record is-gcd (x y : Nat) (gcd : Nat) : Type where field gcd-β£l : gcd β£ x gcd-β£r : gcd β£ y greatest : β {g'} β g' β£ x β g' β£ y β g' β£ gcd open is-gcd is-gcd-is-prop : β {x y z} β is-prop (is-gcd x y z) is-gcd-is-prop p q i .gcd-β£l = β£-is-prop _ _ (p .gcd-β£l) (q .gcd-β£l) i is-gcd-is-prop p q i .gcd-β£r = β£-is-prop _ _ (p .gcd-β£r) (q .gcd-β£r) i is-gcd-is-prop p q i .greatest a b = β£-is-prop _ _ (p .greatest a b) (q .greatest a b) i instance H-Level-is-gcd : β {x y z n} β H-Level (is-gcd x y z) (suc n) H-Level-is-gcd = prop-instance is-gcd-is-prop GCD : Nat β Nat β Type GCD a b = Ξ£ _ (is-gcd a b) GCD-is-prop : β {a b} β is-prop (GCD a b) GCD-is-prop (_ , p) (_ , q) = Ξ£-prop-path! $ β£-antisym (q .greatest (p .gcd-β£l) (p .gcd-β£r)) (p .greatest (gcd-β£l q) (gcd-β£r q)) GCD-magnitude : β {x y g : Nat} β is-gcd x y (suc g) β β {g'} β g' β£ x β g' β£ y β g' β€ suc g GCD-magnitude gcd Ξ± Ξ² = mβ£snβmβ€sn (gcd .greatest Ξ± Ξ²)
The following observations may seem trivial, but it will come in super handy later: If , then the GCD of and is ; Similarly, the GCD function is βsymmetricβ, meaning .
dividesβGCD : β {x y} β x β£ y β GCD x y dividesβGCD {x} {y} w = x , gcd where gcd : is-gcd x y x gcd .gcd-β£l = β£-refl gcd .gcd-β£r = w gcd .greatest g'β£x g'β£y = g'β£x GCD-sym : β {x y} β GCD x y β GCD y x GCD-sym w .fst = w .fst GCD-sym w .snd .gcd-β£l = w .snd .gcd-β£r GCD-sym w .snd .gcd-β£r = w .snd .gcd-β£l GCD-sym w .snd .greatest g'β£y g'β£x = w .snd .greatest g'β£x g'β£y
Euclidβs algorithmπ
To compute greatest common divisors, we use the familiar (recursive) Euclidean algorithm: and otherwise. The difficulty comes in when we want not only to compute the numbers (in which case the definition above is perfectly cromulent), but also to show that they are greatest common divisors.
module Euclid where private variable n m k d d' : Nat
The base case can be established using our existing functions:
is-gcd-0 : is-gcd n 0 n is-gcd-0 .gcd-β£l = β£-refl is-gcd-0 .gcd-β£r = β£-zero is-gcd-0 .greatest g'β£n gβ£n = g'β£n
The inductive step is a bit more complicated. We first have to establish (using the most tedious of arithmetic properties) the following properties of the divisibility relation:
- If is nonzero, , and , then ;
- If is nonzero, , and , then .
These two facts, incarnated in the following helper lemmas lemβ
and lemβ
, will allow us to compute
the inductive step for GCD.
private lemβ : fibre (_* d) (suc n) β fibre (_* d) (m % suc n) β fibre (_* d) m lemβ {d = d} {n} {m} (cβ , p) (cβ , q) = dm.q * cβ + cβ , r where module dm = DivMod (divide-pos m (suc n)) renaming (quot to q ; rem to r) r = (dm.q * cβ + cβ) * d β‘β¨ *-distrib-+r (dm.q * cβ) _ _ β©β‘ dm.q * cβ * d + β cβ * d β β‘β¨ ap! q β©β‘ β dm.q * cβ * d β + dm.r β‘β¨ ap! (*-associative dm.q cβ d) β©β‘ dm.q * β cβ * d β + dm.r β‘β¨ ap! p β©β‘ dm.q * suc n + dm.r β‘Λβ¨ is-divmod m (suc n) β©β‘Λ m β lemβ : fibre (_* d) (suc n) β fibre (_* d) m β fibre (_* d) (m % suc n) lemβ {d = d} {n} {m} (cβ , p) (cβ , q) = cβ - dm.q * cβ , r where module dm = DivMod (divide-pos m (suc n)) renaming (quot to q ; rem to r) r = (cβ - dm.q * cβ) * d β‘β¨ monus-distribr cβ (dm.q * cβ) d β©β‘ cβ * d - β dm.q * cβ * d β β‘β¨ ap! (*-associative dm.q cβ d β ap (dm.q *_) p) β©β‘ β cβ * d β - dm.q * suc n β‘β¨ ap! (q β is-divmod m (suc n)) β©β‘ β dm.q * suc n + dm.r β - dm.q * suc n β‘β¨ ap! (+-commutative (dm.q * suc n) dm.r) β©β‘ (dm.r + dm.q * suc n) - dm.q * suc n β‘β¨ monus-cancelr dm.r 0 (dm.q * suc n) β©β‘ dm.r β
That step is put together here: It says that, if ( is nonzero and) , then , too β so that it will suffice to compute the former if we want the latter.
is-gcd-step : is-gcd (suc n) (m % suc n) d β is-gcd m (suc n) d is-gcd-step x .gcd-β£l = fibreββ£ (lemβ (β£βfibre (x .gcd-β£l)) (β£βfibre (x .gcd-β£r))) is-gcd-step x .gcd-β£r = x .gcd-β£l is-gcd-step x .greatest g'β£m g'β£sucn = x .greatest g'β£sucn (fibreββ£ (lemβ (β£βfibre g'β£sucn) (β£βfibre g'β£m)))
Actually putting this together is a bit indirect, since we are not performing structural induction: Agda canβt see that the argument is less than . We can, however, turn this into well-founded recursion, which demands only that the recursive call be less (in the sense of ) than the original input.
euclid-< : β y x β x < y β GCD y x euclid-< = Wf-induction _ <-wf _ Ξ» where x rec zero p β x , is-gcd-0 x rec (suc y) p β let (d , step) = rec (suc y) p (x % suc y) (x%y<y x (suc y)) in d , is-gcd-step step
With a handy wrapper to put the arguments in the order our induction
worker euclid-<
expects, we can
create the euclid
function, together with
its numerical component gcd
, and a proof that the graph of
gcd
is indeed is-gcd
.
euclid : β x y β GCD x y euclid x y with β€-split y x ... | inl y<x = euclid-< x y y<x ... | inr (inl x<y) = GCD-sym (euclid-< y x x<y) ... | inr (inr y=x) = dividesβGCD (subst (x β£_) (sym y=x) β£-refl) gcd : Nat β Nat β Nat gcd x y = Euclid.euclid x y .fst is-gcd-graphs-gcd : β {m n d} β is-gcd m n d β (gcd m n β‘ d) is-gcd-graphs-gcd {m = m} {n} {d} = prop-ext! (Ξ» x β ap fst $ GCD-is-prop (gcd m n , Euclid.euclid m n .snd) (d , x)) (Ξ» p β subst (is-gcd m n) p (Euclid.euclid m n .snd))