import Mathlib.Algebra.Field.Basic
import Mathlib.Data.Real.Basic
import Mathlib.Tactic.FieldSimp
import Mathlib.Tactic.Ring
import Mathlib.Tactic.Linarith
import Mathlib.Tactic.Positivity

/-! Steffensen (1933): a derivative-free root step for quadratic functions.
We prove its exact error, denominator safety on a positive branch, and
interval preservation. This is a specialization, not general convergence. -/
namespace PaperToLean.Steffensen
noncomputable section
def f (a x : ℝ) : ℝ := x ^ 2 - a
def slope (a x : ℝ) : ℝ := 2 * x + f a x
def step (a x : ℝ) : ℝ := x - f a x / slope a x
def rawStep (a x : ℝ) : ℝ := x - (f a x) ^ 2 / (f a (x + f a x) - f a x)

theorem divided_difference (a x : ℝ) :
    f a (x + f a x) - f a x = f a x * slope a x := by
  simp only [slope, f]
  ring

theorem raw_equivalence (a x : ℝ) (hd : slope a x ≠ 0) :
    rawStep a x = step a x := by
  unfold rawStep step
  rw [divided_difference]
  field_simp

theorem fixed_root (r : ℝ) : step (r ^ 2) r = r := by
  simp [step, f]

theorem exact_error (r x : ℝ) (hd : slope (r ^ 2) x ≠ 0) :
    step (r ^ 2) x - r = (x - r) ^ 2 * (1 + x + r) / slope (r ^ 2) x := by
  unfold step
  apply (eq_div_iff hd).mpr
  field_simp
  unfold slope f
  ring

theorem safe_denominator (r x : ℝ) (hr : 0 < r) (hx : r ≤ x) :
    0 < slope (r ^ 2) x := by
  unfold slope f
  nlinarith

theorem bracket (r x : ℝ) (hr : 0 < r) (hx : r ≤ x) :
    r ≤ step (r ^ 2) x ∧ step (r ^ 2) x ≤ x := by
  have hd := safe_denominator r x hr hx
  constructor
  · rw [← sub_nonneg, exact_error r x (ne_of_gt hd)]
    have hx0 : 0 < x := lt_of_lt_of_le hr hx
    positivity
  · unfold step
    have hf : 0 ≤ f (r ^ 2) x := by unfold f; nlinarith
    exact sub_le_self _ (div_nonneg hf (le_of_lt hd))

theorem quadratic_bound (r x R : ℝ) (hr : 0 < r) (hx : r ≤ x) (hR : x ≤ R) :
    step (r ^ 2) x - r ≤ (1 + R + r) / (2 * r) * (x - r) ^ 2 := by
  have hd := safe_denominator r x hr hx
  rw [exact_error r x (ne_of_gt hd)]
  have hden : 2 * r ≤ slope (r ^ 2) x := by unfold slope f; nlinarith
  have hx0 : 0 < x := lt_of_lt_of_le hr hx
  calc
    _ ≤ (x - r) ^ 2 * (1 + x + r) / (2 * r) :=
      div_le_div_of_nonneg_left (by positivity) (by positivity) hden
    _ ≤ (x - r) ^ 2 * (1 + R + r) / (2 * r) := by
      apply (div_le_div_iff_of_pos_right (by positivity : 0 < 2 * r)).mpr
      exact mul_le_mul_of_nonneg_left (by linarith) (sq_nonneg _)
    _ = _ := by ring

#print axioms divided_difference
#print axioms raw_equivalence
#print axioms fixed_root
#print axioms exact_error
#print axioms safe_denominator
#print axioms bracket
#print axioms quadratic_bound
end
end PaperToLean.Steffensen