-- Real analysis: Differentiation
-- 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.Differentiation
import Continuity

-- Section: Differentiation
def constantRealFunction (L : Set) : Set :=
  constantFunction(ℝ, ℝ, L)

theorem constantRealFunction.function (L : Real)
  shows constantRealFunction(L) : ℝ → ℝ
proof
  have realL : L ∈ ℝ by simp [Real]
  unfold constantRealFunction
  exact constantFunction.function with realL

theorem constantRealFunction.value (L X : Real)
  shows constantRealFunction(L)(X) = L
proof
  have realL : L ∈ ℝ by simp [Real]
  have realX : X ∈ ℝ by simp [Real]
  unfold constantRealFunction
  exact constantFunction.apply with realX, realL

theorem constantRealFunction.continuous_at (L A : Real)
  shows ContinuousAt(constantRealFunction(L), A)
proof
  have function := constantRealFunction.function at L
  constructor
  . exact function
  . intros E realE positive
    witness E
    constructor
    . exact realE
    . constructor
      . exact positive
      . intros X realX near
        have typedX : Real(X) by exact realX
        have valueX := constantRealFunction.value at L, X with
          assumption, typedX
        have valueA := constantRealFunction.value at L, A with
          *{2}
        rw [valueX, valueA]
        exact openBall.center with positive

theorem constantRealFunction.continuous (L : Real)
  shows Continuous(constantRealFunction(L))
proof
  constructor
  . exact constantRealFunction.function at L
  . intros A realA
    have typedA : Real(A) by exact realA
    exact constantRealFunction.continuous_at with assumption, typedA

-- Translation is represented by the graph x ↦ A + x.
def translateBy (A : Real) : Set :=
  pointwiseAdd(constantRealFunction(A), identity(ℝ))

theorem translateBy.function (A : Real)
  shows translateBy(A) : ℝ → ℝ
proof
  unfold translateBy
  exact pointwiseAdd.function at constantRealFunction(A), identity(ℝ) with
    constantRealFunction.function at A with assumption,
    identity.function at ℝ

theorem translateBy.value (A X : Real)
  shows translateBy(A)(X) = A + X
proof
  have constantFunction := constantRealFunction.function at A
  have identityFunction := identity.function at ℝ
  have constantValue := constantRealFunction.value at A, X with
    *{2}
  have realX : X ∈ ℝ by simp [Real]
  have mappedX := identity.relates with realX
  have identityValue := FunctionBetween.apply_value with identityFunction, mappedX
  have sumValue := pointwiseAdd.value at
    constantRealFunction(A), identity(ℝ), X with
      assumption, constantFunction, identityFunction, *{2}
  unfold translateBy
  rw [constantValue, identityValue] at sumValue
  exact sumValue

theorem translateBy.continuous (A : Real)
  shows Continuous(translateBy(A))
proof
  unfold translateBy
  exact Continuous.add at constantRealFunction(A), identity(ℝ) with
    constantRealFunction.continuous at A with assumption,
    identity.continuous

-- The increment h ↦ f(A+h)-f(A) is useful both in the definition of the
-- difference quotient and in the proof that differentiability is continuous.
def incrementAt (f : Set) (A : Real) : Set :=
  pointwiseSub(
    compose(f, translateBy(A), ℝ, ℝ, ℝ),
    constantRealFunction(f(A)))

theorem incrementAt.function (f : Set) (A : Real)
  requires function : f : ℝ → ℝ
  shows incrementAt(f, A) : ℝ → ℝ
proof
  have realA : A ∈ ℝ by simp [Real]
  have fAInReals := FunctionBetween.apply_closed with function, realA
  have typedFA : Real(f(A)) by exact fAInReals
  have translation := translateBy.function at A
  have shifted := compose.function with translation, function
  have constant := constantRealFunction.function with typedFA
  unfold incrementAt
  exact pointwiseSub.function with shifted, constant

theorem incrementAt.value (f : Set) (A H : Real)
  requires function : f : ℝ → ℝ
  shows incrementAt(f, A)(H) =
    f(A + H) - f(A)
proof
  have realH : H ∈ ℝ by simp [Real]
  have realA : A ∈ ℝ by simp [Real]
  have sumReal := realAddValue.real at A, H with *
  have sumInReals : A + H ∈ ℝ by exact sumReal
  have fShiftInReals := FunctionBetween.apply_closed with function, sumInReals
  have fAInReals := FunctionBetween.apply_closed with function, realA
  have typedFShift : Real(f(A + H)) by exact fShiftInReals
  have typedFA : Real(f(A)) by exact fAInReals
  have translation := translateBy.function at A
  have translationValue := translateBy.value at A, H with
    *{2}
  have shifted := compose.function with translation, function
  have shiftedValue := compose.apply with translation, function, realH
  rw translationValue at shiftedValue
  have constant := constantRealFunction.function with typedFA
  have constantValue := constantRealFunction.value at f(A), H with
    typedFA, assumption
  have shiftedAtHReal : Real(compose(f, translateBy(A), ℝ, ℝ, ℝ)(H)) by rwa shiftedValue
  have constantAtHReal : Real(constantRealFunction(f(A))(H)) by rwa constantValue
  have value := pointwiseSub.value at
    compose(f, translateBy(A), ℝ, ℝ, ℝ),
    constantRealFunction(f(A)), H with
      assumption, shifted, constant, shiftedAtHReal, constantAtHReal
  rw [shiftedValue, constantValue] at value
  simp [incrementAt]

