import ProofExplorer.Jacobian.MarkusEight
import ProofExplorer.Jacobian.CalculusLibrary

namespace ProofExplorer.Jacobian.MarkusEight
open MvPolynomial
noncomputable section
set_option maxHeartbeats 0
set_option maxRecDepth 10000
variable {K : Type*} [Field K] [CharZero K]

def polynomialMap : Fin 8 → MvPolynomial (Fin 8) K :=
  ![(((C (-1)) * (X 0)) + ((C (1/2)) * (X 7) * ((X 5) ^ 3)) + ((C (-2)) * ((X 5) ^ 3) * ((X 6) ^ 2))),
    ((X 3) + ((C (-1)) * (X 1)) + ((C (15)) * ((X 5) ^ 2) * ((X 6) ^ 3)) + ((C (-6)) * (X 6) * (X 7) * ((X 5) ^ 2))),
    ((X 4) + ((C (-1)) * (X 2)) + ((C (9)) * ((X 5) ^ 2) * ((X 6) ^ 4)) + ((C (-3)) * (X 7) * ((X 5) ^ 2) * ((X 6) ^ 2))),
    (((C (-1)) * (X 3)) + ((C (12)) * ((X 5) ^ 3) * ((X 6) ^ 4)) + ((C (-3)) * (X 7) * ((X 5) ^ 3) * ((X 6) ^ 2))),
    (((C (-1)) * (X 4)) + ((C (4)) * ((X 5) ^ 3) * ((X 6) ^ 5)) + ((C (-1)) * (X 7) * ((X 5) ^ 3) * ((X 6) ^ 3))),
    ((X 0) + ((C (-1)) * (X 5)) + ((C (3/2)) * (X 6) * ((X 5) ^ 2))),
    ((X 1) + ((C (-1)) * (X 6)) + ((C (-3)) * (X 5) * (X 7))),
    ((C (-1/4)) + (X 2) + ((C (-1)) * (X 7)) + ((C (5)) * (X 5) * ((X 6) ^ 3)) + ((C (-3)) * (X 5) * (X 6) * (X 7)))]

def polynomialJacobian : Matrix (Fin 8) (Fin 8) (MvPolynomial (Fin 8) K) :=
  !![(C (-1)), (C (0)), (C (0)), (C (0)), (C (0)), (((C (-6)) * ((X 5) ^ 2) * ((X 6) ^ 2)) + ((C (3/2)) * (X 7) * ((X 5) ^ 2))), ((C (-4)) * (X 6) * ((X 5) ^ 3)), ((C (1/2)) * ((X 5) ^ 3));
    (C (0)), (C (-1)), (C (0)), (C (1)), (C (0)), (((C (30)) * (X 5) * ((X 6) ^ 3)) + ((C (-12)) * (X 5) * (X 6) * (X 7))), (((C (-6)) * (X 7) * ((X 5) ^ 2)) + ((C (45)) * ((X 5) ^ 2) * ((X 6) ^ 2))), ((C (-6)) * (X 6) * ((X 5) ^ 2));
    (C (0)), (C (0)), (C (-1)), (C (0)), (C (1)), (((C (18)) * (X 5) * ((X 6) ^ 4)) + ((C (-6)) * (X 5) * (X 7) * ((X 6) ^ 2))), (((C (36)) * ((X 5) ^ 2) * ((X 6) ^ 3)) + ((C (-6)) * (X 6) * (X 7) * ((X 5) ^ 2))), ((C (-3)) * ((X 5) ^ 2) * ((X 6) ^ 2));
    (C (0)), (C (0)), (C (0)), (C (-1)), (C (0)), (((C (36)) * ((X 5) ^ 2) * ((X 6) ^ 4)) + ((C (-9)) * (X 7) * ((X 5) ^ 2) * ((X 6) ^ 2))), (((C (48)) * ((X 5) ^ 3) * ((X 6) ^ 3)) + ((C (-6)) * (X 6) * (X 7) * ((X 5) ^ 3))), ((C (-3)) * ((X 5) ^ 3) * ((X 6) ^ 2));
    (C (0)), (C (0)), (C (0)), (C (0)), (C (-1)), (((C (12)) * ((X 5) ^ 2) * ((X 6) ^ 5)) + ((C (-3)) * (X 7) * ((X 5) ^ 2) * ((X 6) ^ 3))), (((C (20)) * ((X 5) ^ 3) * ((X 6) ^ 4)) + ((C (-3)) * (X 7) * ((X 5) ^ 3) * ((X 6) ^ 2))), ((C (-1)) * ((X 5) ^ 3) * ((X 6) ^ 3));
    (C (1)), (C (0)), (C (0)), (C (0)), (C (0)), ((C (-1)) + ((C (3)) * (X 5) * (X 6))), ((C (3/2)) * ((X 5) ^ 2)), (C (0));
    (C (0)), (C (1)), (C (0)), (C (0)), (C (0)), ((C (-3)) * (X 7)), (C (-1)), ((C (-3)) * (X 5));
    (C (0)), (C (0)), (C (1)), (C (0)), (C (0)), (((C (5)) * ((X 6) ^ 3)) + ((C (-3)) * (X 6) * (X 7))), (((C (-3)) * (X 5) * (X 7)) + ((C (15)) * (X 5) * ((X 6) ^ 2))), ((C (-1)) + ((C (-3)) * (X 5) * (X 6)))]

set_option maxHeartbeats 0 in
theorem polynomial_jacobian_entries (i j : Fin 8) :
    pderiv j (polynomialMap (K := K) i) = polynomialJacobian i j := by
  apply MvPolynomial.funext
  intro u
  fin_cases i <;> fin_cases j <;>
    simp [polynomialMap, polynomialJacobian, Derivation.leibniz, smul_eq_mul,
      Library.mv_pderiv_numeral, map_ofNat] <;> ring

theorem polynomialMap_evaluation (p : Fin 8 → K) (i : Fin 8) :
    eval p (polynomialMap i) = field (p 0) (p 1) (p 2) (p 3) (p 4) (p 5) (p 6) (p 7) i := by
  fin_cases i <;> simp [polynomialMap, field] <;> ring

theorem polynomialJacobian_evaluation (p : Fin 8 → K) (i j : Fin 8) :
    eval p (polynomialJacobian i j) = jacobian (p 5) (p 6) (p 7) i j := by
  fin_cases i <;> fin_cases j <;> simp [polynomialJacobian, jacobian] <;> ring

/-- Analytic coordinate derivatives of the actual real polynomial field. -/
theorem analytic_jacobian (p : Fin 8 → ℝ) (i j : Fin 8) :
    HasDerivAt (fun t => let q := Function.update p j t
      field (q 0) (q 1) (q 2) (q 3) (q 4) (q 5) (q 6) (q 7) i)
      (jacobian (p 5) (p 6) (p 7) i j) (p j) := by
  have h := CalculusLibrary.eval_hasDerivAt (polynomialMap i) p j
  simpa only [polynomial_jacobian_entries, polynomialJacobian_evaluation,
    polynomialMap_evaluation] using h

end
end ProofExplorer.Jacobian.MarkusEight