-- Archimedean limits: reciprocals, successor ratios, and powers below one
-- 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.ArchimedeanLimits
import Series
import HarmonicSeries

-- Section: NullSequences
-- A sequence whose absolute values are eventually below every positive
-- radius converges to zero.
theorem ConvergesTo.zero_of_abs_bounds (u : Set)
  requires sequence : RealSequence(u)
  requires small : ∀ E ∈ ℝ, 0 < E →
    ∃ N ∈ ℕ, ∀ n ∈ ℕ, N ≤ n → abs(u[n]) < E
  shows ConvergesTo(u, (0:R))
proof
  constructor
  . exact sequence
  . intros E realE positive
    have typedE : Real(E) by exact realE
    obtain ⟨N, naturalN, close⟩ := small with E, realE, positive
    witness N
    constructor
    . exact naturalN
    . intros n naturalIndex after
      have typedN : Nat(n) by exact nat.omega.backward at n with naturalIndex
      have typedValue := RealSequence.value at u, n with sequence
      have bound := close with n, naturalIndex, after
      have zero := realSubValue.zero_right with typedValue
      apply openBall.membership.backward
      constructor
      . exact typedValue
      . unfold dist
        rwa zero

theorem realInvValue.involutive (A : Real)
  requires nonzero : A ≠ 0
  shows realInvValue(realInvValue(A)) = A
proof
  have inverseNonzero := realInvValue.nonzero with nonzero
  have law := realInvValue.mul_left with nonzero
  have equal := realInvValue.eq_of_mul_right with inverseNonzero, law
  symm at equal
  exact equal

-- Section: Reciprocals
-- The reciprocals 1/(n + 1) tend to zero: this is the Archimedean property.
theorem harmonicTerms.converges_zero
  shows ConvergesTo(harmonicTerms, (0:R))
proof
  have sequence := harmonicTerms.sequence
  apply ConvergesTo.zero_of_abs_bounds with sequence
  intros E realE positive
  have typedE : Real(E) by exact realE
  have nonzeroE : E ≠ 0
    intro zero
    rw zero at positive
    linarith
  have scaleReal := realInvValue.real at E
  have scalePositive := realInvValue.positive with positive
  have unscaled := realInvValue.involutive with nonzeroE
  obtain ⟨N, typedN, above⟩ := Real.below_natural at realInvValue(E)
  have naturalN := nat.omega.forward with typedN
  witness N
  constructor
  . exact naturalN
  . intros n naturalIndex after
    have typedIndex : Nat(n) by exact nat.omega.backward at n with naturalIndex
    have countReal := realNat.real at n
    have denominatorReal : (realNat(n):R) + 1 ∈ ℝ by typecheck
    have monotone := realNat.monotone with after
    have one := realZero.lt_one
    have denominatorAbove : realInvValue(E) < (realNat(n):R) + 1 by linarith
    have termReal := RealSequence.value at harmonicTerms, n with sequence
    have value := harmonicTerms.value at n
    have termPositive := harmonicTerms.positive at n
    have termNonnegative := realLess.implies_le with termPositive
    have absolute := abs.of_nonnegative with termNonnegative
    have inverted := realInvValue.antitone_positive at
      realInvValue(E), (realNat(n):R) + 1 with scalePositive, denominatorAbove
    rw unscaled at inverted
    rw [absolute, value]
    exact inverted

-- The harmonic series is the standard example of a series that fails to sum.
theorem harmonicSeries.not_summable
  shows ¬Summable(harmonicTerms)
proof
  intro summable
  obtain ⟨L, realL, sums⟩ := summable
  have converges := SumsTo.converges with sums
  have limit : ∃ L ∈ ℝ, ConvergesTo(prefixSums(harmonicTerms), L)
    witness L
    constructor
    . exact realL
    . exact converges
  exact harmonicSeries.diverges with limit

-- Section: SuccessorRatios
-- The ratios n/(n + 1).
defgraph successorRatios (n ∈ ℕ) : ℝ :=
  (realNat(n):R) / ((realNat(n):R) + 1)
proof
  have typedN := nat.omega.backward with member
  have countReal := realNat.real at n
  exact realDiv.real at (realNat(n):R), (realNat(n):R) + 1 with *

theorem successorRatios.sequence
  shows RealSequence(successorRatios)
proof
  exact successorRatios.function

-- n/(n + 1) = 1 - 1/(n + 1).
theorem successorRatios.sub_harmonic (n : Nat)
  shows successorRatios[n] = (1:R) - harmonicTerms(n)
proof
  have countReal := realNat.real at n
  have denominatorReal : (realNat(n):R) + 1 ∈ ℝ by typecheck
  have nonzero := realNat.successor_nonzero at n
  have inverseReal := realInvValue.real at (realNat(n):R) + 1
  have law := realInvValue.mul_left with nonzero
  have product : ((1:R) - realInvValue((realNat(n):R) + 1)) *
      ((realNat(n):R) + 1) = (realNat(n):R)
    have expand : ((1:R) - realInvValue((realNat(n):R) + 1)) *
        ((realNat(n):R) + 1) =
      ((realNat(n):R) + 1) -
        realInvValue((realNat(n):R) + 1) * ((realNat(n):R) + 1) by ring
    rw [expand, law]
    ring
  have quotient := realDiv.of_mul_right with nonzero, product
  have value := successorRatios.value at n
  have termValue := harmonicTerms.value at n
  rw [value, termValue]
  symm at quotient
  exact quotient

