-- Real analysis: IntegralOrderBound
-- Part of Cumulus, the Grace17 mathematics library.
-- Copyright (c) 2026 Grant Passmore.
-- Released under the Apache License 2.0. See LICENSE in the repository root.
-- Authors: Grant Passmore

module Grace.Cumulus.Real.Analysis.Integration.IntegralOrderBound
import Set.Choice
import RiemannSums
import CommonRefinementComparison
import Squares

-- Section: IntegralOrderBound
theorem Nat.not_less_zero (i : Nat)
  shows ¬(i < 0)
proof
  intro before
  have weak := Grace.Core.Nat.Order.LessEqual.zero at i
  have typedZero : Nat(0) by typecheck
  have cycle := Less.transitive_LessEqual at i, 0, i with
    assumption, typedZero, assumption, before, weak
  have impossible := Less.irreflexive at i
  exact impossible with cycle

-- A singleton interval has the empty tagged partition at every mesh size.
theorem FineTaggedPartitionsExist.point (A : Real)
  shows FineTaggedPartitionsExist(A, A)
proof
  intros D realD positive
  have sequence := constantSequence.sequence at A
  have atZero := constantSequence.value at A, 0 with *
  have partition : IntervalPartition(constantSequence(A), A, A, 0)
    constructor
    . exact sequence
    . constructor
      . exact atZero
      . constructor
        . exact atZero
        . intros i naturalI before
          exfalso
          have typedI := nat.omega.backward with naturalI
          have impossible := Nat.not_less_zero with typedI
          exact impossible with before
  have tagged : TaggedPartition(constantSequence(A), constantSequence(A), A, A, 0)
    constructor
    . exact partition
    . constructor
      . exact sequence
      . intros i naturalI before
        exfalso
        have typedI := nat.omega.backward with naturalI
        have impossible := Nat.not_less_zero with typedI
        exact impossible with before
  have fine : FinePartition(constantSequence(A), A, A, 0, D)
    constructor
    . exact partition
    . constructor
      . exact positive
      . intros i naturalI before
        exfalso
        have typedI := nat.omega.backward with naturalI
        have impossible := Nat.not_less_zero with typedI
        exact impossible with before
  witness constantSequence(A)
  witness constantSequence(A)
  witness 0
  constructor
  . exact nat.omega.forward at 0 with Naturals.constructor.zero
  . refine ⟨tagged, fine⟩

-- IntegralUniqueness
theorem HasRiemannIntegral.unique (f : Set) (A B I J : Real)
  requires cofinal : FineTaggedPartitionsExist(A, B)
  requires first : HasRiemannIntegral(f, A, B, I)
  requires second : HasRiemannIntegral(f, A, B, J)
  shows I = J
proof
  by_contra unequal
  have distanceReal := dist.real at I, J with *
  have radiusReal := half.real with distanceReal
  have distancePositive := dist.positive_of_ne with unequal
  have radiusPositive := half.positive with distancePositive
  rcases first with ⟨function, ordered, controlI⟩
  rcases second with ⟨functionJ, orderedJ, controlJ⟩
  obtain ⟨D, realD, positiveD, estimateI⟩ :=
    controlI with half(dist(I, J)), radiusPositive
  obtain ⟨E, realE, positiveE, estimateJ⟩ :=
    controlJ with half(dist(I, J)), radiusPositive
  obtain ⟨p, t, N, naturalN, tagged, fineD, fineE⟩ := FineTaggedPartitionsExist.common with
    cofinal, positiveD, positiveE
  have nearI := estimateI with p, t, N, naturalN, tagged, fineD
  have nearJ := estimateJ with p, t, N, naturalN, tagged, fineE
  have typedN := nat.omega.backward with naturalN
  have points : RealSequence(p)
    rcases tagged with ⟨partition, tags, cells⟩
    rcases partition with ⟨sequence, left, right, increasing⟩
    exact sequence
  have tags : RealSequence(t)
    rcases tagged with ⟨partition, sequence, cells⟩
    exact sequence
  have sumReal := riemannSum.real with function, points, tags
  have sumMember : riemannSum(f, p, t, N) ∈ ℝ by exact sumReal
  have ballI : riemannSum(f, p, t, N) ∈ openBall(I, half(dist(I, J)))
    apply openBall.membership.backward
    refine ⟨sumMember, nearI⟩
  have ballJ : riemannSum(f, p, t, N) ∈ openBall(J, half(dist(I, J)))
    apply openBall.membership.backward
    refine ⟨sumMember, nearJ⟩
  have equal := openBall.half_radius_unique with ballI, ballJ
  exact unequal with equal

theorem HasRiemannIntegral.unique_point (f : Set) (A I J : Real)
  requires first : HasRiemannIntegral(f, A, A, I)
  requires second : HasRiemannIntegral(f, A, A, J)
  shows I = J
proof
  exact HasRiemannIntegral.unique at f, A, A, I, J with
    FineTaggedPartitionsExist.point at A with assumption, first, second

-- PointIntegral
theorem HasRiemannIntegral.point (f : Set) (A : Real)
  requires function : f : ℝ → ℝ
  shows HasRiemannIntegral(f, A, A, (0:R))
