module UInt

import Base
import Nat
import UIntDefs
import UIntToFrom
import UIntLess
import UIntAdd
import UIntMonus
import UIntMult

/*
  UInt division, modulo, and divisibility.

  Division and modulo are defined by repeated subtraction over UInt.
  The proofs below establish their Nat correspondence, divisor
  bounds, inverse facts, and divisibility lemmas.
*/

// Quotient by repeated subtraction, returning zero for division by zero.
recfun operator /(n : UInt, m : UInt) -> UInt
  measure n of UInt
{
  if n < m then 0
  else if m = 0 then 0
  else 1 + ((n ∸ m) / m)
}
terminates {
  arbitrary n:UInt, m:UInt
  assume cond: not (n < m) and not (m = 0)
  suffices m + (n ∸ m) < m + n by uint_add_both_sides_of_less[m,n∸m,n]
  suffices n < m + n by {
    have m_n: m ≤ n by apply uint_not_less_implies_less_equal to conjunct 0 of cond
    replace apply uint_monus_add_identity[n,m] to m_n.
  }
  have m_pos: 0 < m by apply uint_not_zero_pos to conjunct 1 of cond
  conclude n < m + n by {
    replace uint_add_commute in
    apply uint_less_add_pos[n, m]
    to expand lit | fromNat in m_pos
  }
}

// Remainder after removing the quotient times the divisor.
fun operator % (n:UInt, m:UInt) {
  n ∸ (n / m) * m
}

// Dividing zero by a positive UInt yields zero.
theorem uint_zero_div: all x:UInt. if 0 < x then 0 / x = 0
proof
  arbitrary x:UInt
  assume zx: 0 < x
  expand operator/
  replace (apply eq_true to zx).
end

// Induction predicate used to prove quotient/remainder existence.
private fun UIntDivPred(n:UInt) {
  all m:UInt. if 0 < m then some r:UInt. (n / m) * m + r = n and r < m
}

// Every positive divisor has a quotient/remainder decomposition.
lemma uint_division: all n:UInt, m:UInt.
  if 0 < m
  then some r:UInt. (n/m)*m + r = n and r < m
proof
  have SI: all nn:UInt. if (all i:UInt. if i < nn then UIntDivPred(i))
                        then UIntDivPred(nn) by {
    arbitrary nn:UInt
    assume prem: all i:UInt. if i < nn then UIntDivPred(i)
    expand UIntDivPred
    arbitrary m:UInt
    assume m_pos: 0 < m
    switch nn < m {
      case true assume n_m {
        suffices some r:UInt. r = nn and r < m by {
          expand operator/
          simplify with n_m.
        }
        choose nn
        conclude nn = nn and nn < m by simplify with n_m.
      }
      case false assume not_n_m {
        have m_ne_z: not (m = 0) by apply uint_pos_not_zero to m_pos
        have m_le_n: m ≤ nn by apply uint_not_less_implies_less_equal to not_n_m
        have nm_n: nn ∸ m < nn by {
          have step1: (nn ∸ m) + m = nn by {
            replace uint_add_commute[nn ∸ m, m]
            apply uint_monus_add_identity[nn, m] to m_le_n
          }
          have step2: nn ∸ m < (nn ∸ m) + m
            by apply uint_less_add_pos[nn ∸ m, m] to expand lit | fromNat in m_pos
          replace step1 in step2
        }
        have IH_a: UIntDivPred(nn ∸ m) by apply prem to nm_n
        have IH_c: some r:UInt. ((nn ∸ m) / m) * m + r = (nn ∸ m) and r < m
          by apply (expand UIntDivPred in IH_a)[m] to m_pos
        obtain r0 where R: ((nn ∸ m) / m) * m + r0 = (nn ∸ m) and r0 < m from IH_c
        choose r0
        have conj0: (nn / m) * m + r0 = nn by {
          expand operator/
          replace (apply eq_false to not_n_m) | (apply eq_false to m_ne_z)
          show (1 + (nn ∸ m) / m) * m + r0 = nn
          replace uint_dist_mult_add_right[1, (nn ∸ m) / m, m]
          show m + ((nn ∸ m) / m) * m + r0 = nn
          replace conjunct 0 of R
          apply uint_monus_add_identity[nn, m] to m_le_n
        }
        conj0, conjunct 1 of R
      }
    }
  }
  arbitrary n:UInt
  have R: UIntDivPred(n) by apply uint_strong_induction[UIntDivPred, n] to SI
  expand UIntDivPred in R
end

// The quotient and remainder reconstruct the dividend.
theorem uint_div_mod: all n:UInt, m:UInt.
  if 0 < m
  then (n / m) * m + (n % m) = n
proof
  arbitrary n:UInt, m:UInt
  assume m_pos: 0 < m
  have ex: some r:UInt. (n/m)*m + r = n and r < m
    by apply uint_division[n, m] to m_pos
  obtain r where R: (n/m)*m + r = n and r < m from ex
  expand operator%
  define a = (n/m) * m
  have ar_n: a + r = n by R
  have a_le_n: a ≤ n by {
    have a_le_a_r: a ≤ a + r by uint_less_equal_add
    have eq_n_ar: n = a + r by symmetric (conjunct 0 of R)
    replace eq_n_ar
    a_le_a_r
  }
  have id: a + (n ∸ a) = n by apply uint_monus_add_identity[n, a] to a_le_n
  id
end

// The remainder is always smaller than a positive divisor.
theorem uint_mod_less_divisor: all n:UInt, m:UInt. if 0 < m then n % m < m
proof
  arbitrary n:UInt, m:UInt
  assume m_pos: 0 < m
  expand operator%
  obtain r where R: (n/m)*m + r = n and r < m
    from apply uint_division[n, m] to m_pos
  define a = (n/m)*m
  have ar_n: a + r = n by R
  have r_na: r = n ∸ a by {
    replace symmetric ar_n.
  }
  replace symmetric r_na
  conjunct 1 of R
end

// If the dividend is smaller than the divisor, the remainder is the dividend.
theorem uint_mod_small: all n:UInt, m:UInt. if n < m then n % m = n
proof
  arbitrary n:UInt, m:UInt
  assume n_m: n < m
  expand operator% | operator/
  replace (apply eq_true to n_m)
  uint_monus_zero[n]
end

// Divisibility: `a` divides `b` when `b` is a multiple of `a`.
fun divides(a : UInt, b : UInt) {
  some k:UInt. a * k = b
}

// Divisibility is reflexive.
theorem uint_divides_refl: all n:UInt. divides(n, n)
proof
  arbitrary n:UInt
  expand divides
  choose 1
  uint_mult_one[n]
end

// Every UInt divides zero.
theorem uint_divides_zero: all n:UInt. divides(n, 0)
proof
  arbitrary n:UInt
  expand divides
  choose 0
  uint_mult_zero[n]
end

// Divisibility is transitive.
theorem uint_divides_trans: all a:UInt, b:UInt, c:UInt.
  if divides(a, b) and divides(b, c) then divides(a, c)
proof
  arbitrary a:UInt, b:UInt, c:UInt
  assume prem
  obtain k where ak_b: a * k = b from expand divides in conjunct 0 of prem
  obtain l where bl_c: b * l = c from expand divides in conjunct 1 of prem
  expand divides
  choose k * l
  have eq1: (a * k) * l = c by replace symmetric ak_b in bl_c
  eq1
end

// Divisibility is closed under addition.
theorem uint_divides_add: all d:UInt, m:UInt, n:UInt.
  if divides(d, m) and divides(d, n) then divides(d, m + n)