theorem successorRatios.converges_one
  shows ConvergesTo(successorRatios, (1:R))
proof
  have sequence := successorRatios.sequence
  have harmonics := harmonicTerms.sequence
  have constants := constantSequence.sequence at (1:R) with realOne.real
  have constantConverges := constantSequence.converges at (1:R) with
    realOne.real
  have limit := ConvergesTo.sub with
    constantConverges, harmonicTerms.converges_zero
  have zero := realSubValue.zero_right with realOne.real
  rw zero at limit
  have equal : ∀ n ∈ ℕ, successorRatios[n] =
      subSequence(constantSequence((1:R)), harmonicTerms)[n]
    intros n naturalN
    have typedN : Nat(n) by exact nat.omega.backward at n with naturalN
    have constantReal := RealSequence.value at
      constantSequence((1:R)), n with constants
    have termReal := RealSequence.value at harmonicTerms, n with harmonics
    have difference := subSequence.value at
      constantSequence((1:R)), harmonicTerms, n with
        constants, harmonics, constantReal, termReal
    have constant := constantSequence.value at (1:R), n with *
    have ratio := successorRatios.sub_harmonic at n
    rw [difference, constant, ratio]
    rfl
  exact ConvergesTo.of_equal_terms with sequence, limit, equal

-- Section: Powers
theorem realPow.nonnegative (A : Real, n : Nat)
  requires nonnegative : 0 ≤ A
  shows 0 ≤ A ^ n
proof
  induction n
  | zero =>
    rw realPow.zero
    exact realLess.implies_le with realZero.lt_one
  | succ(k) using ih =>
    have power := realPow.real at A, k
    rw realPow.succ
    exact realMulValue.nonnegative with nonnegative, ih

theorem realPow.le_one (A : Real, n : Nat)
  requires nonnegative : 0 ≤ A
  requires bounded : A ≤ 1
  shows A ^ n ≤ 1
proof
  induction n
  | zero =>
    rw realPow.zero
    exact realLe.reflexive with realOne.real
  | succ(k) using ih =>
    have power := realPow.real at A, k
    have scaled := realLe.mul_left_of_nonnegative with ih, nonnegative
    have one := realMulValue.one_right at A
    rw one at scaled
    rw realPow.succ
    exact realLe.transitive with scaled, bounded

theorem realPow.abs (A : Real, n : Nat)
  shows abs(A ^ n) = abs(A) ^ n
proof
  induction n
  | zero =>
    rw [realPow.zero, realPow.zero]
    have oneNonnegative := realLess.implies_le with realZero.lt_one
    exact abs.of_nonnegative with oneNonnegative
  | succ(k) using ih =>
    have power := realPow.real at A, k
    have product := abs.mul at A, A ^ k with assumption, power
    rw [realPow.succ, realPow.succ, product, ih]
    rfl

-- Bernoulli's inequality.
theorem realPow.bernoulli (H : Real, n : Nat)
  requires nonnegative : 0 ≤ H
  shows (1:R) + (realNat(n):R) * H ≤ ((1:R) + H) ^ n
proof
  induction n
  | zero =>
    rw [realPow.zero, realNat.zero]
    have zeroProduct : (1:R) + (0:R) * H = (1:R) by ring
    rw zeroProduct
    exact realLe.reflexive with realOne.real
  | succ(k) using ih =>
    have countReal := realNat.real at k
    have baseReal : (1:R) + H ∈ ℝ by typecheck
    have one := realZero.lt_one
    have baseNonnegative : 0 ≤ (1:R) + H by linarith
    have power := realPow.real at (1:R) + H, k
    have scaled := realLe.mul_left_of_nonnegative with ih, baseNonnegative
    have countNonnegative := realNat.nonnegative at k
    have square := realMulValue.square_nonnegative at H
    have cross := realMulValue.nonnegative with countNonnegative, square
    have expand : ((1:R) + H) * ((1:R) + (realNat(k):R) * H) =
        (1:R) + ((realNat(k):R) + 1) * H + (realNat(k):R) * (H * H) by ring
    rw [realPow.succ, realNat.succ_add]
    rw expand at scaled
    linarith

-- For 0 ≤ A < 1 the powers decay at least like 1/(1 + n(1 - A)).
theorem realPow.decay_bound (A : Real, n : Nat)
  requires nonnegative : 0 ≤ A
  requires below : A < 1
  shows A ^ n * ((1:R) + (realNat(n):R) * ((1:R) - A)) ≤ 1