proof
  constructor
  . exact function
  . constructor
    . exact realLe.reflexive at A
    . intros E realE positiveE
      witness E
      constructor
      . exact realE
      . constructor
        . exact positiveE
        . intros p t N naturalN tagged fine
          have typedN := nat.omega.backward with naturalN
          rcases tagged with ⟨partition, tags, cells⟩
          have empty := IntervalPartition.point_empty with partition
          rw [empty, riemannSum.zero]
          have same := dist.self with realZero.real
          rwa same

theorem RiemannIntegrable.point (f : Set) (A : Real)
  requires function : f : ℝ → ℝ
  shows RiemannIntegrable(f, A, A)
proof
  witness realZero
  constructor
  . exact realZero.real
  . exact HasRiemannIntegral.point at f, A with function

-- PositiveIntegralUniqueness
theorem FineTaggedPartitionsExist.ordered_interval (A B : Real)
  requires ordered : A ≤ B
  shows FineTaggedPartitionsExist(A, B)
proof
  have alternatives : A = B ∨ ¬(A = B) by tauto
  cases alternatives with equal | unequal
  . rw ← equal
    exact FineTaggedPartitionsExist.point with assumption
  . have le := realLe.characterization.forward with ordered
    rcases le with ⟨realA, realB, subset⟩
    have proper : A ⊂ B by refine ⟨subset, unequal⟩
    have strict : A < B := realLess.characterization.backward with ⟨realA, realB, proper⟩
    exact FineTaggedPartitionsExist.positive_interval with strict

theorem HasRiemannIntegral.unique_positive (f : Set) (A B I J : Real)
  requires ordered : A < B
  requires first : HasRiemannIntegral(f, A, B, I)
  requires second : HasRiemannIntegral(f, A, B, J)
  shows I = J
proof
  have cofinal := FineTaggedPartitionsExist.positive_interval with ordered
  exact HasRiemannIntegral.unique with cofinal, first, second

theorem HasRiemannIntegral.unique_on_interval (f : Set) (A B I J : Real)
  requires ordered : A ≤ B
  requires first : HasRiemannIntegral(f, A, B, I)
  requires second : HasRiemannIntegral(f, A, B, J)
  shows I = J
proof
  have cofinal := FineTaggedPartitionsExist.ordered_interval with ordered
  exact HasRiemannIntegral.unique with cofinal, first, second

-- IntegralOrderBound
theorem HasRiemannIntegral.upper_bound (f : Set) (A B I U : Real)
  requires cofinal : FineTaggedPartitionsExist(A, B)
  requires integral : HasRiemannIntegral(f, A, B, I)
  requires bounds : ∀ p, ∀ t, ∀ N ∈ ℕ,
    TaggedPartition(p, t, A, B, N) → (riemannSum(f, p, t, N):R) ≤ U
  shows I ≤ U
proof
  by_contra notLe
  have strict := realLess.of_not_le with notLe
  have positiveDifference : 0 < I - U by linarith
  rcases integral with ⟨function, ordered, control⟩
  obtain ⟨D, realD, positiveD, estimate⟩ :=
    control with I - U, positiveDifference
  obtain ⟨p, t, N, naturalN, tagged, fine⟩ := cofinal with D, realD, positiveD
  have typedN := nat.omega.backward with naturalN
  have close := estimate with p, t, N, naturalN, tagged, fine
  have upper := bounds with p, t, N, naturalN, tagged
  have points : RealSequence(p)
    rcases tagged with ⟨partition, tags, cells⟩
    rcases partition with ⟨points, start, finish, increasing⟩
    exact points
  have tags : RealSequence(t)
    rcases tagged with ⟨partition, tags, cells⟩
    exact tags
  have sumReal := riemannSum.real with function, points, tags
  have data : riemannSum(f, p, t, N) ∈ ℝ ∧ dist(riemannSum(f, p, t, N), I) < I - U by refine ⟨sumReal, close⟩
  have ball := openBall.membership.backward at I, I - U, riemannSum(f, p, t, N)
    with data
  obtain ⟨lower, upperBall⟩ := openBall.order_characterization.forward with ball
  have cancel : I - (I - U) = U by ring
  rw cancel at lower
  have impossible : U < riemannSum(f, p, t, N) by exact lower
  have cycle := realLe.transitive_less with upper, impossible
  have contradiction := realLess.characterization.forward with cycle
  exact RealLess.irreflexive with contradiction

theorem HasRiemannIntegral.lower_bound (f : Set) (A B I L : Real)
  requires cofinal : FineTaggedPartitionsExist(A, B)
  requires integral : HasRiemannIntegral(f, A, B, I)
  requires bounds : ∀ p, ∀ t, ∀ N ∈ ℕ,
    TaggedPartition(p, t, A, B, N) → L ≤ (riemannSum(f, p, t, N):R)
  shows L ≤ I
