-- The reciprocals of the triangular numbers sum to 2n/(n+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.TriangularSum
import FiniteSum
import RootTwo

-- A quotient is determined by multiplying through its denominator.
theorem realDiv.of_mul_right (X Y D : Real)
  requires nonzero : D ≠ 0
  requires product : X * D = Y
  shows X = Y / D
proof
  have inverse := realInvValue.mul_right with *
  have quotient := realDiv.eq_mul_inv at Y, D with *
  rw quotient
  linear_combination realInvValue(D) * product - X * inverse

theorem realNat.succ_add (k : Nat)
  shows realNat(succ(k)) = (realNat(k):R) + 1
proof
  rw realNat.succ
  unfold realAddValue
  unfold call2
  rfl

theorem realNat.successor_positive (k : Nat)
  shows 0 < (realNat(k):R) + 1
proof
  have nonnegative := realNat.nonnegative with *
  have one := realZero.lt_one
  linarith

theorem realNat.successor_nonzero (k : Nat)
  shows (realNat(k):R) + 1 ≠ 0
proof
  have positive := realNat.successor_positive with *
  intro zero
  rw zero at positive
  exact realLess.not_self with *, positive

-- The k-th term is 1/T_{k+1} = 2/((k+1)(k+2)).
defgraph triangularReciprocals (k ∈ ℕ) : ℝ :=
  (2:R) / ((realNat(k) + 1) * (realNat(k) + 2))
proof
  have typedK := nat.omega.backward with member
  exact realDiv.real at (2:R), (realNat(k) + 1) * (realNat(k) + 2) with *

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

-- The sequence whose adjacent differences are the reciprocals: -2/(k+1).
defgraph triangularTails (k ∈ ℕ) : ℝ := -((2:R) / (realNat(k) + 1))
proof
  have typedK := nat.omega.backward with member
  have quotient := realDiv.real at (2:R), (realNat(k):R) + 1 with *
  exact realNegValue.real with quotient

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

-- 2/(k+1) - 2/(k+2) = 2/((k+1)(k+2)).
theorem triangularReciprocals.difference (k : Nat)
  shows triangularReciprocals[k] = adjacentDifferences(triangularTails)[k]
proof
  have sequence := triangularTails.sequence
  have term := triangularReciprocals.value at k with *
  have difference := adjacentDifferences.value with sequence
  have tailK := triangularTails.value at k with *
  have tailNext := triangularTails.value at succ(k) with *
  rw [term, difference, tailNext, tailK]
  have next := realNat.succ_add with *
  rw next
  have firstNonzero := realNat.successor_nonzero with *
  have secondPositive : 0 < (realNat(k):R) + 2
    have nonnegative := realNat.nonnegative with *
    have one := realZero.lt_one
    linarith
  have secondNonzero : (realNat(k):R) + 2 ≠ 0
    intro zero
    rw zero at secondPositive
    exact realLess.not_self with *, secondPositive
  have productPositive : 0 < (realNat(k) + 1) * (realNat(k) + 2)
    have firstPositive := realNat.successor_positive with *
    exact realMulValue.positive with *
  have productNonzero : (realNat(k) + 1) * (realNat(k) + 2) ≠ 0
    intro zero
    rw zero at productPositive
    exact realLess.not_self with *, productPositive
  have shifted : (realNat(k):R) + 1 + 1 = (realNat(k):R) + 2 by ring
  rw shifted
  have firstQuotient := realDiv.mul_right at (2:R), (realNat(k):R) + 1 with *
  have secondQuotient := realDiv.mul_right at (2:R), (realNat(k):R) + 2 with *
  have product : (-((2:R) / (realNat(k) + 2)) - -((2:R) / (realNat(k) + 1))) * ((realNat(k) + 1) * (realNat(k) + 2)) = (2:R)
    linear_combination (realNat(k) + 2) * firstQuotient - (realNat(k) + 1) * secondQuotient
  have quotient := realDiv.of_mul_right with product
  symm at quotient
  exact quotient

theorem triangularTails.zero
  shows triangularTails[0] = -(2:R)
proof
  have value := triangularTails.value at 0 with *
  rw [value, realNat.zero]
  have sum : (0:R) + 1 = (1:R) by ring
  rw sum
  have oneNonzero : (1:R) ≠ 0
    intro zero
    have one := realZero.lt_one
    rw zero at one
    exact realLess.not_self with *, one
  have quotient := realDiv.mul_right at (2:R), (1:R) with *
  have simplified : (2:R) / (1:R) = (2:R) by linear_combination quotient
  rw simplified
  rfl

-- Σ_{k=1}^{n} 1/T_k = 2n/(n+1).
theorem realPartialSum.triangular_reciprocals (n : Nat)
  shows ∑ i < n, triangularReciprocals(i) = (2:R) * (realNat(n):R) / (realNat(n) + 1)
proof
  have sequence := triangularReciprocals.sequence
  have tails := triangularTails.sequence
  have differences := adjacentDifferences.sequence with tails
  have agree : ∀ k ∈ ℕ, k < n → triangularReciprocals[k] = adjacentDifferences(triangularTails)[k]
    intros k kMember less
    have kNat := nat.omega.backward with kMember
    exact triangularReciprocals.difference at k with *
  have same := realPartialSum.eq_prefix at triangularReciprocals, adjacentDifferences(triangularTails), n with *
  have telescoped := realPartialSum.telescoping at triangularTails, n with *
  have tailN := triangularTails.value at n with *
  have tailZero := triangularTails.zero
  rw [same, telescoped, tailN, tailZero]
  have nonzero := realNat.successor_nonzero at n with *
  have quotient := realDiv.mul_right at (2:R), (realNat(n):R) + 1 with *
  have product : (-((2:R) / (realNat(n) + 1)) - -(2:R)) * (realNat(n) + 1) = (2:R) * (realNat(n):R)
    have expand : (-((2:R) / ((realNat(n):R) + 1)) - -(2:R)) * ((realNat(n):R) + 1) = (2:R) * ((realNat(n):R) + 1) - ((2:R) / ((realNat(n):R) + 1)) * ((realNat(n):R) + 1)
      ring
    rw [expand, quotient]
    ring
  exact realDiv.of_mul_right at -((2:R) / (realNat(n) + 1)) - -(2:R), (2:R) * (realNat(n):R), (realNat(n):R) + 1 with *