Rat.half

-- The HALVING law for the harmonic bounds: 1/(2n+2) + 1/(2n+2) is
-- 1/(n+1). This is the one arithmetic fact that makes addition on ℝ
-- possible — Bishop's regular sequences carry their modulus in the
-- type, so a sum has to be SAMPLED at doubled indices, and the
-- doubling is paid for exactly here.
--
-- `dbl n` is the index whose bound halves n's: qInvNat k denotes
-- 1/(k+1), so qInvNat (dbl n) is 1/(2n+2).
--
-- The proof never touches a difference-pair representative. Both
-- fractions are 1/·, so cross-multiplication reduces to a RING
-- identity in Int — ((g+g)+(g+g))·g ≡ (g+g)·(g+g), two
-- distributivities — once the denominators are pushed through
-- nzToInt.

import Natural (+, *, sucPlus, plusSucId)
import Int (Int, intOne)
import Int.add (+)
import Int.mul (*, intMulOneL, intMulDistribL, intMulDistribR)
import Int.order (intOfNat, intOfNatPlus)
import Rat.frac (NZ, nzPos, nzMul, nzToInt, Rat, mkRat, num, den, ratAdd, intScale, intScaleToInt)
import Rat (Q, qcls, +, qZero, qAddCls, RatR, dInt, clsEqOfRel, nzToIntMul)
import Rat.order (≤)
import Rat.bound (leQSelfAdd, leQZeroOfPos)
import Real (qInvNat, qInvNatPos)
import Core.equality (trans, sym, cong, transport)

dbl : ℕ → ℕ
dbl = λn. S (n + n)

-- ===== the Int identity behind it =====
-- 4g² two ways: distribute right, then collect left
intDblSquare : (g : Int) → (g + g + (g + g)) * g ≡ (g + g) * (g + g) using (Int.Int.unfold)
intDblSquare =
  λg. (g + g + (g + g)) * g
    ≡⟨ intMulDistribR (g + g) (g + g) g ⟩ (g + g) * g + (g + g) * g
    ≡⟨ sym _ _ (intMulDistribL (g + g) g g) ⟩ (g + g) * (g + g)

-- ===== pushing the doubled index through nzToInt =====
nzPosInt : (k : ℕ) → nzToInt (nzPos k) ≡ intOfNat (S k)
  using (Int.order.intOfNat.eq, Int.Int.unfold, Rat.frac.nzPos.eq, Rat.frac.nzToInt.eq)
nzPosInt = λk. ⋆

-- S (dbl n) is 2·(n+1) — one successor moved across a sum
sucDbl : (n : ℕ) → S (dbl n) ≡ S n + S n using (sucPlus.rw, plusSucId.rw, Rat.half.dbl.eq)
sucDbl = λn. ⋆

-- the doubled denominator IS the original one, doubled in Int
nzPosDbl : (n : ℕ) → nzToInt (nzPos (dbl n)) ≡ nzToInt (nzPos n) + nzToInt (nzPos n)
  using (Int.Int.unfold)
nzPosDbl =
  λn. nzToInt (nzPos (dbl n))
    ≡⟨ nzPosInt (dbl n) ⟩ intOfNat (S (dbl n))
    ≡⟨ cong (λw. Int) (λw. intOfNat w) (sucDbl n) ⟩ intOfNat (S n + S n)
    ≡⟨ sym _ _ (intOfNatPlus (S n) (S n)) ⟩ intOfNat (S n) + intOfNat (S n)
    ≡⟨ sym _ _ (cong (λw. Int) (λw. w + w) (nzPosInt n)) ⟩ nzToInt (nzPos n) + nzToInt (nzPos n)

-- ===== the cross-multiplication =====
-- the unit fraction at index k, as a raw fraction
unitRat : ℕ → Rat using (Rat.frac.Rat.unfold)
unitRat = λk. mkRat intOne (nzPos k)

-- scaling 1 by d is just d
scaleOne : {d : NZ} → intScale d intOne ≡ nzToInt d using (Int.intOne.eq, Rat.frac.NZ.unfold)
scaleOne = λd. intScaleToInt

-- the numerator of 1/d + 1/d is d + d
unitAddNum : {d : NZ} → num (ratAdd (mkRat intOne d) (mkRat intOne d)) ≡ nzToInt d + nzToInt d
  using (Int.intNeg.eq,
    Int.intOne.eq,
    Int.add.+.eq,
    Rat.frac.NZ.unfold,
    Rat.frac.den.eq,
    Rat.frac.intScale.eq,
    Rat.frac.intScaleN.eq,
    Rat.frac.mkRat.eq,
    Rat.frac.num.eq,
    Rat.frac.ratAdd.eq)
unitAddNum = λd. cong (λw. Int) (λw. w + w) {intScale d intOne} {nzToInt d} scaleOne

-- ...and its denominator is d·d
unitAddDen : {d : NZ} → dInt (ratAdd (mkRat intOne d) (mkRat intOne d)) ≡ nzToInt d * nzToInt d
  using (Int.intOne.eq,
    Int.add.+.eq,
    Rat.frac.NZ.unfold,
    Rat.frac.den.eq,
    Rat.frac.intScale.eq,
    Rat.frac.mkRat.eq,
    Rat.frac.num.eq,
    Rat.frac.nzMul.eq,
    Rat.frac.nzNeg.eq,
    Rat.frac.nzPos.eq,
    Rat.frac.nzToInt.eq,
    Rat.frac.ratAdd.eq,
    Rat.dInt.eq)
unitAddDen = λd. nzToIntMul d d

