Natural.div

-- EUCLIDEAN DIVISION on ℕ, by a successor (so the divisor is never
-- zero and no side condition is needed). The recursion is on FUEL —
-- `a ∸ S b` is smaller than `a` but not structurally so — and `a`
-- itself is always enough fuel.
--
-- This is the arithmetic ℚ needs in order to have a CEILING, and the
-- ceiling is what ℝ needs in order to have multiplication: the product
-- of two regular sequences must be sampled at a depth proportional to
-- a bound on both, and a bound has to be a natural, computed from a
-- rational, by a function — i.e. by something the quotient respects.
-- ===== casing on a natural, with the equation in hand =====

import Natural (+, *, plusZeroId, plusSucId, zeroPlusId, sucPlus, plusComm, plusAssoc, sucInj, multZeroId, multSucId, zeroMult, sucMult, multComm, multAssoc)
import Natural.order (≤, leRefl, leOfEq, leZero, leZeroInv, leSucMono, leTrans, leTotal, leMultMonoR, leStep)
import Natural.more (∸, pred, plusCancelL, multDistribR, monusZeroR, monusSucR, zeroMonus, sucMonusSuc, plusMonus, monusPlus, monusSelf, monusZeroOfLe, leMonusPlus)
import Natural.eq (zNotS)
import Core.equality (trans, sym, cong, transport)
import Core.id (Id, idToEq, eqToId)

natCase : {C : 𝕌} → C → (ℕ → C) → ℕ → C
natCase = λC z s n. ℕ-elim z (m rec. s m) n

-- the eliminator for it: each branch may use the equation that put it
-- there, which is what lets a division step know a > b. Two copies,
-- because a motive into Ω and a motive into 𝕌 are different rules and
-- the division spec needs one of each (the quotient equation is a
-- a prop, the remainder bound is data).
natCaseElimEqD : (C : 𝕌)
  (P : C → 𝕌)
  (z : C)
  (s : ℕ → C)
  (n : ℕ)
  → ((n ≡ Z) → P z) → ((m : ℕ) → (n ≡ S m) → P (s m)) → P (natCase z s n)
  using (Natural.div.natCase.eq)
natCaseElimEqD = λC P z s n. ℕ-elim (λh0 h1. h0 ⋆) (k' rec. λh0 h1. h1 k' ⋆) n

natCaseElimEq : (C : 𝕌)
  (P : C → Ω)
  (z : C)
  (s : ℕ → C)
  (n : ℕ)
  → ((n ≡ Z) → P z) → ((m : ℕ) → (n ≡ S m) → P (s m)) → P (natCase z s n)
  using (Natural.div.natCase.eq)
natCaseElimEq = λC P z s n. ℕ-elim (λh0 h1. h0 ⋆) (k' rec. λh0 h1. h1 k' ⋆) n

-- ===== comparison, read off monus =====
leOfMonusZero : {a b : ℕ} → (a ∸ b ≡ Z) → a ≤ b
  using (Natural.more.monusZeroR, Natural.order.≤.unfold)
