Real.recip

-- THE RECIPROCAL of a POSITIVE real.
--
--   realRecip : (x : Real) → (PosR x) → Real
--
-- This is where the MODULUS is unavoidable. √ could make do with a
-- nonnegative presentation, because its estimate only needed the
-- samples to be ≥ 0. 1/x needs them bounded AWAY from zero — a rate —
-- and a rate is exactly what PosR carries and what (LtR 0 x)
-- cannot give up.
--
-- With modulus k, sample at
--
--   ψ n = mulIdx C n,   C = mulIdx (dbl k) (dbl k),  so ψ n + 1 = (2k+2)²(n+1)
--
-- Two things land exactly, and they are the same two facts:
--
--   * rBound (ψ m) (ψ n) = 1/(2k+2)² · rBound m n, so multiplying by
--     the bound (2k+2) twice — once for each reciprocal in the
--     difference identity — returns rBound m n on the nose.
--   * rBound (ψ n) (ψ n) ≤ 1/(2k+2), which is what makes every sample
--     at least 1/(2k+2) and hence invertible.

import Natural (+, *)
import Rat (Q, +, *, qNeg, qZero, qOne, qAddComm, qAddAssoc, qAddNegR, qAddZeroR, qMulComm, qMulAssoc, qMulOneL, qMulOneR, qDistribL, qDistribR)
import Rat.order (Sign, sPos, sgnQ, ≤, leQTrans, leQMulMono, leQPlusMonoL)
import Rat.bound (Bnd, bndEq, bndEqB, bndSubSym, bndVia, bndWeaken, bndOfBothLe, leQAdd, leQNegFlip, leQZeroOfNonNeg, qNegZeroQ)
import Rat.half (dbl, qInvHalf, leQZeroInvNat)
import Rat.nat (qOfNat, qOfNatAdd, qOfNatMul)
import Rat.ceil (leQZeroNat)
import Rat.arch (leQSubShift, leQAddShift, leQUnsquash)
import Rat.abs (bndMul, leQZeroMul)
import Rat.sqrt (qMulInvNat, qInvNatMulIdx, qMulSwap3, qMulLeftComm)
import Rat.algInv (qInv)
import Rat.invOrder (qInvDiff, leQInvBound, leQZeroInv, sgnQPosOfLe)
import Real (Real, RSeq, REq, Regular, qInvNat, qInvNatPos, rBound, leQZeroBound, realOfQ, realZero)
import Real.neg (seqOf, regOf, regBnd, bndReg, reqOf, realEqOfREq)
import Core.bracket (br, brElim)
import Real.order (≤)
import Real.mul (mulIdx, sucMulIdx, *)
import Real.pos (leQOneInv, PosR, PosPayload, posRep, posK, posBoundRep, posRepEq, reqOfClassEq)
import Real.ring (leQInvIdx, reqOfClose, bndDeep, leQZeroAdd)
import Core.equality (trans, sym, cong, transport, transportP)

recipC : ℕ → ℕ
recipC = λk. mulIdx (dbl k) (dbl k)

recipIdx : ℕ → ℕ → ℕ
recipIdx = λk n. mulIdx (recipC k) n

recipM : ℕ → Q
recipM = λk. qOfNat (S (dbl k))

recipLo : ℕ → Q
recipLo = λk. qInvNat (dbl k)

recipLoPos : (k : ℕ) → sgnQ (recipLo k) ≡ sPos using (recipLo.eq)
recipLoPos = λk. qInvNatPos (dbl k)

recipLoM : (k : ℕ) → recipLo k * recipM k ≡ qOne using (recipLo.eq, recipM.eq)
recipLoM =
  λk. trans
    qInvNat (dbl k) * qOfNat (S (dbl k))
    _
    _
    qMulComm (qInvNat (dbl k)) (qOfNat (S (dbl k)))
    qMulInvNat