proof
  arbitrary d:UInt, m:UInt, n:UInt
  assume prem
  obtain k where dk_m: d * k = m from expand divides in conjunct 0 of prem
  obtain l where dl_n: d * l = n from expand divides in conjunct 1 of prem
  expand divides
  choose k + l
  replace uint_dist_mult_add[d, k, l]
  replace dk_m | dl_n.
end

// Divisibility is closed under multiplying the divided value on the right.
theorem uint_divides_mult_right: all d:UInt, n:UInt, m:UInt.
  if divides(d, n) then divides(d, n * m)
proof
  arbitrary d:UInt, n:UInt, m:UInt
  assume dn
  obtain k where dk_n: d * k = n from expand divides in dn
  expand divides
  choose k * m
  replace dk_n.
end

// Divisibility is closed under multiplying the divided value on the left.
theorem uint_divides_mult_left: all d:UInt, n:UInt, m:UInt.
  if divides(d, n) then divides(d, m * n)
proof
  arbitrary d:UInt, n:UInt, m:UInt
  assume dn
  replace uint_mult_commute[m, n]
  apply uint_divides_mult_right[d, n, m] to dn
end

// A common divisor of two values also divides their truncated difference.
theorem uint_divides_monus: all d:UInt, m:UInt, n:UInt.
  if divides(d, m) and divides(d, n) then divides(d, m ∸ n)
proof
  arbitrary d:UInt, m:UInt, n:UInt
  assume prem
  obtain k where dk_m: d * k = m from expand divides in conjunct 0 of prem
  obtain l where dl_n: d * l = n from expand divides in conjunct 1 of prem
  expand divides
  choose k ∸ l
  replace uint_dist_mult_monus[d, k, l]
  replace dk_m | dl_n.
end

// If a divisor divides both `m` and `n`, it divides `m % n`.
theorem uint_divides_mod_of_divides: all d:UInt, m:UInt, n:UInt.
  if divides(d, m) and divides(d, n) then divides(d, m % n)
proof
  arbitrary d:UInt, m:UInt, n:UInt
  assume prem
  expand operator%
  have d_m: divides(d, m) by conjunct 0 of prem
  have d_n: divides(d, n) by conjunct 1 of prem
  have d_qn: divides(d, (m / n) * n)
    by apply uint_divides_mult_left[d, n, m / n] to d_n
  have both: divides(d, m) and divides(d, (m / n) * n) by d_m, d_qn
  apply uint_divides_monus[d, m, (m / n) * n] to both
end

// Euclidean algorithm for the greatest common divisor.
recfun gcd(a : UInt, b : UInt) -> UInt
  measure b of UInt
{
  if b = 0 then a
  else gcd(b, a % b)
}
terminates {
  arbitrary a:UInt, b:UInt
  assume bnz: not (b = 0)
  have b_pos: 0 < b by apply or_not to uint_zero_or_positive[b], bnz
  conclude a % b < b by apply uint_mod_less_divisor[a,b] to b_pos
}

theorem uint_divides_mod: all d:UInt, m:UInt, n:UInt.
  if divides(d, n) and divides(d, m % n) and 0 < n then divides(d, m)
proof
  arbitrary d:UInt, m:UInt, n:UInt
  assume prem
  obtain k1 where dk1_n: d*k1 = n from expand divides in conjunct 0 of prem
  obtain k2 where dk2_mn: d*k2 = m % n from expand divides in conjunct 1 of prem
  have n_pos: 0 < n by prem
  have eq1: (m / n) * n + m % n = m by apply uint_div_mod[m, n] to n_pos
  expand divides
  have eq2: (m / n) * n + d*k2 = m by replace symmetric dk2_mn in eq1
  define X = m / n
  have eq3: (m / n)*d*k1 + d*k2 = m by expand X in replace symmetric dk1_n in eq2
  have eq4: d*(m/n)*k1 + d*k2 = m by replace uint_mult_commute[m/n, d] in eq3
  have eq5: d*((m/n)*k1 + k2) = m by replace symmetric uint_dist_mult_add[d, (m/n)*k1, k2] in eq4
  choose (m/n)*k1 + k2
  eq5
end

theorem uint_gcd_divides: all b:UInt, a:UInt. divides(gcd(a,b), a) and divides(gcd(a,b), b)
proof
  define P = fun b':UInt {all a:UInt. divides(gcd(a,b'), a) and divides(gcd(a,b'), b')}
  have X: all j:UInt. (if (all i:UInt. (if i < j then P(i))) then P(j)) by {
    arbitrary j:UInt
    assume IH: all i:UInt. (if i < j then P(i))
    expand P
    cases uint_zero_or_positive[j]
      case j_z {
        arbitrary a:UInt
        have A: divides(gcd(a, j), a) by {
          replace j_z
          expand divides | gcd
          choose 1
          conclude a * 1 = a by uint_mult_one
        }
        have B: divides(gcd(a, j), j) by {
          replace j_z
          expand divides | gcd
          choose 0
          conclude a * 0 = 0 by uint_mult_zero
        }
        A, B
      }
      case j_pos {
        arbitrary a:UInt
        have smaller: a % j < j
          by apply uint_mod_less_divisor[a,j] to j_pos
        have div_j_div_aj: divides(gcd(j, a % j), j) and divides(gcd(j, a % j), a % j)
          by (expand P in apply IH[a%j] to smaller)[j]
        have j_ne_z: not (j = 0) by apply uint_pos_not_zero to j_pos
        have A: divides(gcd(a, j), a) by {
          expand gcd
          replace apply eq_false to j_ne_z
          conclude divides(gcd(j, a % j), a)
            by apply uint_divides_mod[gcd(j, a % j), a, j] to div_j_div_aj, j_pos
        }
        have B: divides(gcd(a, j), j) by {
          expand gcd
          replace apply eq_false to j_ne_z
          conclude divides(gcd(j, a % j), j) by div_j_div_aj
        }
        A, B
      }
  }
  arbitrary b:UInt
  expand P in apply uint_strong_induction[P,b] to X
end

theorem uint_gcd_divides_left: all a:UInt, b:UInt. divides(gcd(a,b), a)
proof
  arbitrary a:UInt, b:UInt
  conjunct 0 of uint_gcd_divides[b, a]
end

theorem uint_gcd_divides_right: all a:UInt, b:UInt. divides(gcd(a,b), b)
proof
  arbitrary a:UInt, b:UInt
  conjunct 1 of uint_gcd_divides[b, a]
end

theorem uint_divides_less_equal: all a:UInt, b:UInt.
  if divides(a, b) and 0 < b then a ≤ b
proof
  arbitrary a:UInt, b:UInt
  assume prem
  obtain k where ak_b: a * k = b from expand divides in conjunct 0 of prem
  have k_ne_z: not (k = 0) by {
    assume k_z
    have a0_b: a * 0 = b by replace k_z in ak_b
    have z_b: 0 = b by a0_b
    have b_z: b = 0 by symmetric z_b
    have b_nz: not (b = 0) by apply uint_pos_not_zero to conjunct 1 of prem
    apply b_nz to b_z
  }
  have k_pos: 0 < k by apply uint_not_zero_pos to k_ne_z
  have one_le_k: 1 ≤ k by apply uint_pos_implies_one_le to k_pos
  have a_le_ak: a * 1 ≤ a * k by apply uint_mult_mono_le[a, 1, k] to one_le_k
  have a_le_ak2: a ≤ a * k by a_le_ak
  replace ak_b in a_le_ak2
end

theorem uint_divides_antisymmetric: all a:UInt, b:UInt.
  if divides(a, b) and divides(b, a) then a = b