proof
  have gapReal := realSubValue.real at (1:R), A with realOne.real, assumption
  have gapPositive := realLess.iff_pos_sub.forward with below
  have gapNonnegative := realLess.implies_le with gapPositive
  have countReal := realNat.real at n
  have power := realPow.real at A, n
  have powerNonnegative := realPow.nonnegative at A, n with nonnegative
  have bernoulli := realPow.bernoulli at (1:R) - A, n with gapNonnegative
  have baseReal : (1:R) + ((1:R) - A) ∈ ℝ by typecheck
  have basePower := realPow.real at (1:R) + ((1:R) - A), n
  have factorReal : (1:R) + (realNat(n):R) * ((1:R) - A) ∈ ℝ by typecheck
  have scaled := realLe.mul_left_of_nonnegative with bernoulli, powerNonnegative
  have combined := realPow.mul at A, (1:R) + ((1:R) - A), n
  have productReal : A * ((1:R) + ((1:R) - A)) ∈ ℝ by typecheck
  have one := realZero.lt_one
  have baseNonnegative : 0 ≤ (1:R) + ((1:R) - A) by linarith
  have productNonnegative := realMulValue.nonnegative with
    nonnegative, baseNonnegative
  have productBounded : A * ((1:R) + ((1:R) - A)) ≤ 1
    have square := realMulValue.square_nonnegative at (1:R) - A
    have expand : A * ((1:R) + ((1:R) - A)) =
        (1:R) - ((1:R) - A) * ((1:R) - A) by ring
    rw expand
    linarith
  have small := realPow.le_one at A * ((1:R) + ((1:R) - A)), n with
    productNonnegative, productBounded
  rw combined at small
  exact realLe.transitive with scaled, small

-- Powers of a number in [0, 1) eventually fall below every positive bound.
theorem realPow.eventually_small (A E : Real)
  requires nonnegative : 0 ≤ A
  requires below : A < 1
  requires positive : 0 < E
  shows ∃ N ∈ ℕ, ∀ n ∈ ℕ, N ≤ n → A ^ n < E
proof
  have gapReal := realSubValue.real at (1:R), A with realOne.real, assumption
  have gapPositive := realLess.iff_pos_sub.forward with below
  have gapNonnegative := realLess.implies_le with gapPositive
  have scaleReal := realMulValue.real at E, (1:R) - A with assumption, gapReal
  have scalePositive := realMulValue.positive with positive, gapPositive
  have scaleNonzero : E * ((1:R) - A) ≠ 0
    intro zero
    rw zero at scalePositive
    linarith
  have thresholdReal := realInvValue.real at E * ((1:R) - A)
  have law := realInvValue.mul_left with scaleNonzero
  obtain ⟨N, typedN, above⟩ :=
    Real.below_natural at realInvValue(E * ((1:R) - A))
  have naturalN := nat.omega.forward with typedN
  witness N
  constructor
  . exact naturalN
  . intros n naturalIndex after
    have typedIndex : Nat(n) by exact nat.omega.backward at n with naturalIndex
    have countReal := realNat.real at n
    have monotone := realNat.monotone with after
    have aboveIndex : realInvValue(E * ((1:R) - A)) < realNat(n) by linarith
    have multiplied := realLess.mul_right_of_positive with
      aboveIndex, scalePositive
    rw law at multiplied
    have power := realPow.real at A, n
    have bound := realPow.decay_bound at A, n with nonnegative, below
    have factorReal : (1:R) + (realNat(n):R) * ((1:R) - A) ∈ ℝ by typecheck
    have countNonnegative := realNat.nonnegative at n
    have cross := realMulValue.nonnegative with countNonnegative, gapNonnegative
    have one := realZero.lt_one
    have factorNonnegative : 0 ≤ (1:R) + (realNat(n):R) * ((1:R) - A) by
      linarith
    by_contra notLess
    have reversed := realLe.of_not_less with notLess
    have scaledReverse := realLe.mul_right_of_nonnegative with
      reversed, factorNonnegative
    have expand : E * ((1:R) + (realNat(n):R) * ((1:R) - A)) =
        E + (realNat(n):R) * (E * ((1:R) - A)) by ring
    rw expand at scaledReverse
    linarith

-- The geometric sequence X^n tends to zero when |X| < 1.
theorem geometricTerms.converges_zero (X : Real)
  requires small : abs(X) < 1
  shows ConvergesTo(geometricTerms(X), (0:R))
proof
  have sequence := geometricTerms.sequence at X
  have absReal := abs.real at X
  have absNonnegative := abs.nonnegative at X
  apply ConvergesTo.zero_of_abs_bounds with sequence
  intros E realE positive
  have typedE : Real(E) by exact realE
  obtain ⟨N, naturalN, eventually⟩ := realPow.eventually_small at
    abs(X), E with absNonnegative, small, positive
  witness N
  constructor
  . exact naturalN
  . intros n naturalIndex after
    have typedN : Nat(n) by exact nat.omega.backward at n with naturalIndex
    have value := geometricTerms.value at X, n with *
    have absolute := realPow.abs at X, n
    have bound := eventually with n, naturalIndex, after
    rw [value, absolute]
    exact bound