-- The harmonic series diverges: doubling the number of terms adds at least a
-- half, so the partial sums pass every real number and have no limit. The
-- Archimedean property of the reals is proved on the way
-- 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.HarmonicSeries
import TriangularSum
import RationalEmbedding

-- Every real lies below some natural number.
theorem Real.below_natural (M : Real)
  shows ∃ n, Nat(n) ∧ M < realNat(n)
proof
  have realM : Real(M) by typecheck
  have cut := Real.characterization.forward with realM
  obtain ⟨u, uInQ, uPositive, uMissing⟩ := DedekindCut.exists_positive_upper with cut
  have uRational : Rational(u) by typecheck
  have vRational : Rational(u + rationalOne) by typecheck
  have uLessV : u < u + rationalOne
    have zeroLtOne := Grace.Cumulus.Rat.OrderedOperations.rationalZero.lt_one
    have step := rationalLess.add_left at rationalZero, rationalOne, u with zeroLtOne
    have zero := rationalAddValue.zero_right at u
    rwa zero at step
  have cutReal := rationalCut.real with vRational
  have subset := DedekindCut.subset_rationals with cut
  have below : RealLess(M, rationalCut(u + rationalOne))
    constructor
    . exact realM
    . constructor
      . exact cutReal
      . constructor
        . intros a ha
          have aLess := DedekindCut.lt_missing with cut, ha, uInQ, uMissing
          have aInQ := subset with a, ha
          have aLessV := rationalLess.transitive with aLess, uLessV
          apply rationalCut.membership.backward
          exact ⟨aInQ, aLessV⟩
        . intro equal
          have uInCut : u ∈ rationalCut(u + rationalOne)
            apply rationalCut.membership.backward
            exact ⟨uInQ, uLessV⟩
          rw equal at uMissing
          exact uMissing with uInCut
  have first := realLess.characterization.backward with below
  obtain ⟨n, nNat, vBelow⟩ := Rational.below_natural at u + rationalOne
  have class := rationalNat.class at n
  have vBelowNat : u + rationalOne < rationalNat(n) by rwa class
  have mono := rationalCut.strict_mono with vBelowNat
  have natural := rationalCut.natural at n
  rw natural at mono
  have second := realLess.characterization.backward with mono
  witness n
  exact realLess.transitive with first, second

-- Reciprocation reverses order on the positive reals.
theorem realInvValue.antitone_le (A B : Real)
  requires positive : 0 < A
  requires le : A ≤ B
  shows realInvValue(B) ≤ realInvValue(A)
proof
  have cases := realLess.trichotomy at A, B
  cases cases with less | rest
  . have strict := realInvValue.antitone_positive with positive, less
    exact realLess.implies_le with strict
  . cases rest with equal | greater
    . rw equal
      linarith
    . exfalso
      linarith

theorem realNat.monotone (m n : Nat)
  requires bound : m ≤ n
  shows realNat(m) ≤ realNat(n)
proof
  have less : m < succ(n) by arith
  have inSucc := membership.of_Less with less
  have split := bootstrap.succ.membership.forward with inSucc
  cases split with below | same
  . have strict := realNat.lt_of_lt with (Less.of_membership with below)
    exact realLess.implies_le with strict
  . rw same
    linarith

-- A constant sequence sums to the count times the constant.
theorem realPartialSum.constant (C : Real) (m : Nat)
  shows ∑ i < m, constantSequence(C)(i) = (realNat(m):R) * C
proof
  induction m
  | zero =>
    rw [realPartialSum.zero, realNat.zero]
    ring
  | succ(k) using ih =>
    have count := realNat.real at k
    have value := constantSequence.value at C, k
    rw [realPartialSum.succ, value, realNat.succ]
    linear_combination ih

-- The harmonic terms 1, 1/2, 1/3, …
defgraph harmonicTerms (n ∈ ℕ) : ℝ := realInvValue((realNat(n):R) + (1:R))
proof
  have typedN := nat.omega.backward with member
  have count := realNat.real with typedN
  exact realInvValue.real at (realNat(n):R) + (1:R)

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

theorem harmonicTerms.positive (n : Nat)
  shows (0:R) < harmonicTerms(n)
proof
  have value := harmonicTerms.value at n
  rw value
  exact realInvValue.positive with (realNat.successor_positive at n)

-- Each of the m terms after the first m is at least 1/(2m).
theorem harmonicTerms.tail_bound (m j : Nat)
  requires below : j < m
  shows realInvValue((2:R) * (realNat(m):R)) ≤ harmonicTerms(add(m, j))