proof
  arbitrary a:UInt, b:UInt
  assume prem
  cases uint_zero_or_positive[a]
  case a_z {
    obtain k where ak_b: a * k = b from expand divides in conjunct 0 of prem
    have z_b: 0 = b by replace a_z in ak_b
    have b_z: b = 0 by symmetric z_b
    replace a_z | b_z.
  }
  case a_pos {
    have b_ne_z: not (b = 0) by {
      assume b_z
      obtain k where bk_a: b * k = a from expand divides in conjunct 1 of prem
      have z_k_a: 0 * k = a by replace b_z in bk_a
      have z_a: 0 = a by z_k_a
      have a_z: a = 0 by symmetric z_a
      have a_ne_z: not (a = 0) by apply uint_pos_not_zero to a_pos
      apply a_ne_z to a_z
    }
    have b_pos: 0 < b by apply uint_not_zero_pos to b_ne_z
    have div_ab: divides(a, b) by conjunct 0 of prem
    have div_ba: divides(b, a) by conjunct 1 of prem
    have ab: divides(a, b) and 0 < b by div_ab, b_pos
    have ba: divides(b, a) and 0 < a by div_ba, a_pos
    have a_le_b: a ≤ b by apply uint_divides_less_equal[a, b] to ab
    have b_le_a: b ≤ a by apply uint_divides_less_equal[b, a] to ba
    apply uint_less_equal_antisymmetric to a_le_b, b_le_a
  }
end

theorem uint_gcd_greatest: all d:UInt, a:UInt, b:UInt.
  if divides(d, a) and divides(d, b) then divides(d, gcd(a,b))
proof
  define P = fun b':UInt {all a:UInt, d:UInt.
    if divides(d, a) and divides(d, b') then divides(d, gcd(a,b'))}
  have X: all j:UInt. (if (all i:UInt. (if i < j then P(i))) then P(j)) by {
    arbitrary j:UInt
    assume IH: all i:UInt. (if i < j then P(i))
    expand P
    cases uint_zero_or_positive[j]
      case j_z {
        arbitrary a:UInt, d:UInt
        assume prem
        replace j_z
        expand gcd
        conjunct 0 of prem
      }
      case j_pos {
        arbitrary a:UInt, d:UInt
        assume prem
        have smaller: a % j < j
          by apply uint_mod_less_divisor[a,j] to j_pos
        have d_aj: divides(d, a % j)
          by apply uint_divides_mod_of_divides[d, a, j] to prem
        have IH_aj: all a0:UInt, d0:UInt.
          if divides(d0, a0) and divides(d0, a % j) then divides(d0, gcd(a0, a % j))
          by expand P in apply IH[a%j] to smaller
        have j_ne_z: not (j = 0) by apply uint_pos_not_zero to j_pos
        expand gcd
        replace apply eq_false to j_ne_z
        have d_j: divides(d, j) by conjunct 1 of prem
        have both: divides(d, j) and divides(d, a % j) by d_j, d_aj
        apply IH_aj[j, d] to both
      }
  }
  arbitrary d:UInt, a:UInt, b:UInt
  assume prem
  have Pb: P(b) by apply uint_strong_induction[P,b] to X
  apply (expand P in Pb)[a, d] to prem
end

// Divisibility is preserved by multiplying both sides by a common factor.
theorem uint_divides_mult_both: all k:UInt, a:UInt, b:UInt.
  if divides(a, b) then divides(k * a, k * b)
proof
  arbitrary k:UInt, a:UInt, b:UInt
  assume ab
  obtain x where ax_b: a * x = b from expand divides in ab
  expand divides
  choose x
  equations
    (k * a) * x = k * (a * x)  by uint_mult_assoc[k, a, x]
            ... = k * b        by replace ax_b.
end

// A positive common factor may be cancelled from a divisibility fact.
theorem uint_divides_mult_cancel: all k:UInt, a:UInt, b:UInt.
  if 0 < k and divides(k * a, k * b) then divides(a, b)
proof
  arbitrary k:UInt, a:UInt, b:UInt
  assume prem
  have k_pos: 0 < k by prem
  obtain x where kax_kb: (k * a) * x = k * b from expand divides in conjunct 1 of prem
  have e1: k * (a * x) = k * b by replace symmetric uint_mult_assoc[k, a, x] in kax_kb
  have e2: a * x = b by apply uint_pos_mult_left_cancel[k, a * x, b] to k_pos, e1
  expand divides
  choose x
  e2
end

// `gcd` distributes over a common multiplicative factor. Proved from the
// universal property (`uint_gcd_greatest`/`uint_gcd_divides_*`) via
// antisymmetry, so it holds for every `k` including zero. This is the key
// ingredient behind the least-common-multiple universal property below.
theorem uint_gcd_mult_distributive: all k:UInt, a:UInt, b:UInt.
  gcd(k * a, k * b) = k * gcd(a, b)
proof
  arbitrary k:UInt, a:UInt, b:UInt
  cases uint_zero_or_positive[k]
  case k_z {
    suffices gcd((0:UInt) * a, 0 * b) = 0 * gcd(a, b)  by replace k_z.
    evaluate
  }
  case k_pos {
    have d1: divides(k * gcd(a, b), gcd(k * a, k * b)) by {
      have kga: divides(k * gcd(a, b), k * a)
        by apply uint_divides_mult_both[k, gcd(a, b), a] to uint_gcd_divides_left[a, b]
      have kgb: divides(k * gcd(a, b), k * b)
        by apply uint_divides_mult_both[k, gcd(a, b), b] to uint_gcd_divides_right[a, b]
      apply uint_gcd_greatest[k * gcd(a, b), k * a, k * b] to kga, kgb
    }
    have d2: divides(gcd(k * a, k * b), k * gcd(a, b)) by {
      have kka: divides(k, k * a) by { expand divides choose a. }
      have kkb: divides(k, k * b) by { expand divides choose b. }
      have k_div_G: divides(k, gcd(k * a, k * b))
        by apply uint_gcd_greatest[k, k * a, k * b] to kka, kkb
      obtain e where ke_G: k * e = gcd(k * a, k * b) from expand divides in k_div_G
      have ke_ka: divides(k * e, k * a)
        by replace symmetric ke_G in uint_gcd_divides_left[k * a, k * b]
      have ke_kb: divides(k * e, k * b)
        by replace symmetric ke_G in uint_gcd_divides_right[k * a, k * b]
      have e_a: divides(e, a) by apply uint_divides_mult_cancel[k, e, a] to k_pos, ke_ka
      have e_b: divides(e, b) by apply uint_divides_mult_cancel[k, e, b] to k_pos, ke_kb
      have e_g: divides(e, gcd(a, b)) by apply uint_gcd_greatest[e, a, b] to e_a, e_b
      have ke_kg: divides(k * e, k * gcd(a, b))
        by apply uint_divides_mult_both[k, e, gcd(a, b)] to e_g
      conclude divides(gcd(k * a, k * b), k * gcd(a, b)) by replace ke_G in ke_kg
    }
    apply uint_divides_antisymmetric[gcd(k * a, k * b), k * gcd(a, b)] to d2, d1
  }
end

theorem uint_div_cancel: all y:UInt. if 0 < y then y / y = 1
proof
  arbitrary y:UInt
  assume y_pos
  have y_ne_zero: not (y = 0) by apply uint_pos_not_zero to y_pos
  expand operator/
  replace (apply eq_false to y_ne_zero)
  show 1 + 0 / y = 1
  replace (apply uint_zero_div to y_pos).
