import Mathlib.Algebra.Field.Basic
import Mathlib.Tactic.FieldSimp
import Mathlib.Tactic.Ring
import Mathlib.Topology.Instances.Real.Lemmas
import Mathlib.Topology.Algebra.Field

/-!
# Aitken's delta-squared transformation

A. C. Aitken, On Bernoulli's Numerical Solution of Algebraic Equations,
Proc. Roy. Soc. Edinburgh 46, 289–305. Usually cited as 1926; Cambridge's
publication record dates the volume 1927. DOI: 10.1017/S0370164600022070.

Selected results: exact removal of one geometric error mode, the exact
residual for two modes, affine equivariance, and an explicit singular example;
then a modern real-sequence proof of asymptotic acceleration and preservation
of the limit. This is not a formalization of the whole paper or a floating-point
error analysis. Algebraic statements use arbitrary fields; limits use real numbers.

Lean / Mathlib 4.33.1. Run in the supplied Lake project:
  lake env lean -DwarningAsError=true Aitken.lean
License: MIT. See the website's prior-work note for the bounded novelty search.
-/

namespace PaperToLean.Aitken

variable {K : Type*} [Field K]

/-- The second forward difference of three successive samples. -/
def delta2 (x y z : K) : K := z - 2 * y + x

/-- Defined algebraically for every input, but its numerical interpretation
requires a nonzero second difference. Lean assigns zero to division by zero. -/
def transform (x y z : K) : K := x - (y - x) ^ 2 / delta2 x y z

theorem determinant_form (x y z : K) (h : delta2 x y z ≠ 0) :
    transform x y z = (x * z - y ^ 2) / delta2 x y z := by
  unfold transform
  field_simp
  simp only [delta2]
  ring

/-- The denominator is explicit, so exactness cannot hide a singular case. -/
theorem geometric_denominator (L a r : K) :
    delta2 (L + a) (L + a * r) (L + a * r ^ 2) = a * (r - 1) ^ 2 := by
  unfold delta2
  ring

/-- Three samples determine the target L when there is one geometric error. -/
theorem geometric_exact (L a r : K) (ha : a ≠ 0) (hr : r ≠ 1) :
    transform (L + a) (L + a * r) (L + a * r ^ 2) = L := by
  have hd : a * (r - 1) ^ 2 ≠ 0 := mul_ne_zero ha (pow_ne_zero _ (sub_ne_zero.mpr hr))
  rw [determinant_form _ _ _ (by rw [geometric_denominator]; exact hd)]
  rw [geometric_denominator]
  apply (div_eq_iff hd).2
  ring

/-- The same result holds at every index of an infinite geometric model. -/
theorem geometric_sequence_exact (L a r : K) (n : Nat)
    (ha : a ≠ 0) (hr0 : r ≠ 0) (hr1 : r ≠ 1) :
    transform (L + a * r ^ n) (L + a * r ^ (n + 1))
      (L + a * r ^ (n + 2)) = L := by
  simpa [pow_add, mul_assoc] using
    geometric_exact L (a * r ^ n) r (mul_ne_zero ha (pow_ne_zero _ hr0)) hr1

/-- Exact two-mode error, with no asymptotic or small-error approximation. -/
theorem two_mode_residual (L a b r q : K)
    (hd : a * (r - 1) ^ 2 + b * (q - 1) ^ 2 ≠ 0) :
    transform (L + a + b) (L + a * r + b * q)
      (L + a * r ^ 2 + b * q ^ 2) - L =
      a * b * (r - q) ^ 2 / (a * (r - 1) ^ 2 + b * (q - 1) ^ 2) := by
  have hden : delta2 (L + a + b) (L + a * r + b * q)
      (L + a * r ^ 2 + b * q ^ 2) = a * (r - 1) ^ 2 + b * (q - 1) ^ 2 := by
    unfold delta2
    ring
  unfold transform
  rw [hden]
  field_simp
  ring

/-- Changing units and origin commutes with the nonsingular transformation. -/
theorem affine_equivariant (x y z scale shift : K)
    (hs : scale ≠ 0) (hd : delta2 x y z ≠ 0) :
    transform (scale * x + shift) (scale * y + shift) (scale * z + shift) =
      scale * transform x y z + shift := by
  have hden : delta2 (scale * x + shift) (scale * y + shift) (scale * z + shift) =
      scale * delta2 x y z := by
    unfold delta2
    ring
  unfold transform
  rw [hden]
  field_simp
  ring

/-- A constant nonzero error has zero second difference; it is not removed.
This documents why r ≠ 1 is a mathematical assumption, not boilerplate. -/
theorem constant_error_is_not_removed (L a : K) (ha : a ≠ 0) :
    delta2 (L + a) (L + a) (L + a) = 0 ∧
      transform (L + a) (L + a) (L + a) ≠ L := by
  have hd : delta2 (L + a) (L + a) (L + a) = 0 := by unfold delta2; ring
  refine ⟨hd, ?_⟩
  simpa [transform, hd] using ha

