-- The geometric series: Σ X^k = 1/(1 - X) for |X| < 1
-- 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.GeometricSeries
import ArchimedeanLimits

-- The partial sums (1 - X^n)/(1 - X), read as a quotient of sequences.
theorem geometricSeries.partial_sums (X : Real, n : Nat)
  requires nonzero : (1:R) - X ≠ 0
  shows prefixSums(geometricTerms(X))[n] =
    divSequence(subSequence(constantSequence((1:R)), geometricTerms(X)),
      constantSequence((1:R) - X))[n]
proof
  have sequence := geometricTerms.sequence at X
  have ones := constantSequence.sequence at (1:R) with realOne.real
  have denominatorReal := realSubValue.real at (1:R), X with
    realOne.real, assumption
  have denominators := constantSequence.sequence at (1:R) - X with
    denominatorReal
  have differences := subSequence.sequence with ones, sequence
  have oneReal := RealSequence.value at constantSequence((1:R)), n with ones
  have termReal := RealSequence.value at geometricTerms(X), n with sequence
  have quotient := divSequence.value at
    subSequence(constantSequence((1:R)), geometricTerms(X)),
    constantSequence((1:R) - X), n with
      differences, denominators
  have difference := subSequence.value at
    constantSequence((1:R)), geometricTerms(X), n with
      ones, sequence, oneReal, termReal
  have one := constantSequence.value at (1:R), n with realOne.real, assumption
  have denominator := constantSequence.value at (1:R) - X, n with
    denominatorReal, assumption
  have term := geometricTerms.value at X, n with *
  have value := prefixSums.value at geometricTerms(X), n with sequence
  have formula := realPartialSum.geometric_quotient at X, n with nonzero
  rw [value, formula, quotient, difference, one, denominator, term]
  rfl

-- Σ_{k=0}^{∞} X^k = 1/(1 - X) for |X| < 1.
theorem geometricSeries.sums_to (X : Real)
  requires small : abs(X) < 1
  shows SumsTo(geometricTerms(X), (1:R) / ((1:R) - X))
proof
  have sequence := geometricTerms.sequence at X
  have prefixes := prefixSums.sequence with sequence
  have absReal := abs.real at X
  have upper := abs.upper at X
  have belowOne : X < 1 by linarith
  have denominatorPositive := realLess.iff_pos_sub.forward with belowOne
  have denominatorReal := realSubValue.real at (1:R), X with
    realOne.real, assumption
  have nonzero : (1:R) - X ≠ 0
    intro zero
    rw zero at denominatorPositive
    linarith
  have ones := constantSequence.converges at (1:R) with realOne.real
  have powers := geometricTerms.converges_zero with small
  have numerators := ConvergesTo.sub with ones, powers
  have denominators := constantSequence.converges at (1:R) - X with
    denominatorReal
  have limit := ConvergesTo.div with numerators, denominators, nonzero
  have zero := realSubValue.zero_right with realOne.real
  rw zero at limit
  have equal : ∀ n ∈ ℕ, prefixSums(geometricTerms(X))[n] =
      divSequence(subSequence(constantSequence((1:R)), geometricTerms(X)),
        constantSequence((1:R) - X))[n]
    intros n naturalN
    have typedN : Nat(n) by exact nat.omega.backward at n with naturalN
    exact geometricSeries.partial_sums at X, n with nonzero
  have converges := ConvergesTo.of_equal_terms with prefixes, limit, equal
  exact SumsTo.of_converges with sequence, converges

-- The classic instance: 1 + 1/2 + 1/4 + ⋯ = 2.
theorem geometricSeries.half_sums_to_two
  shows SumsTo(geometricTerms(half((1:R))), (2:R))
proof
  have halfReal := half.real with realOne.real
  have halfPositive := half.positive with realZero.lt_one
  have halfBelow := half.lt_self with realOne.real, realZero.lt_one
  have halfNonnegative := realLess.implies_le with halfPositive
  have absolute := abs.of_nonnegative with halfNonnegative
  have small : abs(half((1:R))) < 1 by rwa absolute
  have sums := geometricSeries.sums_to with small
  have doubled := half.add_self with realOne.real
  have gapReal := realSubValue.real at (1:R), half((1:R)) with
    realOne.real, halfReal
  have gapNonzero : (1:R) - half((1:R)) ≠ 0
    intro zero
    linarith
  have product : (2:R) * ((1:R) - half((1:R))) = (1:R)
    have expand : (2:R) * ((1:R) - half((1:R))) =
        (2:R) - (half((1:R)) + half((1:R))) by ring
    rw [expand, doubled]
    ring
  have quotient := realDiv.of_mul_right with gapNonzero, product
  rw ← quotient at sums
  exact sums