-- The total graph h ↦ (f(A+h)-f(A))/h. Its value at zero is harmless because
-- differentiation uses the punctured limit.
def differenceQuotient (f : Set) (A : Real) : Set :=
  pointwiseDiv(incrementAt(f, A), identity(ℝ))

theorem differenceQuotient.function (f : Set) (A : Real)
  requires function : f : ℝ → ℝ
  shows differenceQuotient(f, A) : ℝ → ℝ
proof
  have numerator := incrementAt.function with function
  unfold differenceQuotient
  exact pointwiseDiv.function at incrementAt(f, A), identity(ℝ) with
    numerator, identity.function at ℝ

theorem differenceQuotient.value (f : Set) (A H : Real)
  requires function : f : ℝ → ℝ
  shows differenceQuotient(f, A)(H) =
    (f(A + H) - f(A)) / H
proof
  have realH : H ∈ ℝ by simp [Real]
  have numerator := incrementAt.function at f, A with function
  have numeratorValue := incrementAt.value at f, A, H with
    *{2}, function
  have identityMapsH := identity.relates with realH
  have identityValue := FunctionBetween.apply_value at
    identity(ℝ), ℝ, ℝ, H, H with
      identity.function at ℝ, realH, identityMapsH
  have numeratorAtHReal := FunctionBetween.apply_closed with numerator, realH
  have typedNumeratorAtH : Real(incrementAt(f, A)(H))
    exact numeratorAtHReal
  have identityAtHReal : Real(identity(ℝ)(H)) by rwa identityValue
  have quotientValue := pointwiseDiv.value at
    incrementAt(f, A), identity(ℝ), H with
      assumption, numerator, identity.function at ℝ,
      typedNumeratorAtH, identityAtHReal
  unfold differenceQuotient
  rw [numeratorValue, identityValue] at quotientValue
  exact quotientValue

predicate HasDerivativeAt (f : Set) (D A : Real) :=
  FunctionBetween(f, ℝ, ℝ) ∧
  TendsToAt(differenceQuotient(f, A), D, realZero)

predicate DifferentiableAt (f : Set) (A : Real) :=
  ∃ D ∈ ℝ, HasDerivativeAt(f, D, A)

theorem HasDerivativeAt.function (f : Set) (D A : Real)
  requires derivative : HasDerivativeAt(f, D, A)
  shows f : ℝ → ℝ
proof
  simp [HasDerivativeAt]

theorem HasDerivativeAt.limit (f : Set) (D A : Real)
  requires derivative : HasDerivativeAt(f, D, A)
  shows TendsToAt(differenceQuotient(f, A), D, realZero)
proof
  simp [HasDerivativeAt]

theorem HasDerivativeAt.differentiable_at (f : Set) (D A : Real)
  requires derivative : HasDerivativeAt(f, D, A)
  shows DifferentiableAt(f, A)
proof
  have realD : D ∈ ℝ by simp [Real]
  unfold DifferentiableAt
  refine ⟨D, realD, derivative⟩

-- In increment coordinates, differentiability says exactly that
-- f(A+h) tends to f(A) as h tends to zero.
theorem HasDerivativeAt.shifted_limit (f : Set) (D A : Real)
  requires derivative : HasDerivativeAt(f, D, A)
  shows TendsToAt(
    compose(f, translateBy(A), ℝ, ℝ, ℝ), f(A), realZero)