/-- Normalize the transformed error by the current error. The denominator
on the right is a normalized second difference, not an assumed limit. -/
theorem normalized_error (L x y z : K)
    (hx : x ≠ 0) (hy : y ≠ 0) (hd : z - 2*y + x ≠ 0) :
    (transform (L+x) (L+y) (L+z) - L) / x =
      ((y/x) * (z/y) - (y/x)^2) / ((y/x) * (z/y) - 2*(y/x) + 1) := by
  have hden : delta2 (L+x) (L+y) (L+z) = z - 2*y + x := by unfold delta2; ring
  have hn : (y/x) * (z/y) - 2*(y/x) + 1 = (z - 2*y + x) / x := by
    field_simp [hx, hy]
  have hd' : x - y*2 + z ≠ 0 := by
    convert hd using 1
    ring
  unfold transform
  rw [hden, hn]
  field_simp [hx, hy, hd, hd']
  ring_nf
  field_simp [hd']
  ring

open Filter Topology

/-- If consecutive nonzero errors have ratio tending to r ≠ 1, Aitken's
new error is negligible compared with the old error. Only eventual
nonvanishing is required; early singular windows do not affect the limit. -/
theorem accelerates (s : Nat → ℝ) (L r : ℝ)
    (he : ∀ᶠ n in atTop, s n - L ≠ 0)
    (hr : Tendsto (fun n => (s (n+1) - L) / (s n - L)) atTop (𝓝 r))
    (hne : r ≠ 1) :
    Tendsto (fun n => (transform (s n) (s (n+1)) (s (n+2)) - L) / (s n - L))
      atTop (𝓝 0) := by
  let err := fun n => s n - L
  let ratio := fun n => err (n+1) / err n
  have hr' : Tendsto ratio atTop (𝓝 r) := hr
  have hnext : Tendsto (fun n => ratio (n+1)) atTop (𝓝 r) :=
    hr'.comp (tendsto_add_atTop_nat 1)
  have hden : Tendsto (fun n => ratio n * ratio (n+1) - 2*ratio n + 1)
      atTop (𝓝 (r*r - 2*r + 1)) :=
    ((hr'.mul hnext).sub (tendsto_const_nhds.mul hr')).add tendsto_const_nhds
  have hlim : r*r - 2*r + 1 ≠ 0 := by
    have hsq := pow_ne_zero 2 (sub_ne_zero.mpr hne)
    convert hsq using 1
    ring
  have hnonzero := hden.eventually_ne hlim
  have he_next : ∀ᶠ n in atTop, err (n+1) ≠ 0 :=
    (tendsto_add_atTop_nat 1).eventually he
  have hquot := ((hr'.mul hnext).sub (hr'.pow 2)).div hden hlim
  have hzero : (r*r - r^2) / (r*r - 2*r + 1) = 0 := by ring
  rw [hzero] at hquot
  apply hquot.congr'
  filter_upwards [he, he_next, hnonzero] with n hx hy hn
  have hnorm : ratio n * ratio (n+1) - 2*ratio n + 1 =
      (err (n+2) - 2*err (n+1) + err n) / err n := by
    dsimp [ratio]
    have hx' : err n ≠ 0 := hx
    have hy' : err (n+1) ≠ 0 := hy
    field_simp [hx', hy', Nat.add_assoc]
  have hd : err (n+2) - 2*err (n+1) + err n ≠ 0 := by
    intro h
    rw [hnorm, h, zero_div] at hn
    exact hn rfl
  have h := normalized_error L (err n) (err (n+1)) (err (n+2)) hx hy hd
  simpa [err, ratio, Nat.add_assoc] using h.symm

/-- With convergence of the original sequence, acceleration preserves its
limit as well as improving its asymptotic relative error. -/
theorem preserves_limit (s : Nat → ℝ) (L r : ℝ)
    (hs : Tendsto s atTop (𝓝 L))
    (he : ∀ᶠ n in atTop, s n - L ≠ 0)
    (hr : Tendsto (fun n => (s (n+1) - L) / (s n - L)) atTop (𝓝 r))
    (hne : r ≠ 1) :
    Tendsto (fun n => transform (s n) (s (n+1)) (s (n+2))) atTop (𝓝 L) := by
  have hacc := accelerates s L r he hr hne
  have hprod := hacc.mul (hs.sub (tendsto_const_nhds (x := L)))
  have hadd := hprod.add (tendsto_const_nhds (x := L))
  simp only [sub_self, mul_zero, zero_add] at hadd
  apply hadd.congr'
  filter_upwards [he] with n hn
  rw [div_mul_cancel₀ _ hn, sub_add_cancel]

#print axioms determinant_form
#print axioms geometric_denominator
#print axioms geometric_exact
#print axioms geometric_sequence_exact
#print axioms two_mode_residual
#print axioms affine_equivariant
#print axioms constant_error_is_not_removed
#print axioms normalized_error
#print axioms accelerates
#print axioms preserves_limit

end PaperToLean.Aitken