leQZeroRecipM : (k : ℕ) → qZero ≤ recipM k using (recipM.eq)
leQZeroRecipM = λk. leQZeroNat (S (dbl k))

-- (2k+2)² is S C, so 1/(C+1) cancels it
recipMSq : (k : ℕ) → recipM k * recipM k ≡ qOfNat (S (recipC k)) using (recipM.eq, recipC.eq)
recipMSq =
  λk. trans
    qOfNat (S (dbl k)) * qOfNat (S (dbl k))
    _
    qOfNat (S (mulIdx (dbl k) (dbl k)))
    qOfNatMul (S (dbl k)) (S (dbl k))
    cong (λw. Q) (λw. qOfNat w) (sym _ _ (sucMulIdx (dbl k) (dbl k)))

-- ===== the scaling, exactly =====
leQMulMonoL : {x y c : Q} → x ≤ y → qZero ≤ c → c * x ≤ c * y
leQMulMonoL =
  λx y c h hc. transport
    λw. w ≤ c * y
    qMulComm x c
    transport (λw. x * c ≤ w) (qMulComm y c) (leQMulMono _ _ h hc)

rBoundRecipEq : (k m n : ℕ) → rBound (recipIdx k m) (recipIdx k n) ≡ qInvNat (recipC k) * rBound m n
  using (Real.rBound.eq, recipIdx.eq)
rBoundRecipEq =
  λk m n. trans
    qInvNat (mulIdx (recipC k) m) + qInvNat (mulIdx (recipC k) n)
    _
    qInvNat (recipC k) * (qInvNat m + qInvNat n)
    trans
      _
      _
      _
      cong
        λw. Q
        λw. w + qInvNat (mulIdx (recipC k) n)
        {qInvNat (mulIdx (recipC k) m)}
        {qInvNat (recipC k) * qInvNat m}
        qInvNatMulIdx
      cong
        λw. Q
        λw. qInvNat (recipC k) * qInvNat m + w
        {qInvNat (mulIdx (recipC k) n)}
        {qInvNat (recipC k) * qInvNat n}
        qInvNatMulIdx
    sym _ _ (qDistribL (qInvNat (recipC k)) (qInvNat m) (qInvNat n))

recipScaleEq : (k m n : ℕ) → rBound (recipIdx k m) (recipIdx k n) * recipM k * recipM k ≡ rBound m n
recipScaleEq =
  λk m n. rBound (recipIdx k m) (recipIdx k n) * recipM k * recipM k
    ≡⟨ cong (λw. Q) (λw. w * recipM k * recipM k) (rBoundRecipEq k m n) ⟩
      qInvNat (recipC k) * rBound m n * recipM k * recipM k
    ≡⟨ qMulAssoc (qInvNat (recipC k) * rBound m n) (recipM k) (recipM k) ⟩
      qInvNat (recipC k) * rBound m n * (recipM k * recipM k)
    ≡⟨ cong (λw. Q) (λw. qInvNat (recipC k) * rBound m n * w) (recipMSq k) ⟩
      qInvNat (recipC k) * rBound m n * qOfNat (S (recipC k))
    ≡⟨ qMulSwap3 (qInvNat (recipC k)) (rBound m n) (qOfNat (S (recipC k))) ⟩
      qInvNat (recipC k) * qOfNat (S (recipC k)) * rBound m n
    ≡⟨ cong
      λw. Q
      λw. w * rBound m n
      trans _ _ qOne (qMulComm (qInvNat (recipC k)) (qOfNat (S (recipC k)))) qMulInvNat ⟩
      qOne * rBound m n
    ≡⟨ qMulOneL (rBound m n) ⟩ rBound m n

-- ===== every sample is at least 1/(2k+2) =====
sqSumRearr : (A N : Q) → A * A * (N + N) ≡ (A + A) * (A * N)
sqSumRearr =
  λA N. A * A * (N + N)
    ≡⟨ qDistribL (A * A) N N ⟩ A * A * N + A * A * N
    ≡⟨ cong (λw. Q) (λw. w + w) (qMulAssoc A A N) ⟩ A * (A * N) + A * (A * N)
    ≡⟨ sym _ _ (qDistribR (A * N) A A) ⟩ (A + A) * (A * N)