proof
  have function := HasDerivativeAt.function with derivative
  have derivativeLimit := HasDerivativeAt.limit with derivative
  have quotientFunction := differenceQuotient.function at f, A with function
  have identityLimit := identity.tends_to_at with realZero.real
  have productLimit := TendsToAt.mul with identityLimit, derivativeLimit
  have zeroProduct := realMulValue.zero_left at D with assumption
  rw zeroProduct at productLimit
  have incrementFunction := incrementAt.function at f, A with function
  have equalIncrement : ∀ H ∈ ℝ, (H ≠ realZero →
      incrementAt(f, A)(H) =
        pointwiseMul(identity(ℝ), differenceQuotient(f, A))(H))
    intros H realH nonzero
    have typedH : Real(H)
      exact realH
    have incrementValue := incrementAt.value at f, A, H with
      assumption, typedH, function
    have quotientValue := differenceQuotient.value at f, A, H with
      assumption, typedH, function
    have numeratorReal := realSubValue.real at
      f(A + H), f(A) with *{2}
    have quotientReal := FunctionBetween.apply_closed with quotientFunction, realH
    have typedQuotient : Real(differenceQuotient(f, A)(H))
      exact quotientReal
    have identityMaps := identity.relates with realH
    have identityValue := FunctionBetween.apply_value at
      identity(ℝ), ℝ, ℝ, H, H with
        identity.function at ℝ, realH, identityMaps
    have typedIdentity : Real(identity(ℝ)(H)) by rwa identityValue
    have productValue := pointwiseMul.value at
      identity(ℝ), differenceQuotient(f, A), H with
        typedH, identity.function at ℝ, quotientFunction,
        typedIdentity, typedQuotient
    have cancellation := realDiv.mul_left with numeratorReal, typedH, nonzero
    calc
      incrementAt(f, A)(H) =
          f(A + H) - f(A) := incrementValue
      _ = H * ((f(A + H) - f(A)) / H) := by cc
      _ = H * differenceQuotient(f, A)(H) :=
        by simp [quotientValue]
      _ = pointwiseMul(identity(ℝ),
          differenceQuotient(f, A))(H) := by cc
  have incrementLimit := TendsToAt.congr with productLimit, incrementFunction, equalIncrement
  have realA : A ∈ ℝ by simp [Real]
  have fAInReals := FunctionBetween.apply_closed with function, realA
  have typedFA : Real(f(A))
    exact fAInReals
  have constantContinuous := constantRealFunction.continuous_at with typedFA, realZero.real
  have constantLimit := ContinuousAt.tends_to with constantContinuous
  have constantValue := constantRealFunction.value with typedFA, realZero.real
  rw constantValue at constantLimit
  have sumLimit := TendsToAt.add with constantLimit, incrementLimit
  have zeroSum := realAddValue.zero_right with typedFA
  rw zeroSum at sumLimit
  have translation := translateBy.function at A with assumption
  have shiftedFunction := compose.function with translation, function
  have equalShift : ∀ H ∈ ℝ, (H ≠ realZero →
      compose(f, translateBy(A), ℝ, ℝ, ℝ)(H) =
        pointwiseAdd(constantRealFunction(f(A)),
          incrementAt(f, A))(H))
    intros H realH nonzero
    have typedH : Real(H)
      exact realH
    have sumReal := realAddValue.real at A, H with assumption, typedH
    have shiftedValue := compose.apply with translation, function, realH
    have translationValue := translateBy.value at A, H with
      assumption, typedH
    rw translationValue at shiftedValue
    have incrementValue := incrementAt.value at f, A, H with
      assumption, typedH, function
    have constantFunction := constantRealFunction.function with typedFA
    have constantAtH := constantRealFunction.value with typedFA, typedH
    have sumInReals : A + H ∈ ℝ
      exact sumReal
    have shiftedInReals := FunctionBetween.apply_closed with function, sumInReals
    have typedShifted : Real(f(A + H))
      exact shiftedInReals
    have incrementAtHReal := FunctionBetween.apply_closed with incrementFunction, realH
    have typedIncrement : Real(incrementAt(f, A)(H))
      exact incrementAtHReal
    have typedConstant : Real(constantRealFunction(f(A))(H)) by rwa constantAtH
    have sumValue := pointwiseAdd.value with
      typedH, constantFunction, incrementFunction, typedConstant, typedIncrement
    have cancel := realSubValue.add_right_cancel with typedShifted, typedFA
    have differenceReal := realSubValue.real with typedShifted, typedFA
    have commute := realAddValue.commutative with typedFA, differenceReal
    have reconstructed : f(A) +
        (f(A + H) - f(A)) = f(A + H)
      calc
        f(A) + (f(A + H) - f(A)) =
            (f(A + H) - f(A)) + f(A) := commute
        _ = f(A + H) := cancel
    rw ← reconstructed at shiftedValue
    rw [constantAtH, incrementValue] at sumValue
    rwa ← sumValue at shiftedValue
  exact TendsToAt.congr with sumLimit, shiftedFunction, equalShift

-- Every differentiable real function is continuous at the point of
-- differentiation.  Translation carries the increment-coordinate result
-- back from zero to A.
theorem HasDerivativeAt.continuous_at (f : Set) (D A : Real)
  requires derivative : HasDerivativeAt(f, D, A)
  shows ContinuousAt(f, A)
proof
  have function := HasDerivativeAt.function with derivative
  have shiftedLimit := HasDerivativeAt.shifted_limit with derivative
  have translation := translateBy.function at A
  have shiftedFunction := compose.function with translation, function
  have realA : A ∈ ℝ by simp [Real]
  have fAInReals := FunctionBetween.apply_closed with function
  have zeroInReals : realZero ∈ ℝ by exact realZero.real
  have shiftedAtZero := compose.apply with translation, function, zeroInReals
  have translationAtZero := translateBy.value at A, realZero with
    assumption, realZero.real
  have addZero := realAddValue.zero_right at A
  rw [translationAtZero, addZero] at shiftedAtZero
  have shiftedContinuous := TendsToAt.continuous_at with shiftedLimit, shiftedAtZero
  have negAReal := realNegValue.real at A
  have backContinuousEverywhere := translateBy.continuous with negAReal
  have backContinuous := Continuous.at at translateBy(-A), A with
    backContinuousEverywhere
  have backFunction := translateBy.function with negAReal
  have backAtA := translateBy.value at -A, A with negAReal, assumption
  have cancelA := realAddValue.neg_left at A
  rw cancelA at backAtA
  have shiftedAtBack : ContinuousAt(
      compose(f, translateBy(A), ℝ, ℝ, ℝ), translateBy(-A)(A))
    rwa backAtA
  have composedContinuous := ContinuousAt.comp with backContinuous, shiftedAtBack
  have equal : ∀ X ∈ ℝ,
      compose(
        compose(f, translateBy(A), ℝ, ℝ, ℝ),
        translateBy(-A), ℝ, ℝ, ℝ)(X) = f(X)
    intros X realX
    have typedX : Real(X) by exact realX
    have backValue := translateBy.value with negAReal, typedX
    have backXInReals := FunctionBetween.apply_closed with backFunction, realX
    have typedBackX : Real(translateBy(-A)(X)) by exact backXInReals
    have outerValue := compose.apply with backFunction, shiftedFunction, realX
    have innerValue := compose.apply with translation, function, backXInReals
    have forwardValue := translateBy.value at
      A, translateBy(-A)(X) with assumption, typedBackX
    rw [forwardValue, backValue] at innerValue
    rw backValue at outerValue
    have associate := realAddValue.associative at A, -A, X with
      assumption, negAReal, typedX
    have cancel := realAddValue.neg_right at A
    have zeroX := realAddValue.zero_left with typedX
    have restored : A + ((-A) + X) = X by rwa [← associate, cancel]
    rw restored at innerValue
    calc
      compose(
        compose(f, translateBy(A), ℝ, ℝ, ℝ),
        translateBy(-A), ℝ, ℝ, ℝ)(X) =
          compose(f, translateBy(A), ℝ, ℝ, ℝ)((-A) + X) := outerValue
      _ = f(X) := innerValue
  exact ContinuousAt.congr with composedContinuous, function, equal