proof
  by_contra notLe
  have strict := realLess.of_not_le with notLe
  have positiveDifference : 0 < L - I by linarith
  rcases integral with ⟨function, ordered, control⟩
  obtain ⟨D, realD, positiveD, estimate⟩ :=
    control with L - I, positiveDifference
  obtain ⟨p, t, N, naturalN, tagged, fine⟩ := cofinal with D, realD, positiveD
  have typedN := nat.omega.backward with naturalN
  have close := estimate with p, t, N, naturalN, tagged, fine
  have lower := bounds with p, t, N, naturalN, tagged
  have points : RealSequence(p)
    rcases tagged with ⟨partition, tags, cells⟩
    rcases partition with ⟨points, start, finish, increasing⟩
    exact points
  have tags : RealSequence(t)
    rcases tagged with ⟨partition, tags, cells⟩
    exact tags
  have sumReal := riemannSum.real with function, points, tags
  have data : riemannSum(f, p, t, N) ∈ ℝ ∧ dist(riemannSum(f, p, t, N), I) < L - I by refine ⟨sumReal, close⟩
  have ball := openBall.membership.backward at I, L - I, riemannSum(f, p, t, N)
    with data
  obtain ⟨lowerBall, upperBall⟩ := openBall.order_characterization.forward with ball
  have cancel : I + (L - I) = L by ring
  rw cancel at upperBall
  have cycle := realLe.transitive_less with lower, upperBall
  have contradiction := realLess.characterization.forward with cycle
  exact RealLess.irreflexive with contradiction

theorem HasRiemannIntegral.between_constants (f : Set) (A B I L U : Real)
  requires ordered : A < B
  requires integral : HasRiemannIntegral(f, A, B, I)
  requires lower : ∀ X ∈ closedInterval(A, B), L ≤ (f(X):R)
  requires upper : ∀ X ∈ closedInterval(A, B), (f(X):R) ≤ U
  shows L * (B - A) ≤ I ∧ I ≤ U * (B - A)
proof
  have function : f : ℝ → ℝ
    rcases integral with ⟨function, ordered, estimates⟩
    exact function
  have cofinal := FineTaggedPartitionsExist.positive_interval with ordered
  have lowerSums : ∀ p, ∀ t, ∀ N ∈ ℕ,
      TaggedPartition(p, t, A, B, N) → L * (B - A) ≤ (riemannSum(f, p, t, N):R)
    intros p t N naturalN tagged
    have typedN := nat.omega.backward with naturalN
    have partition : IntervalPartition(p, A, B, N)
      rcases tagged with ⟨partition, tags, cells⟩
      exact partition
    have lowerTags : ∀ i ∈ ℕ, i < N → L ≤ (f(t[i]):R)
      intros i naturalI before
      have typedI := nat.omega.backward with naturalI
      have cell := TaggedPartition.cell with tagged, before
      have subset := IntervalPartition.cell_subset with partition, before
      have global := subset with t[i], cell
      exact lower with t[i], global
    have upperTags : ∀ i ∈ ℕ, i < N → (f(t[i]):R) ≤ U
      intros i naturalI before
      have typedI := nat.omega.backward with naturalI
      have cell := TaggedPartition.cell with tagged, before
      have subset := IntervalPartition.cell_subset with partition, before
      have global := subset with t[i], cell
      exact upper with t[i], global
    obtain ⟨low, high⟩ := riemannSum.between_constants with function, tagged, lowerTags, upperTags
    exact low
  have upperSums : ∀ p, ∀ t, ∀ N ∈ ℕ,
      TaggedPartition(p, t, A, B, N) → (riemannSum(f, p, t, N):R) ≤ U * (B - A)
    intros p t N naturalN tagged
    have typedN := nat.omega.backward with naturalN
    have partition : IntervalPartition(p, A, B, N)
      rcases tagged with ⟨partition, tags, cells⟩
      exact partition
    have lowerTags : ∀ i ∈ ℕ, i < N → L ≤ (f(t[i]):R)
      intros i naturalI before
      have typedI := nat.omega.backward with naturalI
      have cell := TaggedPartition.cell with tagged, before
      have subset := IntervalPartition.cell_subset with partition, before
      have global := subset with t[i], cell
      exact lower with t[i], global
    have upperTags : ∀ i ∈ ℕ, i < N → (f(t[i]):R) ≤ U
      intros i naturalI before
      have typedI := nat.omega.backward with naturalI
      have cell := TaggedPartition.cell with tagged, before
      have subset := IntervalPartition.cell_subset with partition, before
      have global := subset with t[i], cell
      exact upper with t[i], global
    obtain ⟨low, high⟩ := riemannSum.between_constants with function, tagged, lowerTags, upperTags
    exact high
  have low := HasRiemannIntegral.lower_bound with cofinal, integral, lowerSums
  have high := HasRiemannIntegral.upper_bound with cofinal, integral, upperSums
  refine ⟨low, high⟩

-- Section: LocalIntegralBounds
theorem realLe.divide_positive (A B C : Real)
  requires bound : A ≤ B
  requires positive : 0 < C
  shows A / C ≤ B / C
proof
  have inverseReal := realInvValue.real at C
  have inversePositive := realInvValue.positive with positive
  have nonnegative := realLess.implies_le with realZero.real, inversePositive
  have scaled := realLe.mul_right_of_nonnegative at A, B, realInvValue(C) with
    *{2}, inverseReal, bound, nonnegative
  have leftShape := realDiv.eq_mul_inv at A, C with *
  have rightShape := realDiv.eq_mul_inv at B, C with *
  rwa [← leftShape, ← rightShape] at scaled