end

theorem uint_zero_mod: all x:UInt. 0 % x = 0
proof
  arbitrary x:UInt
  expand operator%
  expand lit | fromNat
  uint_bzero_monus
end

theorem uint_mod_self_zero: all y:UInt. y % y = 0
proof
  arbitrary y:UInt
  have y_z_p: y = 0 or 0 < y by uint_zero_or_positive[y]
  cases y_z_p
  case y_z {
    replace y_z
    uint_zero_mod
  }
  case y_pos {
    expand operator%
    have yyc: y/y = 1 by apply uint_div_cancel to y_pos
    replace yyc.
  }
end

theorem uint_mod_one: all n:UInt. n % 1 = 0
proof
  arbitrary n:UInt
  have one_pos: 0 < 1 by .
  have nm_lt_1: n % 1 < 1 by apply uint_mod_less_divisor[n, 1] to one_pos
  have nm_le_0: n % 1 ≤ 0
    by apply uint_less_add_one_implies_less_equal[n % 1, 0] to nm_lt_1
  apply uint_less_equal_zero to nm_le_0
end

theorem uint_div_one: all n:UInt. n / 1 = n
proof
  arbitrary n:UInt
  have one_pos: 0 < 1 by .
  have eq1: (n / 1) * 1 + (n % 1) = n by apply uint_div_mod[n, 1] to one_pos
  have eq2: (n / 1) + (n % 1) = n by eq1
  replace uint_mod_one in eq2
end

theorem uint_add_div_one: all n:UInt, m:UInt.
  if 0 < m
  then (n + m) / m = 1 + n / m
proof
  arbitrary n:UInt, m:UInt
  assume m_pos
  have m_nz: not (m = 0) by apply uint_pos_not_zero to m_pos
  have m_le_nm: m ≤ n + m by {
    have h: m ≤ m + n by uint_less_equal_add[m, n]
    replace uint_add_commute[m, n] in h
  }
  have not_nm_m: not (n + m < m) by {
    assume nm_m: n + m < m
    have nm_le_m: n + m ≤ m by apply uint_less_implies_less_equal to nm_m
    have eq_nm: n + m = m by apply uint_less_equal_antisymmetric to nm_le_m, m_le_nm
    have m_lt_m: m < m by replace eq_nm in nm_m
    apply uint_less_irreflexive to m_lt_m
  }
  equations
          (n + m) / m
        = 1 + ((n + m) ∸ m) / m  by {
            expand operator/
            replace (apply eq_false to not_nm_m)
                  | (apply eq_false to m_nz).
          }
    ... = 1 + n / m              by {
            replace uint_add_commute[n, m].
          }
end

theorem uint_mult_div_inverse: all n:UInt, m:UInt.
  (if 0 < m then (n * m) / m = n)
proof
  define P = fun n:UInt { all m:UInt. if 0 < m then (n * m) / m = n }
  have base_case: P(0) by {
    expand P
    arbitrary m:UInt
    assume m_pos
    show (0 * m) / m = 0
    apply uint_zero_div to m_pos
  }
  have ind: all n:UInt. if P(n) then P(1 + n) by {
    arbitrary n:UInt
    expand P
    assume IH: all m:UInt. if 0 < m then (n * m) / m = n
    arbitrary m:UInt
    assume m_pos
    show ((1 + n) * m) / m = 1 + n
    have eq1: (1 + n) * m = n * m + m by {
      replace uint_dist_mult_add_right[1, n, m]
      replace uint_add_commute[m, n * m].
    }
    replace eq1
    show (n * m + m) / m = 1 + n
    equations
      (n * m + m) / m = 1 + (n * m) / m  by apply uint_add_div_one[n * m, m] to m_pos
                  ... = 1 + n             by replace apply IH[m] to m_pos.
  }
  expand P in apply uint_induction[P] to base_case, ind
end

theorem uint_mult_div_left_inverse: all n:UInt, m:UInt.
  (if 0 < m then (m * n) / m = n)
proof
  arbitrary n:UInt, m:UInt
  assume m_pos: 0 < m
  replace uint_mult_commute[m, n]
  apply uint_mult_div_inverse[n, m] to m_pos
end

theorem uint_mult_mod_right_zero: all n:UInt, m:UInt.
  if 0 < m then (n * m) % m = 0
proof
  arbitrary n:UInt, m:UInt
  assume m_pos: 0 < m
  expand operator%
  replace (apply uint_mult_div_inverse[n, m] to m_pos)
  uint_monus_cancel[n * m]
end

theorem uint_mult_mod_left_zero: all n:UInt, m:UInt.
  if 0 < m then (m * n) % m = 0
proof
  arbitrary n:UInt, m:UInt
  assume m_pos: 0 < m
  replace uint_mult_commute[m, n]
  apply uint_mult_mod_right_zero[n, m] to m_pos
end

theorem uint_mult_add_div: all k:UInt, n:UInt, m:UInt.
  if 0 < m then (k * m + n) / m = k + n / m
proof
  define P = fun k:UInt { all n:UInt, m:UInt.
    if 0 < m then (k * m + n) / m = k + n / m }
  have base_case: P(0) by {
    expand P.
  }
  have ind: all k:UInt. if P(k) then P(1 + k) by {
    arbitrary k:UInt
    assume IH: P(k)
    expand P
    arbitrary n:UInt, m:UInt
    assume m_pos: 0 < m
    have step_arg: (1 + k) * m + n = (k * m + n) + m by {
      suffices toNat((1 + k) * m + n) = toNat((k * m + n) + m)
        by uint_toNat_injective
      replace toNat_add
      replace toNat_mult[1 + k, m] | toNat_mult[k, m]
      replace toNat_add[1, k] | toNat_add[n, m]
      have one_toNat: toNat(1) = ℕ1 by evaluate
      replace one_toNat
      replace dist_mult_add_right[ℕ1, toNat(k), toNat(m)]
      replace add_commute[toNat(m), toNat(k) * toNat(m)]
      replace add_commute[toNat(m), toNat(n)].
    }
    equations
      ((1 + k) * m + n) / m
          = ((k * m + n) + m) / m  by replace step_arg.
      ... = 1 + (k * m + n) / m    by apply uint_add_div_one[k * m + n, m] to m_pos
      ... = 1 + (k + n / m)        by {
              have IH_all: all n0:UInt, m0:UInt.
                if 0 < m0 then (k * m0 + n0) / m0 = k + n0 / m0
                by expand P in IH
              replace apply IH_all[n, m] to m_pos.
            }
      ... = (1 + k) + n / m        by .
  }
  expand P in apply uint_induction[P] to base_case, ind
end

theorem uint_add_mult_div: all n:UInt, k:UInt, m:UInt.
  if 0 < m then (n + k * m) / m = k + n / m
proof
  arbitrary n:UInt, k:UInt, m:UInt
  assume m_pos: 0 < m
  replace uint_add_commute[n, k * m]
  apply uint_mult_add_div[k, n, m] to m_pos
end

theorem uint_mult_add_mod: all k:UInt, n:UInt, m:UInt.
  if 0 < m then (k * m + n) % m = n % m
proof
  arbitrary k:UInt, n:UInt, m:UInt
  assume m_pos: 0 < m
  expand operator%
  replace (apply uint_mult_add_div[k, n, m] to m_pos)
  replace uint_dist_mult_add_right[k, n / m, m]
  uint_add_both_monus[k * m, n, (n / m) * m]
end

theorem uint_add_mult_mod: all n:UInt, k:UInt, m:UInt.
  if 0 < m then (n + k * m) % m = n % m