-- Difference quotients commute with negation.
theorem differenceQuotient.neg (f : Set) (A H : Real)
  requires function : f : ℝ → ℝ
  shows differenceQuotient(pointwiseNeg(f), A)(H) =
    -differenceQuotient(f, A)(H)
proof
  have negFunction := pointwiseNeg.function with function
  have sumReal := realAddValue.real at A, H with *
  have realA : A ∈ ℝ by simp [Real]
  have sumInReals : A + H ∈ ℝ by exact sumReal
  have fShiftInReals := FunctionBetween.apply_closed with function, sumInReals
  have fAInReals := FunctionBetween.apply_closed with function, realA
  have typedFShift : Real(f(A + H)) by exact fShiftInReals
  have typedFA : Real(f(A)) by exact fAInReals
  have shiftedNeg := pointwiseNeg.value with sumReal, function, typedFShift
  have centerNeg := pointwiseNeg.value at f, A with
    assumption, function, typedFA
  have negQuotient := differenceQuotient.value at
    pointwiseNeg(f), A, H with *{2}, negFunction
  have quotient := differenceQuotient.value at f, A, H with
    *{2}, function
  have negNumerator := realSubValue.neg_neg with typedFShift, typedFA
  have numeratorReal := realSubValue.real with typedFShift, typedFA
  have negDivision := realDiv.neg_left at
    f(A + H) - f(A), H with
      numeratorReal, assumption
  calc
    differenceQuotient(pointwiseNeg(f), A)(H) =
        (pointwiseNeg(f)(A + H) -
          pointwiseNeg(f)(A)) / H := negQuotient
    _ = ((-f(A + H)) - (-f(A))) / H :=
      by simp [shiftedNeg, centerNeg]
    _ = (-(f(A + H) - f(A))) / H :=
      by simp [negNumerator]
    _ = -((f(A + H) - f(A)) / H) := negDivision
    _ = -differenceQuotient(f, A)(H) := by simp [quotient]

-- The quotient of a pointwise sum is the sum of the quotients.
theorem differenceQuotient.add (f g : Set) (A H : Real)
  requires fFunction : f : ℝ → ℝ
  requires gFunction : g : ℝ → ℝ
  shows differenceQuotient(pointwiseAdd(f, g), A)(H) =
    pointwiseAdd(differenceQuotient(f, A),
      differenceQuotient(g, A))(H)
proof
  have sumFunction := pointwiseAdd.function with fFunction, gFunction
  have quotientF := differenceQuotient.function at f, A with fFunction
  have quotientG := differenceQuotient.function at g, A with gFunction
  have sumReal := realAddValue.real at A, H with *
  have realA : A ∈ ℝ by simp [Real]
  have realH : H ∈ ℝ by simp [Real]
  have sumInReals : A + H ∈ ℝ by exact sumReal
  have fShiftInReals := FunctionBetween.apply_closed with fFunction, sumInReals
  have gShiftInReals := FunctionBetween.apply_closed with gFunction, sumInReals
  have fAInReals := FunctionBetween.apply_closed with fFunction, realA
  have gAInReals := FunctionBetween.apply_closed with gFunction, realA
  have typedFShift : Real(f(A + H)) by exact fShiftInReals
  have typedGShift : Real(g(A + H)) by exact gShiftInReals
  have typedFA : Real(f(A)) by exact fAInReals
  have typedGA : Real(g(A)) by exact gAInReals
  have shiftedSum := pointwiseAdd.value with sumReal, fFunction, gFunction, typedFShift, typedGShift
  have centerSum := pointwiseAdd.value at f, g, A with
    assumption, fFunction, gFunction, typedFA, typedGA
  have sumQuotient := differenceQuotient.value at
    pointwiseAdd(f, g), A, H with *{2}, sumFunction
  have fQuotient := differenceQuotient.value at f, A, H with
    *{2}, fFunction
  have gQuotient := differenceQuotient.value at g, A, H with
    *{2}, gFunction
  have fQuotientInReals := FunctionBetween.apply_closed with quotientF, realH
  have gQuotientInReals := FunctionBetween.apply_closed with quotientG, realH
  have typedFQuotient : Real(differenceQuotient(f, A)(H)) by exact fQuotientInReals
  have typedGQuotient : Real(differenceQuotient(g, A)(H)) by exact gQuotientInReals
  have pointwiseValue := pointwiseAdd.value at
    differenceQuotient(f, A), differenceQuotient(g, A), H with
      assumption, quotientF, quotientG, typedFQuotient, typedGQuotient
  have rearranged := realSubValue.add_add with typedFShift, typedGShift, typedFA, typedGA
  have fDifferenceReal := realSubValue.real with typedFShift, typedFA
  have gDifferenceReal := realSubValue.real with typedGShift, typedGA
  have dividedSum := realDiv.add_left at
    f(A + H) - f(A),
    g(A + H) - g(A), H with
      fDifferenceReal, gDifferenceReal, assumption
  calc
    differenceQuotient(pointwiseAdd(f, g), A)(H) =
        (pointwiseAdd(f, g)(A + H) -
          pointwiseAdd(f, g)(A)) / H := sumQuotient
    _ = ((f(A + H) + g(A + H)) -
        (f(A) + g(A))) / H :=
      by simp [shiftedSum, centerSum]
    _ = ((f(A + H) - f(A)) +
        (g(A + H) - g(A))) / H := by simp [rearranged]
    _ = (f(A + H) - f(A)) / H +
        (g(A + H) - g(A)) / H := dividedSum
    _ = differenceQuotient(f, A)(H) +
        differenceQuotient(g, A)(H) :=
      by simp [fQuotient, gQuotient]
    _ = pointwiseAdd(differenceQuotient(f, A),
        differenceQuotient(g, A))(H) := by simp [pointwiseValue]

