import Mathlib.LinearAlgebra.BilinearForm.Properties
import Mathlib.Data.Real.Basic
import Mathlib.Tactic.FieldSimp
import Mathlib.Tactic.Ring
import Mathlib.Tactic.Linarith
import Mathlib.Tactic.Positivity

/-!
Fletcher and Powell (1963), DOI 10.1093/comjnl/6.2.163, sections 3--4.
DFP update in bilinear-form notation: B(x,z) = xᵀHz and l(x) = sᵀx.
The curvature condition l(y) > 0 is explicit. No line search, finite
termination, or convergence of the optimization algorithm is formalized.
-/
namespace PaperToLean.FletcherPowell
noncomputable section
variable {V : Type*} [AddCommGroup V] [Module ℝ V]
open LinearMap (BilinForm)

def update (B : BilinForm ℝ V) (l : V →ₗ[ℝ] ℝ) (y x z : V) : ℝ :=
  B x z + l x * l z / l y - B x y * B y z / B y y

theorem symmetric (B : BilinForm ℝ V) (l : V →ₗ[ℝ] ℝ) (y : V)
    (hB : B.IsSymm) (x z : V) : update B l y x z = update B l y z x := by
  unfold update
  rw [hB.eq z x, hB.eq y z, hB.eq y x]
  ring

theorem add_left (B : BilinForm ℝ V) (l : V →ₗ[ℝ] ℝ) (y x z w : V) :
    update B l y (x + z) w = update B l y x w + update B l y z w := by
  simp only [update, map_add, LinearMap.add_apply, div_eq_mul_inv]
  ring

theorem smul_left (B : BilinForm ℝ V) (l : V →ₗ[ℝ] ℝ)
    (y x z : V) (c : ℝ) : update B l y (c • x) z = c * update B l y x z := by
  simp only [update, map_smul, LinearMap.smul_apply, smul_eq_mul, div_eq_mul_inv]
  ring

theorem secant (B : BilinForm ℝ V) (l : V →ₗ[ℝ] ℝ) (y x : V)
    (hl : l y ≠ 0) (hy : B y y ≠ 0) : update B l y y x = l x := by
  unfold update
  field_simp
  ring

theorem completed_square (B : BilinForm ℝ V) (l : V →ₗ[ℝ] ℝ) (y x : V)
    (hB : B.IsSymm) (hy : B y y ≠ 0) :
    update B l y x x =
      B (x - (B y x / B y y) • y) (x - (B y x / B y y) • y) + (l x) ^ 2 / l y := by
  simp only [update, BilinForm.sub_left, BilinForm.sub_right,
    BilinForm.smul_left, BilinForm.smul_right]
  rw [hB.eq x y]
  field_simp
  ring

/-- Strict positivity is preserved; the null direction in the first term
is controlled by the positive curvature term. -/
theorem positive_definite (B : BilinForm ℝ V) (l : V →ₗ[ℝ] ℝ) (y : V)
    (hB : B.IsSymm) (hpos : ∀ x : V, x ≠ 0 → 0 < B x x)
    (hcurv : 0 < l y) : ∀ x : V, x ≠ 0 → 0 < update B l y x x := by
  have hy0 : y ≠ 0 := by
    intro h
    simp [h] at hcurv
  have hy := hpos y hy0
  intro x hx
  rw [completed_square B l y x hB (ne_of_gt hy)]
  let z := x - (B y x / B y y) • y
  change 0 < B z z + (l x) ^ 2 / l y
  have hterm : 0 ≤ (l x) ^ 2 / l y := by positivity
  by_cases hz : z = 0
  · have hxeq : x = (B y x / B y y) • y := sub_eq_zero.mp hz
    have hlx : l x ≠ 0 := by
      intro hzero
      have hcoeff : B y x / B y y = 0 := by
        have heq := congrArg l hxeq
        simp only [map_smul, smul_eq_mul, hzero] at heq
        exact (mul_eq_zero.mp heq.symm).resolve_right (ne_of_gt hcurv)
      apply hx
      rw [hxeq, hcoeff, zero_smul]
    have hstrict : 0 < (l x) ^ 2 / l y := by positivity
    simpa [hz] using hstrict
  · exact add_pos_of_pos_of_nonneg (hpos z hz) hterm

#print axioms symmetric
#print axioms add_left
#print axioms smul_left
#print axioms secant
#print axioms completed_square
#print axioms positive_definite
end
end PaperToLean.FletcherPowell