proof
  arbitrary n:UInt, k:UInt, m:UInt
  assume m_pos: 0 < m
  replace uint_add_commute[n, k * m]
  apply uint_mult_add_mod[k, n, m] to m_pos
end

theorem uint_mod_mod: all n:UInt, m:UInt.
  if 0 < m then (n % m) % m = n % m
proof
  arbitrary n:UInt, m:UInt
  assume m_pos: 0 < m
  have nm_lt_m: n % m < m by apply uint_mod_less_divisor[n, m] to m_pos
  apply uint_mod_small[n % m, m] to nm_lt_m
end

theorem uint_add_mod: all a:UInt, b:UInt, m:UInt.
  if 0 < m then (a + b) % m = ((a % m) + (b % m)) % m
proof
  arbitrary a:UInt, b:UInt, m:UInt
  assume m_pos: 0 < m
  have a_div: (a / m) * m + (a % m) = a
    by apply uint_div_mod[a, m] to m_pos
  have b_div: (b / m) * m + (b % m) = b
    by apply uint_div_mod[b, m] to m_pos
  have decomp:
    a + b = ((a / m) + (b / m)) * m + ((a % m) + (b % m)) by {
    suffices toNat(a + b) =
             toNat(((a / m) + (b / m)) * m + ((a % m) + (b % m)))
      by uint_toNat_injective
    have a_nat: toNat((a / m) * m + (a % m)) = toNat(a)
      by replace a_div.
    have b_nat: toNat((b / m) * m + (b % m)) = toNat(b)
      by replace b_div.
    have a_nat2: toNat(a / m) * toNat(m) + toNat(a % m) = toNat(a)
      by replace toNat_add | toNat_mult in a_nat
    have b_nat2: toNat(b / m) * toNat(m) + toNat(b % m) = toNat(b)
      by replace toNat_add | toNat_mult in b_nat
    replace toNat_add | toNat_mult
    replace symmetric a_nat2 | symmetric b_nat2
    replace toNat_add[a / m, b / m] | toNat_add[a % m, b % m]
    replace add_commute[toNat(a % m), toNat(b / m) * toNat(m)]
    replace symmetric dist_mult_add_right[toNat(a / m), toNat(b / m), toNat(m)].
  }
  replace decomp
  apply uint_mult_add_mod[(a / m) + (b / m), (a % m) + (b % m), m] to m_pos
end

theorem uint_mult_mod: all a:UInt, b:UInt, m:UInt.
  if 0 < m then (a * b) % m = ((a % m) * (b % m)) % m
proof
  arbitrary a:UInt, b:UInt, m:UInt
  assume m_pos: 0 < m
  have a_div: (a / m) * m + (a % m) = a
    by apply uint_div_mod[a, m] to m_pos
  have b_div: (b / m) * m + (b % m) = b
    by apply uint_div_mod[b, m] to m_pos
  have a_prod: a * b = ((a / m) * b) * m + (a % m) * b by {
    suffices toNat(a * b) = toNat(((a / m) * b) * m + (a % m) * b)
      by uint_toNat_injective
    have a_nat: toNat((a / m) * m + (a % m)) = toNat(a)
      by replace a_div.
    have a_nat2: toNat(a / m) * toNat(m) + toNat(a % m) = toNat(a)
      by replace toNat_add | toNat_mult in a_nat
    replace toNat_add | toNat_mult
    replace symmetric a_nat2
    replace toNat_mult[b, m]
    replace dist_mult_add_right[toNat(a / m) * toNat(m), toNat(a % m), toNat(b)]
    replace mult_commute[toNat(m), toNat(b)].
  }
  have a_zero: (((a / m) * b) * m) % m = 0
    by apply uint_mult_mod_right_zero[(a / m) * b, m] to m_pos
  have reduce_a: (a * b) % m = ((a % m) * b) % m by {
    equations
      (a * b) % m
          = ((((a / m) * b) * m + (a % m) * b) % m) by replace a_prod.
      ... = (((((a / m) * b) * m) % m) + (((a % m) * b) % m)) % m
            by apply uint_add_mod[((a / m) * b) * m, (a % m) * b, m] to m_pos
      ... = (0 + (((a % m) * b) % m)) % m by replace a_zero.
      ... = (((a % m) * b) % m) % m by .
      ... = ((a % m) * b) % m by apply uint_mod_mod[(a % m) * b, m] to m_pos
  }
  have b_prod: (a % m) * b = ((a % m) * (b / m)) * m + (a % m) * (b % m) by {
    suffices toNat((a % m) * b) = toNat(((a % m) * (b / m)) * m + (a % m) * (b % m))
      by uint_toNat_injective
    have b_nat: toNat((b / m) * m + (b % m)) = toNat(b)
      by replace b_div.
    have b_nat2: toNat(b / m) * toNat(m) + toNat(b % m) = toNat(b)
      by replace toNat_add | toNat_mult in b_nat
    replace toNat_add | toNat_mult
    replace symmetric b_nat2
    replace toNat_mult[b / m, m]
    replace dist_mult_add[toNat(a % m), toNat(b / m) * toNat(m), toNat(b % m)].
  }
  have b_zero: (((a % m) * (b / m)) * m) % m = 0
    by apply uint_mult_mod_right_zero[(a % m) * (b / m), m] to m_pos
  have reduce_b: ((a % m) * b) % m = ((a % m) * (b % m)) % m by {
    equations
      ((a % m) * b) % m
          = ((((a % m) * (b / m)) * m + (a % m) * (b % m)) % m)
            by replace b_prod.
      ... = (((((a % m) * (b / m)) * m) % m) + (((a % m) * (b % m)) % m)) % m
            by apply uint_add_mod[((a % m) * (b / m)) * m, (a % m) * (b % m), m] to m_pos
      ... = (0 + (((a % m) * (b % m)) % m)) % m by replace b_zero.
      ... = (((a % m) * (b % m)) % m) % m by .
      ... = ((a % m) * (b % m)) % m by apply uint_mod_mod[(a % m) * (b % m), m] to m_pos
  }
  equations
    (a * b) % m = ((a % m) * b) % m by reduce_a
    ... = ((a % m) * (b % m)) % m by reduce_b
end

theorem uint_div_less_equal: all n:UInt, m:UInt.
  if 0 < m then n / m ≤ n
proof
  arbitrary n:UInt, m:UInt
  assume m_pos: 0 < m
  have one_le_m: 1 ≤ m by apply uint_pos_implies_one_le to m_pos
  have q_le_qm: (n / m) * 1 ≤ (n / m) * m
    by apply uint_mult_mono_le[n / m, 1, m] to one_le_m
  have qm_le_n: (n / m) * m ≤ n by {
    have qmr: (n / m) * m ≤ (n / m) * m + (n % m) by uint_less_equal_add
    replace apply uint_div_mod[n, m] to m_pos in qmr
  }
  apply uint_less_equal_trans to q_le_qm, qm_le_n
end

theorem uint_div_less: all n:UInt, m:UInt.
  if 0 < n and 1 < m then n / m < n