rBoundSelfEq : (n : ℕ) → rBound n n ≡ qInvNat n + qInvNat n using (Real.rBound.eq)
rBoundSelfEq = λn. ⋆

recipCInv : (k : ℕ) → qInvNat (recipC k) ≡ qInvNat (dbl k) * qInvNat (dbl k) using (recipC.eq)
recipCInv = λk. qInvNatMulIdx

rBoundRecipSelf : (k n : ℕ)
  → rBound (recipIdx k n) (recipIdx k n) ≡ qInvNat k * (recipLo k * qInvNat n)
  using (recipLo.eq)
rBoundRecipSelf =
  λk n. rBound (recipIdx k n) (recipIdx k n)
    ≡⟨ rBoundRecipEq k n n ⟩ qInvNat (recipC k) * rBound n n
    ≡⟨ cong (λw. Q) (λw. qInvNat (recipC k) * w) (rBoundSelfEq n) ⟩
      qInvNat (recipC k) * (qInvNat n + qInvNat n)
    ≡⟨ cong (λw. Q) (λw. w * (qInvNat n + qInvNat n)) (recipCInv k) ⟩
      qInvNat (dbl k) * qInvNat (dbl k) * (qInvNat n + qInvNat n)
    ≡⟨ sqSumRearr (qInvNat (dbl k)) (qInvNat n) ⟩
      (qInvNat (dbl k) + qInvNat (dbl k)) * (qInvNat (dbl k) * qInvNat n)
    ≡⟨ cong (λw. Q) (λw. w * (qInvNat (dbl k) * qInvNat n)) (qInvHalf k) ⟩
      qInvNat k * (qInvNat (dbl k) * qInvNat n)

leQRBoundRecip : (k n : ℕ) → rBound (recipIdx k n) (recipIdx k n) ≤ recipLo k using (recipLo.eq)
leQRBoundRecip =
  λk n. transport
    λw. w ≤ qInvNat (dbl k)
    sym _ (qInvNat k * (qInvNat (dbl k) * qInvNat n)) (rBoundRecipSelf k n)
    leQTrans
      _
      _
      _
      leQMulMonoL
        transport
          λw. qInvNat (dbl k) * qInvNat n ≤ w
          qMulOneR (qInvNat (dbl k))
          leQMulMonoL (leQOneInv n) (leQZeroInvNat (dbl k))
        leQZeroInvNat k
      transport
        λw. qInvNat k * qInvNat (dbl k) ≤ w
        qMulOneL (qInvNat (dbl k))
        leQMulMono _ _ (leQOneInv k) (leQZeroInvNat (dbl k))

-- ===== the modulus, spent =====
--
-- 1/(k+1) is TWICE 1/(2k+2), and the harmonic bound at the sampled
-- index is at most 1/(2k+2), so the sample keeps at least the other
-- half. That half is the whole point of the modulus.
halfSubEq : (A : Q) → A + A + qNeg A ≡ A
halfSubEq =
  λA. trans
    _
    _
    _
    qAddAssoc A A (qNeg A)
    trans _ _ _ (cong (λw. Q) (λw. A + w) (qAddNegR A)) (qAddZeroR A)

posSampleLo : {p : RSeq}
  {k : ℕ}
  → realOfQ (qInvNat k) ≤ class p → (n : ℕ) → recipLo k ≤ seqOf p (recipIdx k n)
  using (Real.Real.unfold,
    Real.RSeq.unfold,
    Real.Regular.unfold,
    Real.realOfQ.eq,
    Real.constReg.eq,
    Real.order.≤.eq,
    Real.order.RLeP.eq,
    Real.neg.seqOf.eq,
    recipLo.eq)