-- The product quotient is the sum of the two textbook terms.  The first
-- factor is evaluated at A+h, which is what makes the identity hold
-- before taking a limit.
theorem HasDerivativeAt.neg (f : Set) (D A : Real)
  requires derivative : HasDerivativeAt(f, D, A)
  shows HasDerivativeAt(pointwiseNeg(f), -D, A)
proof
  have function := HasDerivativeAt.function with derivative
  have limit := HasDerivativeAt.limit with derivative
  have negFunction := pointwiseNeg.function with function
  have quotientFunction := differenceQuotient.function at
    pointwiseNeg(f), A with negFunction
  have negLimit := TendsToAt.neg with limit
  have equalAway : ∀ H ∈ ℝ, (H ≠ realZero →
      differenceQuotient(pointwiseNeg(f), A)(H) =
        pointwiseNeg(differenceQuotient(f, A))(H))
    intros H realH nonzero
    have typedH : Real(H) by exact realH
    have quotientF := differenceQuotient.function at f, A with function
    have valueReal := FunctionBetween.apply_closed with quotientF, realH
    have typedValue : Real(differenceQuotient(f, A)(H)) by exact valueReal
    have negValue := pointwiseNeg.value with quotientF
    have quotientNeg := differenceQuotient.neg at f, A, H with
      assumption, typedH, function
    calc
      differenceQuotient(pointwiseNeg(f), A)(H) =
          -differenceQuotient(f, A)(H) := quotientNeg
      _ = pointwiseNeg(differenceQuotient(f, A))(H) :=
        by simp [negValue]
  have quotientLimit := TendsToAt.congr with negLimit, quotientFunction, equalAway
  unfold HasDerivativeAt
  refine ⟨negFunction, quotientLimit⟩

theorem HasDerivativeAt.add (f g : Set) (D E A : Real)
  requires derivativeF : HasDerivativeAt(f, D, A)
  requires derivativeG : HasDerivativeAt(g, E, A)
  shows HasDerivativeAt(pointwiseAdd(f, g), D + E, A)
proof
  have fFunction := HasDerivativeAt.function with derivativeF
  have gFunction := HasDerivativeAt.function with derivativeG
  have limitF := HasDerivativeAt.limit with derivativeF
  have limitG := HasDerivativeAt.limit with derivativeG
  have sumFunction := pointwiseAdd.function with fFunction, gFunction
  have quotientFunction := differenceQuotient.function at
    pointwiseAdd(f, g), A with sumFunction
  have sumLimit := TendsToAt.add with limitF, limitG
  have equalAway : ∀ H ∈ ℝ, (H ≠ realZero →
      differenceQuotient(pointwiseAdd(f, g), A)(H) =
        pointwiseAdd(differenceQuotient(f, A),
          differenceQuotient(g, A))(H))
    intros H realH nonzero
    exact differenceQuotient.add with fFunction, gFunction
  have quotientLimit := TendsToAt.congr with sumLimit, quotientFunction, equalAway
  unfold HasDerivativeAt
  refine ⟨sumFunction, quotientLimit⟩

-- The product rule in its familiar handwritten form.
theorem HasDerivativeAt.sub (f g : Set) (D E A : Real)
  requires derivativeF : HasDerivativeAt(f, D, A)
  requires derivativeG : HasDerivativeAt(g, E, A)
  shows HasDerivativeAt(pointwiseSub(f, g), D - E, A)
proof
  have negDerivative := HasDerivativeAt.neg with derivativeG
  have sumDerivative := HasDerivativeAt.add with derivativeF, negDerivative
  unfold pointwiseSub
  unfold realSubValue
  exact sumDerivative

theorem DifferentiableAt.neg (f : Set) (A : Real)
  requires differentiable : DifferentiableAt(f, A)
  shows DifferentiableAt(pointwiseNeg(f), A)
proof
  obtain ⟨D, realD, derivative⟩ := differentiable
  have typedD : Real(D) by exact realD
  have negDReal := realNegValue.real with typedD
  have negDerivative := HasDerivativeAt.neg with derivative
  witness -D
  constructor
    exact negDReal
  . exact negDerivative

theorem DifferentiableAt.add (f g : Set) (A : Real)
  requires differentiableF : DifferentiableAt(f, A)
  requires differentiableG : DifferentiableAt(g, A)
  shows DifferentiableAt(pointwiseAdd(f, g), A)
proof
  obtain ⟨D, realD, derivativeF⟩ := differentiableF
  obtain ⟨E, realE, derivativeG⟩ := differentiableG
  have typedD : Real(D) by exact realD
  have typedE : Real(E) by exact realE
  have sumReal := realAddValue.real with typedD, typedE
  have sumDerivative := HasDerivativeAt.add with derivativeF, derivativeG
  witness D + E
  constructor
    exact sumReal
  . exact sumDerivative

theorem DifferentiableAt.sub (f g : Set) (A : Real)
  requires differentiableF : DifferentiableAt(f, A)
  requires differentiableG : DifferentiableAt(g, A)
  shows DifferentiableAt(pointwiseSub(f, g), A)
