-- Real analysis: Taylor's theorem with the Lagrange remainder
-- 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.Taylor
import TaylorPolynomial

-- Section: Taylor
-- Taylor's theorem with the Lagrange remainder: when f is n + 1 times
-- differentiable on [A, B], some C strictly between A and B has
-- f(B) = ∑_{k ≤ n} f⁽ᵏ⁾(A)(B − A)ᵏ/k! + f⁽ⁿ⁺¹⁾(C)(B − A)ⁿ⁺¹/(n + 1)!.
-- The proof applies the Cauchy mean value theorem to the approximant
-- x ↦ ∑_{k ≤ n} f⁽ᵏ⁾(x)(B − x)ᵏ/k! and to x ↦ (B − x)ⁿ⁺¹.
theorem Taylor.theorem (f : ℝ → ℝ) (A B : Real, n : Nat)
  requires ordered : A < B
  requires smooth : ∀ k ≤ n,
    ∀ X ∈ closedInterval(A, B), DifferentiableAt(f⁽ᵏ⁾, X)
  shows ∃ C ∈ ℝ, A < C ∧ C < B ∧
    f(B) = ∑ k ≤ n, f⁽ᵏ⁾(A) * (B - A) ^ k / k! +
      f⁽ⁿ⁺¹⁾(C) * (B - A) ^ (n + 1) / (n + 1)!
proof
  have orderedWeak := realLess.implies_le with ordered
  have memberA := closedInterval.left_mem with orderedWeak
  have function : f : ℝ → ℝ by assumption
  have approximantFunction := taylorApproximant.function at f, B, n with
    function
  have powerFunction := basePower.function at B, n + 1
  have intervalReal := closedInterval.subset_reals at A, B
  have pointwise : ∀ X ∈ closedInterval(A, B), ∀ k, Nat(k) → k ≤ n →
      DifferentiableAt(iteratedDerivative(f, k), X)
    intros X memberX k naturalK below
    have onInterval := smooth with k, naturalK, below
    exact onInterval with X, memberX
  have approximantDifferentiable : ∀ X ∈ closedInterval(A, B),
      DifferentiableAt(taylorApproximant(f, B, n), X)
    intros X memberX
    have realX := intervalReal with X, memberX
    have typedX : Real(X) by exact realX
    have atX := pointwise with X, memberX
    have implication := taylorApproximant.derivative at f, B, X, n with function
    have hasDerivative := implication with atX
    exact HasDerivativeAt.differentiable_at with hasDerivative
  have powerDifferentiable : ∀ X ∈ closedInterval(A, B),
      DifferentiableAt(basePower(B, n + 1), X)
    intros X memberX
    have realX := intervalReal with X, memberX
    have typedX : Real(X) by exact realX
    have hasDerivative := basePower.derivative at B, X, n
    exact HasDerivativeAt.differentiable_at with hasDerivative
  have continuousApproximant := ContinuousOn.of_differentiable with
    approximantFunction, intervalReal, approximantDifferentiable
  have continuousPower := ContinuousOn.of_differentiable with
    powerFunction, intervalReal, powerDifferentiable
  have interior : ∀ C ∈ ℝ, A < C → C < B → C ∈ closedInterval(A, B)
    intros C realC left right
    have lower := realLess.implies_le with left
    have upper := realLess.implies_le with right
    have data : C ∈ ℝ ∧ A ≤ C ∧ C ≤ B by refine ⟨realC, lower, upper⟩
    exact closedInterval.membership.backward with data
  have interiorApproximant : ∀ C ∈ ℝ, A < C → C < B →
      DifferentiableAt(taylorApproximant(f, B, n), C)
    intros C realC left right
    have member := interior with C, realC, left, right
    exact approximantDifferentiable with C, member
  have interiorPower : ∀ C ∈ ℝ, A < C → C < B →
      DifferentiableAt(basePower(B, n + 1), C)
    intros C realC left right
    have member := interior with C, realC, left, right
    exact powerDifferentiable with C, member
  obtain ⟨C, realC, left, right, equation⟩ := CauchyMeanValue.theorem with
    approximantFunction, powerFunction, ordered, continuousApproximant,
    continuousPower, interiorApproximant, interiorPower
  have typedC : Real(C) by exact realC
  have memberC := interior with C, realC, left, right
  have atC := pointwise with C, memberC
  have implication := taylorApproximant.derivative at f, B, C, n with function
  have approximantDerivative := implication with atC
  have approximantValue := derivative.of_has_derivative with
    approximantDerivative
  have powerDerivative := basePower.derivative at B, C, n
  have powerDerivativeValue := derivative.of_has_derivative with powerDerivative
  have powerAtB : basePower(B, n + 1)(B) = (0:R)
    have value := basePower.value at B, B, n + 1
    have self := realSubValue.self at B
    have one : n + 1 = succ(n) := add.right_one with *
    rw [value, self, one]
    exact realPow.zero_succ at n
  have powerAtA := basePower.value at B, A, n + 1
  have approximantAtB := taylorApproximant.at_endpoint at f, B, n with function
  have approximantAtA := taylorApproximant.value at f, B, A, n with function
  rw [approximantValue, powerDerivativeValue, powerAtB, powerAtA,
    approximantAtB, approximantAtA] at equation
  have realB : B ∈ ℝ by simp [Real]
  have typedFB : Real(f(B))
    exact FunctionBetween.apply_closed with function, realB
  have nextFunction := iteratedDerivative.function at f, n + 1 with function
  have typedNext : Real(iteratedDerivative(f, n + 1)(C))
    exact FunctionBetween.apply_closed with nextFunction, realC
  have typedWidth := realSubValue.real at B, A
  have typedGap := realSubValue.real at B, C
  have gapNonzero : B - C ≠ 0
    intro zero
    linarith
  have powerNonzero := realPow.nonzero at B - C, n with typedGap, gapNonzero
  have typedPowerGap := realPow.real at B - C, n
  have typedPowerWidth := realPow.real at B - A, n + 1
  have typedNat := realNat.real at n
  have successorNonzero := realNat.succ_nonzero at n
  have factorialNat := factorial.closed at n
  have factorialSuccNat := factorial.closed at n + 1
  have typedFactorial := realNat.real at n!
  have typedFactorialSucc := realNat.real at (n + 1)!
  have factorialNonzero := realNat.factorial_nonzero at n
  have factorialSplit := realNat.factorial_succ at n
  have inverseSplit : realInvValue(realNat((n + 1)!)) =
      realInvValue((realNat(n):R) + 1) * realInvValue(realNat(n!))
    rw factorialSplit
    exact realInvValue.mul with successorNonzero, factorialNonzero
  have unitGap := realInvValue.mul_right with powerNonzero
  have unitSuccessor := realInvValue.mul_right with successorNonzero
  have typedGapInverse := realInvValue.real with typedPowerGap
  have typedSuccessorInverse := realInvValue.real at (realNat(n):R) + 1
  have typedFactorialInverse := realInvValue.real with typedFactorial
  have typedFactorialSuccInverse := realInvValue.real with typedFactorialSucc
  have sequence : RealSequence(taylorTerms(f, A, B))
    unfold RealSequence
    exact taylorTerms.function at f, A, B with function
  have typedSum := realPartialSum.real at taylorTerms(f, A, B), n + 1 with
    sequence
  have typedProductGap := realMulValue.real with typedNext, typedPowerGap
  have typedProductWidth := realMulValue.real with typedNext, typedPowerWidth
  have divisionGap := realDiv.eq_mul_inv with typedProductGap, typedFactorial
  have divisionWidth := realDiv.eq_mul_inv with
    typedProductWidth, typedFactorialSucc
  rw divisionGap at equation
  witness C
  refine ⟨realC, left, right, _⟩
  rw divisionWidth
  symm
  linear_combination ((f(B) - (∑ k < n + 1, taylorTerms(f, A, B)(k):R)) *
        ((realNat(n):R) + 1) * realInvValue((realNat(n):R) + 1) -
      iteratedDerivative(f, n + 1)(C) * realInvValue(realNat(n!)) *
        (B - A) ^ (n + 1) * realInvValue((realNat(n):R) + 1)) * unitGap +
    (f(B) - (∑ k < n + 1, taylorTerms(f, A, B)(k):R)) * unitSuccessor +
    (iteratedDerivative(f, n + 1)(C) * (B - A) ^ (n + 1)) * inverseSplit -
    (realInvValue((B - C) ^ n) * realInvValue((realNat(n):R) + 1)) * equation