theorem realDiv.neg_both (A B : Real)
  requires nonzero : B ≠ realZero
  shows (-A) / (-B) = A / B
proof
  have negAReal := realNegValue.real at A
  have negBReal := realNegValue.real at B
  have negNonzero : -B ≠ realZero
    intro equal
    have involutive := realNegValue.involutive at B
    rw [equal, realNegValue.zero] at involutive
    apply nonzero
    cc
  have inverseLaw := realInvValue.mul_right with nonzero
  have negLaw : (-B) * (-realInvValue(B)) = 1
    have shape : (-B) * (-realInvValue(B)) = B * realInvValue(B) by ring
    rwa shape
  have negInverse := realInvValue.eq_of_mul_right with negNonzero, negLaw
  have leftShape := realDiv.eq_mul_inv with negAReal, negBReal
  have rightShape := realDiv.eq_mul_inv at A, B with *
  rw [leftShape, rightShape, ← negInverse]
  ring

-- LocalIntegralBounds
theorem ContinuousOn.local_image_bounds (f : Set) (A B X E : Real)
  requires continuous : ContinuousOn(f, closedInterval(A, B))
  requires inside : X ∈ closedInterval(A, B)
  requires positive : 0 < E
  shows ∃ D ∈ ℝ, 0 < D ∧
    (∀ Z ∈ closedInterval(A, B), dist(X, Z) < D →
      f(X) - E < (f(Z):R) ∧
        (f(Z):R) < f(X) + E)
proof
  rcases continuous with ⟨function, subset, atPoints⟩
  obtain ⟨functionX, subsetX, control⟩ := atPoints with X, inside
  obtain ⟨D, realD, positiveD, maps⟩ := control with E, positive
  witness D
  constructor
  . exact realD
  . constructor
    . exact positiveD
    . intros Z member near
      have realZ : Z ∈ ℝ
        obtain ⟨realZ, lower, upper⟩ := closedInterval.membership.forward with member
        exact realZ
      have realX : X ∈ ℝ
        obtain ⟨realX, lower, upper⟩ := closedInterval.membership.forward with inside
        exact realX
      have symmetry := dist.symmetric at X, Z with *
      rw symmetry at near
      have data : Z ∈ ℝ ∧ dist(Z, X) < D by refine ⟨realZ, near⟩
      have inputBall := openBall.membership.backward with data
      have outputBall := maps with Z, member, inputBall
      exact openBall.order_characterization.forward with outputBall

theorem ContinuousOn.right_integral_bounds (f : Set) (A B X E D Y I : Real)
  requires continuous : ContinuousOn(f, closedInterval(A, B))
  requires inside : X ∈ closedInterval(A, B)
  requires control : ∀ Z ∈ closedInterval(A, B), dist(X, Z) < D →
    f(X) - E < (f(Z):R) ∧
      (f(Z):R) < f(X) + E
  requires ordered : X < Y
  requires right : Y ≤ B
  requires small : Y - X < D
  requires integral : HasRiemannIntegral(f, X, Y, I)
  shows (f(X) - E) * (Y - X) ≤ I ∧
    I ≤ (f(X) + E) * (Y - X)
proof
  obtain ⟨realX, leftBound, rightBound⟩ := closedInterval.membership.forward with inside
  have function : f : ℝ → ℝ
    rcases continuous with ⟨function, subset, points⟩
    exact function
  have valueX := FunctionBetween.apply_closed with function, realX
  have weak := realLess.implies_le with ordered
  have subset := closedInterval.nested with leftBound, right
  have lowerValues : ∀ Z ∈ closedInterval(X, Y),
      f(X) - E ≤ (f(Z):R)
    intros Z local
    have global := subset with Z, local
    obtain ⟨realZ, xLeZ, zLeY⟩ := closedInterval.membership.forward with local
    have distance := dist.le_interval_width with xLeZ, zLeY
    have near := realLe.transitive_less with distance, small
    obtain ⟨lowValue, highValue⟩ := control with Z, global, near
    have valueZ := FunctionBetween.apply_closed with function, realZ
    exact realLess.implies_le with lowValue
  have upperValues : ∀ W ∈ closedInterval(X, Y),
      (f(W):R) ≤ f(X) + E
    intros W local
    have global := subset with W, local
    obtain ⟨realW, xLeW, wLeY⟩ := closedInterval.membership.forward with local
    have distance := dist.le_interval_width with xLeW, wLeY
    have near := realLe.transitive_less with distance, small
    obtain ⟨lowValue, highValue⟩ := control with W, global, near
    have valueW := FunctionBetween.apply_closed with function, realW
    exact realLess.implies_le with highValue
  exact HasRiemannIntegral.between_constants with ordered, integral, lowerValues, upperValues

theorem ContinuousOn.right_integral_average (f : Set) (A B X E D Y I : Real)
  requires continuous : ContinuousOn(f, closedInterval(A, B))
  requires inside : X ∈ closedInterval(A, B)
  requires control : ∀ Z ∈ closedInterval(A, B), dist(X, Z) < D →
    f(X) - E < (f(Z):R) ∧
      (f(Z):R) < f(X) + E
  requires ordered : X < Y
  requires right : Y ≤ B
  requires small : Y - X < D
  requires integral : HasRiemannIntegral(f, X, Y, I)
  shows f(X) - E ≤ I / (Y - X) ∧
    I / (Y - X) ≤ f(X) + E