leOfMonusZero =
  λa. ℕ-elim
    λb h. leZero _
    a' ih. λb. ℕ-elim
      λh. 𝟘-elim (zNotS (sym (S a') Z h))
      b' rec. λh. leSucMono (ih b' (trans _ _ _ (sym _ _ (sucMonusSuc a' b')) h))
      b
    a

ltOfMonusSuc : {a b : ℕ} (m : ℕ) → (a ∸ b ≡ S m) → S b ≤ a using (Natural.order.≤.unfold)
ltOfMonusSuc =
  λa. ℕ-elim
    λb m h. 𝟘-elim (zNotS (trans _ _ _ (sym _ _ (zeroMonus b)) h))
    a' ih. λb. ℕ-elim
      λm h. leSucMono (leZero a')
      b' rec. λm h. leSucMono (ih b' m (trans _ _ _ (sym _ _ (sucMonusSuc a' b')) h))
      b
    a

-- ===== the fuel decreases =====
-- if a ≤ S f' and b < a then a − (b+1) ≤ f' : peel the successor off
-- both sides of (b+1) + (a − (b+1)) + k ≡ S f'
fuelStep : {a b f' : ℕ} → a ≤ S f' → S b ≤ a → a ∸ S b ≤ f'
  using (Core.id.Id.unfold, Natural.order.≤.unfold)
fuelStep =
  λa b f' le lt. (,)
    b + le .π₁
    eqToId
      _
      _
      sucInj
        a ∸ S b + (b + le .π₁)
        f'
        S (a ∸ S b + (b + le .π₁))
          ≡⟨ sym
            _
            _
            trans
              _
              _
              _
              sucPlus (a ∸ S b + b) (le .π₁)
              cong (λu. ℕ) (λw. S w) (plusAssoc (a ∸ S b) b (le .π₁)) ⟩
            S (a ∸ S b + b) + le .π₁
          ≡⟨ cong
            λu. ℕ
            λw. w + le .π₁
            trans
              _
              _
              _
              cong (λu. ℕ) (λw. S w) (plusComm b (a ∸ S b))
              sym _ _ (sucPlus b (a ∸ S b)) ⟩
            S b + (a ∸ S b) + le .π₁
          ≡⟨ cong (λu. ℕ) (λw. w + le .π₁) (leMonusPlus lt) ⟩ a + le .π₁
          ≡⟨ idToEq _ _ _ (le .π₂) ⟩ S f'

-- ===== division =====
-- dmAux f a b = (q , r) with a ≡ q·(b+1) + r and r ≤ b, whenever a ≤ f
dmAux : ℕ → ℕ → ℕ → ℕ × ℕ
dmAux =
  λf. ℕ-elim
    λa b. Z, a
    f' ih. λa b. natCase (Z, a) (λm. S (ih (a ∸ S b) b .π₁), ih (a ∸ S b) b .π₂) (a ∸ b)
    f

divN : ℕ → ℕ → ℕ
divN = λa b. dmAux a a b .π₁

modN : ℕ → ℕ → ℕ
modN = λa b. dmAux a a b .π₂

-- ===== the specification =====
dmAuxEq : {f a b : ℕ} → a ≤ f → a ≡ dmAux f a b .π₁ * S b + dmAux f a b .π₂
  using (hyp.rw,
    Core.id.Id.unfold,
    Natural.div.dmAux.eq,
    Natural.div.natCase.eq,
    Natural.order.≤.unfold,
    zeroMult,
    zeroMult.rw,
    zeroPlusId,
    zeroPlusId.rw)
dmAuxEq =
  λf. ℕ-elim
    λa b le. ⋆
    f' ih. λa b le. natCaseElimEq
      ℕ × ℕ
      λp. a ≡ p .π₁ * S b + p .π₂
      Z, a
      λm. S (dmAux f' (a ∸ S b) b .π₁), dmAux f' (a ∸ S b) b .π₂
      a ∸ b
      λh0. ⋆
      λm hm. trans
        _
        _
        _
        sym _ _ (leMonusPlus (ltOfMonusSuc _ hm))
        trans
          _
          _
          _
          cong (λu. ℕ) (λw. S b + w) (ih (a ∸ S b) b (fuelStep le (ltOfMonusSuc _ hm)))
          trans
            _
            _
            _
            sym _ _ (plusAssoc (S b) (dmAux f' (a ∸ S b) b .π₁ * S b) (dmAux f' (a ∸ S b) b .π₂))
            cong
              λu. ℕ
              λw. w + dmAux f' (a ∸ S b) b .π₂
              sym _ _ (sucMult (dmAux f' (a ∸ S b) b .π₁) (S b))
    f

dmAuxLe : {f a b : ℕ} → a ≤ f → dmAux f a b .π₂ ≤ b
  using (Core.id.Id.eq,
    Core.id.Id.unfold,
    Natural.div.dmAux.eq,
    Natural.div.natCase.eq,
    Natural.order.≤.unfold)
dmAuxLe =
  λf. ℕ-elim
    λa b le. transport (λw. w ≤ b) (sym _ _ (leZeroInv le)) (leZero b)
    f' ih. λa b le. natCaseElimEqD
      ℕ × ℕ
      λp. p .π₂ ≤ b
      Z, a
      λm. S (dmAux f' (a ∸ S b) b .π₁), dmAux f' (a ∸ S b) b .π₂
      a ∸ b
      λh0. leOfMonusZero h0
      λm hm. ih (a ∸ S b) b (fuelStep le (ltOfMonusSuc _ hm))
    f

-- ===== division, specified =====
-- a itself is always enough fuel, so the two halves of the Euclidean
-- specification hold unconditionally
divModEq : (a b : ℕ) → a ≡ divN a b * S b + modN a b
  using (Natural.div.divN.eq, Natural.div.dmAux.eq, Natural.div.modN.eq)
divModEq = λa b. dmAuxEq (leRefl a)

modLe : (a b : ℕ) → modN a b ≤ b
  using (Core.id.Id.eq, Natural.div.dmAux.eq, Natural.div.modN.eq, Natural.order.≤.unfold)
modLe = λa b. dmAuxLe (leRefl a)

-- and it computes: 7 = 3·2 + 1
divTest : divN 7 (S Z) ≡ 3
  using (Natural.div.divN.eq,
    Natural.div.dmAux.eq,
    Natural.div.natCase.eq,
    Natural.more.pred.eq,
    Natural.more.∸.eq)
divTest = ⋆

modTest : modN 7 (S Z) ≡ S Z
  using (Natural.div.dmAux.eq,
    Natural.div.modN.eq,
    Natural.div.natCase.eq,
    Natural.more.pred.eq,
    Natural.more.∸.eq)
modTest = ⋆

-- ===== uniqueness of the quotient =====
-- no natural is below its own successor's floor
leSucNotSelf : (a : ℕ) → S a ≤ a → 𝟘
  using (Core.id.Id.unfold, Natural.plusSucId, Natural.order.≤.unfold)
leSucNotSelf =
  λa le. zNotS
    sym
      _
      _
      plusCancelL
        _
        _
        _
        trans
          _
          _
          _
          trans
            _
            _
            _
            sym _ _ (trans _ _ (a + S (le .π₁)) (sucPlus a (le .π₁)) ⋆)
            idToEq _ _ _ (le .π₂)
          sym _ _ (plusZeroId a)

-- d·(d'+1) is strictly below (d+1)·(d'+1)
ltProdSuc : (d d' : ℕ) → S (d * S d') ≤ S d * S d' using (Natural.plusSucId, Natural.order.≤.unfold)
ltProdSuc =
  λd d'. (,)
    d'
    eqToId
      _
      _
      trans
        _
        _
        _
        trans _ _ (d * S d' + S d') (sucPlus (d * S d') d') ⋆
        trans _ _ _ (plusComm (S d') (d * S d')) (sym _ _ (sucMult d (S d')))

-- the two quotients' common part cancels, leaving the remainders
-- against a whole multiple of (d+1)(d'+1)
divCancel : (q : ℕ)
  {r k r' d d' : ℕ}
  → ((q * S d + r) * S d' ≡ ((q + k) * S d' + r') * S d) → r * S d' ≡ k * (S d' * S d) + r' * S d
divCancel =
  λq r k r' d d' h. plusCancelL
    _
    _
    _
    trans
      _
      _
      _
      sym
        _
        _
        trans
          _
          _
          _
          multDistribR (S d') (q * S d) r
          cong
            λu. ℕ
            λw. w + r * S d'
            trans
              _
              _
              _
              multAssoc q (S d) (S d')
              cong (λu. ℕ) (λw. q * w) (multComm (S d') (S d))
      trans
        _
        _
        _
        h
        trans
          _
          _
          _
          multDistribR (S d) ((q + k) * S d') r'
          trans
            _
            _
            _
            cong
              λu. ℕ
              λw. w + r' * S d
              trans _ _ _ (multAssoc (q + k) (S d') (S d)) (multDistribR (S d' * S d) q k)
            plusAssoc (q * (S d' * S d)) (k * (S d' * S d)) (r' * S d)

-- ...and that multiple has to be zero: one more copy of (d+1)(d'+1)
-- would already exceed r·(d'+1), which r ≤ d keeps strictly below it
divKZero : (r d r' d' : ℕ) {k : ℕ} → r ≤ d → (r * S d' ≡ k * (S d' * S d) + r' * S d) → k ≡ Z
  using (Natural.order.≤.unfold)
divKZero =
  λr d r' d' k. ℕ-elim
    λb h. ⋆
    k' rec. λb h. 𝟘-elim
      leSucNotSelf
        _
        leTrans
          _
          _
          _
          ltProdSuc d d'
          leTrans
            S d * S d'
            _
            _
            (,)
              k' * (S d' * S d) + r' * S d
              eqToId
                _
                _
                sym
                  _
                  _
                  trans
                    _
                    _
                    _
                    h
                    trans
                      _
                      _
                      _
                      cong (λu. ℕ) (λw. w + r' * S d) (sucMult k' (S d' * S d))
                      trans
                        _
                        _
                        _
                        plusAssoc (S d' * S d) (k' * (S d' * S d)) (r' * S d)
                        cong
                          λu. ℕ
                          λw. w + (k' * (S d' * S d) + r' * S d)
                          multComm (S d) (S d')
            leMultMonoR (S d') b
    k

-- ===== equal rationals have equal floors =====
divNUniqueLe : {m d m' d' : ℕ}
  → (m * S d' ≡ m' * S d) → divN m d ≤ divN m' d' → divN m d ≡ divN m' d'
  using (Natural.order.≤.unfold)
divNUniqueLe =
  λm d m' d' h le. trans
    _
    _
    _
    sym
      _
      _
      trans
        _
        _
        _
        cong
          λu. ℕ
          λw. divN m d + w
          divKZero
            _
            _
            _
            _
            modLe m d
            divCancel
              _
              trans
                _
                _
                _
                trans _ _ _ (cong (λu. ℕ) (λw. w * S d') (sym _ _ (divModEq m d))) h
                cong
                  λu. ℕ
                  λw. w * S d
                  trans
                    _
                    _
                    _
                    divModEq m' d'
                    cong (λu. ℕ) (λw. w * S d' + modN m' d') (sym _ _ (idToEq _ _ _ (le .π₂)))
        plusZeroId (divN m d)
    idToEq _ _ _ (le .π₂)

divNUnique : {m d m' d' : ℕ} → (m * S d' ≡ m' * S d) → divN m d ≡ divN m' d'
divNUnique =
  λm d m' d' h. ⊎-elim
    le. divNUniqueLe h le
    le. sym _ _ (divNUniqueLe (sym _ _ h) le)
    leTotal (divN m d) (divN m' d')

-- a is strictly below (a div (b+1) + 1)·(b+1): the remainder is at
-- most b, so one more divisor overshoots
leDivSuc : (m b : ℕ) → m ≤ S (divN m b) * S b using (Natural.order.≤.unfold)
leDivSuc =
  λm b. (,)
    S b ∸ modN m b
    eqToId
      _
      _
      trans
        _
        _
        _
        trans
          _
          _
          _
          cong (λu. ℕ) (λw. w + (S b ∸ modN m b)) (divModEq m b)
          plusAssoc (divN m b * S b) (modN m b) (S b ∸ modN m b)
        trans
          _
          _
          _
          cong (λu. ℕ) (λw. divN m b * S b + w) (leMonusPlus (leStep (modLe m b)))
          trans _ _ _ (plusComm (S b) (divN m b * S b)) (sym _ _ (sucMult (divN m b) (S b)))