proof
  have value := harmonicTerms.value at add(m, j)
  have split := realNat.add at m, j
  have succLe : succ(j) ≤ m by arith
  have mono := realNat.monotone with succLe
  have succValue := realNat.succ_add at j
  rw succValue at mono
  have jNonneg := realNat.nonnegative at j
  have argPositive : (0:R) < (realNat(add(m, j)):R) + (1:R)
    rw split
    linarith
  have two : (2:R) * (realNat(m):R) = (realNat(m):R) + (realNat(m):R) by ring
  have argLe : (realNat(add(m, j)):R) + (1:R) ≤ (2:R) * (realNat(m):R)
    rw [split, two]
    linarith
  rw value
  exact realInvValue.antitone_le with argPositive, argLe

-- Doubling the number of terms adds at least one half to the sum.
theorem harmonicSeries.doubling (m : Nat)
  requires positive : 0 < m
  shows (2:R) * (∑ i < m, harmonicTerms(i):R) + (1:R) ≤ (2:R) * (realPartialSum(harmonicTerms, add(m, m)):R)
proof
  have sequence := harmonicTerms.sequence
  have split := realPartialSum.split_at with sequence
  have shifted := shiftedRealSequence.sequence with sequence
  have mReal := realNat.real at m
  have mPositive : (0:R) < realNat(m)
    have strict := realNat.lt_of_lt with positive
    rwa realNat.zero at strict
  have two : (2:R) * (realNat(m):R) = (realNat(m):R) + (realNat(m):R) by ring
  have twoM : (0:R) < (2:R) * (realNat(m):R)
    rw two
    linarith
  have invReal := realInvValue.real at (2:R) * (realNat(m):R)
  have constant := constantSequence.sequence at realInvValue((2:R) * (realNat(m):R))
  have bounded : ∀ k ∈ ℕ, k < m → (constantSequence(realInvValue((2:R) * (realNat(m):R)))(k):R) ≤ (shiftedRealSequence(harmonicTerms, m)(k):R)
    intros k hk hkm
    have kNat : Nat(k) by typecheck
    have cValue := constantSequence.value at realInvValue((2:R) * (realNat(m):R)), k
    have sValue := shiftedRealSequence.value at harmonicTerms, m, k with sequence
    rw [cValue, sValue]
    exact harmonicTerms.tail_bound with hkm
  have mono := realPartialSum.monotone_prefix with constant, shifted, bounded
  have constSum := realPartialSum.constant at realInvValue((2:R) * (realNat(m):R)), m
  rw constSum at mono
  have nonzero : (2:R) * (realNat(m):R) ≠ (0:R)
    intro zero
    rw zero at twoM
    linarith
  have cancel := realInvValue.mul_right with nonzero
  have expand : (2:R) * ((realNat(m):R) * realInvValue((2:R) * (realNat(m):R))) = (1:R)
    have rearrange : (2:R) * ((realNat(m):R) * realInvValue((2:R) * (realNat(m):R))) = ((2:R) * (realNat(m):R)) * realInvValue((2:R) * (realNat(m):R))
      ring
    rw [rearrange, cancel]
    rfl
  have doubleTerm : (2:R) * ((realNat(m):R) * realInvValue((2:R) * (realNat(m):R))) = (realNat(m):R) * realInvValue((2:R) * (realNat(m):R)) + (realNat(m):R) * realInvValue((2:R) * (realNat(m):R))
    ring
  rw doubleTerm at expand
  have doubleSum : (2:R) * (∑ i < m, harmonicTerms(i):R) = (∑ i < m, harmonicTerms(i):R) + (∑ i < m, harmonicTerms(i):R)
    ring
  have doubleTotal : (2:R) * ((∑ i < m, harmonicTerms(i):R) + (∑ i < m, shiftedRealSequence(harmonicTerms, m)(i):R)) = (∑ i < m, harmonicTerms(i):R) + (∑ i < m, harmonicTerms(i):R) + (∑ i < m, shiftedRealSequence(harmonicTerms, m)(i):R) + (∑ i < m, shiftedRealSequence(harmonicTerms, m)(i):R)
    ring
  rw [split, doubleSum, doubleTotal]
  linarith

-- Some partial sum reaches 1 + k/2, for every k.
theorem harmonicSeries.lower_bound (k : Nat)
  shows ∃ m, Nat(m) ∧ 0 < m ∧ (2:R) + (realNat(k):R) ≤ (2:R) * (∑ l < m, harmonicTerms(l):R)
