theory Prove_Less
imports
  Affine_Code
begin

text {* TODO: taken from ODEs *}
definition "split_aform_largest_uncond X =
    (let (i, x) = max_pdev (snd X) in split_aform X i)"

definition "split_aform_largest p t X =
  (let (a, b) =
      split_aform_largest_uncond (fst X, summarize_threshold p t (degree_aform X) (snd X))
    in [a, b])"

subsection {* Positivity by splitting *}

approximate_affine schwefel "\<lambda>(x1, x2, x3).
  (FloatR 58806 0 * inverse (FloatR 6103515625 14) +
    (x1 + - (x2*x2))*(x1 + - (x2*x2)) +
    (x2 + - 1)*(x2 + - 1) +
    (x1 + - (x3*x3))*(x1 + - (x3*x3)) +
    (x3 + - 1)*(x3 + - 1))*\<^sub>R1::real"
term schwefel
value "Inf_aform' 80 (the (schwefel 80 (FloatR 1 (-10)) (aform_of_ivl (0,0,0) (0,0,0)) []))"
term split_aform
primrec prove_pos::"nat \<Rightarrow> nat \<Rightarrow> real \<Rightarrow> real \<Rightarrow> (nat \<Rightarrow> real \<Rightarrow> 'a::executable_euclidean_space aform \<Rightarrow> _ \<Rightarrow> real aform option) \<Rightarrow> 'a aform \<Rightarrow> bool" where
  "prove_pos 0 p t st F X = False"
| "prove_pos (Suc i) p t st F X = (if 0 < Inf_aform' p (the (F p t X [])) then True else
    case split_aform_largest p st X of
      [a, b] \<Rightarrow> prove_pos i p t st F a \<and> prove_pos i p t st F b
    | _ \<Rightarrow> False)"

value "prove_pos 80 20 (FloatR  1 (-1)) (FloatR 1 1) schwefel (aform_of_ivl (-10,-10,-10) (10,10,10))"

approximate_affine delta6 "\<lambda>(x1, x2, x3, x4, x5, x6).
  (x1 * x4 * (-x1 + x2 + x3 + -x4 + x5 + x6) +
   x2 * x5 * (x1 + -x2 + x3 + x4 + - x5 + x6) +
   x3 * x6 * (x1 + x2 + - x3 + x4 + x5 + -x6) +
   - x2 * x3 * x4 +
   - x1 * x3 * x5 +
   - x1 * x2 * x6 +
   - x4 * x5 * x6)*\<^sub>R1::real"

definition "all6 x = (x, x, x, x, x, x)"

value "prove_pos 20 10 (FloatR 1 (-1)) (FloatR 1 1) delta6 (aform_of_ivl (all6 4) (all6 (FloatR 104045 (-14))))"

definition "all7 x = (x, x, x, x, x, x, x)"

approximate_affine magnetism "\<lambda>(x1, x2, x3, x4, x5, x6, x7).
  (FloatR 409616384 (-14) + x1 * x1 + 2 * x2*x2 + 2 * x3 * x3 + 2 * x4 * x4 + 2 * x5 * x5 + 2* x6 * x6 + 2 * x7 * x7 + - x1)*\<^sub>R1::real"

value "prove_pos 20 10 (FloatR 1 (-1)) (FloatR 1 1) magnetism (aform_of_ivl (all7 (-1)) (all7 1))"

definition "all4 x = (x, x, x, x)"
approximate_affine caprasse "\<lambda>(x1, x2, x3, x4).
  (FloatR 7160948587500457 (-51) + - x1 * x3 * x3 * x3 + 4 * x2 * x3 * x3 * x4 + 4 * x1 * x3 * x4 * x4 + 2 * x2 * x4 * x4 * x4
    + 4 * x1*x3 + 4*x3*x3 + - 10*x2*x4 + -10*x4*x4 + 2)*\<^sub>R1::real"

value "prove_pos 50 40 (FloatR 1 (-3)) (FloatR 1 (-1)) caprasse (aform_of_ivl (all4 (-FloatR 1 (-1))) (all4 (FloatR 1 (-1))))"

approximate_affine heart "\<lambda>(x1, x2, x3, x4, x5, x6, x7, x8).
  (FloatR 7852025950320460 (-52) + - x1 * x6*x6*x6 + 3 * x1 * x6 * x7 * x7 + - x3 * x7 * x7 * x7 + 3 * x3 * x7 * x6 * x6 + - x2 * x5 * x5 * x5 +
    3*x2*x5*x8*x8 + - x4 * x8*x8*x8 + 3 * x4 * x8 * x5 * x5 + FloatR (-8613992673435050) (-53))*\<^sub>R1::real"

definition "p1 = FloatR 7205759403792794 (-56)"

value "prove_pos 80 20 (FloatR 1 (-2)) (FloatR 1 1) heart (aform_of_ivl (-p1,4*p1,-7*p1,-7*p1,p1,-p1,-3*p1,-1 - p1) (4*p1,1,-4*p1,4*p1,2*p1,2*p1,1 + p1,-3*p1))"

end