proof
  have bounds := ContinuousOn.right_integral_bounds with
    continuous, inside, control, ordered, right, small, integral
  rcases bounds with ⟨lower, upper⟩
  have positive : 0 < Y - X by linarith
  have width := realSubValue.real at Y, X with *
  have nonzero : Y - X ≠ realZero
    intro equal
    have strict := realLess.characterization.forward with positive
    rw equal at strict
    exact RealLess.irreflexive with strict
  have function : f : ℝ → ℝ
    rcases continuous with ⟨function, subset, points⟩
    exact function
  have realX : X ∈ ℝ by simp [Real]
  have valueX := FunctionBetween.apply_closed with function, realX
  have typedFX : Real(f(X)) by exact valueX
  have lowReal := realSubValue.real at f(X), E with typedFX, assumption
  have highReal := realAddValue.real at f(X), E with typedFX, assumption
  have lowProduct := realMulValue.real with lowReal, width
  have highProduct := realMulValue.real with highReal, width
  have lowDiv := realLe.divide_positive at (f(X) - E) * (Y - X), I,
    Y - X with lowProduct, assumption, width, lower, positive
  have highDiv := realLe.divide_positive at I, (f(X) + E) * (Y - X),
    Y - X with assumption, highProduct, width, upper, positive
  have lowShape := realDiv.mul_numerator_left with lowReal, width, width
  have highShape := realDiv.mul_numerator_left with highReal, width, width
  have cancel := realDiv.self with nonzero
  rw lowShape at lowDiv
  rw highShape at highDiv
  rw cancel at lowDiv
  rw cancel at highDiv
  have oneLow := realMulValue.one_right with lowReal
  have oneHigh := realMulValue.one_right with highReal
  rw oneLow at lowDiv
  rw oneHigh at highDiv
  refine ⟨lowDiv, highDiv⟩

theorem ContinuousOn.left_integral_bounds (f : Set) (A B X E D Y I : Real)
  requires continuous : ContinuousOn(f, closedInterval(A, B))
  requires inside : X ∈ closedInterval(A, B)
  requires control : ∀ Z ∈ closedInterval(A, B), dist(X, Z) < D →
    f(X) - E < (f(Z):R) ∧
      (f(Z):R) < f(X) + E
  requires ordered : Y < X
  requires left : A ≤ Y
  requires small : X - Y < D
  requires integral : HasRiemannIntegral(f, Y, X, I)
  shows (f(X) - E) * (X - Y) ≤ I ∧
    I ≤ (f(X) + E) * (X - Y)
proof
  obtain ⟨realX, leftBound, rightBound⟩ := closedInterval.membership.forward with inside
  have function : f : ℝ → ℝ
    rcases continuous with ⟨function, subset, points⟩
    exact function
  have valueX := FunctionBetween.apply_closed with function, realX
  have subset := closedInterval.nested with left, rightBound
  have lowerValues : ∀ Z ∈ closedInterval(Y, X),
      f(X) - E ≤ (f(Z):R)
    intros Z local
    have global := subset with Z, local
    obtain ⟨realZ, yLeZ, zLeX⟩ := closedInterval.membership.forward with local
    have distance := dist.le_interval_width at Z, X, X with
      zLeX, realLe.reflexive at X with assumption
    have difference : X - Z ≤ X - Y by linarith
    have distanceBound := realLe.transitive with distance, difference
    have symmetry := dist.symmetric at X, Z with *
    rw ← symmetry at distanceBound
    have near := realLe.transitive_less with distanceBound, small
    obtain ⟨lowValue, highValue⟩ := control with Z, global, near
    have valueZ := FunctionBetween.apply_closed with function, realZ
    exact realLess.implies_le with lowValue
  have upperValues : ∀ W ∈ closedInterval(Y, X),
      (f(W):R) ≤ f(X) + E
    intros W local
    have global := subset with W, local
    obtain ⟨realW, yLeW, wLeX⟩ := closedInterval.membership.forward with local
    have distance := dist.le_interval_width at W, X, X with
      wLeX, realLe.reflexive at X with assumption
    have difference : X - W ≤ X - Y by linarith
    have distanceBound := realLe.transitive with distance, difference
    have symmetry := dist.symmetric at X, W with *
    rw ← symmetry at distanceBound
    have near := realLe.transitive_less with distanceBound, small
    obtain ⟨lowValue, highValue⟩ := control with W, global, near
    have valueW := FunctionBetween.apply_closed with function, realW
    exact realLess.implies_le with highValue
  exact HasRiemannIntegral.between_constants with ordered, integral, lowerValues, upperValues

theorem ContinuousOn.left_integral_average (f : Set) (A B X E D Y I : Real)
  requires continuous : ContinuousOn(f, closedInterval(A, B))
  requires inside : X ∈ closedInterval(A, B)
  requires control : ∀ Z ∈ closedInterval(A, B), dist(X, Z) < D →
    f(X) - E < (f(Z):R) ∧
      (f(Z):R) < f(X) + E
  requires ordered : Y < X
  requires left : A ≤ Y
  requires small : X - Y < D
  requires integral : HasRiemannIntegral(f, Y, X, I)
  shows f(X) - E ≤ I / (X - Y) ∧
    I / (X - Y) ≤ f(X) + E