proof
  induction k
  | zero =>
    have sequence := harmonicTerms.sequence
    have step := realPartialSum.succ at harmonicTerms, 0
    have zeroSum := realPartialSum.zero at harmonicTerms
    have value := harmonicTerms.value at 0
    have oneNonzero : (1:R) ≠ (0:R)
      intro equal
      linarith
    have cancel := realInvValue.mul_right at (1:R) with oneNonzero
    have one := realMulValue.one_left at realInvValue((1:R))
    rw one at cancel
    have shift := realAddValue.zero_left at (1:R)
    rw [realNat.zero, shift, cancel] at value
    have sumReal := realPartialSum.real at harmonicTerms, 0 with sequence
    have termReal := RealSequence.value at harmonicTerms, 0 with sequence
    have total : realPartialSum(harmonicTerms, succ(0)) = (1:R)
      rw step
      linarith
    have positive : 0 < succ(0) by arith
    have twoReal : (2:R) ∈ ℝ by typecheck
    have twoOne := realMulValue.one_right at (2:R)
    have twoZero := realAddValue.zero_right at (2:R)
    witness succ(0)
    constructor
    . exact positive
    . rw [total, realNat.zero, twoOne, twoZero]
      exact realLe.reflexive at (2:R)
  | succ(j) using ih =>
    obtain ⟨m, mNat, mPositive, bound⟩ := ih
    have doubling := harmonicSeries.doubling at m with mPositive
    have jReal := realNat.real at j
    have sequence := harmonicTerms.sequence
    have sumReal := realPartialSum.real at harmonicTerms, m with sequence
    have doubleReal := realPartialSum.real at harmonicTerms, add(m, m) with sequence
    have positive : 0 < add(m, m) by arith
    witness add(m, m)
    constructor
    . exact positive
    . rw realNat.succ
      linarith

-- The partial sums pass every real number.
theorem harmonicSeries.unbounded (M : Real)
  shows ∃ N, Nat(N) ∧ M < ∑ i < N, harmonicTerms(i)
proof
  obtain ⟨k, kNat, above⟩ := Real.below_natural at (2:R) * M
  obtain ⟨m, mNat, mPositive, bound⟩ := harmonicSeries.lower_bound at k
  have doubleM : (2:R) * M = M + M by ring
  have doubleSum : (2:R) * (∑ i < m, harmonicTerms(i):R) = (∑ i < m, harmonicTerms(i):R) + (∑ i < m, harmonicTerms(i):R)
    ring
  rw doubleM at above
  rw doubleSum at bound
  witness m
  linarith

-- The harmonic series has no limit.
theorem harmonicSeries.diverges
  shows ¬∃ L ∈ ℝ, ConvergesTo(prefixSums(harmonicTerms), L)
proof
  intro converges
  obtain ⟨L, LReal, limit⟩ := converges
  have sequence := harmonicTerms.sequence
  obtain ⟨N, NIn, close⟩ := ConvergesTo.eventually with limit, realZero.lt_one
  have NNat : Nat(N) by typecheck
  have shiftReal : L + (1:R) ∈ ℝ by typecheck
  obtain ⟨m, mNat, above⟩ := harmonicSeries.unbounded at L + (1:R)
  have nNat : Nat(add(m, N)) by typecheck
  have nIn := nat.omega.forward with nNat
  have le : N ≤ add(m, N) by arith
  have inBall := close with nIn, le
  have parts := openBall.membership.forward with inBall
  have value := prefixSums.value at harmonicTerms, add(m, N) with *
  have split := realPartialSum.split_at at harmonicTerms, m, N with sequence
  have shifted := shiftedRealSequence.sequence at harmonicTerms, m with sequence
  have sumReal := realPartialSum.real at harmonicTerms, m with sequence
  have tailReal := realPartialSum.real at shiftedRealSequence(harmonicTerms, m), N with shifted
  have tailNonneg : (0:R) ≤ (∑ i < N, shiftedRealSequence(harmonicTerms, m)(i):R)
    have constant := constantSequence.sequence at (0:R)
    have bounded : ∀ k ∈ ℕ, k < N → (constantSequence((0:R))(k):R) ≤ (shiftedRealSequence(harmonicTerms, m)(k):R)
      intros k hk hkN
      have kNat : Nat(k) by typecheck
      have cValue := constantSequence.value at (0:R), k
      have sValue := shiftedRealSequence.value at harmonicTerms, m, k with sequence
      rw [cValue, sValue]
      exact realLess.implies_le with (harmonicTerms.positive at add(m, k))
    have mono := realPartialSum.monotone_prefix with constant, shifted, bounded
    have constSum := realPartialSum.constant at (0:R), N
    have zero := realMulValue.zero_right at realNat(N)
    rwa [constSum, zero] at mono
  have distance := parts.right
  rw [value, split] at distance
  have same : dist((∑ i < m, harmonicTerms(i):R) + (∑ i < N, shiftedRealSequence(harmonicTerms, m)(i):R), L) = abs((∑ i < m, harmonicTerms(i):R) + (∑ i < N, shiftedRealSequence(harmonicTerms, m)(i):R) - L)
    unfold dist
    rfl
  rw same at distance
  have upper := abs.upper at (∑ i < m, harmonicTerms(i):R) + (∑ i < N, shiftedRealSequence(harmonicTerms, m)(i):R) - L
  linarith