-- the unit fraction's numerator and denominator, computed
unitNum : (k : ℕ) → num (unitRat k) ≡ intOne
  using (Int.Int.unfold, Rat.half.unitRat.eq, Rat.frac.mkRat.eq, Rat.frac.num.eq)
unitNum = λk. ⋆

unitDen : (k : ℕ) → dInt (unitRat k) ≡ nzToInt (nzPos k)
  using (Int.Int.unfold, Rat.half.unitRat.eq, Rat.frac.den.eq, Rat.frac.mkRat.eq, Rat.dInt.eq)
unitDen = λk. ⋆

-- unitAddNum/unitAddDen, respelled at unitRat
unitAddNumU : (k : ℕ) → num (ratAdd (unitRat k) (unitRat k)) ≡ nzToInt (nzPos k) + nzToInt (nzPos k)
  using (Rat.half.unitRat.eq)
unitAddNumU = λk. unitAddNum

unitAddDenU : (k : ℕ)
  → dInt (ratAdd (unitRat k) (unitRat k)) ≡ nzToInt (nzPos k) * nzToInt (nzPos k)
  using (Rat.half.unitRat.eq)
unitAddDenU = λk. unitAddDen

-- the cross-multiplication itself, as one chain in folded vocabulary
halfCross : {n : ℕ}
  → num (ratAdd (unitRat (dbl n)) (unitRat (dbl n))) * dInt (unitRat n)
    ≡ num (unitRat n) * dInt (ratAdd (unitRat (dbl n)) (unitRat (dbl n)))
  using (Int.Int.unfold)
halfCross =
  λn. num (ratAdd (unitRat (dbl n)) (unitRat (dbl n))) * dInt (unitRat n)
    ≡⟨ cong (λw. Int) (λw. num (ratAdd (unitRat (dbl n)) (unitRat (dbl n))) * w) (unitDen n) ⟩
      num (ratAdd (unitRat (dbl n)) (unitRat (dbl n))) * nzToInt (nzPos n)
    ≡⟨ cong (λw. Int) (λw. w * nzToInt (nzPos n)) (unitAddNumU (dbl n)) ⟩
      (nzToInt (nzPos (dbl n)) + nzToInt (nzPos (dbl n))) * nzToInt (nzPos n)
    ≡⟨ cong (λw. Int) (λw. (w + w) * nzToInt (nzPos n)) (nzPosDbl n) ⟩
      (nzToInt (nzPos n) + nzToInt (nzPos n) + (nzToInt (nzPos n) + nzToInt (nzPos n)))
        * nzToInt (nzPos n)
    ≡⟨ intDblSquare (nzToInt (nzPos n)) ⟩
      (nzToInt (nzPos n) + nzToInt (nzPos n)) * (nzToInt (nzPos n) + nzToInt (nzPos n))
    ≡⟨ sym _ _ (cong (λw. Int) (λw. w * w) (nzPosDbl n)) ⟩
      nzToInt (nzPos (dbl n)) * nzToInt (nzPos (dbl n))
    ≡⟨ sym _ _ (unitAddDenU (dbl n)) ⟩ dInt (ratAdd (unitRat (dbl n)) (unitRat (dbl n)))
    ≡⟨ sym _ _ (intMulOneL (dInt (ratAdd (unitRat (dbl n)) (unitRat (dbl n))))) ⟩
      intOne * dInt (ratAdd (unitRat (dbl n)) (unitRat (dbl n)))
    ≡⟨ cong (λw. Int) (λw. w * dInt (ratAdd (unitRat (dbl n)) (unitRat (dbl n)))) (unitNum n) ⟩
      num (unitRat n) * dInt (ratAdd (unitRat (dbl n)) (unitRat (dbl n)))

-- 1/(2n+2) + 1/(2n+2) ~ 1/(n+1): the relation is halfCross, decoded
halfRel : (n : ℕ) → RatR (ratAdd (unitRat (dbl n)) (unitRat (dbl n))) (unitRat n)
  using (Rat.RatR.eq, Rat.RatR.unfold)
halfRel = λn. halfCross

-- ===== the halving law =====
qInvHalf : (n : ℕ) → qInvNat (dbl n) + qInvNat (dbl n) ≡ qInvNat n
  using (Int.intOne.eq,
    Rat.half.dbl.eq,
    Rat.half.unitRat.eq,
    Rat.frac.mkRat.eq,
    Rat.frac.nzPos.eq,
    Rat.frac.ratAdd.eq,
    Rat.+.eq,
    Rat.qcls.eq,
    Real.qInvNat.eq)
qInvHalf =
  λn. trans
    _
    qcls (ratAdd (unitRat (dbl n)) (unitRat (dbl n)))
    _
    qAddCls (unitRat (dbl n)) (unitRat (dbl n))
    clsEqOfRel _ _ (halfRel n)

-- ===== the halved bound, as an inequality =====
leQZeroInvNat : (n : ℕ) → qZero ≤ qInvNat n using (Rat.order.≤.unfold, Rat.order.NonNegS.unfold)
leQZeroInvNat = λn. leQZeroOfPos _ (qInvNatPos n)

-- 1/(2n+2) ≤ 1/(n+1): the halved bound is at most the original,
-- because the other half is nonnegative
leQInvDbl : (n : ℕ) → qInvNat (dbl n) ≤ qInvNat n
  using (Rat.order.≤.unfold, Rat.order.NonNegS.unfold)
leQInvDbl =
  λn. transport
    λw. qInvNat (dbl n) ≤ w
    qInvHalf n
    leQSelfAdd (qInvNat (dbl n)) (leQZeroInvNat (dbl n))