proof
  have bounds := ContinuousOn.left_integral_bounds with
    continuous, inside, control, ordered, left, small, integral
  rcases bounds with ⟨lower, upper⟩
  have positive : 0 < X - Y by linarith
  have width := realSubValue.real at X, Y with *
  have nonzero : X - Y ≠ realZero
    intro equal
    have strict := realLess.characterization.forward with positive
    rw equal at strict
    exact RealLess.irreflexive with strict
  have function : f : ℝ → ℝ
    rcases continuous with ⟨function, subset, points⟩
    exact function
  have realX : X ∈ ℝ by simp [Real]
  have valueX := FunctionBetween.apply_closed with function, realX
  have typedFX : Real(f(X)) by exact valueX
  have lowReal := realSubValue.real at f(X), E with typedFX, assumption
  have highReal := realAddValue.real at f(X), E with typedFX, assumption
  have lowProduct := realMulValue.real with lowReal, width
  have highProduct := realMulValue.real with highReal, width
  have lowDiv := realLe.divide_positive at (f(X) - E) * (X - Y), I,
    X - Y with lowProduct, assumption, width, lower, positive
  have highDiv := realLe.divide_positive at I, (f(X) + E) * (X - Y),
    X - Y with assumption, highProduct, width, upper, positive
  have lowShape := realDiv.mul_numerator_left with lowReal, width, width
  have highShape := realDiv.mul_numerator_left with highReal, width, width
  have cancel := realDiv.self with nonzero
  rw lowShape at lowDiv
  rw highShape at highDiv
  rw cancel at lowDiv
  rw cancel at highDiv
  have oneLow := realMulValue.one_right with lowReal
  have oneHigh := realMulValue.one_right with highReal
  rw oneLow at lowDiv
  rw oneHigh at highDiv
  refine ⟨lowDiv, highDiv⟩

-- Section: CellMeanValueChoice
defgraph functionIncrements (f p : Set) (n ∈ ℕ) : ℝ :=
  realSubValue(f(p[succ(n)]), f(p[n]))
requires function : f : ℝ → ℝ
requires points : RealSequence(p)
proof
  have pointN := RealSequence.value with points
  have pointNext := RealSequence.value at p, succ(n) with points
  unfold Real at pointN
  unfold Real at pointNext
  have valueN := FunctionBetween.apply_closed with function, pointN
  have valueNext := FunctionBetween.apply_closed with function, pointNext
  have realN : Real(f(p[n])) by exact valueN
  have realNext : Real(f(p[succ(n)])) by exact valueNext
  exact realSubValue.real with realNext, realN

theorem functionIncrements.sequence (f p : Set)
  requires function : f : ℝ → ℝ
  requires points : RealSequence(p)
  shows RealSequence(functionIncrements(f, p))
proof
  unfold RealSequence
  exact functionIncrements.function with function, points

theorem functionIncrements.telescoping (f p : Set) (N : Nat)
  requires function : f : ℝ → ℝ
  requires points : RealSequence(p)
  shows realPartialSum(functionIncrements(f, p), N) =
    f(p[N]) - f(p[0])
proof
  induction N
  | zero =>
    rw realPartialSum.zero
    have point0 := RealSequence.value at p, 0 with points
    unfold Real at point0
    have value0 := FunctionBetween.apply_closed with function, point0
    ring
  | succ(k) using ih =>
    rw [realPartialSum.succ, ih]
    have value := functionIncrements.value with function, points
    rw value
    have point0 := RealSequence.value with points
    have pointK := RealSequence.value with points
    have pointNext := RealSequence.value with points
    unfold Real at point0
    unfold Real at pointK
    unfold Real at pointNext
    have value0 := FunctionBetween.apply_closed with function, point0
    have valueK := FunctionBetween.apply_closed with function, pointK
    have valueNext := FunctionBetween.apply_closed with function, pointNext
    ring

theorem IntervalPartition.sum_function_increments (f p : Set) (A B : Real, N : Nat)
  requires function : f : ℝ → ℝ
  requires partition : IntervalPartition(p, A, B, N)
  shows realPartialSum(functionIncrements(f, p), N) =
    f(B) - f(A)
proof
  rcases partition with ⟨points, left, right, increasing⟩
  have sum := functionIncrements.telescoping with function, points
  rwa [left, right] at sum

-- This is the finite algebraic part of FTC. A mean value theorem will later
-- produce tags satisfying the displayed cell equation.
theorem TaggedPartition.riemann_sum_of_cell_increments (f F p t : Set) (A B : Real, N : Nat)
  requires functionF : F : ℝ → ℝ
  requires functionf : f : ℝ → ℝ
  requires tagged : TaggedPartition(p, t, A, B, N)
  requires cells : ∀ i ∈ ℕ, i < N →
    realMulValue(f(t[i]), realSubValue(p[succ(i)], p[i])) =
      realSubValue(F(p[succ(i)]), F(p[i]))
  shows riemannSum(f, p, t, N) = F(B) - F(A)