posSampleLo =
  λp k h n. leQUnsquash
    qInvNat (dbl k)
    squash-elim
      h
      u. ⋆
        transport
          λw. w ≤ seqOf p (recipIdx k n)
          halfSubEq (qInvNat (dbl k))
          transport
            λw. w + qNeg (qInvNat (dbl k)) ≤ seqOf p (recipIdx k n)
            sym _ _ (qInvHalf k)
            leQSubShift
              _
              _
              transport
                λw. qInvNat k ≤ w
                qAddComm (seqOf p (recipIdx k n)) (qInvNat (dbl k))
                leQTrans
                  qInvNat k
                  seqOf p (recipIdx k n) + rBound (recipIdx k n) (recipIdx k n)
                  _
                  leQAddShift (qInvNat k) (u (recipIdx k n))
                  leQPlusMonoL _ (qInvNat (dbl k)) (seqOf p (recipIdx k n)) (leQRBoundRecip k n)

posSampleSgn : (p : RSeq)
  (k : ℕ)
  (h : realOfQ (qInvNat k) ≤ class p)
  (n : ℕ)
  → sgnQ (seqOf p (recipIdx k n)) ≡ sPos
  using (Real.Real.unfold, Real.RSeq.unfold, Real.Regular.unfold)
posSampleSgn = λp k h n. sgnQPosOfLe _ (recipLoPos k) (posSampleLo h n)

-- ===== the estimate =====
bndOfNonNegLe : {M u : Q} → qZero ≤ M → qZero ≤ u → u ≤ M → Bnd M u using (Rat.bound.Bnd.unfold)
bndOfNonNegLe =
  λM u hM0 h0 hM. leQTrans _ _ _ (transport (λw. qNeg M ≤ w) qNegZeroQ (leQNegFlip _ _ hM0)) h0, hM

-- |1/a − 1/b| ≤ d·Ma·Mb, when a ≥ ca, b ≥ cb and each c inverts its M.
-- The two bounds are kept SEPARATE because the modulus-independence
-- proof compares samples taken with different moduli.
bndRecipDiff : {a b : Q}
  (ca cb : Q)
  {Ma Mb d : Q}
  → (sgnQ ca ≡ sPos)
    → (sgnQ cb ≡ sPos)
      → ca ≤ a
        → cb ≤ b
          → (ca * Ma ≡ qOne)
            → (cb * Mb ≡ qOne)
              → qZero ≤ Ma
                → qZero ≤ Mb
                  → qZero ≤ d → Bnd d (b + qNeg a) → Bnd (d * Ma * Mb) (qInv a + qNeg (qInv b))
bndRecipDiff =
  λa b ca cb Ma Mb d hca hcb ha hb hcMa hcMb hMa0 hMb0 hd0 gap. bndEq
    _
    _
    qInvDiff (sgnQPosOfLe _ hca ha) (sgnQPosOfLe _ hcb hb)
    bndMul
      leQZeroMul hd0 hMa0
      hMb0
      bndMul
        hd0
        hMa0
        gap
        bndOfNonNegLe hMa0 (leQZeroInv (sgnQPosOfLe _ hca ha)) (leQInvBound _ hca ha hcMa)
      bndOfNonNegLe hMb0 (leQZeroInv (sgnQPosOfLe _ hcb hb)) (leQInvBound _ hcb hb hcMb)

-- ===== the sequence, and its regularity =====
recipSeq : RSeq → ℕ → ℕ → Q
recipSeq = λp k n. qInv (seqOf p (recipIdx k n))

recipRegAt : (p : RSeq)
  (k : ℕ)
  (h : realOfQ (qInvNat k) ≤ class p)
  (m n : ℕ)
  → Bnd (rBound m n) (recipSeq p k m + qNeg (recipSeq p k n))
  using (Real.Real.unfold, Real.RSeq.unfold, Real.Regular.unfold, recipSeq.eq)
