import Mathlib.Algebra.Field.Basic
import Mathlib.Tactic.FieldSimp
import Mathlib.Tactic.Ring

/-!
# Shanks's second-order sequence transformation

Daniel Shanks, Non-linear Transformations of Divergent and Slowly Convergent
Sequences, Journal of Mathematics and Physics 34 (1955), 1–42.
DOI: 10.1002/sapm19553411.

Scope: equation (2) at order k = 2 and the two-mode instance of the
finite-transient exactness argument in equations (17)–(23). We use a window
starting at n instead of Shanks's centered index. The determinants are
written explicitly, including all the difference rows. Additional lemmas
factor the denominator and establish exactness at every sequence index.

No claim about general-order transforms, floating-point stability, arbitrary
convergence, or Shanks's rules for singular quotients is made.
Lean / Mathlib 4.33.1. Run: lake env lean -DwarningAsError=true Shanks.lean
License: MIT. See the website's bounded prior-work search before claiming novelty.
-/

namespace PaperToLean.Shanks

variable {K : Type*} [Field K]

/-- The ordinary 3 by 3 determinant, expanded to six products. -/
def det3 (a b c d e f g h i : K) : K :=
  a * (e * i - f * h) - b * (d * i - f * g) + c * (d * h - e * g)

/-- Bottom rows are successive first differences from a five-sample window. -/
def denominator (s0 s1 s2 s3 s4 : K) : K :=
  det3 1 1 1 (s1 - s0) (s2 - s1) (s3 - s2)
    (s2 - s1) (s3 - s2) (s4 - s3)

def numerator (s0 s1 s2 s3 s4 : K) : K :=
  det3 s0 s1 s2 (s1 - s0) (s2 - s1) (s3 - s2)
    (s2 - s1) (s3 - s2) (s4 - s3)

/-- Shanks's order-two determinant quotient, on nonsingular windows. -/
def transform (s0 s1 s2 s3 s4 : K) : K :=
  numerator s0 s1 s2 s3 s4 / denominator s0 s1 s2 s3 s4

/-- Samples of a target plus two independent geometric error modes. -/
def sample (L a b r q : K) (n : Nat) : K := L + a * r ^ n + b * q ^ n

/-- Explicit factorization identifies all degeneracies in this window. -/
theorem denominator_factor (L a b r q : K) :
    denominator (L + a + b) (L + a*r + b*q) (L + a*r^2 + b*q^2)
      (L + a*r^3 + b*q^3) (L + a*r^4 + b*q^4) =
      a * b * (r - 1)^2 * (q - 1)^2 * (r - q)^2 := by
  unfold denominator det3
  ring

/-- Cancellation in the numerator is an identity, before any division. -/
theorem numerator_factor (L a b r q : K) :
    numerator (L + a + b) (L + a*r + b*q) (L + a*r^2 + b*q^2)
      (L + a*r^3 + b*q^3) (L + a*r^4 + b*q^4) =
      L * (a * b * (r - 1)^2 * (q - 1)^2 * (r - q)^2) := by
  unfold numerator det3
  ring

/-- Five samples recover L exactly, even for growing or oscillating modes. -/
theorem two_mode_exact (L a b r q : K)
    (ha : a ≠ 0) (hb : b ≠ 0) (hr : r ≠ 1) (hq : q ≠ 1) (hrq : r ≠ q) :
    transform (L + a + b) (L + a*r + b*q) (L + a*r^2 + b*q^2)
      (L + a*r^3 + b*q^3) (L + a*r^4 + b*q^4) = L := by
  have hd : a * b * (r - 1)^2 * (q - 1)^2 * (r - q)^2 ≠ 0 :=
    mul_ne_zero (mul_ne_zero (mul_ne_zero (mul_ne_zero ha hb)
      (pow_ne_zero _ (sub_ne_zero.mpr hr))) (pow_ne_zero _ (sub_ne_zero.mpr hq)))
      (pow_ne_zero _ (sub_ne_zero.mpr hrq))
  unfold transform
  rw [numerator_factor, denominator_factor]
  exact mul_div_cancel_right₀ L hd

/-- A shift of the window only rescales the two error amplitudes. -/
theorem sample_shift (L a b r q : K) (n k : Nat) :
    sample L a b r q (n + k) = sample L (a * r^n) (b * q^n) r q k := by
  unfold sample
  simp only [pow_add]
  ring

/-- A complete sequence-level statement, for every natural-number index. -/
theorem sequence_exact (L a b r q : K)
    (ha : a ≠ 0) (hb : b ≠ 0) (hr0 : r ≠ 0) (hq0 : q ≠ 0)
    (hr1 : r ≠ 1) (hq1 : q ≠ 1) (hrq : r ≠ q) (n : Nat) :
    transform (sample L a b r q n) (sample L a b r q (n+1))
      (sample L a b r q (n+2)) (sample L a b r q (n+3))
      (sample L a b r q (n+4)) = L := by
  have h := two_mode_exact L (a*r^n) (b*q^n) r q
    (mul_ne_zero ha (pow_ne_zero _ hr0)) (mul_ne_zero hb (pow_ne_zero _ hq0))
    hr1 hq1 hrq
  simpa [sample, pow_add, mul_assoc] using h

/-- A single geometric mode makes the order-two denominator singular.
The raw quotient cannot replace Shanks's separate order-reduction rule. -/
theorem single_mode_singular (L a r : K) :
    denominator (L+a) (L+a*r) (L+a*r^2) (L+a*r^3) (L+a*r^4) = 0 := by
  simpa using denominator_factor L a 0 r 0

#print axioms denominator_factor
#print axioms numerator_factor
#print axioms two_mode_exact
#print axioms sample_shift
#print axioms sequence_exact
#print axioms single_mode_singular

end PaperToLean.Shanks