proof
  rcases tagged with ⟨partition, tags, tagCells⟩
  have points : RealSequence(p)
    rcases partition with ⟨sequence, left, right, increasing⟩
    exact sequence
  have fTerms : RealSequence(riemannTerms(f, p, t))
    unfold RealSequence
    exact riemannTerms.function with functionf, points, tags
  have increments := functionIncrements.sequence with functionF, points
  have equalTerms : ∀ i ∈ ℕ, i < N →
      riemannTerms(f, p, t)[i] = functionIncrements(F, p)[i]
    intros i naturalI before
    have riemannValue := riemannTerms.value at f, p, t, i with
      functionf, points, tags
    have incrementValue := functionIncrements.value at F, p, i with
      functionF, points
    rw [riemannValue, incrementValue]
    exact cells with i, naturalI, before
  have sameSum := realPartialSum.eq_prefix with fTerms, increments, equalTerms
  have endpoint := IntervalPartition.sum_function_increments with functionF, partition
  simp [riemannSum]

-- CellMeanValueBridge
predicate CellMeanValueTagsExist (f F : Set) (A B : Real) :=
  ∀ p, ∀ N ∈ ℕ, IntervalPartition(p, A, B, N) →
    ∃ t, TaggedPartition(p, t, A, B, N) ∧
      (∀ i ∈ ℕ, i < N →
        realMulValue(f(t[i]), realSubValue(p[succ(i)], p[i])) =
          realSubValue(F(p[succ(i)]), F(p[i])))

-- For every epsilon, take a fine partition, then its mean-value
-- tags. Its Riemann sum telescopes exactly, so the proposed integral cannot
-- differ from the endpoint increment.
theorem HasRiemannIntegral.eq_endpoint_of_cell_tags (f F : Set) (A B I : Real)
  requires functionF : F : ℝ → ℝ
  requires integral : HasRiemannIntegral(f, A, B, I)
  requires meanTags : CellMeanValueTagsExist(f, F, A, B)
  shows I = F(B) - F(A)
proof
  by_contra unequal
  have functionf : f : ℝ → ℝ
    rcases integral with ⟨function, ordered, control⟩
    exact function
  have realA : A ∈ ℝ by simp [Real]
  have realB : B ∈ ℝ by simp [Real]
  have valueA := FunctionBetween.apply_closed with functionF, realA
  have valueB := FunctionBetween.apply_closed with functionF, realB
  have distinct : F(B) - F(A) ≠ I
    intro same
    apply unequal
    rw same
    rfl
  have positive := dist.positive_of_ne with distinct
  rcases integral with ⟨functionAgain, ordered, control⟩
  obtain ⟨D, realD, positiveD, estimate⟩ :=
    control with dist(F(B) - F(A), I), positive
  have cofinal := FineTaggedPartitionsExist.ordered_interval with ordered
  obtain ⟨p, preliminaryTags, N, naturalN, preliminary, fine⟩ :=
    cofinal with D, realD, positiveD
  have partition : IntervalPartition(p, A, B, N)
    rcases preliminary with ⟨points, tags, cells⟩
    exact points
  obtain ⟨t, tagged, cells⟩ := meanTags with p, N, naturalN, partition
  have typedN := nat.omega.backward with naturalN
  have sum := TaggedPartition.riemann_sum_of_cell_increments with functionF, functionf, tagged, cells
  have near := estimate with p, t, N, naturalN, tagged, fine
  rw sum at near
  have contradiction := realLess.characterization.forward with near
  exact RealLess.irreflexive with contradiction

-- CellMeanValueChoice
predicate CellMeanValueWitnesses (f F p : Set) (N : Nat) :=
  ∀ i ∈ ℕ, i < N →
    ∃ X ∈ closedInterval(p[i], p[succ(i)]),
      realMulValue(f(X), realSubValue(p[succ(i)], p[i])) =
        realSubValue(F(p[succ(i)]), F(p[i]))

-- Outside the finite prefix a tag is the left endpoint, which makes the
-- sequence total without a condition on unused cells.
defchoose cellMeanValueTag (f F p : Set, N i : Nat) (X ∈ ℝ) :=
  (i < N → X ∈ closedInterval(p[i], p[succ(i)]) ∧
    realMulValue(f(X), realSubValue(p[succ(i)], p[i])) =
      realSubValue(F(p[succ(i)]), F(p[i]))) ∧
  (¬(i < N) → X = p[i])
requires points : RealSequence(p)
requires witnesses : CellMeanValueWitnesses(f, F, p, N)
proof
  have alternatives : i < N ∨ ¬(i < N)
    tauto
  cases alternatives with before | beyond
  . unfold CellMeanValueWitnesses at witnesses
    have naturalI := nat.omega.forward at i with assumption
    obtain ⟨X, inCell, equation⟩ := witnesses with i, naturalI, before
    have pointI := RealSequence.value at p, i with points
    have pointNext := RealSequence.value at p, succ(i) with points
    have data := closedInterval.membership.forward at p[i], p[succ(i)], X
      with inCell
    rcases data with ⟨realX, lower, upper⟩
    witness X
    constructor
    . exact realX
    . constructor
      . intro beforeAgain
        refine ⟨inCell, equation⟩
      . intro notBefore
        exfalso
        exact notBefore with before
  . have pointI := RealSequence.value at p, i with points
    witness p[i]
    constructor
    . unfold Real at pointI
      exact pointI
    . constructor
      . intro before
        exfalso
        exact beyond with before
      . intro notBefore
        rfl