proof
  have negG := DifferentiableAt.neg with differentiableG
  have sum := DifferentiableAt.add with differentiableF, negG
  unfold pointwiseSub
  exact sum

theorem constantRealFunction.has_derivative_at (C A : Real)
  shows HasDerivativeAt(constantRealFunction(C), realZero, A)
proof
  have function := constantRealFunction.function at C
  have quotientFunction := differenceQuotient.function at
    constantRealFunction(C), A with function
  have zeroContinuous := constantRealFunction.continuous_at with realZero.real, realZero.real
  have zeroLimit := ContinuousAt.tends_to with zeroContinuous
  have zeroValue := constantRealFunction.value with realZero.real, realZero.real
  rw zeroValue at zeroLimit
  have equalAway : ∀ H ∈ ℝ, (H ≠ realZero →
      differenceQuotient(constantRealFunction(C), A)(H) =
        constantRealFunction(realZero)(H))
    intros H realH nonzero
    have typedH : Real(H) by exact realH
    have sumReal := realAddValue.real at A, H with assumption, typedH
    have quotientValue := differenceQuotient.value at
      constantRealFunction(C), A, H with
        assumption, typedH, function
    have shiftedValue := constantRealFunction.value at C, A + H with
      assumption, sumReal
    have centerValue := constantRealFunction.value at C, A with
      *{2}
    have numeratorZero := realSubValue.self at C
    have quotientZero := realDiv.zero_left with typedH
    have constantZero := constantRealFunction.value with realZero.real, typedH
    calc
      differenceQuotient(constantRealFunction(C), A)(H) =
          (constantRealFunction(C)(A + H) -
            constantRealFunction(C)(A)) / H := quotientValue
      _ = (C - constantRealFunction(C)(A)) / H := by simp [shiftedValue]
      _ = (C - C) / H := by simp [centerValue]
      _ = 0 / H := by simp [numeratorZero]
      _ = 0 := quotientZero
      _ = constantRealFunction(realZero)(H) := by simp [constantZero]
  have quotientLimit := TendsToAt.congr with zeroLimit, quotientFunction, equalAway
  unfold HasDerivativeAt
  refine ⟨function, quotientLimit⟩

theorem identity.has_derivative_at (A : Real)
  shows HasDerivativeAt(identity(ℝ), realOne, A)
proof
  have function := identity.function at ℝ
  have quotientFunction := differenceQuotient.function with function
  have oneContinuous := constantRealFunction.continuous_at with realOne.real, realZero.real
  have oneLimit := ContinuousAt.tends_to with oneContinuous
  have oneValue := constantRealFunction.value with realOne.real, realZero.real
  rw oneValue at oneLimit
  have equalAway : ∀ H ∈ ℝ, (H ≠ realZero →
      differenceQuotient(identity(ℝ), A)(H) =
        constantRealFunction(realOne)(H))
    intros H realH nonzero
    have typedH : Real(H) by exact realH
    have sumReal := realAddValue.real at A, H with assumption, typedH
    have sumInReals : A + H ∈ ℝ by exact sumReal
    have quotientValue := differenceQuotient.value at
      identity(ℝ), A, H with assumption, typedH, function
    have shiftedMaps := identity.relates with sumInReals
    have shiftedValue := FunctionBetween.apply_value with function, shiftedMaps
    have realA : A ∈ ℝ by simp [Real]
    have centerMaps := identity.relates with realA
    have centerValue := FunctionBetween.apply_value with function, centerMaps
    have numerator := realSubValue.add_cancel_left at A, H with
      assumption, typedH
    have quotientOne := realDiv.self with nonzero
    have constantOne := constantRealFunction.value with realOne.real, typedH
    calc
      differenceQuotient(identity(ℝ), A)(H) =
          (identity(ℝ)(A + H) - identity(ℝ)(A)) / H :=
            quotientValue
      _ = ((A + H) - identity(ℝ)(A)) / H := by simp [shiftedValue]
      _ = ((A + H) - A) / H := by simp [centerValue]
      _ = H / H := by simp [numerator]
      _ = 1 := quotientOne
      _ = constantRealFunction(realOne)(H) := by simp [constantOne]
  have quotientLimit := TendsToAt.congr with oneLimit, quotientFunction, equalAway
  unfold HasDerivativeAt
  refine ⟨function, quotientLimit⟩

theorem HasDerivativeAt.unique (f : Set) (D E A : Real)
  requires first : HasDerivativeAt(f, D, A)
  requires second : HasDerivativeAt(f, E, A)
  shows D = E
proof
  have firstLimit := HasDerivativeAt.limit with first
  have secondLimit := HasDerivativeAt.limit with second
  exact TendsToAt.unique with firstLimit, secondLimit

-- Section: ProductRule
theorem differenceQuotient.mul (f g : Set) (A H : Real)
  requires fFunction : f : ℝ → ℝ
  requires gFunction : g : ℝ → ℝ
  shows differenceQuotient(pointwiseMul(f, g), A)(H) =
    pointwiseAdd(
      pointwiseMul(differenceQuotient(f, A),
        constantRealFunction(g(A))),
      pointwiseMul(compose(f, translateBy(A), ℝ, ℝ, ℝ),
        differenceQuotient(g, A)))(H)
