import Mathlib.Algebra.Ring.Basic
import Mathlib.Tactic.NoncommRing

/-! Schulz (1933): algebraic error squaring for inverse iteration.
Works in any unital ring, including square matrices. No commutativity,
invertibility, convergence, or floating-point stability is assumed. -/
namespace PaperToLean.Schulz
variable {R : Type*} [Ring R]

def step (a x : R) : R := x * (2 - a * x)
def iterate (a x : R) : ℕ → R
  | 0 => x
  | n + 1 => step a (iterate a x n)

theorem right_error_square (a x : R) :
    1 - a * step a x = (1 - a * x) ^ 2 := by
  simp only [step, ← one_add_one_eq_two, mul_sub, mul_add, mul_one]
  noncomm_ring

theorem left_error_square (a x : R) :
    1 - step a x * a = (1 - x * a) ^ 2 := by
  simp only [step, ← one_add_one_eq_two, mul_sub, mul_add, mul_one,
    sub_mul, add_mul]
  noncomm_ring

theorem right_error_exact (a x : R) (n : ℕ) :
    1 - a * iterate a x n = (1 - a * x) ^ (2 ^ n) := by
  induction n with
  | zero => simp [iterate]
  | succ n ih =>
    rw [iterate, right_error_square, ih, ← pow_mul, pow_succ]

theorem left_error_exact (a x : R) (n : ℕ) :
    1 - iterate a x n * a = (1 - x * a) ^ (2 ^ n) := by
  induction n with
  | zero => simp [iterate]
  | succ n ih =>
    rw [iterate, left_error_square, ih, ← pow_mul, pow_succ]

theorem nilpotent_right_exact (a x : R) (n : ℕ)
    (h : (1 - a * x) ^ (2 ^ n) = 0) : a * iterate a x n = 1 := by
  have he := right_error_exact a x n
  rw [h, sub_eq_zero] at he
  exact he.symm

theorem inverse_fixed (a x : R) (h : a * x = 1) : step a x = x := by
  simp [step, h, ← one_add_one_eq_two]

#print axioms right_error_square
#print axioms left_error_square
#print axioms right_error_exact
#print axioms left_error_exact
#print axioms nilpotent_right_exact
#print axioms inverse_fixed
end PaperToLean.Schulz