defgraph cellMeanValueTags (f F p : Set, N : Nat) (i ∈ ℕ) : ℝ :=
  cellMeanValueTag(f, F, p, N, i)
requires points : RealSequence(p)
requires witnesses : CellMeanValueWitnesses(f, F, p, N)
proof
  have selected := cellMeanValueTag.selected at f, F, p, N, i with
    assumption, points, witnesses
  rcases selected with ⟨realTag, properties⟩
  exact realTag

theorem CellMeanValueWitnesses.tags (f F p : Set) (A B : Real, N : Nat)
  requires partition : IntervalPartition(p, A, B, N)
  requires witnesses : CellMeanValueWitnesses(f, F, p, N)
  shows TaggedPartition(p, cellMeanValueTags(f, F, p, N), A, B, N) ∧
    (∀ i ∈ ℕ, i < N →
      realMulValue(f(cellMeanValueTags(f, F, p, N)[i]),
        realSubValue(p[succ(i)], p[i])) =
          realSubValue(F(p[succ(i)]), F(p[i])))
proof
  have points : RealSequence(p)
    rcases partition with ⟨sequence, left, right, increasing⟩
    exact sequence
  have sequence : RealSequence(cellMeanValueTags(f, F, p, N))
    unfold RealSequence
    exact cellMeanValueTags.function with points, witnesses
  have cell : ∀ i ∈ ℕ, i < N → cellMeanValueTags(f, F, p, N)[i] ∈ closedInterval(p[i], p[succ(i)]) ∧ realMulValue(f(cellMeanValueTags(f, F, p, N)[i]), realSubValue(p[succ(i)], p[i])) = realSubValue(F(p[succ(i)]), F(p[i]))
    intros i naturalI before
    have typedI := nat.omega.backward with naturalI
    have value := cellMeanValueTags.value at f, F, p, N, i with
      points, witnesses
    have selected := cellMeanValueTag.selected at f, F, p, N, i with
      points, witnesses
    rcases selected with ⟨realTag, valid, unused⟩
    have selectedCell := valid with before
    rwa value
  constructor
    constructor
    . exact partition
    . constructor
      . exact sequence
      . intros j naturalJ before
        obtain ⟨member, equation⟩ := cell with j, naturalJ, before
        exact member
  . intros j naturalJ before
    obtain ⟨member, equation⟩ := cell with j, naturalJ, before
    exact equation

theorem CellMeanValueWitnesses.all_tags (f F : Set) (A B : Real)
  requires local : ∀ p, ∀ N ∈ ℕ,
    IntervalPartition(p, A, B, N) → CellMeanValueWitnesses(f, F, p, N)
  shows CellMeanValueTagsExist(f, F, A, B)
proof
  intros p N naturalN partition
  have witnesses := local with p, N, naturalN, partition
  have typedN := nat.omega.backward with naturalN
  obtain ⟨tagged, equations⟩ := CellMeanValueWitnesses.tags with partition, witnesses
  witness cellMeanValueTags(f, F, p, N)
  refine ⟨tagged, equations⟩

-- This is the interval form of the mean value theorem needed by FTC. The
-- choice construction above turns its individual cell witnesses into a
-- single tagged partition without an additional choice assumption.
predicate CellMeanValueProperty (f F : Set) :=
  ∀ X ∈ ℝ, ∀ Y ∈ ℝ, (X:R) < (Y:R) →
    ∃ T ∈ closedInterval(X, Y),
      realMulValue(f(T), realSubValue(Y, X)) =
        realSubValue(F(Y), F(X))

theorem CellMeanValueProperty.partition_witnesses (f F p : Set) (A B : Real, N : Nat)
  requires partition : IntervalPartition(p, A, B, N)
  requires meanValue : CellMeanValueProperty(f, F)
  shows CellMeanValueWitnesses(f, F, p, N)
proof
  intros i naturalI before
  have typedI := nat.omega.backward with naturalI
  have points : RealSequence(p)
    rcases partition with ⟨sequence, left, right, increasing⟩
    exact sequence
  have pointI := RealSequence.value at p, i with points
  have pointNext := RealSequence.value at p, succ(i) with points
  have increasing := IntervalPartition.cell with partition, before
  exact meanValue with p[i], pointI, p[succ(i)], pointNext, increasing

theorem HasRiemannIntegral.eq_endpoint_of_mean_value (f F : Set) (A B I : Real)
  requires functionF : F : ℝ → ℝ
  requires integral : HasRiemannIntegral(f, A, B, I)
  requires meanValue : CellMeanValueProperty(f, F)
  shows I = F(B) - F(A)
proof
  have local : ∀ p, ∀ N ∈ ℕ,
      IntervalPartition(p, A, B, N) → CellMeanValueWitnesses(f, F, p, N)
    intros p N naturalN partition
    have typedN := nat.omega.backward with naturalN
    exact CellMeanValueProperty.partition_witnesses with partition, meanValue
  have tags := CellMeanValueWitnesses.all_tags with local
  exact HasRiemannIntegral.eq_endpoint_of_cell_tags with functionF, integral, tags