proof
  arbitrary n:UInt, m:UInt
  assume prem
  have n_pos: 0 < n by prem
  have one_lt_m: 1 < m by prem
  have zero_lt_one: 0 < 1 by .
  have m_pos: 0 < m by apply uint_less_trans to zero_lt_one, one_lt_m
  have q_le_n: n / m ≤ n by apply uint_div_less_equal[n, m] to m_pos
  cases expand operator≤ in q_le_n
  case q_lt_n: n / m < n {
    q_lt_n
  }
  case q_eq_n: n / m = n {
    have n_lt_nm: n < n * m by {
      have q_pos: 0 < n / m by replace symmetric q_eq_n in n_pos
      have q_lt_qm: (n / m) * 1 < (n / m) * m
        by apply uint_pos_mult_both_sides_of_less[n / m, 1, m] to q_pos, one_lt_m
      replace q_eq_n in q_lt_qm
    }
    have nm_le_n: n * m ≤ n by {
      have qm_le_n: (n / m) * m ≤ n by {
        have qmr: (n / m) * m ≤ (n / m) * m + (n % m) by uint_less_equal_add
        replace apply uint_div_mod[n, m] to m_pos in qmr
      }
      replace q_eq_n in qm_le_n
    }
    have nat_n_lt_nm: toNat(n) < toNat(n * m) by apply toNat_less to n_lt_nm
    have nat_nm_le_n: toNat(n * m) ≤ toNat(n) by apply toNat_less_equal to nm_le_n
    have nat_n_lt_n: toNat(n) < toNat(n)
      by apply less_le_trans[toNat(n), toNat(n * m), toNat(n)] to nat_n_lt_nm, nat_nm_le_n
    conclude false by apply less_irreflexive to nat_n_lt_n
  }
end

theorem uint_div_zero: all n:UInt. n / 0 = 0
proof
  arbitrary n:UInt
  expand operator/
  replace (apply eq_false to uint_not_less_zero[n]).
end

theorem fromNat_div: all x:Nat, y:Nat. fromNat(x) / fromNat(y) = fromNat(x / y)
proof
  have lem: all y:Nat. all x:Nat. fromNat(x) / fromNat(y) = fromNat(x / y) by {
    arbitrary y:Nat
    cases nat_zero_or_positive[y]
    case y_z {
      arbitrary x:Nat
      replace y_z
      have xz: x / ℕ0 = ℕ0 by expand operator/ | lit.
      replace xz
      uint_div_zero[fromNat(x)]
    }
    case y_pos {
      have y_pos_nat: zero < y by expand lit in y_pos
      have y_ne_zero: not (y = zero) by {
        assume y_z
        replace y_z in y_pos_nat
      }
      have fy_pos: 0 < fromNat(y) by replace from_zero in apply less_fromNat to y_pos
      have fy_ne_zero: not (fromNat(y) = 0) by apply uint_pos_not_zero to fy_pos
      define P = fun x:Nat { fromNat(x) / fromNat(y) = fromNat(x / y) }
      have SI: all i:Nat. if (all j:Nat. if j < i then P(j)) then P(i) by {
        arbitrary i:Nat
        assume IH: all j:Nat. if j < i then P(j)
        expand P
        switch i < y {
          case true assume i_y_t {
            have i_less_y: i < y by simplify with i_y_t.
            have fi_less_fy: fromNat(i) < fromNat(y) by apply less_fromNat to i_less_y
            have lhs_zero: fromNat(i) / fromNat(y) = 0 by {
              expand operator/
              replace (apply eq_true to fi_less_fy).
            }
            have rhs_zero: i / y = zero by {
              expand operator/
              replace (apply eq_true to i_less_y).
            }
            replace lhs_zero | rhs_zero
            expand lit.
          }
          case false assume i_y_f {
            have not_i_y: not (i < y) by i_y_f
            have y_le_i: y ≤ i by apply not_less_implies_less_equal to not_i_y
            have not_fi_fy: not (fromNat(i) < fromNat(y)) by {
              assume fi_fy
              have A: toNat(fromNat(i)) < toNat(fromNat(y)) by apply toNat_less to fi_fy
              have i_less_y: i < y by replace uint_toNat_fromNat in A
              apply not_i_y to i_less_y
            }
            have monus_eq: fromNat(i) ∸ fromNat(y) = fromNat(i ∸ y) by {
              suffices toNat(fromNat(i) ∸ fromNat(y)) = toNat(fromNat(i ∸ y))
                by uint_toNat_injective
              replace toNat_monus | uint_toNat_fromNat.
            }
            have im_less_i: i ∸ y < i by {
              suffices y + (i ∸ y) < y + i by add_both_sides_of_less[y, i ∸ y, i]
              suffices i < y + i by replace apply monus_add_identity[i, y] to y_le_i.
              replace add_commute[y, i]
              apply nat_less_add_pos[i, y] to y_pos
            }
            have ih_app: P(i ∸ y) by apply IH to im_less_i
            have ih_eq: fromNat(i ∸ y) / fromNat(y) = fromNat((i ∸ y) / y)
              by expand P in ih_app
            have rhs_step: i / y = suc(zero) + (i ∸ y) / y by {
              equations
                i / y
                    = suc(zero) + (i ∸ y) / y by {
                        expand operator/
                        replace (apply eq_false to not_i_y)
                              | (apply eq_false to y_ne_zero).
                      }
            }
            equations
              fromNat(i) / fromNat(y)
                  = 1 + (fromNat(i) ∸ fromNat(y)) / fromNat(y)   by {
                      expand operator/
                      replace (apply eq_false to not_fi_fy)
                            | (apply eq_false to fy_ne_zero).
                    }
              ... = 1 + fromNat(i ∸ y) / fromNat(y)              by replace monus_eq.
              ... = 1 + fromNat((i ∸ y) / y)                     by replace ih_eq.
              ... = fromNat(lit(suc(zero)) + (i ∸ y) / y)        by symmetric fromNat_add[lit(suc(zero)), (i ∸ y) / y]
              ... = fromNat(suc(zero) + (i ∸ y) / y)             by expand lit.
              ... = # fromNat(i / y) #                           by replace rhs_step.
          }
        }
      }
      arbitrary x:Nat
      have hx: P(x) by apply strong_induction[P, x] to SI
      expand P in hx
    }
  }
  arbitrary x:Nat, y:Nat
  lem[y, x]
end

theorem fromNat_mod: all x:Nat, y:Nat. fromNat(x) % fromNat(y) = fromNat(x % y)
proof
  arbitrary x:Nat, y:Nat
  suffices toNat(fromNat(x) % fromNat(y)) = toNat(fromNat(x % y))
    by uint_toNat_injective
  expand operator%
  replace toNat_monus | toNat_mult | fromNat_div | uint_toNat_fromNat.
end

theorem toNat_mod: all x:UInt, y:UInt. toNat(x % y) = toNat(x) % toNat(y)
proof
  arbitrary x:UInt, y:UInt
  replace symmetric uint_fromNat_toNat[x]
  replace symmetric uint_fromNat_toNat[y]
  replace fromNat_mod[toNat(x), toNat(y)]
  replace uint_toNat_fromNat[toNat(x)]
  replace uint_toNat_fromNat[toNat(y)]
  uint_toNat_fromNat[toNat(x) % toNat(y)]
end

theorem uint_lit_div: all x:Nat, y:Nat. (fromNat(lit(x)) / fromNat(lit(y))) = fromNat(lit(x) / lit(y))
proof
  arbitrary x:Nat, y:Nat
  fromNat_div[lit(x), lit(y)]
end

auto uint_lit_div

// Least common multiple: zero if either operand is zero, otherwise
// `(a * b) / gcd(a, b)`. The product-form definition keeps `lcm`
// expressible without a second recursive function and makes the
// `gcd * lcm = a * b` identity an algebraic consequence of the
// quotient/remainder spec.
fun lcm(a : UInt, b : UInt) {
  if a = 0 or b = 0 then 0
  else (a * b) / gcd(a, b)
}

