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

/-! Halley (1694), square-root specialization of the rational iteration.
Exact arithmetic, positive target root. Cubic error identities are proved;
the general smooth-function convergence theorem is outside this file. -/
namespace PaperToLean.Halley
noncomputable section
def step (a x : ℝ) : ℝ := x * (x ^ 2 + 3 * a) / (3 * x ^ 2 + a)
def ratio (r x : ℝ) : ℝ := (x - r) / (x + r)
def iterate (a x : ℝ) : ℕ → ℝ
  | 0 => x
  | n + 1 => step a (iterate a x n)

theorem denominator_pos (r x : ℝ) (hr : 0 < r) : 0 < 3 * x ^ 2 + r ^ 2 := by
  positivity

theorem positive_step (r x : ℝ) (hr : 0 < r) (hx : 0 < x) : 0 < step (r ^ 2) x := by
  unfold step
  positivity

theorem cubic_error (r x : ℝ) (hr : 0 < r) :
    step (r ^ 2) x - r = (x - r) ^ 3 / (3 * x ^ 2 + r ^ 2) := by
  have hd := ne_of_gt (denominator_pos r x hr)
  unfold step
  field_simp
  ring

theorem bracket (r x : ℝ) (hr : 0 < r) (hx : r ≤ x) :
    r ≤ step (r ^ 2) x ∧ step (r ^ 2) x ≤ x := by
  have hd := denominator_pos r x hr
  constructor
  · have he := cubic_error r x hr
    have hp : 0 ≤ (x - r) ^ 3 / (3 * x ^ 2 + r ^ 2) := by positivity
    linarith
  · unfold step
    apply (div_le_iff₀ hd).mpr
    have hx0 : 0 ≤ x := le_trans (le_of_lt hr) hx
    nlinarith [mul_nonneg hx0 (show 0 ≤ x ^ 2 - r ^ 2 by nlinarith)]

theorem cubic_bound (r x : ℝ) (hr : 0 < r) (hx : r ≤ x) :
    step (r ^ 2) x - r ≤ (x - r) ^ 3 / (4 * r ^ 2) := by
  rw [cubic_error r x hr]
  apply div_le_div_of_nonneg_left (by positivity) (by positivity)
  nlinarith

theorem ratio_cubed (r x : ℝ) (hr : 0 < r) (hx : 0 < x) :
    ratio r (step (r ^ 2) x) = (ratio r x) ^ 3 := by
  have hd := ne_of_gt (denominator_pos r x hr)
  have hs : x + r ≠ 0 := ne_of_gt (add_pos hx hr)
  have hc : x ^ 3 + 3 * x ^ 2 * r + 3 * x * r ^ 2 + r ^ 3 ≠ 0 := by positivity
  unfold ratio step
  field_simp
  ring

theorem iterate_positive (r x : ℝ) (hr : 0 < r) (hx : 0 < x) (n : ℕ) :
    0 < iterate (r ^ 2) x n := by
  induction n with
  | zero => exact hx
  | succ n ih => exact positive_step r _ hr ih

theorem ratio_at_every_step (r x : ℝ) (hr : 0 < r) (hx : 0 < x) (n : ℕ) :
    ratio r (iterate (r ^ 2) x n) = (ratio r x) ^ (3 ^ n) := by
  induction n with
  | zero => simp [iterate]
  | succ n ih =>
    rw [iterate, ratio_cubed r _ hr (iterate_positive r x hr hx n), ih,
      ← pow_mul, pow_succ]

#print axioms denominator_pos
#print axioms positive_step
#print axioms cubic_error
#print axioms bracket
#print axioms cubic_bound
#print axioms ratio_cubed
#print axioms iterate_positive
#print axioms ratio_at_every_step
end
end PaperToLean.Halley