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 gcd⁑(a,b)\gcd(a,b) 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 gcd⁑\gcd. 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 x∣yx | y, then the GCD of xx and yy is xx; Similarly, the GCD function is β€œsymmetric”, meaning gcd⁑(a,b)=gcd⁑(b,a)\gcd(a,b) = \gcd(b,a).

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: gcd⁑(a,0)=a\gcd(a,0) = a and gcd⁑(a,b)=gcd⁑(b,a%b)\gcd(a, b) = \gcd(b, a \% b) 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 nn is nonzero, d∣nd|n, and d∣(m%n)d|(m \% n), then d∣md|m;
  • If nn is nonzero, d∣nd|n, and d∣md|m, then d∣(m%n)d|(m \% n).

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 (nn is nonzero and) gcd⁑(n,m%n)=d\gcd(n, m \% n) = d, then gcd⁑(m,n)=d\gcd(m, n) = d, 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 m%nm \% n is less than nn. 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))