// `lcm(0, b) = 0`: definitional, since the `or` short-circuits.
theorem uint_lcm_zero_left: all b:UInt. lcm(0, b) = 0
proof
  arbitrary b:UInt
  expand lcm.
end

auto uint_lcm_zero_left

// `lcm(a, 0) = 0`: the `or` clause `0 = 0` collapses the conditional.
theorem uint_lcm_zero_right: all a:UInt. lcm(a, 0) = 0
proof
  arbitrary a:UInt
  expand lcm.
end

auto uint_lcm_zero_right

// A positive first operand forces a positive gcd, since `gcd(a, b)`
// divides `a` and only zero divides into zero.
theorem uint_gcd_pos: all a:UInt, b:UInt. if 0 < a then 0 < gcd(a, b)
proof
  arbitrary a:UInt, b:UInt
  assume a_pos
  have g_div_a: divides(gcd(a, b), a) by uint_gcd_divides_left[a, b]
  have g_ne_z: not (gcd(a, b) = 0) by {
    assume g_z
    have div0a: divides(0, a) by replace g_z in g_div_a
    obtain j where z_a: 0 = a from expand divides in div0a
    have a_ne: not (a = 0) by apply uint_pos_not_zero to a_pos
    apply a_ne to (symmetric z_a)
  }
  apply uint_not_zero_pos to g_ne_z
end

// `a` divides its least common multiple with `b`. In the interesting
// case `lcm(a, b) = (a * b) / gcd(a, b)`; since `gcd(a, b)` divides `b`
// as `gcd(a, b) * k = b`, the quotient collapses to `a * k`.
theorem uint_divides_lcm_left: all a:UInt, b:UInt. divides(a, lcm(a, b))
proof
  arbitrary a:UInt, b:UInt
  cases uint_zero_or_positive[a]
  case a_z {
    replace a_z
    uint_divides_zero[0]
  }
  case a_pos {
    cases uint_zero_or_positive[b]
    case b_z {
      replace b_z
      uint_divides_zero[a]
    }
    case b_pos {
      have g_pos: 0 < gcd(a, b) by apply uint_gcd_pos[a, b] to a_pos
      have g_div_b: divides(gcd(a, b), b) by uint_gcd_divides_right[a, b]
      obtain k where gk_b: gcd(a, b) * k = b from expand divides in g_div_b
      have a_ne: not (a = 0) by apply uint_pos_not_zero to a_pos
      have b_ne: not (b = 0) by apply uint_pos_not_zero to b_pos
      have cond_false: not (a = 0 or b = 0) by {
        assume c
        cases c
        case l { apply a_ne to l }
        case r { apply b_ne to r }
      }
      have lcm_eq: lcm(a, b) = a * k by {
        expand lcm
        replace (apply eq_false to cond_false)
        have prod_eq: a * b = gcd(a, b) * (a * k) by {
          have s1: a * b = a * (gcd(a, b) * k) by replace gk_b.
          replace s1
          replace uint_mult_commute[a, gcd(a, b)].
        }
        replace prod_eq
        apply uint_mult_div_left_inverse[a * k, gcd(a, b)] to g_pos
      }
      expand divides
      choose k
      symmetric lcm_eq
    }
  }
end

// `b` divides its least common multiple with `a`. Mirrors
// `uint_divides_lcm_left`, using that `gcd(a, b)` divides `a`.
theorem uint_divides_lcm_right: all a:UInt, b:UInt. divides(b, lcm(a, b))
proof
  arbitrary a:UInt, b:UInt
  cases uint_zero_or_positive[a]
  case a_z {
    replace a_z
    uint_divides_zero[b]
  }
  case a_pos {
    cases uint_zero_or_positive[b]
    case b_z {
      replace b_z
      uint_divides_zero[0]
    }
    case b_pos {
      have g_pos: 0 < gcd(a, b) by apply uint_gcd_pos[a, b] to a_pos
      have g_div_a: divides(gcd(a, b), a) by uint_gcd_divides_left[a, b]
      obtain j where gj_a: gcd(a, b) * j = a from expand divides in g_div_a
      have a_ne: not (a = 0) by apply uint_pos_not_zero to a_pos
      have b_ne: not (b = 0) by apply uint_pos_not_zero to b_pos
      have cond_false: not (a = 0 or b = 0) by {
        assume c
        cases c
        case l { apply a_ne to l }
        case r { apply b_ne to r }
      }
      have lcm_eq: lcm(a, b) = j * b by {
        expand lcm
        replace (apply eq_false to cond_false)
        have prod_eq: a * b = gcd(a, b) * (j * b) by {
          have s1: a * b = (gcd(a, b) * j) * b by replace gj_a.
          replace s1.
        }
        replace prod_eq
        apply uint_mult_div_left_inverse[j * b, gcd(a, b)] to g_pos
      }
      expand divides
      choose j
      replace lcm_eq
      uint_mult_commute[b, j]
    }
  }
end

// Universal property of the least common multiple: every common multiple of
// `a` and `b` is a multiple of `lcm(a, b)`. Combined with
// `uint_divides_lcm_left`/`uint_divides_lcm_right`, this characterizes
// `lcm(a, b)` as the least common multiple up to divisibility. The proof
// runs entirely through `uint_gcd_mult_distributive`: applied with factor
// `m` it gives `gcd(m*a, m*b) = m*gcd(a,b)`, and applied with factor `a*b`
// (after rewriting `m*a`/`m*b` as `(a*b)*t`/`(a*b)*s`) it gives
// `gcd(m*a, m*b) = (a*b)*gcd(t,s)`; equating and cancelling `gcd(a,b)` shows
// `m = lcm(a,b) * gcd(t,s)`.
theorem uint_lcm_least: all a:UInt, b:UInt, m:UInt.
  if divides(a, m) and divides(b, m) then divides(lcm(a, b), m)