proof
  have productFunction := pointwiseMul.function with fFunction, gFunction
  have quotientF := differenceQuotient.function at f, A with fFunction
  have quotientG := differenceQuotient.function at g, A with gFunction
  have realA : A ∈ ℝ by simp [Real]
  have realH : H ∈ ℝ by simp [Real]
  have sumReal := realAddValue.real at A, H with *
  have sumInReals : A + H ∈ ℝ by exact sumReal
  have fShiftInReals := FunctionBetween.apply_closed with fFunction, sumInReals
  have gShiftInReals := FunctionBetween.apply_closed with gFunction, sumInReals
  have fAInReals := FunctionBetween.apply_closed with fFunction, realA
  have gAInReals := FunctionBetween.apply_closed with gFunction, realA
  have typedFShift : Real(f(A + H)) by exact fShiftInReals
  have typedGShift : Real(g(A + H)) by exact gShiftInReals
  have typedFA : Real(f(A)) by exact fAInReals
  have typedGA : Real(g(A)) by exact gAInReals
  have shiftedProduct := pointwiseMul.value with sumReal, fFunction, gFunction, typedFShift, typedGShift
  have centerProduct := pointwiseMul.value at f, g, A with
    assumption, fFunction, gFunction, typedFA, typedGA
  have productQuotient := differenceQuotient.value at
    pointwiseMul(f, g), A, H with *{2}, productFunction
  have fQuotient := differenceQuotient.value at f, A, H with
    *{2}, fFunction
  have gQuotient := differenceQuotient.value at g, A, H with
    *{2}, gFunction
  have fDifferenceReal := realSubValue.real with typedFShift, typedFA
  have gDifferenceReal := realSubValue.real with typedGShift, typedGA
  have firstProductReal := realMulValue.real with typedFShift, gDifferenceReal
  have secondProductReal := realMulValue.real with typedGA, fDifferenceReal
  have rearranged := realSubValue.mul_mul with typedFShift, typedGShift, typedFA, typedGA
  have dividedSum := realDiv.add_left at
    f(A + H) * (g(A + H) - g(A)),
    g(A) * (f(A + H) - f(A)), H with
      firstProductReal, secondProductReal, assumption
  have factorFirst := realDiv.mul_numerator_left at
    f(A + H), g(A + H) - g(A), H with
      typedFShift, gDifferenceReal, assumption
  have factorSecond := realDiv.mul_numerator_left at
    g(A), f(A + H) - f(A), H with
      typedGA, fDifferenceReal, assumption
  have fQuotientInReals := FunctionBetween.apply_closed with quotientF, realH
  have gQuotientInReals := FunctionBetween.apply_closed with quotientG, realH
  have typedFQuotient : Real(differenceQuotient(f, A)(H)) by exact fQuotientInReals
  have typedGQuotient : Real(differenceQuotient(g, A)(H)) by exact gQuotientInReals
  have constantGA := constantRealFunction.function with typedGA
  have constantValue := constantRealFunction.value at g(A), H with
    typedGA, assumption
  have typedConstant : Real(constantRealFunction(g(A))(H)) by rwa constantValue
  have translation := translateBy.function at A
  have shiftedF := compose.function with translation, fFunction
  have translationValue := translateBy.value at A, H with
    *{2}
  have shiftedValue := compose.apply with translation, fFunction, realH
  rw translationValue at shiftedValue
  have typedShifted : Real(compose(f, translateBy(A), ℝ, ℝ, ℝ)(H)) by rwa shiftedValue
  have firstTermFunction := pointwiseMul.function with quotientF, constantGA
  have secondTermFunction := pointwiseMul.function with shiftedF, quotientG
  have firstTermValue := pointwiseMul.value at
    differenceQuotient(f, A), constantRealFunction(g(A)), H with
      assumption, quotientF, constantGA, typedFQuotient, typedConstant
  have secondTermValue := pointwiseMul.value at
    compose(f, translateBy(A), ℝ, ℝ, ℝ), differenceQuotient(g, A), H with
      assumption, shiftedF, quotientG, typedShifted, typedGQuotient
  have typedFirstTerm := realMulValue.real with typedFQuotient, typedConstant
  have typedSecondTerm := realMulValue.real with typedShifted, typedGQuotient
  have typedFirstTermValue : Real(pointwiseMul(differenceQuotient(f, A),
        constantRealFunction(g(A)))(H))
    rwa firstTermValue
  have typedSecondTermValue : Real(pointwiseMul(compose(f, translateBy(A), ℝ, ℝ, ℝ),
        differenceQuotient(g, A))(H))
    rwa secondTermValue
  have sumValue := pointwiseAdd.value at
    pointwiseMul(differenceQuotient(f, A),
      constantRealFunction(g(A))),
    pointwiseMul(compose(f, translateBy(A), ℝ, ℝ, ℝ),
      differenceQuotient(g, A)), H with
        assumption, firstTermFunction, secondTermFunction,
        typedFirstTermValue, typedSecondTermValue
  have commuteFirst := realMulValue.commutative with typedGA, typedFQuotient
  have firstRawReal := realMulValue.real with typedFShift, typedGQuotient
  have secondRawReal := realMulValue.real with typedGA, typedFQuotient
  have commuteSum := realAddValue.commutative with firstRawReal, secondRawReal
  calc
    differenceQuotient(pointwiseMul(f, g), A)(H) =
        (pointwiseMul(f, g)(A + H) -
          pointwiseMul(f, g)(A)) / H := productQuotient
    _ = (f(A + H) * g(A + H) -
        f(A) * g(A)) / H :=
      by simp [shiftedProduct, centerProduct]
    _ = (f(A + H) * (g(A + H) - g(A)) +
        g(A) * (f(A + H) - f(A))) / H :=
      by simp [rearranged]
    _ = (f(A + H) * (g(A + H) - g(A))) / H +
        (g(A) * (f(A + H) - f(A))) / H := dividedSum
    _ = f(A + H) *
          ((g(A + H) - g(A)) / H) +
        g(A) * ((f(A + H) - f(A)) / H) :=
      by simp [factorFirst, factorSecond]
    _ = f(A + H) * differenceQuotient(g, A)(H) +
        g(A) * differenceQuotient(f, A)(H) :=
      by simp [fQuotient, gQuotient]
    _ = differenceQuotient(f, A)(H) * g(A) +
        f(A + H) * differenceQuotient(g, A)(H) := by cc
    _ = pointwiseMul(differenceQuotient(f, A),
          constantRealFunction(g(A)))(H) +
        pointwiseMul(compose(f, translateBy(A), ℝ, ℝ, ℝ),
          differenceQuotient(g, A))(H) :=
      by simp [firstTermValue, secondTermValue, constantValue, shiftedValue]
    _ = pointwiseAdd(
        pointwiseMul(differenceQuotient(f, A),
          constantRealFunction(g(A))),
        pointwiseMul(compose(f, translateBy(A), ℝ, ℝ, ℝ),
          differenceQuotient(g, A)))(H) := by cc