-- With n = 0 the sum is the single term f(A) and the remainder is
-- f′(C)(B − A): Taylor's theorem specializes to the mean value theorem.
theorem Taylor.mean_value (f : ℝ → ℝ) (A B : Real)
  requires ordered : A < B
  requires differentiable : ∀ X ∈ closedInterval(A, B), DifferentiableAt(f, X)
  shows ∃ C ∈ ℝ, A < C ∧ C < B ∧ f(B) = f(A) + f′(C) * (B - A)
proof
  have function : f : ℝ → ℝ by assumption
  have smooth : ∀ k ≤ (0:N),
      ∀ X ∈ closedInterval(A, B), DifferentiableAt(f⁽ᵏ⁾, X)
    intros k naturalK below
    have split := natural.zero_or_successor with naturalK
    cases split with zero | successor
    . intros X memberX
      rw [zero, iteratedDerivative.zero]
      exact differentiable with X, memberX
    . obtain ⟨j, naturalJ, equal⟩ := successor
      rw equal at below
      exfalso
      arith
  obtain ⟨C, realC, left, right, expansion⟩ := Taylor.theorem at f, A, B, (0:N) with
    function, ordered, smooth
  have typedC : Real(C) by exact realC
  have one : (0:N) + 1 = succ(0) := add.right_one at (0:N)
  rw one at expansion
  rw [realPartialSum.succ, realPartialSum.zero] at expansion
  have termValue := taylorTerms.value at f, A, B, (0:N) with function
  rw [termValue, iteratedDerivative.zero, realPow.zero] at expansion
  have factorialZero : (0:N)! = 1 by simp [factorial]
  have factorialOne : succ(0)! = 1 by simp [factorial]
  rw [factorialZero, factorialOne, realNat.one] at expansion
  rw [← one, iteratedDerivative.succ, iteratedDerivative.zero, one] at expansion
  have realA : A ∈ ℝ by simp [Real]
  have typedFA : Real(f(A)) by exact FunctionBetween.apply_closed with function, realA
  have typedDerivative := derivative.real at f, C
  have typedWidth := realSubValue.real at B, A
  have powerOne : (B - A) ^ succ(0) = B - A
    rw [realPow.succ, realPow.zero]
    exact realMulValue.one_right with typedWidth
  have unitA := realMulValue.one_right with typedFA
  have typedProduct := realMulValue.real with typedDerivative, typedWidth
  rw [powerOne, unitA, realDiv.one_right, realDiv.one_right] at expansion
  have zeroLeft := realAddValue.zero_left with typedFA
  unfold call2 at zeroLeft
  unfold realAddValue at expansion
  rw zeroLeft at expansion
  witness C
  refine ⟨realC, left, right, _⟩
  exact expansion