proof
  arbitrary a:UInt, b:UInt, m:UInt
  assume prem
  have a_div_m: divides(a, m) by conjunct 0 of prem
  have b_div_m: divides(b, m) by conjunct 1 of prem
  cases uint_zero_or_positive[a]
  case a_z {
    obtain k where zk_m: 0 * k = m from expand divides in (replace a_z in a_div_m)
    have m_z: m = 0 by symmetric zk_m
    replace m_z
    uint_divides_zero[lcm(a, b)]
  }
  case a_pos {
    cases uint_zero_or_positive[b]
    case b_z {
      obtain k where zk_m: 0 * k = m from expand divides in (replace b_z in b_div_m)
      have m_z: m = 0 by symmetric zk_m
      replace m_z
      uint_divides_zero[lcm(a, b)]
    }
    case b_pos {
      have g_pos: 0 < gcd(a, b) by apply uint_gcd_pos[a, b] to a_pos
      have a_ne: not (a = 0) by apply uint_pos_not_zero to a_pos
      have b_ne: not (b = 0) by apply uint_pos_not_zero to b_pos
      have cond_false: not (a = 0 or b = 0) by {
        assume c
        cases c
        case l { apply a_ne to l }
        case r { apply b_ne to r }
      }
      have lcm_def: lcm(a, b) = (a * b) / gcd(a, b) by {
        expand lcm
        replace (apply eq_false to cond_false).
      }
      obtain s where as_m: a * s = m from expand divides in a_div_m
      obtain t where bt_m: b * t = m from expand divides in b_div_m
      // Rewrite the two products m*a and m*b as multiples of a*b.
      have ma_eq: m * a = (a * b) * t by {
        replace symmetric bt_m
        replace uint_mult_commute[t, a] | uint_mult_commute[b, a * t]
              | uint_mult_commute[t, b].
      }
      have mb_eq: m * b = (a * b) * s by {
        replace symmetric as_m
        replace uint_mult_commute[s, b].
      }
      have E1: gcd(m * a, m * b) = m * gcd(a, b) by uint_gcd_mult_distributive[m, a, b]
      have E2: gcd((a * b) * t, (a * b) * s) = (a * b) * gcd(t, s)
        by uint_gcd_mult_distributive[a * b, t, s]
      have E3: m * gcd(a, b) = (a * b) * gcd(t, s) by {
        equations
          m * gcd(a, b) = gcd(m * a, m * b)              by symmetric E1
                    ... = gcd((a * b) * t, (a * b) * s)  by replace ma_eq | mb_eq.
                    ... = (a * b) * gcd(t, s)            by E2
      }
      // a*b = gcd(a,b) * lcm(a,b), and lcm(a,b) = w below.
      have g_div_ab: divides(gcd(a, b), a * b)
        by apply uint_divides_mult_right[gcd(a, b), a, b] to uint_gcd_divides_left[a, b]
      obtain w where gw_ab: gcd(a, b) * w = a * b from expand divides in g_div_ab
      have lcm_eq_w: lcm(a, b) = w by {
        replace lcm_def
        have ab_eq: a * b = gcd(a, b) * w by symmetric gw_ab
        replace ab_eq
        apply uint_mult_div_left_inverse[w, gcd(a, b)] to g_pos
      }
      have final: m = lcm(a, b) * gcd(t, s) by {
        have step: gcd(a, b) * m = gcd(a, b) * (w * gcd(t, s)) by {
          equations
            gcd(a, b) * m = m * gcd(a, b)                by uint_mult_commute[gcd(a, b), m]
                      ... = (a * b) * gcd(t, s)          by E3
                      ... = (gcd(a, b) * w) * gcd(t, s)  by replace symmetric gw_ab.
                      ... = gcd(a, b) * (w * gcd(t, s))  by uint_mult_assoc[gcd(a, b), w, gcd(t, s)]
        }
        have m_eq: m = w * gcd(t, s)
          by apply uint_pos_mult_left_cancel[gcd(a, b), m, w * gcd(t, s)] to g_pos, step
        replace lcm_eq_w
        m_eq
      }
      expand divides
      choose gcd(t, s)
      symmetric final
    }
  }
end

// `gcd(a, 0) = a`: definitional, since the `b = 0` guard selects `a`.
theorem uint_gcd_zero_right: all a:UInt. gcd(a, 0) = a
proof
  arbitrary a:UInt
  expand gcd.
end

// Commutativity of `gcd`. Proved through the universal property rather
// than the (asymmetric) Euclidean recursion: each of `gcd(a,b)` and
// `gcd(b,a)` divides the other by `uint_gcd_greatest`, so they are
// equal by `uint_divides_antisymmetric`.
theorem uint_gcd_commutative: all a:UInt, b:UInt. gcd(a, b) = gcd(b, a)
proof
  arbitrary a:UInt, b:UInt
  have ab_a: divides(gcd(a, b), a) by uint_gcd_divides_left[a, b]
  have ab_b: divides(gcd(a, b), b) by uint_gcd_divides_right[a, b]
  have ba_a: divides(gcd(b, a), a) by uint_gcd_divides_right[b, a]
  have ba_b: divides(gcd(b, a), b) by uint_gcd_divides_left[b, a]
  have d1: divides(gcd(a, b), gcd(b, a))
    by apply uint_gcd_greatest[gcd(a, b), b, a] to ab_b, ab_a
  have d2: divides(gcd(b, a), gcd(a, b))
    by apply uint_gcd_greatest[gcd(b, a), a, b] to ba_a, ba_b
  apply uint_divides_antisymmetric[gcd(a, b), gcd(b, a)] to d1, d2
end

// `gcd(0, b) = b`, via commutativity and `uint_gcd_zero_right`.
theorem uint_gcd_zero_left: all b:UInt. gcd(0, b) = b
proof
  arbitrary b:UInt
  replace uint_gcd_commutative[0, b]
  uint_gcd_zero_right[b]
end

// `gcd(a, a) = a`, again through the universal property.
theorem uint_gcd_self: all a:UInt. gcd(a, a) = a
proof
  arbitrary a:UInt
  have a_a: divides(a, a) by uint_divides_refl[a]
  have d1: divides(a, gcd(a, a))
    by apply uint_gcd_greatest[a, a, a] to a_a, a_a
  have d2: divides(gcd(a, a), a) by uint_gcd_divides_left[a, a]
  apply uint_divides_antisymmetric[gcd(a, a), a] to d2, d1
end

// Commutativity of `lcm`. An algebraic consequence of the product-form
// definition: `a*b = b*a` and `gcd(a,b) = gcd(b,a)`.
theorem uint_lcm_commutative: all a:UInt, b:UInt. lcm(a, b) = lcm(b, a)
proof
  arbitrary a:UInt, b:UInt
  cases uint_zero_or_positive[a]
  case a_z {
    replace a_z.
  }
  case a_pos {
    cases uint_zero_or_positive[b]
    case b_z {
      replace b_z.
    }
    case b_pos {
      have a_ne: not (a = 0) by apply uint_pos_not_zero to a_pos
      have b_ne: not (b = 0) by apply uint_pos_not_zero to b_pos
      have cond_ab: not (a = 0 or b = 0) by {
        assume c
        cases c
        case l { apply a_ne to l }
        case r { apply b_ne to r }
      }
      have cond_ba: not (b = 0 or a = 0) by {
        assume c
        cases c
        case l { apply b_ne to l }
        case r { apply a_ne to r }
      }
      expand lcm
      replace (apply eq_false to cond_ab) | (apply eq_false to cond_ba)
      replace uint_mult_commute[a, b] | uint_gcd_commutative[a, b].
    }
  }
end

// The product of the gcd and lcm recovers the product of the arguments:
// `gcd(a, b) * lcm(a, b) = a * b`. In the interesting case
// `lcm(a, b) = (a * b) / gcd(a, b)`; since `gcd(a, b)` divides `a` (hence
// `a * b`), the quotient is exact, so multiplying it back by `gcd(a, b)`
// recovers `a * b`.
theorem uint_gcd_lcm_product: all a:UInt, b:UInt. gcd(a, b) * lcm(a, b) = a * b
proof
  arbitrary a:UInt, b:UInt
  cases uint_zero_or_positive[a]
  case a_z {
    replace a_z.
  }
  case a_pos {
    cases uint_zero_or_positive[b]
    case b_z {
      replace b_z.
    }
    case b_pos {
      have g_pos: 0 < gcd(a, b) by apply uint_gcd_pos[a, b] to a_pos
      have g_div_a: divides(gcd(a, b), a) by uint_gcd_divides_left[a, b]
      obtain k where gk_a: gcd(a, b) * k = a from expand divides in g_div_a
      have a_ne: not (a = 0) by apply uint_pos_not_zero to a_pos
      have b_ne: not (b = 0) by apply uint_pos_not_zero to b_pos
      have cond_false: not (a = 0 or b = 0) by {
        assume c
        cases c
        case l { apply a_ne to l }
        case r { apply b_ne to r }
      }
      have prod_eq: a * b = gcd(a, b) * (k * b) by {
        have s1: a * b = (gcd(a, b) * k) * b by replace gk_a.
        replace s1.
      }
      expand lcm
      replace (apply eq_false to cond_false)
      replace prod_eq
      replace (apply uint_mult_div_left_inverse[k * b, gcd(a, b)] to g_pos).
    }
  }
end