theorem HasDerivativeAt.mul (f g : Set) (D E A : Real)
  requires derivativeF : HasDerivativeAt(f, D, A)
  requires derivativeG : HasDerivativeAt(g, E, A)
  shows HasDerivativeAt(pointwiseMul(f, g),
    realAddValue(
      realMulValue(D, g(A)),
      realMulValue(f(A), E)), A)
proof
  have fFunction := HasDerivativeAt.function with derivativeF
  have gFunction := HasDerivativeAt.function with derivativeG
  have limitF := HasDerivativeAt.limit with derivativeF
  have limitG := HasDerivativeAt.limit with derivativeG
  have productFunction := pointwiseMul.function with fFunction, gFunction
  have quotientFunction := differenceQuotient.function at
    pointwiseMul(f, g), A with productFunction
  have realA : A ∈ ℝ by simp [Real]
  have fAInReals := FunctionBetween.apply_closed with fFunction, realA
  have gAInReals := FunctionBetween.apply_closed with gFunction, realA
  have typedFA : Real(f(A)) by exact fAInReals
  have typedGA : Real(g(A)) by exact gAInReals
  have firstDerivativeReal := realMulValue.real at
    D, g(A) with assumption, typedGA
  have secondDerivativeReal := realMulValue.real at
    f(A), E with typedFA, assumption
  have productDerivativeReal := realAddValue.real with firstDerivativeReal, secondDerivativeReal
  have constantContinuous := constantRealFunction.continuous_at with typedGA, realZero.real
  have constantLimit := ContinuousAt.tends_to with constantContinuous
  have constantValue := constantRealFunction.value with typedGA, realZero.real
  rw constantValue at constantLimit
  have firstLimit := TendsToAt.mul with limitF, constantLimit
  have shiftedLimit := HasDerivativeAt.shifted_limit with derivativeF
  have secondLimit := TendsToAt.mul with shiftedLimit, limitG
  have sumLimit := TendsToAt.add with firstLimit, secondLimit
  have equalAway : ∀ H ∈ ℝ, (H ≠ realZero →
      differenceQuotient(pointwiseMul(f, g), A)(H) =
        pointwiseAdd(
          pointwiseMul(differenceQuotient(f, A),
            constantRealFunction(g(A))),
          pointwiseMul(compose(f, translateBy(A), ℝ, ℝ, ℝ),
            differenceQuotient(g, A)))(H))
    intros H realH nonzero
    exact differenceQuotient.mul with fFunction, gFunction
  have quotientLimit := TendsToAt.congr at
    pointwiseAdd(
      pointwiseMul(differenceQuotient(f, A),
        constantRealFunction(g(A))),
      pointwiseMul(compose(f, translateBy(A), ℝ, ℝ, ℝ),
        differenceQuotient(g, A))),
    differenceQuotient(pointwiseMul(f, g), A),
    D * g(A) + f(A) * E, realZero with
      sumLimit, quotientFunction, equalAway
  unfold realAddValue
  unfold realMulValue
  unfold call2 at productDerivativeReal
  unfold call2 at quotientLimit
  unfold HasDerivativeAt
  refine ⟨productFunction, quotientLimit⟩

theorem DifferentiableAt.mul (f g : Set) (A : Real)
  requires differentiableF : DifferentiableAt(f, A)
  requires differentiableG : DifferentiableAt(g, A)
  shows DifferentiableAt(pointwiseMul(f, g), A)
proof
  obtain ⟨D, realD, derivativeF⟩ := differentiableF
  obtain ⟨E, realE, derivativeG⟩ := differentiableG
  have typedD : Real(D) by exact realD
  have typedE : Real(E) by exact realE
  have fFunction := HasDerivativeAt.function with derivativeF
  have gFunction := HasDerivativeAt.function with derivativeG
  have realA : A ∈ ℝ by simp [Real]
  have fAInReals := FunctionBetween.apply_closed with fFunction, realA
  have gAInReals := FunctionBetween.apply_closed with gFunction, realA
  have typedFA : Real(f(A)) by exact fAInReals
  have typedGA : Real(g(A)) by exact gAInReals
  have firstReal := realMulValue.real with typedD, typedGA
  have secondReal := realMulValue.real with typedFA, typedE
  have derivativeReal := realAddValue.real with firstReal, secondReal
  have productDerivative := HasDerivativeAt.mul with derivativeF, derivativeG
  unfold realAddValue at productDerivative
  unfold realMulValue at productDerivative
  witness D * g(A) + f(A) * E
  constructor
    exact derivativeReal
  . unfold call2
    exact productDerivative