open import 1Lab.Prelude open import Data.Wellfounded.Properties open import Data.Wellfounded.Base open import Data.Nat.Properties open import Data.Nat.Divisible open import Data.Nat.DivMod open import Data.Nat.Order open import Data.Nat.Base open import Data.Sum.Base 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))