recipRegAt =
  λp k h m n. bndEqB
    _
    _
    qInv (seqOf p (recipIdx k m)) + qNeg (qInv (seqOf p (recipIdx k n)))
    recipScaleEq k m n
    bndRecipDiff
      _
      _
      recipLoPos k
      recipLoPos k
      posSampleLo h m
      posSampleLo h n
      recipLoM k
      recipLoM k
      leQZeroRecipM k
      leQZeroRecipM k
      leQZeroBound (recipIdx k m) (recipIdx k n)
      bndSubSym _ _ _ (regBnd _ (regOf p) (recipIdx k m) (recipIdx k n))

rRecip : (p : RSeq) (k : ℕ) → realOfQ (qInvNat k) ≤ class p → RSeq
  using (Real.Real.unfold, Real.RSeq.unfold, Real.Regular.unfold)
rRecip = λp k h. recipSeq p k, bndReg (recipRegAt _ _ h)

-- ===== well-definedness and independence, in one =====
--
-- Two payloads for the same x differ in BOTH representative and
-- modulus. The samples then sit at different depths, so the gap is
-- bounded by twice rBound n n rather than once, and the constant that
-- comes out is 2·(2k+2)(2k'+2). reqOfClose absorbs any constant, which
-- is exactly why the exact scaling is needed only for regularity.
recipK : ℕ → ℕ → ℕ
recipK = λk k'. S (dbl k) * S (dbl k') + S (dbl k) * S (dbl k')

recipMM : (k k' : ℕ) → recipM k * recipM k' ≡ qOfNat (S (dbl k) * S (dbl k')) using (recipM.eq)
recipMM = λk k'. qOfNatMul (S (dbl k)) (S (dbl k'))

recipKEq : (k k' : ℕ)
  → qOfNat (S (dbl k) * S (dbl k')) + qOfNat (S (dbl k) * S (dbl k')) ≡ qOfNat (recipK k k')
  using (recipK.eq)
recipKEq = λk k'. qOfNatAdd (S (dbl k) * S (dbl k')) (S (dbl k) * S (dbl k'))

recipBoundEq : (k k' : ℕ) (R : Q) → (R + R) * recipM k * recipM k' ≡ qOfNat (recipK k k') * R
recipBoundEq =
  λk k' R. (R + R) * recipM k * recipM k'
    ≡⟨ qMulAssoc (R + R) (recipM k) (recipM k') ⟩ (R + R) * (recipM k * recipM k')
    ≡⟨ cong (λw. Q) (λw. (R + R) * w) (recipMM k k') ⟩ (R + R) * qOfNat (S (dbl k) * S (dbl k'))
    ≡⟨ qDistribR (qOfNat (S (dbl k) * S (dbl k'))) R R ⟩
      R * qOfNat (S (dbl k) * S (dbl k')) + R * qOfNat (S (dbl k) * S (dbl k'))
    ≡⟨ cong (λw. Q) (λw. w + w) (qMulComm R (qOfNat (S (dbl k) * S (dbl k')))) ⟩
      qOfNat (S (dbl k) * S (dbl k')) * R + qOfNat (S (dbl k) * S (dbl k')) * R
    ≡⟨ sym _ _ (qDistribR R (qOfNat (S (dbl k) * S (dbl k'))) (qOfNat (S (dbl k) * S (dbl k')))) ⟩
      (qOfNat (S (dbl k) * S (dbl k')) + qOfNat (S (dbl k) * S (dbl k'))) * R
    ≡⟨ cong (λw. Q) (λw. w * R) (recipKEq k k') ⟩ qOfNat (recipK k k') * R

recipCloseAt : (p p' : RSeq)
  (k k' n : ℕ)
  → recipLo k ≤ seqOf p (recipIdx k n)
    → recipLo k' ≤ seqOf p' (recipIdx k' n)
      → Bnd (rBound n n + rBound n n) (seqOf p' (recipIdx k' n) + qNeg (seqOf p (recipIdx k n)))
        → Bnd (qOfNat (recipK k k') * rBound n n) (recipSeq p k n + qNeg (recipSeq p' k' n))
  using (recipSeq.eq)
recipCloseAt =
  λp p' k k' n lo lo' gap. bndEqB
    _
    _
    qInv (seqOf p (recipIdx k n)) + qNeg (qInv (seqOf p' (recipIdx k' n)))
    recipBoundEq k k' (rBound n n)
    bndRecipDiff
      _
      _
      recipLoPos k
      recipLoPos k'
      lo
      lo'
      recipLoM k
      recipLoM k'
      leQZeroRecipM k
      leQZeroRecipM k'
      leQZeroAdd (leQZeroBound n n) (leQZeroBound n n)
      gap

deepR : (k n : ℕ) → qInvNat (recipIdx k n) ≤ qInvNat n using (recipIdx.eq)
deepR = λk n. leQInvIdx (recipC k) n

recipWDGen : {p p' : RSeq}
  (k k' : ℕ)
  (h : realOfQ (qInvNat k) ≤ class p)
  (h' : realOfQ (qInvNat k') ≤ class p')
  → REq p p' → REq (rRecip _ _ h) (rRecip _ _ h')
  using (rRecip.eq,
    Rat.bound.Bnd.unfold,
    Real.REq.unfold,
    Real.RSeq.unfold,
    Real.Real.unfold,
    Real.Regular.unfold,
    Real.rBound.eq,
    Real.neg.seqOf.eq,
    recipIdx.eq,
    recipSeq.eq)
recipWDGen =
  λp p' k k' h h' he. squash-elim
    he
    u. reqOfClose
      recipK k k'
      λn. recipCloseAt
        _
        _
        _
        _
        _
        posSampleLo h n
        posSampleLo h' n
        bndVia
          _
          _
          _
          _
          _
          bndDeep (seqOf p') (regBnd _ (regOf p')) _ _ _ (deepR k' n) (deepR k n)
          bndWeaken
            rBound (recipIdx k n) (recipIdx k n)
            rBound n n
            _
            leQAdd _ _ _ _ (deepR k n) (deepR k n)
            bndSubSym
              rBound (recipIdx k n) (recipIdx k n)
              seqOf p (recipIdx k n)
              seqOf p' (recipIdx k n)
              u (recipIdx k n)

-- ===== the descent, and division =====
recipFromPayload : (x : Real) → PosPayload x → Real
  using (Real.Real.unfold, Real.RSeq.unfold, Real.Regular.unfold)
recipFromPayload = λx d. class (rRecip _ _ (posBoundRep _ d))

recipPayloadConst : (x : Real) (d d' : PosPayload x) → recipFromPayload _ d ≡ recipFromPayload _ d'
  using (Real.Real.unfold, Real.RSeq.unfold, Real.Regular.unfold, recipFromPayload.eq)
recipPayloadConst =
  λx d d'. realEqOfREq
    _
    _
    recipWDGen
      _
      _
      posBoundRep _ d
      posBoundRep _ d'
      reqOfClassEq (trans _ _ _ (posRepEq _ d) (sym _ _ (posRepEq _ d')))

realRecip : (x : Real) → PosR x → Real using (PosR.eq)
realRecip = λx w. brElim (recipFromPayload x) (recipPayloadConst x) w

realRecipRep : (x : Real)
  (d : PosPayload x)
  → realRecip _ (br d) ≡ class (rRecip _ _ (posBoundRep _ d))
  using (Real.Real.unfold,
    Real.RSeq.unfold,
    Real.Regular.unfold,
    PosR.eq,
    Core.bracket.br.eq,
    Core.bracket.brElim.eq,
    recipFromPayload.eq,
    Real.recip.realRecip.eq)
realRecipRep = λx d. ⋆

-- ...and division is multiplication by the reciprocal
realDiv : (u x : Real) → PosR x → Real
realDiv = λu x w. u * realRecip _ w