renamed Sum_Of_Squares to Sum_of_Squares;
authorwenzelm
Sat Jan 08 17:39:51 2011 +0100 (2011-01-08)
changeset 4147460d091240485
parent 41473 3717fc42ebe9
child 41475 fe4f0d9f9dbb
renamed Sum_Of_Squares to Sum_of_Squares;
etc/components
src/HOL/IsaMakefile
src/HOL/Library/Library.thy
src/HOL/Library/Sum_Of_Squares.thy
src/HOL/Library/Sum_Of_Squares/etc/settings
src/HOL/Library/Sum_Of_Squares/neos_csdp_client
src/HOL/Library/Sum_Of_Squares/positivstellensatz_tools.ML
src/HOL/Library/Sum_Of_Squares/sos_wrapper.ML
src/HOL/Library/Sum_Of_Squares/sum_of_squares.ML
src/HOL/Library/Sum_of_Squares.thy
src/HOL/Library/Sum_of_Squares/etc/settings
src/HOL/Library/Sum_of_Squares/neos_csdp_client
src/HOL/Library/Sum_of_Squares/positivstellensatz_tools.ML
src/HOL/Library/Sum_of_Squares/sos_wrapper.ML
src/HOL/Library/Sum_of_Squares/sum_of_squares.ML
     1.1 --- a/etc/components	Sat Jan 08 17:30:05 2011 +0100
     1.2 +++ b/etc/components	Sat Jan 08 17:39:51 2011 +0100
     1.3 @@ -14,7 +14,7 @@
     1.4  src/Tools/WWW_Find
     1.5  src/HOL/Tools/ATP
     1.6  src/HOL/Mirabelle
     1.7 -src/HOL/Library/Sum_Of_Squares
     1.8 +src/HOL/Library/Sum_of_Squares
     1.9  src/HOL/Tools/SMT
    1.10  src/HOL/Tools/Predicate_Compile
    1.11  src/HOL/Mutabelle
     2.1 --- a/src/HOL/IsaMakefile	Sat Jan 08 17:30:05 2011 +0100
     2.2 +++ b/src/HOL/IsaMakefile	Sat Jan 08 17:39:51 2011 +0100
     2.3 @@ -460,9 +460,9 @@
     2.4    Library/RBT.thy Library/RBT_Impl.thy Library/README.html		\
     2.5    Library/Set_Algebras.thy Library/State_Monad.thy Library/Ramsey.thy	\
     2.6    Library/Reflection.thy Library/SML_Quickcheck.thy 			\
     2.7 -  Library/Sublist_Order.thy Library/Sum_Of_Squares.thy			\
     2.8 -  Library/Sum_Of_Squares/sos_wrapper.ML					\
     2.9 -  Library/Sum_Of_Squares/sum_of_squares.ML				\
    2.10 +  Library/Sublist_Order.thy Library/Sum_of_Squares.thy			\
    2.11 +  Library/Sum_of_Squares/sos_wrapper.ML					\
    2.12 +  Library/Sum_of_Squares/sum_of_squares.ML				\
    2.13    Library/Transitive_Closure_Table.thy Library/Univ_Poly.thy		\
    2.14    Library/While_Combinator.thy Library/Zorn.thy				\
    2.15    $(SRC)/Tools/adhoc_overloading.ML Library/positivstellensatz.ML	\
     3.1 --- a/src/HOL/Library/Library.thy	Sat Jan 08 17:30:05 2011 +0100
     3.2 +++ b/src/HOL/Library/Library.thy	Sat Jan 08 17:39:51 2011 +0100
     3.3 @@ -56,7 +56,7 @@
     3.4    Set_Algebras
     3.5    SML_Quickcheck
     3.6    State_Monad
     3.7 -  Sum_Of_Squares
     3.8 +  Sum_of_Squares
     3.9    Transitive_Closure_Table
    3.10    Univ_Poly
    3.11    While_Combinator
     4.1 --- a/src/HOL/Library/Sum_Of_Squares.thy	Sat Jan 08 17:30:05 2011 +0100
     4.2 +++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
     4.3 @@ -1,159 +0,0 @@
     4.4 -(*  Title:      HOL/Library/Sum_Of_Squares.thy
     4.5 -    Author:     Amine Chaieb, University of Cambridge
     4.6 -    Author:     Philipp Meyer, TU Muenchen
     4.7 -*)
     4.8 -
     4.9 -header {* A decision method for universal multivariate real arithmetic with addition, 
    4.10 -  multiplication and ordering using semidefinite programming *}
    4.11 -
    4.12 -theory Sum_Of_Squares
    4.13 -imports Complex_Main
    4.14 -uses
    4.15 -  "positivstellensatz.ML"
    4.16 -  "Sum_Of_Squares/sum_of_squares.ML"
    4.17 -  "Sum_Of_Squares/positivstellensatz_tools.ML"
    4.18 -  "Sum_Of_Squares/sos_wrapper.ML"
    4.19 -begin
    4.20 -
    4.21 -text {*
    4.22 -  In order to use the method sos, call it with @{text "(sos
    4.23 -  remote_csdp)"} to use the remote solver.  Or install CSDP
    4.24 -  (https://projects.coin-or.org/Csdp), configure the Isabelle setting
    4.25 -  @{text CSDP_EXE}, and call it with @{text "(sos csdp)"}.  By
    4.26 -  default, sos calls @{text remote_csdp}.  This can take of the order
    4.27 -  of a minute for one sos call, because sos calls CSDP repeatedly.  If
    4.28 -  you install CSDP locally, sos calls typically takes only a few
    4.29 -  seconds.
    4.30 -  sos generates a certificate which can be used to repeat the proof
    4.31 -  without calling an external prover.
    4.32 -*}
    4.33 -
    4.34 -setup Sum_Of_Squares.setup
    4.35 -setup SOS_Wrapper.setup
    4.36 -
    4.37 -text {* Tests *}
    4.38 -
    4.39 -lemma "(3::real) * x + 7 * a < 4 & 3 < 2 * x \<Longrightarrow> a < 0"
    4.40 -by (sos_cert "((R<1 + (((A<1 * R<1) * (R<2 * [1]^2)) + (((A<0 * R<1) * (R<3 * [1]^2)) + ((A<=0 * R<1) * (R<14 * [1]^2))))))")
    4.41 -
    4.42 -lemma "a1 >= 0 & a2 >= 0 \<and> (a1 * a1 + a2 * a2 = b1 * b1 + b2 * b2 + 2) \<and> (a1 * b1 + a2 * b2 = 0) --> a1 * a2 - b1 * b2 >= (0::real)"
    4.43 -by (sos_cert "(((A<0 * R<1) + (([~1/2*a1*b2 + ~1/2*a2*b1] * A=0) + (([~1/2*a1*a2 + 1/2*b1*b2] * A=1) + (((A<0 * R<1) * ((R<1/2 * [b2]^2) + (R<1/2 * [b1]^2))) + ((A<=0 * (A<=1 * R<1)) * ((R<1/2 * [b2]^2) + ((R<1/2 * [b1]^2) + ((R<1/2 * [a2]^2) + (R<1/2 * [a1]^2))))))))))")
    4.44 -
    4.45 -lemma "(3::real) * x + 7 * a < 4 & 3 < 2 * x --> a < 0"
    4.46 -by (sos_cert "((R<1 + (((A<1 * R<1) * (R<2 * [1]^2)) + (((A<0 * R<1) * (R<3 * [1]^2)) + ((A<=0 * R<1) * (R<14 * [1]^2))))))")
    4.47 -
    4.48 -lemma "(0::real) <= x & x <= 1 & 0 <= y & y <= 1  --> x^2 + y^2 < 1 |(x - 1)^2 + y^2 < 1 | x^2 + (y - 1)^2 < 1 | (x - 1)^2 + (y - 1)^2 < 1"
    4.49 -by (sos_cert "((R<1 + (((A<=3 * (A<=4 * R<1)) * (R<1 * [1]^2)) + (((A<=2 * (A<=7 * R<1)) * (R<1 * [1]^2)) + (((A<=1 * (A<=6 * R<1)) * (R<1 * [1]^2)) + ((A<=0 * (A<=5 * R<1)) * (R<1 * [1]^2)))))))")
    4.50 -
    4.51 -lemma "(0::real) <= x & 0 <= y & 0 <= z & x + y + z <= 3 --> x * y + x * z + y * z >= 3 * x * y * z"
    4.52 -by (sos_cert "(((A<0 * R<1) + (((A<0 * R<1) * (R<1/2 * [1]^2)) + (((A<=2 * R<1) * (R<1/2 * [~1*x + y]^2)) + (((A<=1 * R<1) * (R<1/2 * [~1*x + z]^2)) + (((A<=1 * (A<=2 * (A<=3 * R<1))) * (R<1/2 * [1]^2)) + (((A<=0 * R<1) * (R<1/2 * [~1*y + z]^2)) + (((A<=0 * (A<=2 * (A<=3 * R<1))) * (R<1/2 * [1]^2)) + ((A<=0 * (A<=1 * (A<=3 * R<1))) * (R<1/2 * [1]^2))))))))))")
    4.53 -
    4.54 -lemma "((x::real)^2 + y^2 + z^2 = 1) --> (x + y + z)^2 <= 3"
    4.55 -by (sos_cert "(((A<0 * R<1) + (([~3] * A=0) + (R<1 * ((R<2 * [~1/2*x + ~1/2*y + z]^2) + (R<3/2 * [~1*x + y]^2))))))")
    4.56 -
    4.57 -lemma "(w^2 + x^2 + y^2 + z^2 = 1) --> (w + x + y + z)^2 <= (4::real)"
    4.58 -by (sos_cert "(((A<0 * R<1) + (([~4] * A=0) + (R<1 * ((R<3 * [~1/3*w + ~1/3*x + ~1/3*y + z]^2) + ((R<8/3 * [~1/2*w + ~1/2*x + y]^2) + (R<2 * [~1*w + x]^2)))))))")
    4.59 -
    4.60 -lemma "(x::real) >= 1 & y >= 1 --> x * y >= x + y - 1"
    4.61 -by (sos_cert "(((A<0 * R<1) + ((A<=0 * (A<=1 * R<1)) * (R<1 * [1]^2))))")
    4.62 -
    4.63 -lemma "(x::real) > 1 & y > 1 --> x * y > x + y - 1"
    4.64 -by (sos_cert "((((A<0 * A<1) * R<1) + ((A<=0 * R<1) * (R<1 * [1]^2))))") 
    4.65 -
    4.66 -lemma "abs(x) <= 1 --> abs(64 * x^7 - 112 * x^5 + 56 * x^3 - 7 * x) <= (1::real)"
    4.67 -by (sos_cert "((((A<0 * R<1) + ((A<=1 * R<1) * (R<1 * [~8*x^3 + ~4*x^2 + 4*x + 1]^2)))) & ((((A<0 * A<1) * R<1) + ((A<=1 * (A<0 * R<1)) * (R<1 * [8*x^3 + ~4*x^2 + ~4*x + 1]^2)))))")
    4.68 -
    4.69 -(* ------------------------------------------------------------------------- *)
    4.70 -(* One component of denominator in dodecahedral example.                     *)
    4.71 -(* ------------------------------------------------------------------------- *)
    4.72 -
    4.73 -lemma "2 <= x & x <= 125841 / 50000 & 2 <= y & y <= 125841 / 50000 & 2 <= z & z <= 125841 / 50000 --> 2 * (x * z + x * y + y * z) - (x * x + y * y + z * z) >= (0::real)"
    4.74 -by (sos_cert "(((A<0 * R<1) + ((R<1 * ((R<5749028157/5000000000 * [~25000/222477*x + ~25000/222477*y + ~25000/222477*z + 1]^2) + ((R<864067/1779816 * [419113/864067*x + 419113/864067*y + z]^2) + ((R<320795/864067 * [419113/1283180*x + y]^2) + (R<1702293/5132720 * [x]^2))))) + (((A<=4 * (A<=5 * R<1)) * (R<3/2 * [1]^2)) + (((A<=3 * (A<=5 * R<1)) * (R<1/2 * [1]^2)) + (((A<=2 * (A<=4 * R<1)) * (R<1 * [1]^2)) + (((A<=2 * (A<=3 * R<1)) * (R<3/2 * [1]^2)) + (((A<=1 * (A<=5 * R<1)) * (R<1/2 * [1]^2)) + (((A<=1 * (A<=3 * R<1)) * (R<1/2 * [1]^2)) + (((A<=0 * (A<=4 * R<1)) * (R<1 * [1]^2)) + (((A<=0 * (A<=2 * R<1)) * (R<1 * [1]^2)) + ((A<=0 * (A<=1 * R<1)) * (R<3/2 * [1]^2)))))))))))))")
    4.75 -
    4.76 -(* ------------------------------------------------------------------------- *)
    4.77 -(* Over a larger but simpler interval.                                       *)
    4.78 -(* ------------------------------------------------------------------------- *)
    4.79 -
    4.80 -lemma "(2::real) <= x & x <= 4 & 2 <= y & y <= 4 & 2 <= z & z <= 4 --> 0 <= 2 * (x * z + x * y + y * z) - (x * x + y * y + z * z)"
    4.81 -by (sos_cert "((R<1 + ((R<1 * ((R<1 * [~1/6*x + ~1/6*y + ~1/6*z + 1]^2) + ((R<1/18 * [~1/2*x + ~1/2*y + z]^2) + (R<1/24 * [~1*x + y]^2)))) + (((A<0 * R<1) * (R<1/12 * [1]^2)) + (((A<=4 * (A<=5 * R<1)) * (R<1/6 * [1]^2)) + (((A<=2 * (A<=4 * R<1)) * (R<1/6 * [1]^2)) + (((A<=2 * (A<=3 * R<1)) * (R<1/6 * [1]^2)) + (((A<=0 * (A<=4 * R<1)) * (R<1/6 * [1]^2)) + (((A<=0 * (A<=2 * R<1)) * (R<1/6 * [1]^2)) + ((A<=0 * (A<=1 * R<1)) * (R<1/6 * [1]^2)))))))))))")
    4.82 -
    4.83 -(* ------------------------------------------------------------------------- *)
    4.84 -(* We can do 12. I think 12 is a sharp bound; see PP's certificate.          *)
    4.85 -(* ------------------------------------------------------------------------- *)
    4.86 -
    4.87 -lemma "2 <= (x::real) & x <= 4 & 2 <= y & y <= 4 & 2 <= z & z <= 4 --> 12 <= 2 * (x * z + x * y + y * z) - (x * x + y * y + z * z)"
    4.88 -by (sos_cert "(((A<0 * R<1) + (((A<=4 * R<1) * (R<2/3 * [1]^2)) + (((A<=4 * (A<=5 * R<1)) * (R<1 * [1]^2)) + (((A<=3 * (A<=4 * R<1)) * (R<1/3 * [1]^2)) + (((A<=2 * R<1) * (R<2/3 * [1]^2)) + (((A<=2 * (A<=5 * R<1)) * (R<1/3 * [1]^2)) + (((A<=2 * (A<=4 * R<1)) * (R<8/3 * [1]^2)) + (((A<=2 * (A<=3 * R<1)) * (R<1 * [1]^2)) + (((A<=1 * (A<=4 * R<1)) * (R<1/3 * [1]^2)) + (((A<=1 * (A<=2 * R<1)) * (R<1/3 * [1]^2)) + (((A<=0 * R<1) * (R<2/3 * [1]^2)) + (((A<=0 * (A<=5 * R<1)) * (R<1/3 * [1]^2)) + (((A<=0 * (A<=4 * R<1)) * (R<8/3 * [1]^2)) + (((A<=0 * (A<=3 * R<1)) * (R<1/3 * [1]^2)) + (((A<=0 * (A<=2 * R<1)) * (R<8/3 * [1]^2)) + ((A<=0 * (A<=1 * R<1)) * (R<1 * [1]^2))))))))))))))))))")
    4.89 -
    4.90 -(* ------------------------------------------------------------------------- *)
    4.91 -(* Inequality from sci.math (see "Leon-Sotelo, por favor").                  *)
    4.92 -(* ------------------------------------------------------------------------- *)
    4.93 -
    4.94 -lemma "0 <= (x::real) & 0 <= y & (x * y = 1) --> x + y <= x^2 + y^2"
    4.95 -by (sos_cert "(((A<0 * R<1) + (([1] * A=0) + (R<1 * ((R<1 * [~1/2*x + ~1/2*y + 1]^2) + (R<3/4 * [~1*x + y]^2))))))") 
    4.96 -
    4.97 -lemma "0 <= (x::real) & 0 <= y & (x * y = 1) --> x * y * (x + y) <= x^2 + y^2"
    4.98 -by (sos_cert "(((A<0 * R<1) + (([~1*x + ~1*y + 1] * A=0) + (R<1 * ((R<1 * [~1/2*x + ~1/2*y + 1]^2) + (R<3/4 * [~1*x + y]^2))))))") 
    4.99 -
   4.100 -lemma "0 <= (x::real) & 0 <= y --> x * y * (x + y)^2 <= (x^2 + y^2)^2"
   4.101 -by (sos_cert "(((A<0 * R<1) + (R<1 * ((R<1 * [~1/2*x^2 + y^2 + ~1/2*x*y]^2) + (R<3/4 * [~1*x^2 + x*y]^2)))))")
   4.102 -
   4.103 -lemma "(0::real) <= a & 0 <= b & 0 <= c & c * (2 * a + b)^3/ 27 <= x \<longrightarrow> c * a^2 * b <= x"
   4.104 -by (sos_cert "(((A<0 * R<1) + (((A<=3 * R<1) * (R<1 * [1]^2)) + (((A<=1 * (A<=2 * R<1)) * (R<1/27 * [~1*a + b]^2)) + ((A<=0 * (A<=2 * R<1)) * (R<8/27 * [~1*a + b]^2))))))")
   4.105 - 
   4.106 -lemma "(0::real) < x --> 0 < 1 + x + x^2"
   4.107 -by (sos_cert "((R<1 + ((R<1 * (R<1 * [x]^2)) + (((A<0 * R<1) * (R<1 * [1]^2)) + ((A<=0 * R<1) * (R<1 * [1]^2))))))")
   4.108 -
   4.109 -lemma "(0::real) <= x --> 0 < 1 + x + x^2"
   4.110 -by (sos_cert "((R<1 + ((R<1 * (R<1 * [x]^2)) + (((A<=1 * R<1) * (R<1 * [1]^2)) + ((A<=0 * R<1) * (R<1 * [1]^2))))))")
   4.111 -
   4.112 -lemma "(0::real) < 1 + x^2"
   4.113 -by (sos_cert "((R<1 + ((R<1 * (R<1 * [x]^2)) + ((A<=0 * R<1) * (R<1 * [1]^2)))))")
   4.114 -
   4.115 -lemma "(0::real) <= 1 + 2 * x + x^2"
   4.116 -by (sos_cert "(((A<0 * R<1) + (R<1 * (R<1 * [x + 1]^2))))")
   4.117 -
   4.118 -lemma "(0::real) < 1 + abs x"
   4.119 -by (sos_cert "((R<1 + (((A<=1 * R<1) * (R<1/2 * [1]^2)) + ((A<=0 * R<1) * (R<1/2 * [1]^2)))))")
   4.120 -
   4.121 -lemma "(0::real) < 1 + (1 + x)^2 * (abs x)"
   4.122 -by (sos_cert "(((R<1 + (((A<=1 * R<1) * (R<1 * [1]^2)) + ((A<=0 * R<1) * (R<1 * [x + 1]^2))))) & ((R<1 + (((A<0 * R<1) * (R<1 * [x + 1]^2)) + ((A<=0 * R<1) * (R<1 * [1]^2))))))")
   4.123 -
   4.124 -
   4.125 -
   4.126 -lemma "abs ((1::real) + x^2) = (1::real) + x^2"
   4.127 -by (sos_cert "(() & (((R<1 + ((R<1 * (R<1 * [x]^2)) + ((A<1 * R<1) * (R<1/2 * [1]^2))))) & ((R<1 + ((R<1 * (R<1 * [x]^2)) + ((A<0 * R<1) * (R<1 * [1]^2)))))))")
   4.128 -lemma "(3::real) * x + 7 * a < 4 \<and> 3 < 2 * x \<longrightarrow> a < 0"
   4.129 -by (sos_cert "((R<1 + (((A<1 * R<1) * (R<2 * [1]^2)) + (((A<0 * R<1) * (R<3 * [1]^2)) + ((A<=0 * R<1) * (R<14 * [1]^2))))))")
   4.130 -
   4.131 -lemma "(0::real) < x --> 1 < y --> y * x <= z --> x < z"
   4.132 -by (sos_cert "((((A<0 * A<1) * R<1) + (((A<=1 * R<1) * (R<1 * [1]^2)) + ((A<=0 * R<1) * (R<1 * [1]^2)))))")
   4.133 -lemma "(1::real) < x --> x^2 < y --> 1 < y"
   4.134 -by (sos_cert "((((A<0 * A<1) * R<1) + ((R<1 * ((R<1/10 * [~2*x + y + 1]^2) + (R<1/10 * [~1*x + y]^2))) + (((A<1 * R<1) * (R<1/2 * [1]^2)) + (((A<0 * R<1) * (R<1 * [x]^2)) + (((A<=0 * R<1) * ((R<1/10 * [x + 1]^2) + (R<1/10 * [x]^2))) + (((A<=0 * (A<1 * R<1)) * (R<1/5 * [1]^2)) + ((A<=0 * (A<0 * R<1)) * (R<1/5 * [1]^2)))))))))")
   4.135 -lemma "(b::real)^2 < 4 * a * c --> ~(a * x^2 + b * x + c = 0)"
   4.136 -by (sos_cert "(((A<0 * R<1) + (R<1 * (R<1 * [2*a*x + b]^2))))")
   4.137 -lemma "(b::real)^2 < 4 * a * c --> ~(a * x^2 + b * x + c = 0)"
   4.138 -by (sos_cert "(((A<0 * R<1) + (R<1 * (R<1 * [2*a*x + b]^2))))")
   4.139 -lemma "((a::real) * x^2 + b * x + c = 0) --> b^2 >= 4 * a * c"
   4.140 -by (sos_cert "(((A<0 * R<1) + (R<1 * (R<1 * [2*a*x + b]^2))))")
   4.141 -lemma "(0::real) <= b & 0 <= c & 0 <= x & 0 <= y & (x^2 = c) & (y^2 = a^2 * c + b) --> a * c <= y * x"
   4.142 -by (sos_cert "(((A<0 * (A<0 * R<1)) + (((A<=2 * (A<=3 * (A<0 * R<1))) * (R<2 * [1]^2)) + ((A<=0 * (A<=1 * R<1)) * (R<1 * [1]^2)))))")
   4.143 -lemma "abs(x - z) <= e & abs(y - z) <= e & 0 <= u & 0 <= v & (u + v = 1) --> abs((u * x + v * y) - z) <= (e::real)"
   4.144 -by (sos_cert "((((A<0 * R<1) + (((A<=3 * (A<=6 * R<1)) * (R<1 * [1]^2)) + ((A<=1 * (A<=5 * R<1)) * (R<1 * [1]^2))))) & ((((A<0 * A<1) * R<1) + (((A<=3 * (A<=5 * (A<0 * R<1))) * (R<1 * [1]^2)) + ((A<=1 * (A<=4 * (A<0 * R<1))) * (R<1 * [1]^2))))))")
   4.145 -
   4.146 -
   4.147 -(* lemma "((x::real) - y - 2 * x^4 = 0) & 0 <= x & x <= 2 & 0 <= y & y <= 3 --> y^2 - 7 * y - 12 * x + 17 >= 0" by sos *) (* Too hard?*)
   4.148 -
   4.149 -lemma "(0::real) <= x --> (1 + x + x^2)/(1 + x^2) <= 1 + x"
   4.150 -by (sos_cert "(((((A<0 * A<1) * R<1) + ((A<=0 * (A<0 * R<1)) * (R<1 * [x]^2)))) & ((R<1 + ((R<1 * (R<1 * [x]^2)) + ((A<0 * R<1) * (R<1 * [1]^2))))))")
   4.151 -
   4.152 -lemma "(0::real) <= x --> 1 - x <= 1 / (1 + x + x^2)"
   4.153 -by (sos_cert "(((R<1 + (([~4/3] * A=0) + ((R<1 * ((R<1/3 * [3/2*x + 1]^2) + (R<7/12 * [x]^2))) + ((A<=0 * R<1) * (R<1/3 * [1]^2)))))) & (((((A<0 * A<1) * R<1) + ((A<=0 * (A<0 * R<1)) * (R<1 * [x]^2)))) & ((R<1 + ((R<1 * (R<1 * [x]^2)) + (((A<0 * R<1) * (R<1 * [1]^2)) + ((A<=0 * R<1) * (R<1 * [1]^2))))))))")
   4.154 -
   4.155 -lemma "(x::real) <= 1 / 2 --> - x - 2 * x^2 <= - x / (1 - x)"
   4.156 -by (sos_cert "((((A<0 * A<1) * R<1) + ((A<=0 * (A<0 * R<1)) * (R<1 * [x]^2))))")
   4.157 -
   4.158 -lemma "4*r^2 = p^2 - 4*q & r >= (0::real) & x^2 + p*x + q = 0 --> 2*(x::real) = - p + 2*r | 2*x = -p - 2*r"
   4.159 -by (sos_cert "((((((A<0 * A<1) * R<1) + ([~4] * A=0))) & ((((A<0 * A<1) * R<1) + ([4] * A=0)))) & (((((A<0 * A<1) * R<1) + ([4] * A=0))) & ((((A<0 * A<1) * R<1) + ([~4] * A=0)))))")
   4.160 -
   4.161 -end
   4.162 -
     5.1 --- a/src/HOL/Library/Sum_Of_Squares/etc/settings	Sat Jan 08 17:30:05 2011 +0100
     5.2 +++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
     5.3 @@ -1,3 +0,0 @@
     5.4 -# -*- shell-script -*- :mode=shellscript:
     5.5 -
     5.6 -ISABELLE_SUM_OF_SQUARES="$COMPONENT"
     6.1 --- a/src/HOL/Library/Sum_Of_Squares/neos_csdp_client	Sat Jan 08 17:30:05 2011 +0100
     6.2 +++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
     6.3 @@ -1,87 +0,0 @@
     6.4 -#!/usr/bin/env python
     6.5 -import sys
     6.6 -import signal
     6.7 -import xmlrpclib
     6.8 -import time
     6.9 -import re
    6.10 -
    6.11 -# Neos server config
    6.12 -NEOS_HOST="neos.mcs.anl.gov"
    6.13 -NEOS_PORT=3332
    6.14 -
    6.15 -neos=xmlrpclib.Server("http://%s:%d" % (NEOS_HOST, NEOS_PORT))
    6.16 -
    6.17 -jobNumber = 0
    6.18 -password = ""
    6.19 -inputfile = None
    6.20 -outputfile = None
    6.21 -# interrupt handler
    6.22 -def cleanup(signal, frame):
    6.23 -  sys.stdout.write("Caught interrupt, cleaning up\n")
    6.24 -  if jobNumber != 0:
    6.25 -    neos.killJob(jobNumber, password)
    6.26 -  if inputfile != None:
    6.27 -    inputfile.close()
    6.28 -  if outputfile != None:
    6.29 -    outputfile.close()
    6.30 -  sys.exit(21)
    6.31 -
    6.32 -signal.signal(signal.SIGHUP, cleanup)
    6.33 -signal.signal(signal.SIGINT, cleanup)
    6.34 -signal.signal(signal.SIGQUIT, cleanup)
    6.35 -signal.signal(signal.SIGTERM, cleanup)
    6.36 -
    6.37 -if len(sys.argv) <> 3:
    6.38 -  sys.stderr.write("Usage: neos_csdp_client <input_filename> <output_filename>\n")
    6.39 -  sys.exit(19)
    6.40 -
    6.41 -xml_pre = "<document>\n<category>sdp</category>\n<solver>csdp</solver>\n<inputMethod>SPARSE_SDPA</inputMethod>\n<dat><![CDATA["
    6.42 -xml_post = "]]></dat>\n</document>\n"
    6.43 -xml = xml_pre
    6.44 -inputfile = open(sys.argv[1],"r")
    6.45 -buffer = 1
    6.46 -while buffer:
    6.47 -  buffer = inputfile.read()
    6.48 -  xml += buffer
    6.49 -inputfile.close()
    6.50 -xml += xml_post
    6.51 -
    6.52 -(jobNumber,password) = neos.submitJob(xml)
    6.53 -
    6.54 -if jobNumber == 0:
    6.55 -  sys.stdout.write("error submitting job: %s" % password)
    6.56 -  sys.exit(20)
    6.57 -else:
    6.58 -  sys.stdout.write("jobNumber = %d\tpassword = %s\n" % (jobNumber,password))
    6.59 -
    6.60 -offset=0
    6.61 -messages = ""
    6.62 -status="Waiting"
    6.63 -while status == "Running" or status=="Waiting":
    6.64 -  time.sleep(1)
    6.65 -  (msg,offset) = neos.getIntermediateResults(jobNumber,password,offset)
    6.66 -  messages += msg.data
    6.67 -  sys.stdout.write(msg.data)
    6.68 -  status = neos.getJobStatus(jobNumber, password)
    6.69 -
    6.70 -msg = neos.getFinalResults(jobNumber, password).data
    6.71 -sys.stdout.write("---------- Begin CSDP Output -------------\n");
    6.72 -sys.stdout.write(msg)
    6.73 -
    6.74 -# extract solution
    6.75 -result = msg.split("Solution:")
    6.76 -if len(result) > 1:
    6.77 -  solution = result[1].strip()
    6.78 -  if solution != "":
    6.79 -    outputfile = open(sys.argv[2],"w")
    6.80 -    outputfile.write(solution)
    6.81 -    outputfile.close()
    6.82 -
    6.83 -# extract return code
    6.84 -p = re.compile(r"^Error: Command exited with non-zero status (\d+)$", re.MULTILINE)
    6.85 -m = p.search(messages)
    6.86 -if m:
    6.87 -  sys.exit(int(m.group(1)))
    6.88 -else:
    6.89 -  sys.exit(0)
    6.90 -
     7.1 --- a/src/HOL/Library/Sum_Of_Squares/positivstellensatz_tools.ML	Sat Jan 08 17:30:05 2011 +0100
     7.2 +++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
     7.3 @@ -1,155 +0,0 @@
     7.4 -(* Title:      HOL/Library/Sum_Of_Squares/positivstellensatz_tools.ML
     7.5 -   Author:     Philipp Meyer, TU Muenchen
     7.6 -
     7.7 -Functions for generating a certificate from a positivstellensatz and vice versa
     7.8 -*)
     7.9 -
    7.10 -signature POSITIVSTELLENSATZ_TOOLS =
    7.11 -sig
    7.12 -  val pss_tree_to_cert : RealArith.pss_tree -> string
    7.13 -
    7.14 -  val cert_to_pss_tree : Proof.context -> string -> RealArith.pss_tree
    7.15 -end
    7.16 -
    7.17 -
    7.18 -structure PositivstellensatzTools : POSITIVSTELLENSATZ_TOOLS =
    7.19 -struct
    7.20 -
    7.21 -(*** certificate generation ***)
    7.22 -
    7.23 -fun string_of_rat r =
    7.24 -  let
    7.25 -    val (nom, den) = Rat.quotient_of_rat r
    7.26 -  in
    7.27 -    if den = 1 then string_of_int nom
    7.28 -    else string_of_int nom ^ "/" ^ string_of_int den
    7.29 -  end
    7.30 -
    7.31 -(* map polynomials to strings *)
    7.32 -
    7.33 -fun string_of_varpow x k =
    7.34 -  let
    7.35 -    val term = term_of x
    7.36 -    val name = case term of
    7.37 -      Free (n, _) => n
    7.38 -    | _ => error "Term in monomial not free variable"
    7.39 -  in
    7.40 -    if k = 1 then name else name ^ "^" ^ string_of_int k 
    7.41 -  end
    7.42 -
    7.43 -fun string_of_monomial m = 
    7.44 - if FuncUtil.Ctermfunc.is_empty m then "1" 
    7.45 - else 
    7.46 -  let 
    7.47 -   val m' = FuncUtil.dest_monomial m
    7.48 -   val vps = fold_rev (fn (x,k) => cons (string_of_varpow x k)) m' [] 
    7.49 -  in foldr1 (fn (s, t) => s ^ "*" ^ t) vps
    7.50 -  end
    7.51 -
    7.52 -fun string_of_cmonomial (m,c) =
    7.53 -  if FuncUtil.Ctermfunc.is_empty m then string_of_rat c
    7.54 -  else if c = Rat.one then string_of_monomial m
    7.55 -  else (string_of_rat c) ^ "*" ^ (string_of_monomial m);
    7.56 -
    7.57 -fun string_of_poly p = 
    7.58 - if FuncUtil.Monomialfunc.is_empty p then "0" 
    7.59 - else
    7.60 -  let 
    7.61 -   val cms = map string_of_cmonomial
    7.62 -     (sort (prod_ord FuncUtil.monomial_order (K EQUAL)) (FuncUtil.Monomialfunc.dest p))
    7.63 -  in foldr1 (fn (t1, t2) => t1 ^ " + " ^ t2) cms
    7.64 -  end;
    7.65 -
    7.66 -fun pss_to_cert (RealArith.Axiom_eq i) = "A=" ^ string_of_int i
    7.67 -  | pss_to_cert (RealArith.Axiom_le i) = "A<=" ^ string_of_int i
    7.68 -  | pss_to_cert (RealArith.Axiom_lt i) = "A<" ^ string_of_int i
    7.69 -  | pss_to_cert (RealArith.Rational_eq r) = "R=" ^ string_of_rat r
    7.70 -  | pss_to_cert (RealArith.Rational_le r) = "R<=" ^ string_of_rat r
    7.71 -  | pss_to_cert (RealArith.Rational_lt r) = "R<" ^ string_of_rat r
    7.72 -  | pss_to_cert (RealArith.Square p) = "[" ^ string_of_poly p ^ "]^2"
    7.73 -  | pss_to_cert (RealArith.Eqmul (p, pss)) = "([" ^ string_of_poly p ^ "] * " ^ pss_to_cert pss ^ ")"
    7.74 -  | pss_to_cert (RealArith.Sum (pss1, pss2)) = "(" ^ pss_to_cert pss1 ^ " + " ^ pss_to_cert pss2 ^ ")"
    7.75 -  | pss_to_cert (RealArith.Product (pss1, pss2)) = "(" ^ pss_to_cert pss1 ^ " * " ^ pss_to_cert pss2 ^ ")"
    7.76 -
    7.77 -fun pss_tree_to_cert RealArith.Trivial = "()"
    7.78 -  | pss_tree_to_cert (RealArith.Cert pss) = "(" ^ pss_to_cert pss ^ ")"
    7.79 -  | pss_tree_to_cert (RealArith.Branch (t1, t2)) = "(" ^ pss_tree_to_cert t1 ^ " & " ^ pss_tree_to_cert t2 ^ ")"
    7.80 -
    7.81 -(*** certificate parsing ***)
    7.82 -
    7.83 -(* basic parser *)
    7.84 -
    7.85 -val str = Scan.this_string
    7.86 -
    7.87 -val number = Scan.repeat1 (Scan.one Symbol.is_ascii_digit >>
    7.88 -  (fn s => ord s - ord "0")) >>
    7.89 -  foldl1 (fn (n, d) => n * 10 + d)
    7.90 -
    7.91 -val nat = number
    7.92 -val int = Scan.optional (str "~" >> K ~1) 1 -- nat >> op *;
    7.93 -val rat = int --| str "/" -- int >> Rat.rat_of_quotient
    7.94 -val rat_int = rat || int >> Rat.rat_of_int
    7.95 -
    7.96 -(* polynomial parser *)
    7.97 -
    7.98 -fun repeat_sep s f = f ::: Scan.repeat (str s |-- f)
    7.99 -
   7.100 -val parse_id = Scan.one Symbol.is_letter ::: Scan.many Symbol.is_letdig >> implode
   7.101 -
   7.102 -fun parse_varpow ctxt = parse_id -- Scan.optional (str "^" |-- nat) 1 >>
   7.103 -  (fn (x, k) => (cterm_of (ProofContext.theory_of ctxt) (Free (x, @{typ real})), k)) 
   7.104 -
   7.105 -fun parse_monomial ctxt = repeat_sep "*" (parse_varpow ctxt) >>
   7.106 -  (fn xs => fold FuncUtil.Ctermfunc.update xs FuncUtil.Ctermfunc.empty)
   7.107 -
   7.108 -fun parse_cmonomial ctxt =
   7.109 -  rat_int --| str "*" -- (parse_monomial ctxt) >> swap ||
   7.110 -  (parse_monomial ctxt) >> (fn m => (m, Rat.one)) ||
   7.111 -  rat_int >> (fn r => (FuncUtil.Ctermfunc.empty, r))
   7.112 -
   7.113 -fun parse_poly ctxt = repeat_sep "+" (parse_cmonomial ctxt) >>
   7.114 -  (fn xs => fold FuncUtil.Monomialfunc.update xs FuncUtil.Monomialfunc.empty)
   7.115 -
   7.116 -(* positivstellensatz parser *)
   7.117 -
   7.118 -val parse_axiom =
   7.119 -  (str "A=" |-- int >> RealArith.Axiom_eq) ||
   7.120 -  (str "A<=" |-- int >> RealArith.Axiom_le) ||
   7.121 -  (str "A<" |-- int >> RealArith.Axiom_lt)
   7.122 -
   7.123 -val parse_rational =
   7.124 -  (str "R=" |-- rat_int >> RealArith.Rational_eq) ||
   7.125 -  (str "R<=" |-- rat_int >> RealArith.Rational_le) ||
   7.126 -  (str "R<" |-- rat_int >> RealArith.Rational_lt)
   7.127 -
   7.128 -fun parse_cert ctxt input =
   7.129 -  let
   7.130 -    val pc = parse_cert ctxt
   7.131 -    val pp = parse_poly ctxt
   7.132 -  in
   7.133 -  (parse_axiom ||
   7.134 -   parse_rational ||
   7.135 -   str "[" |-- pp --| str "]^2" >> RealArith.Square ||
   7.136 -   str "([" |-- pp --| str "]*" -- pc --| str ")" >> RealArith.Eqmul ||
   7.137 -   str "(" |-- pc --| str "*" -- pc --| str ")" >> RealArith.Product ||
   7.138 -   str "(" |-- pc --| str "+" -- pc --| str ")" >> RealArith.Sum) input
   7.139 -  end
   7.140 -
   7.141 -fun parse_cert_tree ctxt input =
   7.142 -  let
   7.143 -    val pc = parse_cert ctxt
   7.144 -    val pt = parse_cert_tree ctxt
   7.145 -  in
   7.146 -  (str "()" >> K RealArith.Trivial ||
   7.147 -   str "(" |-- pc --| str ")" >> RealArith.Cert ||
   7.148 -   str "(" |-- pt --| str "&" -- pt --| str ")" >> RealArith.Branch) input
   7.149 -  end
   7.150 -
   7.151 -(* scanner *)
   7.152 -
   7.153 -fun cert_to_pss_tree ctxt input_str = Symbol.scanner "bad certificate" (parse_cert_tree ctxt)
   7.154 -  (filter_out Symbol.is_blank (Symbol.explode input_str))
   7.155 -
   7.156 -end
   7.157 -
   7.158 -
     8.1 --- a/src/HOL/Library/Sum_Of_Squares/sos_wrapper.ML	Sat Jan 08 17:30:05 2011 +0100
     8.2 +++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
     8.3 @@ -1,159 +0,0 @@
     8.4 -(*  Title:      HOL/Library/Sum_Of_Squares/sos_wrapper.ML
     8.5 -    Author:     Philipp Meyer, TU Muenchen
     8.6 -
     8.7 -Added functionality for sums of squares, e.g. calling a remote prover.
     8.8 -*)
     8.9 -
    8.10 -signature SOS_WRAPPER =
    8.11 -sig
    8.12 -  datatype prover_result = Success | Failure | Error
    8.13 -  val setup: theory -> theory
    8.14 -  val dest_dir: string Config.T
    8.15 -  val prover_name: string Config.T
    8.16 -end
    8.17 -
    8.18 -structure SOS_Wrapper: SOS_WRAPPER =
    8.19 -struct
    8.20 -
    8.21 -datatype prover_result = Success | Failure | Error
    8.22 -
    8.23 -fun str_of_result Success = "Success"
    8.24 -  | str_of_result Failure = "Failure"
    8.25 -  | str_of_result Error = "Error"
    8.26 -
    8.27 -
    8.28 -(*** calling provers ***)
    8.29 -
    8.30 -val (dest_dir, setup_dest_dir) = Attrib.config_string "sos_dest_dir" (K "")
    8.31 -
    8.32 -fun filename dir name =
    8.33 -  let
    8.34 -    val probfile = Path.basic (name ^ serial_string ())
    8.35 -    val dir_path = Path.explode dir
    8.36 -  in
    8.37 -    if dir = "" then
    8.38 -      File.tmp_path probfile
    8.39 -    else if File.exists dir_path then
    8.40 -      Path.append dir_path probfile
    8.41 -    else error ("No such directory: " ^ dir)
    8.42 -  end
    8.43 -
    8.44 -fun run_solver ctxt name cmd find_failure input =
    8.45 -  let
    8.46 -    val _ = warning ("Calling solver: " ^ name)
    8.47 -
    8.48 -    (* create input file *)
    8.49 -    val dir = Config.get ctxt dest_dir
    8.50 -    val input_file = filename dir "sos_in"
    8.51 -    val _ = File.write input_file input
    8.52 -
    8.53 -    (* call solver *)
    8.54 -    val output_file = filename dir "sos_out"
    8.55 -    val (output, rv) =
    8.56 -      bash_output
    8.57 -       (if File.exists cmd then
    8.58 -          space_implode " "
    8.59 -            [File.shell_path cmd, File.shell_path input_file, File.shell_path output_file]
    8.60 -        else error ("Bad executable: " ^ File.platform_path cmd))
    8.61 -
    8.62 -    (* read and analyze output *)
    8.63 -    val (res, res_msg) = find_failure rv
    8.64 -    val result = if File.exists output_file then File.read output_file else ""
    8.65 -
    8.66 -    (* remove temporary files *)
    8.67 -    val _ =
    8.68 -      if dir = "" then
    8.69 -        (File.rm input_file; if File.exists output_file then File.rm output_file else ())
    8.70 -      else ()
    8.71 -
    8.72 -    val _ =
    8.73 -      if Config.get ctxt Sum_Of_Squares.trace
    8.74 -      then writeln ("Solver output:\n" ^ output)
    8.75 -      else ()
    8.76 -
    8.77 -    val _ = warning (str_of_result res ^ ": " ^ res_msg)
    8.78 -  in
    8.79 -    (case res of
    8.80 -      Success => result
    8.81 -    | Failure => raise Sum_Of_Squares.Failure res_msg
    8.82 -    | Error => error ("Prover failed: " ^ res_msg))
    8.83 -  end
    8.84 -
    8.85 -
    8.86 -(*** various provers ***)
    8.87 -
    8.88 -(* local csdp client *)
    8.89 -
    8.90 -fun find_csdp_failure rv =
    8.91 -  case rv of
    8.92 -    0 => (Success, "SDP solved")
    8.93 -  | 1 => (Failure, "SDP is primal infeasible")
    8.94 -  | 2 => (Failure, "SDP is dual infeasible")
    8.95 -  | 3 => (Success, "SDP solved with reduced accuracy")
    8.96 -  | 4 => (Failure, "Maximum iterations reached")
    8.97 -  | 5 => (Failure, "Stuck at edge of primal feasibility")
    8.98 -  | 6 => (Failure, "Stuck at edge of dual infeasibility")
    8.99 -  | 7 => (Failure, "Lack of progress")
   8.100 -  | 8 => (Failure, "X, Z, or O was singular")
   8.101 -  | 9 => (Failure, "Detected NaN or Inf values")
   8.102 -  | _ => (Error, "return code is " ^ string_of_int rv)
   8.103 -
   8.104 -val csdp = ("$CSDP_EXE", find_csdp_failure)
   8.105 -
   8.106 -
   8.107 -(* remote neos server *)
   8.108 -
   8.109 -fun find_neos_failure rv =
   8.110 -  case rv of
   8.111 -    20 => (Error, "error submitting job")
   8.112 -  | 21 => (Error, "interrupt")
   8.113 -  |  _ => find_csdp_failure rv
   8.114 -
   8.115 -val neos_csdp = ("$ISABELLE_SUM_OF_SQUARES/neos_csdp_client", find_neos_failure)
   8.116 -
   8.117 -
   8.118 -(* named provers *)
   8.119 -
   8.120 -fun prover "remote_csdp" = neos_csdp
   8.121 -  | prover "csdp" = csdp
   8.122 -  | prover name = error ("Unknown prover: " ^ name)
   8.123 -
   8.124 -val (prover_name, setup_prover_name) = Attrib.config_string "sos_prover_name" (K "remote_csdp")
   8.125 -
   8.126 -fun call_solver ctxt opt_name =
   8.127 -  let
   8.128 -    val name = the_default (Config.get ctxt prover_name) opt_name
   8.129 -    val (cmd, find_failure) = prover name
   8.130 -  in run_solver ctxt name (Path.explode cmd) find_failure end
   8.131 -
   8.132 -
   8.133 -(* certificate output *)
   8.134 -
   8.135 -fun output_line cert =
   8.136 -  "To repeat this proof with a certifiate use this command:\n" ^
   8.137 -    Markup.markup Markup.sendback ("by (sos_cert \"" ^ cert ^ "\")")
   8.138 -
   8.139 -val print_cert = warning o output_line o PositivstellensatzTools.pss_tree_to_cert
   8.140 -
   8.141 -
   8.142 -(* method setup *)
   8.143 -
   8.144 -fun sos_solver print method = SIMPLE_METHOD' o Sum_Of_Squares.sos_tac print method
   8.145 -
   8.146 -val setup =
   8.147 -  setup_dest_dir #>
   8.148 -  setup_prover_name #>
   8.149 -  Method.setup @{binding sos}
   8.150 -    (Scan.lift (Scan.option Parse.xname)
   8.151 -      >> (fn opt_name => fn ctxt =>
   8.152 -        sos_solver print_cert
   8.153 -          (Sum_Of_Squares.Prover (call_solver ctxt opt_name)) ctxt))
   8.154 -    "prove universal problems over the reals using sums of squares" #>
   8.155 -  Method.setup @{binding sos_cert}
   8.156 -    (Scan.lift Parse.string
   8.157 -      >> (fn cert => fn ctxt =>
   8.158 -        sos_solver ignore
   8.159 -          (Sum_Of_Squares.Certificate (PositivstellensatzTools.cert_to_pss_tree ctxt cert)) ctxt))
   8.160 -    "prove universal problems over the reals using sums of squares with certificates"
   8.161 -
   8.162 -end
     9.1 --- a/src/HOL/Library/Sum_Of_Squares/sum_of_squares.ML	Sat Jan 08 17:30:05 2011 +0100
     9.2 +++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
     9.3 @@ -1,1435 +0,0 @@
     9.4 -(*  Title:      HOL/Library/Sum_Of_Squares/sum_of_squares.ML
     9.5 -    Author:     Amine Chaieb, University of Cambridge
     9.6 -    Author:     Philipp Meyer, TU Muenchen
     9.7 -
     9.8 -A tactic for proving nonlinear inequalities.
     9.9 -*)
    9.10 -
    9.11 -signature SUM_OF_SQUARES =
    9.12 -sig
    9.13 -  datatype proof_method = Certificate of RealArith.pss_tree | Prover of string -> string
    9.14 -  val sos_tac: (RealArith.pss_tree -> unit) -> proof_method -> Proof.context -> int -> tactic
    9.15 -  val trace: bool Config.T
    9.16 -  val setup: theory -> theory
    9.17 -  exception Failure of string;
    9.18 -end
    9.19 -
    9.20 -structure Sum_Of_Squares: SUM_OF_SQUARES =
    9.21 -struct
    9.22 -
    9.23 -val rat_0 = Rat.zero;
    9.24 -val rat_1 = Rat.one;
    9.25 -val rat_2 = Rat.two;
    9.26 -val rat_10 = Rat.rat_of_int 10;
    9.27 -val rat_1_2 = rat_1 // rat_2;
    9.28 -val max = Integer.max;
    9.29 -val min = Integer.min;
    9.30 -
    9.31 -val denominator_rat = Rat.quotient_of_rat #> snd #> Rat.rat_of_int;
    9.32 -val numerator_rat = Rat.quotient_of_rat #> fst #> Rat.rat_of_int;
    9.33 -fun int_of_rat a =
    9.34 -    case Rat.quotient_of_rat a of (i,1) => i | _ => error "int_of_rat: not an int";
    9.35 -fun lcm_rat x y = Rat.rat_of_int (Integer.lcm (int_of_rat x) (int_of_rat y));
    9.36 -
    9.37 -fun rat_pow r i =
    9.38 - let fun pow r i =
    9.39 -   if i = 0 then rat_1 else
    9.40 -   let val d = pow r (i div 2)
    9.41 -   in d */ d */ (if i mod 2 = 0 then rat_1 else r)
    9.42 -   end
    9.43 - in if i < 0 then pow (Rat.inv r) (~ i) else pow r i end;
    9.44 -
    9.45 -fun round_rat r =
    9.46 - let val (a,b) = Rat.quotient_of_rat (Rat.abs r)
    9.47 -     val d = a div b
    9.48 -     val s = if r </ rat_0 then (Rat.neg o Rat.rat_of_int) else Rat.rat_of_int
    9.49 -     val x2 = 2 * (a - (b * d))
    9.50 - in s (if x2 >= b then d + 1 else d) end
    9.51 -
    9.52 -val abs_rat = Rat.abs;
    9.53 -val pow2 = rat_pow rat_2;
    9.54 -val pow10 = rat_pow rat_10;
    9.55 -
    9.56 -val (trace, setup_trace) = Attrib.config_bool "sos_trace" (K false);
    9.57 -val setup = setup_trace;
    9.58 -
    9.59 -exception Sanity;
    9.60 -
    9.61 -exception Unsolvable;
    9.62 -
    9.63 -exception Failure of string;
    9.64 -
    9.65 -datatype proof_method =
    9.66 -    Certificate of RealArith.pss_tree
    9.67 -  | Prover of (string -> string)
    9.68 -
    9.69 -(* Turn a rational into a decimal string with d sig digits.                  *)
    9.70 -
    9.71 -local
    9.72 -fun normalize y =
    9.73 -  if abs_rat y </ (rat_1 // rat_10) then normalize (rat_10 */ y) - 1
    9.74 -  else if abs_rat y >=/ rat_1 then normalize (y // rat_10) + 1
    9.75 -  else 0
    9.76 - in
    9.77 -fun decimalize d x =
    9.78 -  if x =/ rat_0 then "0.0" else
    9.79 -  let
    9.80 -   val y = Rat.abs x
    9.81 -   val e = normalize y
    9.82 -   val z = pow10(~ e) */ y +/ rat_1
    9.83 -   val k = int_of_rat (round_rat(pow10 d */ z))
    9.84 -  in (if x </ rat_0 then "-0." else "0.") ^
    9.85 -     implode(tl(raw_explode(string_of_int k))) ^
    9.86 -     (if e = 0 then "" else "e"^string_of_int e)
    9.87 -  end
    9.88 -end;
    9.89 -
    9.90 -(* Iterations over numbers, and lists indexed by numbers.                    *)
    9.91 -
    9.92 -fun itern k l f a =
    9.93 -  case l of
    9.94 -    [] => a
    9.95 -  | h::t => itern (k + 1) t f (f h k a);
    9.96 -
    9.97 -fun iter (m,n) f a =
    9.98 -  if n < m then a
    9.99 -  else iter (m+1,n) f (f m a);
   9.100 -
   9.101 -(* The main types.                                                           *)
   9.102 -
   9.103 -type vector = int* Rat.rat FuncUtil.Intfunc.table;
   9.104 -
   9.105 -type matrix = (int*int)*(Rat.rat FuncUtil.Intpairfunc.table);
   9.106 -
   9.107 -fun iszero (k,r) = r =/ rat_0;
   9.108 -
   9.109 -
   9.110 -(* Vectors. Conventionally indexed 1..n.                                     *)
   9.111 -
   9.112 -fun vector_0 n = (n,FuncUtil.Intfunc.empty):vector;
   9.113 -
   9.114 -fun dim (v:vector) = fst v;
   9.115 -
   9.116 -fun vector_const c n =
   9.117 -  if c =/ rat_0 then vector_0 n
   9.118 -  else (n,fold_rev (fn k => FuncUtil.Intfunc.update (k,c)) (1 upto n) FuncUtil.Intfunc.empty) :vector;
   9.119 -
   9.120 -val vector_1 = vector_const rat_1;
   9.121 -
   9.122 -fun vector_cmul c (v:vector) =
   9.123 - let val n = dim v
   9.124 - in if c =/ rat_0 then vector_0 n
   9.125 -    else (n,FuncUtil.Intfunc.map (fn _ => fn x => c */ x) (snd v))
   9.126 - end;
   9.127 -
   9.128 -fun vector_neg (v:vector) = (fst v,FuncUtil.Intfunc.map (K Rat.neg) (snd v)) :vector;
   9.129 -
   9.130 -fun vector_add (v1:vector) (v2:vector) =
   9.131 - let val m = dim v1
   9.132 -     val n = dim v2
   9.133 - in if m <> n then error "vector_add: incompatible dimensions"
   9.134 -    else (n,FuncUtil.Intfunc.combine (curry op +/) (fn x => x =/ rat_0) (snd v1) (snd v2)) :vector
   9.135 - end;
   9.136 -
   9.137 -fun vector_sub v1 v2 = vector_add v1 (vector_neg v2);
   9.138 -
   9.139 -fun vector_dot (v1:vector) (v2:vector) =
   9.140 - let val m = dim v1
   9.141 -     val n = dim v2
   9.142 - in if m <> n then error "vector_dot: incompatible dimensions"
   9.143 -    else FuncUtil.Intfunc.fold (fn (i,x) => fn a => x +/ a)
   9.144 -        (FuncUtil.Intfunc.combine (curry op */) (fn x => x =/ rat_0) (snd v1) (snd v2)) rat_0
   9.145 - end;
   9.146 -
   9.147 -fun vector_of_list l =
   9.148 - let val n = length l
   9.149 - in (n,fold_rev2 (curry FuncUtil.Intfunc.update) (1 upto n) l FuncUtil.Intfunc.empty) :vector
   9.150 - end;
   9.151 -
   9.152 -(* Matrices; again rows and columns indexed from 1.                          *)
   9.153 -
   9.154 -fun matrix_0 (m,n) = ((m,n),FuncUtil.Intpairfunc.empty):matrix;
   9.155 -
   9.156 -fun dimensions (m:matrix) = fst m;
   9.157 -
   9.158 -fun matrix_const c (mn as (m,n)) =
   9.159 -  if m <> n then error "matrix_const: needs to be square"
   9.160 -  else if c =/ rat_0 then matrix_0 mn
   9.161 -  else (mn,fold_rev (fn k => FuncUtil.Intpairfunc.update ((k,k), c)) (1 upto n) FuncUtil.Intpairfunc.empty) :matrix;;
   9.162 -
   9.163 -val matrix_1 = matrix_const rat_1;
   9.164 -
   9.165 -fun matrix_cmul c (m:matrix) =
   9.166 - let val (i,j) = dimensions m
   9.167 - in if c =/ rat_0 then matrix_0 (i,j)
   9.168 -    else ((i,j),FuncUtil.Intpairfunc.map (fn _ => fn x => c */ x) (snd m))
   9.169 - end;
   9.170 -
   9.171 -fun matrix_neg (m:matrix) =
   9.172 -  (dimensions m, FuncUtil.Intpairfunc.map (K Rat.neg) (snd m)) :matrix;
   9.173 -
   9.174 -fun matrix_add (m1:matrix) (m2:matrix) =
   9.175 - let val d1 = dimensions m1
   9.176 -     val d2 = dimensions m2
   9.177 - in if d1 <> d2
   9.178 -     then error "matrix_add: incompatible dimensions"
   9.179 -    else (d1,FuncUtil.Intpairfunc.combine (curry op +/) (fn x => x =/ rat_0) (snd m1) (snd m2)) :matrix
   9.180 - end;;
   9.181 -
   9.182 -fun matrix_sub m1 m2 = matrix_add m1 (matrix_neg m2);
   9.183 -
   9.184 -fun row k (m:matrix) =
   9.185 - let val (i,j) = dimensions m
   9.186 - in (j,
   9.187 -   FuncUtil.Intpairfunc.fold (fn ((i,j), c) => fn a => if i = k then FuncUtil.Intfunc.update (j,c) a else a) (snd m) FuncUtil.Intfunc.empty ) : vector
   9.188 - end;
   9.189 -
   9.190 -fun column k (m:matrix) =
   9.191 -  let val (i,j) = dimensions m
   9.192 -  in (i,
   9.193 -   FuncUtil.Intpairfunc.fold (fn ((i,j), c) => fn a => if j = k then FuncUtil.Intfunc.update (i,c) a else a) (snd m)  FuncUtil.Intfunc.empty)
   9.194 -   : vector
   9.195 - end;
   9.196 -
   9.197 -fun transp (m:matrix) =
   9.198 -  let val (i,j) = dimensions m
   9.199 -  in
   9.200 -  ((j,i),FuncUtil.Intpairfunc.fold (fn ((i,j), c) => fn a => FuncUtil.Intpairfunc.update ((j,i), c) a) (snd m) FuncUtil.Intpairfunc.empty) :matrix
   9.201 - end;
   9.202 -
   9.203 -fun diagonal (v:vector) =
   9.204 - let val n = dim v
   9.205 - in ((n,n),FuncUtil.Intfunc.fold (fn (i, c) => fn a => FuncUtil.Intpairfunc.update ((i,i), c) a) (snd v) FuncUtil.Intpairfunc.empty) : matrix
   9.206 - end;
   9.207 -
   9.208 -fun matrix_of_list l =
   9.209 - let val m = length l
   9.210 - in if m = 0 then matrix_0 (0,0) else
   9.211 -   let val n = length (hd l)
   9.212 -   in ((m,n),itern 1 l (fn v => fn i => itern 1 v (fn c => fn j => FuncUtil.Intpairfunc.update ((i,j), c))) FuncUtil.Intpairfunc.empty)
   9.213 -   end
   9.214 - end;
   9.215 -
   9.216 -(* Monomials.                                                                *)
   9.217 -
   9.218 -fun monomial_eval assig m =
   9.219 -  FuncUtil.Ctermfunc.fold (fn (x, k) => fn a => a */ rat_pow (FuncUtil.Ctermfunc.apply assig x) k)
   9.220 -        m rat_1;
   9.221 -val monomial_1 = FuncUtil.Ctermfunc.empty;
   9.222 -
   9.223 -fun monomial_var x = FuncUtil.Ctermfunc.onefunc (x, 1);
   9.224 -
   9.225 -val monomial_mul =
   9.226 -  FuncUtil.Ctermfunc.combine Integer.add (K false);
   9.227 -
   9.228 -fun monomial_pow m k =
   9.229 -  if k = 0 then monomial_1
   9.230 -  else FuncUtil.Ctermfunc.map (fn _ => fn x => k * x) m;
   9.231 -
   9.232 -fun monomial_divides m1 m2 =
   9.233 -  FuncUtil.Ctermfunc.fold (fn (x, k) => fn a => FuncUtil.Ctermfunc.tryapplyd m2 x 0 >= k andalso a) m1 true;;
   9.234 -
   9.235 -fun monomial_div m1 m2 =
   9.236 - let val m = FuncUtil.Ctermfunc.combine Integer.add
   9.237 -   (fn x => x = 0) m1 (FuncUtil.Ctermfunc.map (fn _ => fn x => ~ x) m2)
   9.238 - in if FuncUtil.Ctermfunc.fold (fn (x, k) => fn a => k >= 0 andalso a) m true then m
   9.239 -  else error "monomial_div: non-divisible"
   9.240 - end;
   9.241 -
   9.242 -fun monomial_degree x m =
   9.243 -  FuncUtil.Ctermfunc.tryapplyd m x 0;;
   9.244 -
   9.245 -fun monomial_lcm m1 m2 =
   9.246 -  fold_rev (fn x => FuncUtil.Ctermfunc.update (x, max (monomial_degree x m1) (monomial_degree x m2)))
   9.247 -          (union (is_equal o FuncUtil.cterm_ord) (FuncUtil.Ctermfunc.dom m1) (FuncUtil.Ctermfunc.dom m2)) (FuncUtil.Ctermfunc.empty);
   9.248 -
   9.249 -fun monomial_multidegree m =
   9.250 - FuncUtil.Ctermfunc.fold (fn (x, k) => fn a => k + a) m 0;;
   9.251 -
   9.252 -fun monomial_variables m = FuncUtil.Ctermfunc.dom m;;
   9.253 -
   9.254 -(* Polynomials.                                                              *)
   9.255 -
   9.256 -fun eval assig p =
   9.257 -  FuncUtil.Monomialfunc.fold (fn (m, c) => fn a => a +/ c */ monomial_eval assig m) p rat_0;
   9.258 -
   9.259 -val poly_0 = FuncUtil.Monomialfunc.empty;
   9.260 -
   9.261 -fun poly_isconst p =
   9.262 -  FuncUtil.Monomialfunc.fold (fn (m, c) => fn a => FuncUtil.Ctermfunc.is_empty m andalso a) p true;
   9.263 -
   9.264 -fun poly_var x = FuncUtil.Monomialfunc.onefunc (monomial_var x,rat_1);
   9.265 -
   9.266 -fun poly_const c =
   9.267 -  if c =/ rat_0 then poly_0 else FuncUtil.Monomialfunc.onefunc(monomial_1, c);
   9.268 -
   9.269 -fun poly_cmul c p =
   9.270 -  if c =/ rat_0 then poly_0
   9.271 -  else FuncUtil.Monomialfunc.map (fn _ => fn x => c */ x) p;
   9.272 -
   9.273 -fun poly_neg p = FuncUtil.Monomialfunc.map (K Rat.neg) p;;
   9.274 -
   9.275 -fun poly_add p1 p2 =
   9.276 -  FuncUtil.Monomialfunc.combine (curry op +/) (fn x => x =/ rat_0) p1 p2;
   9.277 -
   9.278 -fun poly_sub p1 p2 = poly_add p1 (poly_neg p2);
   9.279 -
   9.280 -fun poly_cmmul (c,m) p =
   9.281 - if c =/ rat_0 then poly_0
   9.282 - else if FuncUtil.Ctermfunc.is_empty m
   9.283 -      then FuncUtil.Monomialfunc.map (fn _ => fn d => c */ d) p
   9.284 -      else FuncUtil.Monomialfunc.fold (fn (m', d) => fn a => (FuncUtil.Monomialfunc.update (monomial_mul m m', c */ d) a)) p poly_0;
   9.285 -
   9.286 -fun poly_mul p1 p2 =
   9.287 -  FuncUtil.Monomialfunc.fold (fn (m, c) => fn a => poly_add (poly_cmmul (c,m) p2) a) p1 poly_0;
   9.288 -
   9.289 -fun poly_div p1 p2 =
   9.290 - if not(poly_isconst p2)
   9.291 - then error "poly_div: non-constant" else
   9.292 - let val c = eval FuncUtil.Ctermfunc.empty p2
   9.293 - in if c =/ rat_0 then error "poly_div: division by zero"
   9.294 -    else poly_cmul (Rat.inv c) p1
   9.295 - end;
   9.296 -
   9.297 -fun poly_square p = poly_mul p p;
   9.298 -
   9.299 -fun poly_pow p k =
   9.300 - if k = 0 then poly_const rat_1
   9.301 - else if k = 1 then p
   9.302 - else let val q = poly_square(poly_pow p (k div 2)) in
   9.303 -      if k mod 2 = 1 then poly_mul p q else q end;
   9.304 -
   9.305 -fun poly_exp p1 p2 =
   9.306 -  if not(poly_isconst p2)
   9.307 -  then error "poly_exp: not a constant"
   9.308 -  else poly_pow p1 (int_of_rat (eval FuncUtil.Ctermfunc.empty p2));
   9.309 -
   9.310 -fun degree x p =
   9.311 - FuncUtil.Monomialfunc.fold (fn (m,c) => fn a => max (monomial_degree x m) a) p 0;
   9.312 -
   9.313 -fun multidegree p =
   9.314 -  FuncUtil.Monomialfunc.fold (fn (m, c) => fn a => max (monomial_multidegree m) a) p 0;
   9.315 -
   9.316 -fun poly_variables p =
   9.317 -  sort FuncUtil.cterm_ord (FuncUtil.Monomialfunc.fold_rev (fn (m, c) => union (is_equal o FuncUtil.cterm_ord) (monomial_variables m)) p []);;
   9.318 -
   9.319 -(* Order monomials for human presentation.                                   *)
   9.320 -
   9.321 -val humanorder_varpow = prod_ord FuncUtil.cterm_ord (rev_order o int_ord);
   9.322 -
   9.323 -local
   9.324 - fun ord (l1,l2) = case (l1,l2) of
   9.325 -  (_,[]) => LESS
   9.326 - | ([],_) => GREATER
   9.327 - | (h1::t1,h2::t2) =>
   9.328 -   (case humanorder_varpow (h1, h2) of
   9.329 -     LESS => LESS
   9.330 -   | EQUAL => ord (t1,t2)
   9.331 -   | GREATER => GREATER)
   9.332 -in fun humanorder_monomial m1 m2 =
   9.333 - ord (sort humanorder_varpow (FuncUtil.Ctermfunc.dest m1),
   9.334 -  sort humanorder_varpow (FuncUtil.Ctermfunc.dest m2))
   9.335 -end;
   9.336 -
   9.337 -(* Conversions to strings.                                                   *)
   9.338 -
   9.339 -fun string_of_vector min_size max_size (v:vector) =
   9.340 - let val n_raw = dim v
   9.341 - in if n_raw = 0 then "[]" else
   9.342 -  let
   9.343 -   val n = max min_size (min n_raw max_size)
   9.344 -   val xs = map (Rat.string_of_rat o (fn i => FuncUtil.Intfunc.tryapplyd (snd v) i rat_0)) (1 upto n)
   9.345 -  in "[" ^ space_implode ", " xs ^
   9.346 -  (if n_raw > max_size then ", ...]" else "]")
   9.347 -  end
   9.348 - end;
   9.349 -
   9.350 -fun string_of_matrix max_size (m:matrix) =
   9.351 - let
   9.352 -  val (i_raw,j_raw) = dimensions m
   9.353 -  val i = min max_size i_raw
   9.354 -  val j = min max_size j_raw
   9.355 -  val rstr = map (fn k => string_of_vector j j (row k m)) (1 upto i)
   9.356 - in "["^ space_implode ";\n " rstr ^
   9.357 -  (if j > max_size then "\n ...]" else "]")
   9.358 - end;
   9.359 -
   9.360 -fun string_of_term t =
   9.361 - case t of
   9.362 -   a$b => "("^(string_of_term a)^" "^(string_of_term b)^")"
   9.363 - | Abs x =>
   9.364 -    let val (xn, b) = Term.dest_abs x
   9.365 -    in "(\\"^xn^"."^(string_of_term b)^")"
   9.366 -    end
   9.367 - | Const(s,_) => s
   9.368 - | Free (s,_) => s
   9.369 - | Var((s,_),_) => s
   9.370 - | _ => error "string_of_term";
   9.371 -
   9.372 -val string_of_cterm = string_of_term o term_of;
   9.373 -
   9.374 -fun string_of_varpow x k =
   9.375 -  if k = 1 then string_of_cterm x
   9.376 -  else string_of_cterm x^"^"^string_of_int k;
   9.377 -
   9.378 -fun string_of_monomial m =
   9.379 - if FuncUtil.Ctermfunc.is_empty m then "1" else
   9.380 - let val vps = fold_rev (fn (x,k) => fn a => string_of_varpow x k :: a)
   9.381 -  (sort humanorder_varpow (FuncUtil.Ctermfunc.dest m)) []
   9.382 - in space_implode "*" vps
   9.383 - end;
   9.384 -
   9.385 -fun string_of_cmonomial (c,m) =
   9.386 - if FuncUtil.Ctermfunc.is_empty m then Rat.string_of_rat c
   9.387 - else if c =/ rat_1 then string_of_monomial m
   9.388 - else Rat.string_of_rat c ^ "*" ^ string_of_monomial m;;
   9.389 -
   9.390 -fun string_of_poly p =
   9.391 - if FuncUtil.Monomialfunc.is_empty p then "<<0>>" else
   9.392 - let
   9.393 -  val cms = sort (fn ((m1,_),(m2,_)) => humanorder_monomial m1  m2) (FuncUtil.Monomialfunc.dest p)
   9.394 -  val s = fold (fn (m,c) => fn a =>
   9.395 -             if c </ rat_0 then a ^ " - " ^ string_of_cmonomial(Rat.neg c,m)
   9.396 -             else a ^ " + " ^ string_of_cmonomial(c,m))
   9.397 -          cms ""
   9.398 -  val s1 = String.substring (s, 0, 3)
   9.399 -  val s2 = String.substring (s, 3, String.size s - 3)
   9.400 - in "<<" ^(if s1 = " + " then s2 else "-"^s2)^">>"
   9.401 - end;
   9.402 -
   9.403 -(* Conversion from HOL term.                                                 *)
   9.404 -
   9.405 -local
   9.406 - val neg_tm = @{cterm "uminus :: real => _"}
   9.407 - val add_tm = @{cterm "op + :: real => _"}
   9.408 - val sub_tm = @{cterm "op - :: real => _"}
   9.409 - val mul_tm = @{cterm "op * :: real => _"}
   9.410 - val inv_tm = @{cterm "inverse :: real => _"}
   9.411 - val div_tm = @{cterm "op / :: real => _"}
   9.412 - val pow_tm = @{cterm "op ^ :: real => _"}
   9.413 - val zero_tm = @{cterm "0:: real"}
   9.414 - val is_numeral = can (HOLogic.dest_number o term_of)
   9.415 - fun is_comb t = case t of _$_ => true | _ => false
   9.416 - fun poly_of_term tm =
   9.417 -  if tm aconvc zero_tm then poly_0
   9.418 -  else if RealArith.is_ratconst tm
   9.419 -       then poly_const(RealArith.dest_ratconst tm)
   9.420 -  else
   9.421 -  (let val (lop,r) = Thm.dest_comb tm
   9.422 -   in if lop aconvc neg_tm then poly_neg(poly_of_term r)
   9.423 -      else if lop aconvc inv_tm then
   9.424 -       let val p = poly_of_term r
   9.425 -       in if poly_isconst p
   9.426 -          then poly_const(Rat.inv (eval FuncUtil.Ctermfunc.empty p))
   9.427 -          else error "poly_of_term: inverse of non-constant polyomial"
   9.428 -       end
   9.429 -   else (let val (opr,l) = Thm.dest_comb lop
   9.430 -         in
   9.431 -          if opr aconvc pow_tm andalso is_numeral r
   9.432 -          then poly_pow (poly_of_term l) ((snd o HOLogic.dest_number o term_of) r)
   9.433 -          else if opr aconvc add_tm
   9.434 -           then poly_add (poly_of_term l) (poly_of_term r)
   9.435 -          else if opr aconvc sub_tm
   9.436 -           then poly_sub (poly_of_term l) (poly_of_term r)
   9.437 -          else if opr aconvc mul_tm
   9.438 -           then poly_mul (poly_of_term l) (poly_of_term r)
   9.439 -          else if opr aconvc div_tm
   9.440 -           then let
   9.441 -                  val p = poly_of_term l
   9.442 -                  val q = poly_of_term r
   9.443 -                in if poly_isconst q then poly_cmul (Rat.inv (eval FuncUtil.Ctermfunc.empty q)) p
   9.444 -                   else error "poly_of_term: division by non-constant polynomial"
   9.445 -                end
   9.446 -          else poly_var tm
   9.447 -
   9.448 -         end
   9.449 -         handle CTERM ("dest_comb",_) => poly_var tm)
   9.450 -   end
   9.451 -   handle CTERM ("dest_comb",_) => poly_var tm)
   9.452 -in
   9.453 -val poly_of_term = fn tm =>
   9.454 - if type_of (term_of tm) = @{typ real} then poly_of_term tm
   9.455 - else error "poly_of_term: term does not have real type"
   9.456 -end;
   9.457 -
   9.458 -(* String of vector (just a list of space-separated numbers).                *)
   9.459 -
   9.460 -fun sdpa_of_vector (v:vector) =
   9.461 - let
   9.462 -  val n = dim v
   9.463 -  val strs = map (decimalize 20 o (fn i => FuncUtil.Intfunc.tryapplyd (snd v) i rat_0)) (1 upto n)
   9.464 - in space_implode " " strs ^ "\n"
   9.465 - end;
   9.466 -
   9.467 -fun triple_int_ord ((a,b,c),(a',b',c')) =
   9.468 - prod_ord int_ord (prod_ord int_ord int_ord)
   9.469 -    ((a,(b,c)),(a',(b',c')));
   9.470 -structure Inttriplefunc = FuncFun(type key = int*int*int val ord = triple_int_ord);
   9.471 -
   9.472 -(* String for block diagonal matrix numbered k.                              *)
   9.473 -
   9.474 -fun sdpa_of_blockdiagonal k m =
   9.475 - let
   9.476 -  val pfx = string_of_int k ^" "
   9.477 -  val ents =
   9.478 -    Inttriplefunc.fold (fn ((b,i,j), c) => fn a => if i > j then a else ((b,i,j),c)::a) m []
   9.479 -  val entss = sort (triple_int_ord o pairself fst) ents
   9.480 - in  fold_rev (fn ((b,i,j),c) => fn a =>
   9.481 -     pfx ^ string_of_int b ^ " " ^ string_of_int i ^ " " ^ string_of_int j ^
   9.482 -     " " ^ decimalize 20 c ^ "\n" ^ a) entss ""
   9.483 - end;
   9.484 -
   9.485 -(* String for a matrix numbered k, in SDPA sparse format.                    *)
   9.486 -
   9.487 -fun sdpa_of_matrix k (m:matrix) =
   9.488 - let
   9.489 -  val pfx = string_of_int k ^ " 1 "
   9.490 -  val ms = FuncUtil.Intpairfunc.fold (fn ((i,j), c) => fn  a => if i > j then a else ((i,j),c)::a) (snd m) []
   9.491 -  val mss = sort ((prod_ord int_ord int_ord) o pairself fst) ms
   9.492 - in fold_rev (fn ((i,j),c) => fn a =>
   9.493 -     pfx ^ string_of_int i ^ " " ^ string_of_int j ^
   9.494 -     " " ^ decimalize 20 c ^ "\n" ^ a) mss ""
   9.495 - end;;
   9.496 -
   9.497 -(* ------------------------------------------------------------------------- *)
   9.498 -(* String in SDPA sparse format for standard SDP problem:                    *)
   9.499 -(*                                                                           *)
   9.500 -(*    X = v_1 * [M_1] + ... + v_m * [M_m] - [M_0] must be PSD                *)
   9.501 -(*    Minimize obj_1 * v_1 + ... obj_m * v_m                                 *)
   9.502 -(* ------------------------------------------------------------------------- *)
   9.503 -
   9.504 -fun sdpa_of_problem obj mats =
   9.505 - let
   9.506 -  val m = length mats - 1
   9.507 -  val (n,_) = dimensions (hd mats)
   9.508 - in
   9.509 -  string_of_int m ^ "\n" ^
   9.510 -  "1\n" ^
   9.511 -  string_of_int n ^ "\n" ^
   9.512 -  sdpa_of_vector obj ^
   9.513 -  fold_rev2 (fn k => fn m => fn a => sdpa_of_matrix (k - 1) m ^ a) (1 upto length mats) mats ""
   9.514 - end;
   9.515 -
   9.516 -fun index_char str chr pos =
   9.517 -  if pos >= String.size str then ~1
   9.518 -  else if String.sub(str,pos) = chr then pos
   9.519 -  else index_char str chr (pos + 1);
   9.520 -fun rat_of_quotient (a,b) = if b = 0 then rat_0 else Rat.rat_of_quotient (a,b);
   9.521 -fun rat_of_string s =
   9.522 - let val n = index_char s #"/" 0 in
   9.523 -  if n = ~1 then s |> Int.fromString |> the |> Rat.rat_of_int
   9.524 -  else
   9.525 -   let val SOME numer = Int.fromString(String.substring(s,0,n))
   9.526 -       val SOME den = Int.fromString (String.substring(s,n+1,String.size s - n - 1))
   9.527 -   in rat_of_quotient(numer, den)
   9.528 -   end
   9.529 - end;
   9.530 -
   9.531 -fun isspace x = (x = " ");
   9.532 -fun isnum x = member (op =) ["0","1","2","3","4","5","6","7","8","9"] x;
   9.533 -
   9.534 -(* More parser basics.                                                       *)
   9.535 -
   9.536 - val word = Scan.this_string
   9.537 - fun token s =
   9.538 -  Scan.repeat ($$ " ") |-- word s --| Scan.repeat ($$ " ")
   9.539 - val numeral = Scan.one isnum
   9.540 - val decimalint = Scan.repeat1 numeral >> (rat_of_string o implode)
   9.541 - val decimalfrac = Scan.repeat1 numeral
   9.542 -    >> (fn s => rat_of_string(implode s) // pow10 (length s))
   9.543 - val decimalsig =
   9.544 -    decimalint -- Scan.option (Scan.$$ "." |-- decimalfrac)
   9.545 -    >> (fn (h,NONE) => h | (h,SOME x) => h +/ x)
   9.546 - fun signed prs =
   9.547 -       $$ "-" |-- prs >> Rat.neg
   9.548 -    || $$ "+" |-- prs
   9.549 -    || prs;
   9.550 -
   9.551 -fun emptyin def xs = if null xs then (def,xs) else Scan.fail xs
   9.552 -
   9.553 - val exponent = ($$ "e" || $$ "E") |-- signed decimalint;
   9.554 -
   9.555 - val decimal = signed decimalsig -- (emptyin rat_0|| exponent)
   9.556 -    >> (fn (h, x) => h */ pow10 (int_of_rat x));
   9.557 -
   9.558 - fun mkparser p s =
   9.559 -  let val (x,rst) = p (raw_explode s)
   9.560 -  in if null rst then x
   9.561 -     else error "mkparser: unparsed input"
   9.562 -  end;;
   9.563 -
   9.564 -(* Parse back csdp output.                                                      *)
   9.565 -
   9.566 - fun ignore inp = ((),[])
   9.567 - fun csdpoutput inp =
   9.568 -   ((decimal -- Scan.repeat (Scan.$$ " " |-- Scan.option decimal) >>
   9.569 -    (fn (h,to) => map_filter I ((SOME h)::to))) --| ignore >> vector_of_list) inp
   9.570 - val parse_csdpoutput = mkparser csdpoutput
   9.571 -
   9.572 -(* Run prover on a problem in linear form.                       *)
   9.573 -
   9.574 -fun run_problem prover obj mats =
   9.575 -  parse_csdpoutput (prover (sdpa_of_problem obj mats))
   9.576 -
   9.577 -(* Try some apparently sensible scaling first. Note that this is purely to   *)
   9.578 -(* get a cleaner translation to floating-point, and doesn't affect any of    *)
   9.579 -(* the results, in principle. In practice it seems a lot better when there   *)
   9.580 -(* are extreme numbers in the original problem.                              *)
   9.581 -
   9.582 -  (* Version for (int*int) keys *)
   9.583 -local
   9.584 -  fun max_rat x y = if x </ y then y else x
   9.585 -  fun common_denominator fld amat acc =
   9.586 -      fld (fn (m,c) => fn a => lcm_rat (denominator_rat c) a) amat acc
   9.587 -  fun maximal_element fld amat acc =
   9.588 -    fld (fn (m,c) => fn maxa => max_rat maxa (abs_rat c)) amat acc
   9.589 -fun float_of_rat x = let val (a,b) = Rat.quotient_of_rat x
   9.590 -                     in Real.fromLargeInt a / Real.fromLargeInt b end;
   9.591 -in
   9.592 -
   9.593 -fun pi_scale_then solver (obj:vector)  mats =
   9.594 - let
   9.595 -  val cd1 = fold_rev (common_denominator FuncUtil.Intpairfunc.fold) mats (rat_1)
   9.596 -  val cd2 = common_denominator FuncUtil.Intfunc.fold (snd obj)  (rat_1)
   9.597 -  val mats' = map (FuncUtil.Intpairfunc.map (fn _ => fn x => cd1 */ x)) mats
   9.598 -  val obj' = vector_cmul cd2 obj
   9.599 -  val max1 = fold_rev (maximal_element FuncUtil.Intpairfunc.fold) mats' (rat_0)
   9.600 -  val max2 = maximal_element FuncUtil.Intfunc.fold (snd obj') (rat_0)
   9.601 -  val scal1 = pow2 (20 - trunc(Math.ln (float_of_rat max1) / Math.ln 2.0))
   9.602 -  val scal2 = pow2 (20 - trunc(Math.ln (float_of_rat max2) / Math.ln 2.0))
   9.603 -  val mats'' = map (FuncUtil.Intpairfunc.map (fn _ => fn x => x */ scal1)) mats'
   9.604 -  val obj'' = vector_cmul scal2 obj'
   9.605 - in solver obj'' mats''
   9.606 -  end
   9.607 -end;
   9.608 -
   9.609 -(* Try some apparently sensible scaling first. Note that this is purely to   *)
   9.610 -(* get a cleaner translation to floating-point, and doesn't affect any of    *)
   9.611 -(* the results, in principle. In practice it seems a lot better when there   *)
   9.612 -(* are extreme numbers in the original problem.                              *)
   9.613 -
   9.614 -  (* Version for (int*int*int) keys *)
   9.615 -local
   9.616 -  fun max_rat x y = if x </ y then y else x
   9.617 -  fun common_denominator fld amat acc =
   9.618 -      fld (fn (m,c) => fn a => lcm_rat (denominator_rat c) a) amat acc
   9.619 -  fun maximal_element fld amat acc =
   9.620 -    fld (fn (m,c) => fn maxa => max_rat maxa (abs_rat c)) amat acc
   9.621 -fun float_of_rat x = let val (a,b) = Rat.quotient_of_rat x
   9.622 -                     in Real.fromLargeInt a / Real.fromLargeInt b end;
   9.623 -fun int_of_float x = (trunc x handle Overflow => 0 | Domain => 0)
   9.624 -in
   9.625 -
   9.626 -fun tri_scale_then solver (obj:vector)  mats =
   9.627 - let
   9.628 -  val cd1 = fold_rev (common_denominator Inttriplefunc.fold) mats (rat_1)
   9.629 -  val cd2 = common_denominator FuncUtil.Intfunc.fold (snd obj)  (rat_1)
   9.630 -  val mats' = map (Inttriplefunc.map (fn _ => fn x => cd1 */ x)) mats
   9.631 -  val obj' = vector_cmul cd2 obj
   9.632 -  val max1 = fold_rev (maximal_element Inttriplefunc.fold) mats' (rat_0)
   9.633 -  val max2 = maximal_element FuncUtil.Intfunc.fold (snd obj') (rat_0)
   9.634 -  val scal1 = pow2 (20 - int_of_float(Math.ln (float_of_rat max1) / Math.ln 2.0))
   9.635 -  val scal2 = pow2 (20 - int_of_float(Math.ln (float_of_rat max2) / Math.ln 2.0))
   9.636 -  val mats'' = map (Inttriplefunc.map (fn _ => fn x => x */ scal1)) mats'
   9.637 -  val obj'' = vector_cmul scal2 obj'
   9.638 - in solver obj'' mats''
   9.639 -  end
   9.640 -end;
   9.641 -
   9.642 -(* Round a vector to "nice" rationals.                                       *)
   9.643 -
   9.644 -fun nice_rational n x = round_rat (n */ x) // n;;
   9.645 -fun nice_vector n ((d,v) : vector) =
   9.646 - (d, FuncUtil.Intfunc.fold (fn (i,c) => fn a =>
   9.647 -   let val y = nice_rational n c
   9.648 -   in if c =/ rat_0 then a
   9.649 -      else FuncUtil.Intfunc.update (i,y) a end) v FuncUtil.Intfunc.empty):vector
   9.650 -
   9.651 -fun dest_ord f x = is_equal (f x);
   9.652 -
   9.653 -(* Stuff for "equations" ((int*int*int)->num functions).                         *)
   9.654 -
   9.655 -fun tri_equation_cmul c eq =
   9.656 -  if c =/ rat_0 then Inttriplefunc.empty else Inttriplefunc.map (fn _ => fn d => c */ d) eq;
   9.657 -
   9.658 -fun tri_equation_add eq1 eq2 = Inttriplefunc.combine (curry op +/) (fn x => x =/ rat_0) eq1 eq2;
   9.659 -
   9.660 -fun tri_equation_eval assig eq =
   9.661 - let fun value v = Inttriplefunc.apply assig v
   9.662 - in Inttriplefunc.fold (fn (v, c) => fn a => a +/ value v */ c) eq rat_0
   9.663 - end;
   9.664 -
   9.665 -(* Eliminate among linear equations: return unconstrained variables and      *)
   9.666 -(* assignments for the others in terms of them. We give one pseudo-variable  *)
   9.667 -(* "one" that's used for a constant term.                                    *)
   9.668 -
   9.669 -local
   9.670 -  fun extract_first p l = case l of  (* FIXME : use find_first instead *)
   9.671 -   [] => error "extract_first"
   9.672 - | h::t => if p h then (h,t) else
   9.673 -          let val (k,s) = extract_first p t in (k,h::s) end
   9.674 -fun eliminate vars dun eqs = case vars of
   9.675 -  [] => if forall Inttriplefunc.is_empty eqs then dun
   9.676 -        else raise Unsolvable
   9.677 - | v::vs =>
   9.678 -  ((let
   9.679 -    val (eq,oeqs) = extract_first (fn e => Inttriplefunc.defined e v) eqs
   9.680 -    val a = Inttriplefunc.apply eq v
   9.681 -    val eq' = tri_equation_cmul ((Rat.neg rat_1) // a) (Inttriplefunc.delete_safe v eq)
   9.682 -    fun elim e =
   9.683 -     let val b = Inttriplefunc.tryapplyd e v rat_0
   9.684 -     in if b =/ rat_0 then e else
   9.685 -        tri_equation_add e (tri_equation_cmul (Rat.neg b // a) eq)
   9.686 -     end
   9.687 -   in eliminate vs (Inttriplefunc.update (v,eq') (Inttriplefunc.map (K elim) dun)) (map elim oeqs)
   9.688 -   end)
   9.689 -  handle Failure _ => eliminate vs dun eqs)
   9.690 -in
   9.691 -fun tri_eliminate_equations one vars eqs =
   9.692 - let
   9.693 -  val assig = eliminate vars Inttriplefunc.empty eqs
   9.694 -  val vs = Inttriplefunc.fold (fn (x, f) => fn a => remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig []
   9.695 -  in (distinct (dest_ord triple_int_ord) vs, assig)
   9.696 -  end
   9.697 -end;
   9.698 -
   9.699 -(* Eliminate all variables, in an essentially arbitrary order.               *)
   9.700 -
   9.701 -fun tri_eliminate_all_equations one =
   9.702 - let
   9.703 -  fun choose_variable eq =
   9.704 -   let val (v,_) = Inttriplefunc.choose eq
   9.705 -   in if is_equal (triple_int_ord(v,one)) then
   9.706 -      let val eq' = Inttriplefunc.delete_safe v eq
   9.707 -      in if Inttriplefunc.is_empty eq' then error "choose_variable"
   9.708 -         else fst (Inttriplefunc.choose eq')
   9.709 -      end
   9.710 -    else v
   9.711 -   end
   9.712 -  fun eliminate dun eqs = case eqs of
   9.713 -    [] => dun
   9.714 -  | eq::oeqs =>
   9.715 -    if Inttriplefunc.is_empty eq then eliminate dun oeqs else
   9.716 -    let val v = choose_variable eq
   9.717 -        val a = Inttriplefunc.apply eq v
   9.718 -        val eq' = tri_equation_cmul ((Rat.rat_of_int ~1) // a)
   9.719 -                   (Inttriplefunc.delete_safe v eq)
   9.720 -        fun elim e =
   9.721 -         let val b = Inttriplefunc.tryapplyd e v rat_0
   9.722 -         in if b =/ rat_0 then e
   9.723 -            else tri_equation_add e (tri_equation_cmul (Rat.neg b // a) eq)
   9.724 -         end
   9.725 -    in eliminate (Inttriplefunc.update(v, eq') (Inttriplefunc.map (K elim) dun))
   9.726 -                 (map elim oeqs)
   9.727 -    end
   9.728 -in fn eqs =>
   9.729 - let
   9.730 -  val assig = eliminate Inttriplefunc.empty eqs
   9.731 -  val vs = Inttriplefunc.fold (fn (x, f) => fn a => remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig []
   9.732 - in (distinct (dest_ord triple_int_ord) vs,assig)
   9.733 - end
   9.734 -end;
   9.735 -
   9.736 -(* Solve equations by assigning arbitrary numbers.                           *)
   9.737 -
   9.738 -fun tri_solve_equations one eqs =
   9.739 - let
   9.740 -  val (vars,assigs) = tri_eliminate_all_equations one eqs
   9.741 -  val vfn = fold_rev (fn v => Inttriplefunc.update(v,rat_0)) vars
   9.742 -            (Inttriplefunc.onefunc(one, Rat.rat_of_int ~1))
   9.743 -  val ass =
   9.744 -    Inttriplefunc.combine (curry op +/) (K false)
   9.745 -    (Inttriplefunc.map (K (tri_equation_eval vfn)) assigs) vfn
   9.746 - in if forall (fn e => tri_equation_eval ass e =/ rat_0) eqs
   9.747 -    then Inttriplefunc.delete_safe one ass else raise Sanity
   9.748 - end;
   9.749 -
   9.750 -(* Multiply equation-parametrized poly by regular poly and add accumulator.  *)
   9.751 -
   9.752 -fun tri_epoly_pmul p q acc =
   9.753 - FuncUtil.Monomialfunc.fold (fn (m1, c) => fn a =>
   9.754 -  FuncUtil.Monomialfunc.fold (fn (m2,e) => fn b =>
   9.755 -   let val m =  monomial_mul m1 m2
   9.756 -       val es = FuncUtil.Monomialfunc.tryapplyd b m Inttriplefunc.empty
   9.757 -   in FuncUtil.Monomialfunc.update (m,tri_equation_add (tri_equation_cmul c e) es) b
   9.758 -   end) q a) p acc ;
   9.759 -
   9.760 -(* Usual operations on equation-parametrized poly.                           *)
   9.761 -
   9.762 -fun tri_epoly_cmul c l =
   9.763 -  if c =/ rat_0 then Inttriplefunc.empty else Inttriplefunc.map (K (tri_equation_cmul c)) l;;
   9.764 -
   9.765 -val tri_epoly_neg = tri_epoly_cmul (Rat.rat_of_int ~1);
   9.766 -
   9.767 -val tri_epoly_add = Inttriplefunc.combine tri_equation_add Inttriplefunc.is_empty;
   9.768 -
   9.769 -fun tri_epoly_sub p q = tri_epoly_add p (tri_epoly_neg q);;
   9.770 -
   9.771 -(* Stuff for "equations" ((int*int)->num functions).                         *)
   9.772 -
   9.773 -fun pi_equation_cmul c eq =
   9.774 -  if c =/ rat_0 then Inttriplefunc.empty else Inttriplefunc.map (fn _ => fn d => c */ d) eq;
   9.775 -
   9.776 -fun pi_equation_add eq1 eq2 = Inttriplefunc.combine (curry op +/) (fn x => x =/ rat_0) eq1 eq2;
   9.777 -
   9.778 -fun pi_equation_eval assig eq =
   9.779 - let fun value v = Inttriplefunc.apply assig v
   9.780 - in Inttriplefunc.fold (fn (v, c) => fn a => a +/ value v */ c) eq rat_0
   9.781 - end;
   9.782 -
   9.783 -(* Eliminate among linear equations: return unconstrained variables and      *)
   9.784 -(* assignments for the others in terms of them. We give one pseudo-variable  *)
   9.785 -(* "one" that's used for a constant term.                                    *)
   9.786 -
   9.787 -local
   9.788 -fun extract_first p l = case l of
   9.789 -   [] => error "extract_first"
   9.790 - | h::t => if p h then (h,t) else
   9.791 -          let val (k,s) = extract_first p t in (k,h::s) end
   9.792 -fun eliminate vars dun eqs = case vars of
   9.793 -  [] => if forall Inttriplefunc.is_empty eqs then dun
   9.794 -        else raise Unsolvable
   9.795 - | v::vs =>
   9.796 -   let
   9.797 -    val (eq,oeqs) = extract_first (fn e => Inttriplefunc.defined e v) eqs
   9.798 -    val a = Inttriplefunc.apply eq v
   9.799 -    val eq' = pi_equation_cmul ((Rat.neg rat_1) // a) (Inttriplefunc.delete_safe v eq)
   9.800 -    fun elim e =
   9.801 -     let val b = Inttriplefunc.tryapplyd e v rat_0
   9.802 -     in if b =/ rat_0 then e else
   9.803 -        pi_equation_add e (pi_equation_cmul (Rat.neg b // a) eq)
   9.804 -     end
   9.805 -   in eliminate vs (Inttriplefunc.update (v,eq') (Inttriplefunc.map (K elim) dun)) (map elim oeqs)
   9.806 -   end
   9.807 -  handle Failure _ => eliminate vs dun eqs
   9.808 -in
   9.809 -fun pi_eliminate_equations one vars eqs =
   9.810 - let
   9.811 -  val assig = eliminate vars Inttriplefunc.empty eqs
   9.812 -  val vs = Inttriplefunc.fold (fn (x, f) => fn a => remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig []
   9.813 -  in (distinct (dest_ord triple_int_ord) vs, assig)
   9.814 -  end
   9.815 -end;
   9.816 -
   9.817 -(* Eliminate all variables, in an essentially arbitrary order.               *)
   9.818 -
   9.819 -fun pi_eliminate_all_equations one =
   9.820 - let
   9.821 -  fun choose_variable eq =
   9.822 -   let val (v,_) = Inttriplefunc.choose eq
   9.823 -   in if is_equal (triple_int_ord(v,one)) then
   9.824 -      let val eq' = Inttriplefunc.delete_safe v eq
   9.825 -      in if Inttriplefunc.is_empty eq' then error "choose_variable"
   9.826 -         else fst (Inttriplefunc.choose eq')
   9.827 -      end
   9.828 -    else v
   9.829 -   end
   9.830 -  fun eliminate dun eqs = case eqs of
   9.831 -    [] => dun
   9.832 -  | eq::oeqs =>
   9.833 -    if Inttriplefunc.is_empty eq then eliminate dun oeqs else
   9.834 -    let val v = choose_variable eq
   9.835 -        val a = Inttriplefunc.apply eq v
   9.836 -        val eq' = pi_equation_cmul ((Rat.rat_of_int ~1) // a)
   9.837 -                   (Inttriplefunc.delete_safe v eq)
   9.838 -        fun elim e =
   9.839 -         let val b = Inttriplefunc.tryapplyd e v rat_0
   9.840 -         in if b =/ rat_0 then e
   9.841 -            else pi_equation_add e (pi_equation_cmul (Rat.neg b // a) eq)
   9.842 -         end
   9.843 -    in eliminate (Inttriplefunc.update(v, eq') (Inttriplefunc.map (K elim) dun))
   9.844 -                 (map elim oeqs)
   9.845 -    end
   9.846 -in fn eqs =>
   9.847 - let
   9.848 -  val assig = eliminate Inttriplefunc.empty eqs
   9.849 -  val vs = Inttriplefunc.fold (fn (x, f) => fn a => remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig []
   9.850 - in (distinct (dest_ord triple_int_ord) vs,assig)
   9.851 - end
   9.852 -end;
   9.853 -
   9.854 -(* Solve equations by assigning arbitrary numbers.                           *)
   9.855 -
   9.856 -fun pi_solve_equations one eqs =
   9.857 - let
   9.858 -  val (vars,assigs) = pi_eliminate_all_equations one eqs
   9.859 -  val vfn = fold_rev (fn v => Inttriplefunc.update(v,rat_0)) vars
   9.860 -            (Inttriplefunc.onefunc(one, Rat.rat_of_int ~1))
   9.861 -  val ass =
   9.862 -    Inttriplefunc.combine (curry op +/) (K false)
   9.863 -    (Inttriplefunc.map (K (pi_equation_eval vfn)) assigs) vfn
   9.864 - in if forall (fn e => pi_equation_eval ass e =/ rat_0) eqs
   9.865 -    then Inttriplefunc.delete_safe one ass else raise Sanity
   9.866 - end;
   9.867 -
   9.868 -(* Multiply equation-parametrized poly by regular poly and add accumulator.  *)
   9.869 -
   9.870 -fun pi_epoly_pmul p q acc =
   9.871 - FuncUtil.Monomialfunc.fold (fn (m1, c) => fn a =>
   9.872 -  FuncUtil.Monomialfunc.fold (fn (m2,e) => fn b =>
   9.873 -   let val m =  monomial_mul m1 m2
   9.874 -       val es = FuncUtil.Monomialfunc.tryapplyd b m Inttriplefunc.empty
   9.875 -   in FuncUtil.Monomialfunc.update (m,pi_equation_add (pi_equation_cmul c e) es) b
   9.876 -   end) q a) p acc ;
   9.877 -
   9.878 -(* Usual operations on equation-parametrized poly.                           *)
   9.879 -
   9.880 -fun pi_epoly_cmul c l =
   9.881 -  if c =/ rat_0 then Inttriplefunc.empty else Inttriplefunc.map (K (pi_equation_cmul c)) l;;
   9.882 -
   9.883 -val pi_epoly_neg = pi_epoly_cmul (Rat.rat_of_int ~1);
   9.884 -
   9.885 -val pi_epoly_add = Inttriplefunc.combine pi_equation_add Inttriplefunc.is_empty;
   9.886 -
   9.887 -fun pi_epoly_sub p q = pi_epoly_add p (pi_epoly_neg q);;
   9.888 -
   9.889 -fun allpairs f l1 l2 =  fold_rev (fn x => (curry (op @)) (map (f x) l2)) l1 [];
   9.890 -
   9.891 -(* Hence produce the "relevant" monomials: those whose squares lie in the    *)
   9.892 -(* Newton polytope of the monomials in the input. (This is enough according  *)
   9.893 -(* to Reznik: "Extremal PSD forms with few terms", Duke Math. Journal,       *)
   9.894 -(* vol 45, pp. 363--374, 1978.                                               *)
   9.895 -(*                                                                           *)
   9.896 -(* These are ordered in sort of decreasing degree. In particular the         *)
   9.897 -(* constant monomial is last; this gives an order in diagonalization of the  *)
   9.898 -(* quadratic form that will tend to display constants.                       *)
   9.899 -
   9.900 -(* Diagonalize (Cholesky/LDU) the matrix corresponding to a quadratic form.  *)
   9.901 -
   9.902 -local
   9.903 -fun diagonalize n i m =
   9.904 - if FuncUtil.Intpairfunc.is_empty (snd m) then []
   9.905 - else
   9.906 -  let val a11 = FuncUtil.Intpairfunc.tryapplyd (snd m) (i,i) rat_0
   9.907 -  in if a11 </ rat_0 then raise Failure "diagonalize: not PSD"
   9.908 -    else if a11 =/ rat_0 then
   9.909 -          if FuncUtil.Intfunc.is_empty (snd (row i m)) then diagonalize n (i + 1) m
   9.910 -          else raise Failure "diagonalize: not PSD ___ "
   9.911 -    else
   9.912 -     let
   9.913 -      val v = row i m
   9.914 -      val v' = (fst v, FuncUtil.Intfunc.fold (fn (i, c) => fn a =>
   9.915 -       let val y = c // a11
   9.916 -       in if y = rat_0 then a else FuncUtil.Intfunc.update (i,y) a
   9.917 -       end)  (snd v) FuncUtil.Intfunc.empty)
   9.918 -      fun upt0 x y a = if y = rat_0 then a else FuncUtil.Intpairfunc.update (x,y) a
   9.919 -      val m' =
   9.920 -      ((n,n),
   9.921 -      iter (i+1,n) (fn j =>
   9.922 -          iter (i+1,n) (fn k =>
   9.923 -              (upt0 (j,k) (FuncUtil.Intpairfunc.tryapplyd (snd m) (j,k) rat_0 -/ FuncUtil.Intfunc.tryapplyd (snd v) j rat_0 */ FuncUtil.Intfunc.tryapplyd (snd v') k rat_0))))
   9.924 -          FuncUtil.Intpairfunc.empty)
   9.925 -     in (a11,v')::diagonalize n (i + 1) m'
   9.926 -     end
   9.927 -  end
   9.928 -in
   9.929 -fun diag m =
   9.930 - let
   9.931 -   val nn = dimensions m
   9.932 -   val n = fst nn
   9.933 - in if snd nn <> n then error "diagonalize: non-square matrix"
   9.934 -    else diagonalize n 1 m
   9.935 - end
   9.936 -end;
   9.937 -
   9.938 -fun gcd_rat a b = Rat.rat_of_int (Integer.gcd (int_of_rat a) (int_of_rat b));
   9.939 -
   9.940 -(* Adjust a diagonalization to collect rationals at the start.               *)
   9.941 -  (* FIXME : Potentially polymorphic keys, but here only: integers!! *)
   9.942 -local
   9.943 - fun upd0 x y a = if y =/ rat_0 then a else FuncUtil.Intfunc.update(x,y) a;
   9.944 - fun mapa f (d,v) =
   9.945 -  (d, FuncUtil.Intfunc.fold (fn (i,c) => fn a => upd0 i (f c) a) v FuncUtil.Intfunc.empty)
   9.946 - fun adj (c,l) =
   9.947 - let val a =
   9.948 -  FuncUtil.Intfunc.fold (fn (i,c) => fn a => lcm_rat a (denominator_rat c))
   9.949 -    (snd l) rat_1 //
   9.950 -  FuncUtil.Intfunc.fold (fn (i,c) => fn a => gcd_rat a (numerator_rat c))
   9.951 -    (snd l) rat_0
   9.952 -  in ((c // (a */ a)),mapa (fn x => a */ x) l)
   9.953 -  end
   9.954 -in
   9.955 -fun deration d = if null d then (rat_0,d) else
   9.956 - let val d' = map adj d
   9.957 -     val a = fold (lcm_rat o denominator_rat o fst) d' rat_1 //
   9.958 -          fold (gcd_rat o numerator_rat o fst) d' rat_0
   9.959 - in ((rat_1 // a),map (fn (c,l) => (a */ c,l)) d')
   9.960 - end
   9.961 -end;
   9.962 -
   9.963 -(* Enumeration of monomials with given multidegree bound.                    *)
   9.964 -
   9.965 -fun enumerate_monomials d vars =
   9.966 - if d < 0 then []
   9.967 - else if d = 0 then [FuncUtil.Ctermfunc.empty]
   9.968 - else if null vars then [monomial_1] else
   9.969 - let val alts =
   9.970 -  map_range (fn k => let val oths = enumerate_monomials (d - k) (tl vars)
   9.971 -               in map (fn ks => if k = 0 then ks else FuncUtil.Ctermfunc.update (hd vars, k) ks) oths end) (d + 1)
   9.972 - in flat alts
   9.973 - end;
   9.974 -
   9.975 -(* Enumerate products of distinct input polys with degree <= d.              *)
   9.976 -(* We ignore any constant input polynomials.                                 *)
   9.977 -(* Give the output polynomial and a record of how it was derived.            *)
   9.978 -
   9.979 -fun enumerate_products d pols =
   9.980 -if d = 0 then [(poly_const rat_1,RealArith.Rational_lt rat_1)]
   9.981 -else if d < 0 then [] else
   9.982 -case pols of
   9.983 -   [] => [(poly_const rat_1,RealArith.Rational_lt rat_1)]
   9.984 - | (p,b)::ps =>
   9.985 -    let val e = multidegree p
   9.986 -    in if e = 0 then enumerate_products d ps else
   9.987 -       enumerate_products d ps @
   9.988 -       map (fn (q,c) => (poly_mul p q,RealArith.Product(b,c)))
   9.989 -         (enumerate_products (d - e) ps)
   9.990 -    end
   9.991 -
   9.992 -(* Convert regular polynomial. Note that we treat (0,0,0) as -1.             *)
   9.993 -
   9.994 -fun epoly_of_poly p =
   9.995 -  FuncUtil.Monomialfunc.fold (fn (m,c) => fn a => FuncUtil.Monomialfunc.update (m, Inttriplefunc.onefunc ((0,0,0), Rat.neg c)) a) p FuncUtil.Monomialfunc.empty;
   9.996 -
   9.997 -(* String for block diagonal matrix numbered k.                              *)
   9.998 -
   9.999 -fun sdpa_of_blockdiagonal k m =
  9.1000 - let
  9.1001 -  val pfx = string_of_int k ^" "
  9.1002 -  val ents =
  9.1003 -    Inttriplefunc.fold
  9.1004 -      (fn ((b,i,j),c) => fn a => if i > j then a else ((b,i,j),c)::a)
  9.1005 -      m []
  9.1006 -  val entss = sort (triple_int_ord o pairself fst) ents
  9.1007 - in fold_rev (fn ((b,i,j),c) => fn a =>
  9.1008 -     pfx ^ string_of_int b ^ " " ^ string_of_int i ^ " " ^ string_of_int j ^
  9.1009 -     " " ^ decimalize 20 c ^ "\n" ^ a) entss ""
  9.1010 - end;
  9.1011 -
  9.1012 -(* SDPA for problem using block diagonal (i.e. multiple SDPs)                *)
  9.1013 -
  9.1014 -fun sdpa_of_blockproblem nblocks blocksizes obj mats =
  9.1015 - let val m = length mats - 1
  9.1016 - in
  9.1017 -  string_of_int m ^ "\n" ^
  9.1018 -  string_of_int nblocks ^ "\n" ^
  9.1019 -  (space_implode " " (map string_of_int blocksizes)) ^
  9.1020 -  "\n" ^
  9.1021 -  sdpa_of_vector obj ^
  9.1022 -  fold_rev2 (fn k => fn m => fn a => sdpa_of_blockdiagonal (k - 1) m ^ a)
  9.1023 -    (1 upto length mats) mats ""
  9.1024 - end;
  9.1025 -
  9.1026 -(* Run prover on a problem in block diagonal form.                       *)
  9.1027 -
  9.1028 -fun run_blockproblem prover nblocks blocksizes obj mats=
  9.1029 -  parse_csdpoutput (prover (sdpa_of_blockproblem nblocks blocksizes obj mats))
  9.1030 -
  9.1031 -(* 3D versions of matrix operations to consider blocks separately.           *)
  9.1032 -
  9.1033 -val bmatrix_add = Inttriplefunc.combine (curry op +/) (fn x => x =/ rat_0);
  9.1034 -fun bmatrix_cmul c bm =
  9.1035 -  if c =/ rat_0 then Inttriplefunc.empty
  9.1036 -  else Inttriplefunc.map (fn _ => fn x => c */ x) bm;
  9.1037 -
  9.1038 -val bmatrix_neg = bmatrix_cmul (Rat.rat_of_int ~1);
  9.1039 -fun bmatrix_sub m1 m2 = bmatrix_add m1 (bmatrix_neg m2);;
  9.1040 -
  9.1041 -(* Smash a block matrix into components.                                     *)
  9.1042 -
  9.1043 -fun blocks blocksizes bm =
  9.1044 - map (fn (bs,b0) =>
  9.1045 -      let val m = Inttriplefunc.fold
  9.1046 -          (fn ((b,i,j),c) => fn a => if b = b0 then FuncUtil.Intpairfunc.update ((i,j),c) a else a) bm FuncUtil.Intpairfunc.empty
  9.1047 -          val d = FuncUtil.Intpairfunc.fold (fn ((i,j),c) => fn a => max a (max i j)) m 0
  9.1048 -      in (((bs,bs),m):matrix) end)
  9.1049 - (blocksizes ~~ (1 upto length blocksizes));;
  9.1050 -
  9.1051 -(* FIXME : Get rid of this !!!*)
  9.1052 -local
  9.1053 -  fun tryfind_with msg f [] = raise Failure msg
  9.1054 -    | tryfind_with msg f (x::xs) = (f x handle Failure s => tryfind_with s f xs);
  9.1055 -in
  9.1056 -  fun tryfind f = tryfind_with "tryfind" f
  9.1057 -end
  9.1058 -
  9.1059 -(* Positiv- and Nullstellensatz. Flag "linf" forces a linear representation. *)
  9.1060 -
  9.1061 -
  9.1062 -fun real_positivnullstellensatz_general ctxt prover linf d eqs leqs pol =
  9.1063 -let
  9.1064 - val vars = fold_rev (union (op aconvc) o poly_variables)
  9.1065 -   (pol :: eqs @ map fst leqs) []
  9.1066 - val monoid = if linf then
  9.1067 -      (poly_const rat_1,RealArith.Rational_lt rat_1)::
  9.1068 -      (filter (fn (p,c) => multidegree p <= d) leqs)
  9.1069 -    else enumerate_products d leqs
  9.1070 - val nblocks = length monoid
  9.1071 - fun mk_idmultiplier k p =
  9.1072 -  let
  9.1073 -   val e = d - multidegree p
  9.1074 -   val mons = enumerate_monomials e vars
  9.1075 -   val nons = mons ~~ (1 upto length mons)
  9.1076 -  in (mons,
  9.1077 -      fold_rev (fn (m,n) => FuncUtil.Monomialfunc.update(m,Inttriplefunc.onefunc((~k,~n,n),rat_1))) nons FuncUtil.Monomialfunc.empty)
  9.1078 -  end
  9.1079 -
  9.1080 - fun mk_sqmultiplier k (p,c) =
  9.1081 -  let
  9.1082 -   val e = (d - multidegree p) div 2
  9.1083 -   val mons = enumerate_monomials e vars
  9.1084 -   val nons = mons ~~ (1 upto length mons)
  9.1085 -  in (mons,
  9.1086 -      fold_rev (fn (m1,n1) =>
  9.1087 -       fold_rev (fn (m2,n2) => fn  a =>
  9.1088 -        let val m = monomial_mul m1 m2
  9.1089 -        in if n1 > n2 then a else
  9.1090 -          let val c = if n1 = n2 then rat_1 else rat_2
  9.1091 -              val e = FuncUtil.Monomialfunc.tryapplyd a m Inttriplefunc.empty
  9.1092 -          in FuncUtil.Monomialfunc.update(m, tri_equation_add (Inttriplefunc.onefunc((k,n1,n2), c)) e) a
  9.1093 -          end
  9.1094 -        end)  nons)
  9.1095 -       nons FuncUtil.Monomialfunc.empty)
  9.1096 -  end
  9.1097 -
  9.1098 -  val (sqmonlist,sqs) = split_list (map2 mk_sqmultiplier (1 upto length monoid) monoid)
  9.1099 -  val (idmonlist,ids) =  split_list(map2 mk_idmultiplier (1 upto length eqs) eqs)
  9.1100 -  val blocksizes = map length sqmonlist
  9.1101 -  val bigsum =
  9.1102 -    fold_rev2 (fn p => fn q => fn a => tri_epoly_pmul p q a) eqs ids
  9.1103 -            (fold_rev2 (fn (p,c) => fn s => fn a => tri_epoly_pmul p s a) monoid sqs
  9.1104 -                     (epoly_of_poly(poly_neg pol)))
  9.1105 -  val eqns = FuncUtil.Monomialfunc.fold (fn (m,e) => fn a => e::a) bigsum []
  9.1106 -  val (pvs,assig) = tri_eliminate_all_equations (0,0,0) eqns
  9.1107 -  val qvars = (0,0,0)::pvs
  9.1108 -  val allassig = fold_rev (fn v => Inttriplefunc.update(v,(Inttriplefunc.onefunc(v,rat_1)))) pvs assig
  9.1109 -  fun mk_matrix v =
  9.1110 -    Inttriplefunc.fold (fn ((b,i,j), ass) => fn m =>
  9.1111 -        if b < 0 then m else
  9.1112 -         let val c = Inttriplefunc.tryapplyd ass v rat_0
  9.1113 -         in if c = rat_0 then m else
  9.1114 -            Inttriplefunc.update ((b,j,i), c) (Inttriplefunc.update ((b,i,j), c) m)
  9.1115 -         end)
  9.1116 -          allassig Inttriplefunc.empty
  9.1117 -  val diagents = Inttriplefunc.fold
  9.1118 -    (fn ((b,i,j), e) => fn a => if b > 0 andalso i = j then tri_equation_add e a else a)
  9.1119 -    allassig Inttriplefunc.empty
  9.1120 -
  9.1121 -  val mats = map mk_matrix qvars
  9.1122 -  val obj = (length pvs,
  9.1123 -            itern 1 pvs (fn v => fn i => FuncUtil.Intfunc.updatep iszero (i,Inttriplefunc.tryapplyd diagents v rat_0))
  9.1124 -                        FuncUtil.Intfunc.empty)
  9.1125 -  val raw_vec = if null pvs then vector_0 0
  9.1126 -                else tri_scale_then (run_blockproblem prover nblocks blocksizes) obj mats
  9.1127 -  fun int_element (d,v) i = FuncUtil.Intfunc.tryapplyd v i rat_0
  9.1128 -
  9.1129 -  fun find_rounding d =
  9.1130 -   let
  9.1131 -    val _ =
  9.1132 -      if Config.get ctxt trace
  9.1133 -      then writeln ("Trying rounding with limit "^Rat.string_of_rat d ^ "\n")
  9.1134 -      else ()
  9.1135 -    val vec = nice_vector d raw_vec
  9.1136 -    val blockmat = iter (1,dim vec)
  9.1137 -     (fn i => fn a => bmatrix_add (bmatrix_cmul (int_element vec i) (nth mats i)) a)
  9.1138 -     (bmatrix_neg (nth mats 0))
  9.1139 -    val allmats = blocks blocksizes blockmat
  9.1140 -   in (vec,map diag allmats)
  9.1141 -   end
  9.1142 -  val (vec,ratdias) =
  9.1143 -    if null pvs then find_rounding rat_1
  9.1144 -    else tryfind find_rounding (map Rat.rat_of_int (1 upto 31) @
  9.1145 -                                map pow2 (5 upto 66))
  9.1146 -  val newassigs =
  9.1147 -    fold_rev (fn k => Inttriplefunc.update (nth pvs (k - 1), int_element vec k))
  9.1148 -           (1 upto dim vec) (Inttriplefunc.onefunc ((0,0,0), Rat.rat_of_int ~1))
  9.1149 -  val finalassigs =
  9.1150 -    Inttriplefunc.fold (fn (v,e) => fn a => Inttriplefunc.update(v, tri_equation_eval newassigs e) a) allassig newassigs
  9.1151 -  fun poly_of_epoly p =
  9.1152 -    FuncUtil.Monomialfunc.fold (fn (v,e) => fn a => FuncUtil.Monomialfunc.updatep iszero (v,tri_equation_eval finalassigs e) a)
  9.1153 -          p FuncUtil.Monomialfunc.empty
  9.1154 -  fun  mk_sos mons =
  9.1155 -   let fun mk_sq (c,m) =
  9.1156 -    (c,fold_rev (fn k=> fn a => FuncUtil.Monomialfunc.updatep iszero (nth mons (k - 1), int_element m k) a)
  9.1157 -                 (1 upto length mons) FuncUtil.Monomialfunc.empty)
  9.1158 -   in map mk_sq
  9.1159 -   end
  9.1160 -  val sqs = map2 mk_sos sqmonlist ratdias
  9.1161 -  val cfs = map poly_of_epoly ids
  9.1162 -  val msq = filter (fn (a,b) => not (null b)) (map2 pair monoid sqs)
  9.1163 -  fun eval_sq sqs = fold_rev (fn (c,q) => poly_add (poly_cmul c (poly_mul q q))) sqs poly_0
  9.1164 -  val sanity =
  9.1165 -    fold_rev (fn ((p,c),s) => poly_add (poly_mul p (eval_sq s))) msq
  9.1166 -           (fold_rev2 (fn p => fn q => poly_add (poly_mul p q)) cfs eqs
  9.1167 -                    (poly_neg pol))
  9.1168 -
  9.1169 -in if not(FuncUtil.Monomialfunc.is_empty sanity) then raise Sanity else
  9.1170 -  (cfs,map (fn (a,b) => (snd a,b)) msq)
  9.1171 - end
  9.1172 -
  9.1173 -
  9.1174 -(* Iterative deepening.                                                      *)
  9.1175 -
  9.1176 -fun deepen f n =
  9.1177 -  (writeln ("Searching with depth limit " ^ string_of_int n);
  9.1178 -    (f n handle Failure s => (writeln ("failed with message: " ^ s); deepen f (n + 1))));
  9.1179 -
  9.1180 -
  9.1181 -(* Map back polynomials and their composites to a positivstellensatz.        *)
  9.1182 -
  9.1183 -fun cterm_of_sqterm (c,p) = RealArith.Product(RealArith.Rational_lt c,RealArith.Square p);
  9.1184 -
  9.1185 -fun cterm_of_sos (pr,sqs) = if null sqs then pr
  9.1186 -  else RealArith.Product(pr,foldr1 RealArith.Sum (map cterm_of_sqterm sqs));
  9.1187 -
  9.1188 -(* Interface to HOL.                                                         *)
  9.1189 -local
  9.1190 -  open Conv
  9.1191 -  val concl = Thm.dest_arg o cprop_of
  9.1192 -  fun simple_cterm_ord t u = Term_Ord.fast_term_ord (term_of t, term_of u) = LESS
  9.1193 -in
  9.1194 -  (* FIXME: Replace tryfind by get_first !! *)
  9.1195 -fun real_nonlinear_prover proof_method ctxt =
  9.1196 - let
  9.1197 -  val {add,mul,neg,pow,sub,main} =  Semiring_Normalizer.semiring_normalizers_ord_wrapper ctxt
  9.1198 -      (the (Semiring_Normalizer.match ctxt @{cterm "(0::real) + 1"}))
  9.1199 -     simple_cterm_ord
  9.1200 -  val (real_poly_add_conv,real_poly_mul_conv,real_poly_neg_conv,
  9.1201 -       real_poly_pow_conv,real_poly_sub_conv,real_poly_conv) = (add,mul,neg,pow,sub,main)
  9.1202 -  fun mainf cert_choice translator (eqs,les,lts) =
  9.1203 -  let
  9.1204 -   val eq0 = map (poly_of_term o Thm.dest_arg1 o concl) eqs
  9.1205 -   val le0 = map (poly_of_term o Thm.dest_arg o concl) les
  9.1206 -   val lt0 = map (poly_of_term o Thm.dest_arg o concl) lts
  9.1207 -   val eqp0 = map_index (fn (i, t) => (t,RealArith.Axiom_eq i)) eq0
  9.1208 -   val lep0 = map_index (fn (i, t) => (t,RealArith.Axiom_le i)) le0
  9.1209 -   val ltp0 = map_index (fn (i, t) => (t,RealArith.Axiom_lt i)) lt0
  9.1210 -   val (keq,eq) = List.partition (fn (p,_) => multidegree p = 0) eqp0
  9.1211 -   val (klep,lep) = List.partition (fn (p,_) => multidegree p = 0) lep0
  9.1212 -   val (kltp,ltp) = List.partition (fn (p,_) => multidegree p = 0) ltp0
  9.1213 -   fun trivial_axiom (p,ax) =
  9.1214 -    case ax of
  9.1215 -       RealArith.Axiom_eq n => if eval FuncUtil.Ctermfunc.empty p <>/ Rat.zero then nth eqs n
  9.1216 -                     else raise Failure "trivial_axiom: Not a trivial axiom"
  9.1217 -     | RealArith.Axiom_le n => if eval FuncUtil.Ctermfunc.empty p </ Rat.zero then nth les n
  9.1218 -                     else raise Failure "trivial_axiom: Not a trivial axiom"
  9.1219 -     | RealArith.Axiom_lt n => if eval FuncUtil.Ctermfunc.empty p <=/ Rat.zero then nth lts n
  9.1220 -                     else raise Failure "trivial_axiom: Not a trivial axiom"
  9.1221 -     | _ => error "trivial_axiom: Not a trivial axiom"
  9.1222 -   in
  9.1223 -  (let val th = tryfind trivial_axiom (keq @ klep @ kltp)
  9.1224 -   in
  9.1225 -    (fconv_rule (arg_conv (arg1_conv real_poly_conv) then_conv Numeral_Simprocs.field_comp_conv) th, RealArith.Trivial)
  9.1226 -   end)
  9.1227 -   handle Failure _ =>
  9.1228 -     (let val proof =
  9.1229 -       (case proof_method of Certificate certs =>
  9.1230 -         (* choose certificate *)
  9.1231 -         let
  9.1232 -           fun chose_cert [] (RealArith.Cert c) = c
  9.1233 -             | chose_cert (RealArith.Left::s) (RealArith.Branch (l, _)) = chose_cert s l
  9.1234 -             | chose_cert (RealArith.Right::s) (RealArith.Branch (_, r)) = chose_cert s r
  9.1235 -             | chose_cert _ _ = error "certificate tree in invalid form"
  9.1236 -         in
  9.1237 -           chose_cert cert_choice certs
  9.1238 -         end
  9.1239 -       | Prover prover =>
  9.1240 -         (* call prover *)
  9.1241 -         let
  9.1242 -          val pol = fold_rev poly_mul (map fst ltp) (poly_const Rat.one)
  9.1243 -          val leq = lep @ ltp
  9.1244 -          fun tryall d =
  9.1245 -           let val e = multidegree pol
  9.1246 -               val k = if e = 0 then 0 else d div e
  9.1247 -               val eq' = map fst eq
  9.1248 -           in tryfind (fn i => (d,i,real_positivnullstellensatz_general ctxt prover false d eq' leq
  9.1249 -                                 (poly_neg(poly_pow pol i))))
  9.1250 -                   (0 upto k)
  9.1251 -           end
  9.1252 -         val (d,i,(cert_ideal,cert_cone)) = deepen tryall 0
  9.1253 -         val proofs_ideal =
  9.1254 -           map2 (fn q => fn (p,ax) => RealArith.Eqmul(q,ax)) cert_ideal eq
  9.1255 -         val proofs_cone = map cterm_of_sos cert_cone
  9.1256 -         val proof_ne = if null ltp then RealArith.Rational_lt Rat.one else
  9.1257 -           let val p = foldr1 RealArith.Product (map snd ltp)
  9.1258 -           in  funpow i (fn q => RealArith.Product(p,q)) (RealArith.Rational_lt Rat.one)
  9.1259 -           end
  9.1260 -         in
  9.1261 -           foldr1 RealArith.Sum (proof_ne :: proofs_ideal @ proofs_cone)
  9.1262 -         end)
  9.1263 -     in
  9.1264 -        (translator (eqs,les,lts) proof, RealArith.Cert proof)
  9.1265 -     end)
  9.1266 -   end
  9.1267 - in mainf end
  9.1268 -end
  9.1269 -
  9.1270 -fun C f x y = f y x;
  9.1271 -  (* FIXME : This is very bad!!!*)
  9.1272 -fun subst_conv eqs t =
  9.1273 - let
  9.1274 -  val t' = fold (Thm.cabs o Thm.lhs_of) eqs t
  9.1275 - in Conv.fconv_rule (Thm.beta_conversion true) (fold (C Thm.combination) eqs (Thm.reflexive t'))
  9.1276 - end
  9.1277 -
  9.1278 -(* A wrapper that tries to substitute away variables first.                  *)
  9.1279 -
  9.1280 -local
  9.1281 - open Conv
  9.1282 -  fun simple_cterm_ord t u = Term_Ord.fast_term_ord (term_of t, term_of u) = LESS
  9.1283 - val concl = Thm.dest_arg o cprop_of
  9.1284 - val shuffle1 =
  9.1285 -   fconv_rule (rewr_conv @{lemma "(a + x == y) == (x == y - (a::real))" by (atomize (full)) (simp add: field_simps) })
  9.1286 - val shuffle2 =
  9.1287 -    fconv_rule (rewr_conv @{lemma "(x + a == y) ==  (x == y - (a::real))" by (atomize (full)) (simp add: field_simps)})
  9.1288 - fun substitutable_monomial fvs tm = case term_of tm of
  9.1289 -    Free(_,@{typ real}) => if not (member (op aconvc) fvs tm) then (Rat.one,tm)
  9.1290 -                           else raise Failure "substitutable_monomial"
  9.1291 -  | @{term "op * :: real => _"}$c$(t as Free _ ) =>
  9.1292 -     if RealArith.is_ratconst (Thm.dest_arg1 tm) andalso not (member (op aconvc) fvs (Thm.dest_arg tm))
  9.1293 -         then (RealArith.dest_ratconst (Thm.dest_arg1 tm),Thm.dest_arg tm) else raise Failure "substitutable_monomial"
  9.1294 -  | @{term "op + :: real => _"}$s$t =>
  9.1295 -       (substitutable_monomial (Thm.add_cterm_frees (Thm.dest_arg tm) fvs) (Thm.dest_arg1 tm)
  9.1296 -        handle Failure _ => substitutable_monomial (Thm.add_cterm_frees (Thm.dest_arg1 tm) fvs) (Thm.dest_arg tm))
  9.1297 -  | _ => raise Failure "substitutable_monomial"
  9.1298 -
  9.1299 -  fun isolate_variable v th =
  9.1300 -   let val w = Thm.dest_arg1 (cprop_of th)
  9.1301 -   in if v aconvc w then th
  9.1302 -      else case term_of w of
  9.1303 -           @{term "op + :: real => _"}$s$t =>
  9.1304 -              if Thm.dest_arg1 w aconvc v then shuffle2 th
  9.1305 -              else isolate_variable v (shuffle1 th)
  9.1306 -          | _ => error "isolate variable : This should not happen?"
  9.1307 -   end
  9.1308 -in
  9.1309 -
  9.1310 -fun real_nonlinear_subst_prover prover ctxt =
  9.1311 - let
  9.1312 -  val {add,mul,neg,pow,sub,main} =  Semiring_Normalizer.semiring_normalizers_ord_wrapper ctxt
  9.1313 -      (the (Semiring_Normalizer.match ctxt @{cterm "(0::real) + 1"}))
  9.1314 -     simple_cterm_ord
  9.1315 -
  9.1316 -  val (real_poly_add_conv,real_poly_mul_conv,real_poly_neg_conv,
  9.1317 -       real_poly_pow_conv,real_poly_sub_conv,real_poly_conv) = (add,mul,neg,pow,sub,main)
  9.1318 -
  9.1319 -  fun make_substitution th =
  9.1320 -   let
  9.1321 -    val (c,v) = substitutable_monomial [] (Thm.dest_arg1(concl th))
  9.1322 -    val th1 = Drule.arg_cong_rule (Thm.capply @{cterm "op * :: real => _"} (RealArith.cterm_of_rat (Rat.inv c))) (mk_meta_eq th)
  9.1323 -    val th2 = fconv_rule (binop_conv real_poly_mul_conv) th1
  9.1324 -   in fconv_rule (arg_conv real_poly_conv) (isolate_variable v th2)
  9.1325 -   end
  9.1326 -   fun oprconv cv ct =
  9.1327 -    let val g = Thm.dest_fun2 ct
  9.1328 -    in if g aconvc @{cterm "op <= :: real => _"}
  9.1329 -         orelse g aconvc @{cterm "op < :: real => _"}
  9.1330 -       then arg_conv cv ct else arg1_conv cv ct
  9.1331 -    end
  9.1332 -  fun mainf cert_choice translator =
  9.1333 -   let
  9.1334 -    fun substfirst(eqs,les,lts) =
  9.1335 -      ((let
  9.1336 -           val eth = tryfind make_substitution eqs
  9.1337 -           val modify = fconv_rule (arg_conv (oprconv(subst_conv [eth] then_conv real_poly_conv)))
  9.1338 -       in  substfirst
  9.1339 -             (filter_out (fn t => (Thm.dest_arg1 o Thm.dest_arg o cprop_of) t
  9.1340 -                                   aconvc @{cterm "0::real"}) (map modify eqs),
  9.1341 -                                   map modify les,map modify lts)
  9.1342 -       end)
  9.1343 -       handle Failure  _ => real_nonlinear_prover prover ctxt cert_choice translator (rev eqs, rev les, rev lts))
  9.1344 -    in substfirst
  9.1345 -   end
  9.1346 -
  9.1347 -
  9.1348 - in mainf
  9.1349 - end
  9.1350 -
  9.1351 -(* Overall function. *)
  9.1352 -
  9.1353 -fun real_sos prover ctxt =
  9.1354 -  RealArith.gen_prover_real_arith ctxt (real_nonlinear_subst_prover prover ctxt)
  9.1355 -end;
  9.1356 -
  9.1357 -val known_sos_constants =
  9.1358 -  [@{term "op ==>"}, @{term "Trueprop"},
  9.1359 -   @{term HOL.implies}, @{term HOL.conj}, @{term HOL.disj},
  9.1360 -   @{term "Not"}, @{term "op = :: bool => _"},
  9.1361 -   @{term "All :: (real => _) => _"}, @{term "Ex :: (real => _) => _"},
  9.1362 -   @{term "op = :: real => _"}, @{term "op < :: real => _"},
  9.1363 -   @{term "op <= :: real => _"},
  9.1364 -   @{term "op + :: real => _"}, @{term "op - :: real => _"},
  9.1365 -   @{term "op * :: real => _"}, @{term "uminus :: real => _"},
  9.1366 -   @{term "op / :: real => _"}, @{term "inverse :: real => _"},
  9.1367 -   @{term "op ^ :: real => _"}, @{term "abs :: real => _"},
  9.1368 -   @{term "min :: real => _"}, @{term "max :: real => _"},
  9.1369 -   @{term "0::real"}, @{term "1::real"}, @{term "number_of :: int => real"},
  9.1370 -   @{term "number_of :: int => nat"},
  9.1371 -   @{term "Int.Bit0"}, @{term "Int.Bit1"},
  9.1372 -   @{term "Int.Pls"}, @{term "Int.Min"}];
  9.1373 -
  9.1374 -fun check_sos kcts ct =
  9.1375 - let
  9.1376 -  val t = term_of ct
  9.1377 -  val _ = if not (null (Term.add_tfrees t [])
  9.1378 -                  andalso null (Term.add_tvars t []))
  9.1379 -          then error "SOS: not sos. Additional type varables" else ()
  9.1380 -  val fs = Term.add_frees t []
  9.1381 -  val _ = if exists (fn ((_,T)) => not (T = @{typ "real"})) fs
  9.1382 -          then error "SOS: not sos. Variables with type not real" else ()
  9.1383 -  val vs = Term.add_vars t []
  9.1384 -  val _ = if exists (fn ((_,T)) => not (T = @{typ "real"})) vs
  9.1385 -          then error "SOS: not sos. Variables with type not real" else ()
  9.1386 -  val ukcs = subtract (fn (t,p) => Const p aconv t) kcts (Term.add_consts t [])
  9.1387 -  val _ = if  null ukcs then ()
  9.1388 -              else error ("SOSO: Unknown constants in Subgoal:" ^ commas (map fst ukcs))
  9.1389 -in () end
  9.1390 -
  9.1391 -fun core_sos_tac print_cert prover = SUBPROOF (fn {concl, context, ...} =>
  9.1392 -  let
  9.1393 -    val _ = check_sos known_sos_constants concl
  9.1394 -    val (ths, certificates) = real_sos prover context (Thm.dest_arg concl)
  9.1395 -    val _ = print_cert certificates
  9.1396 -  in rtac ths 1 end)
  9.1397 -
  9.1398 -fun default_SOME f NONE v = SOME v
  9.1399 -  | default_SOME f (SOME v) _ = SOME v;
  9.1400 -
  9.1401 -fun lift_SOME f NONE a = f a
  9.1402 -  | lift_SOME f (SOME a) _ = SOME a;
  9.1403 -
  9.1404 -
  9.1405 -local
  9.1406 - val is_numeral = can (HOLogic.dest_number o term_of)
  9.1407 -in
  9.1408 -fun get_denom b ct = case term_of ct of
  9.1409 -  @{term "op / :: real => _"} $ _ $ _ =>
  9.1410 -     if is_numeral (Thm.dest_arg ct) then get_denom b (Thm.dest_arg1 ct)
  9.1411 -     else default_SOME (get_denom b) (get_denom b (Thm.dest_arg ct))   (Thm.dest_arg ct, b)
  9.1412 - | @{term "op < :: real => _"} $ _ $ _ => lift_SOME (get_denom true) (get_denom true (Thm.dest_arg ct)) (Thm.dest_arg1 ct)
  9.1413 - | @{term "op <= :: real => _"} $ _ $ _ => lift_SOME (get_denom true) (get_denom true (Thm.dest_arg ct)) (Thm.dest_arg1 ct)
  9.1414 - | _ $ _ => lift_SOME (get_denom b) (get_denom b (Thm.dest_fun ct)) (Thm.dest_arg ct)
  9.1415 - | _ => NONE
  9.1416 -end;
  9.1417 -
  9.1418 -fun elim_one_denom_tac ctxt =
  9.1419 -CSUBGOAL (fn (P,i) =>
  9.1420 - case get_denom false P of
  9.1421 -   NONE => no_tac
  9.1422 - | SOME (d,ord) =>
  9.1423 -     let
  9.1424 -      val ss = simpset_of ctxt addsimps @{thms field_simps}
  9.1425 -               addsimps [@{thm nonzero_power_divide}, @{thm power_divide}]
  9.1426 -      val th = instantiate' [] [SOME d, SOME (Thm.dest_arg P)]
  9.1427 -         (if ord then @{lemma "(d=0 --> P) & (d>0 --> P) & (d<(0::real) --> P) ==> P" by auto}
  9.1428 -          else @{lemma "(d=0 --> P) & (d ~= (0::real) --> P) ==> P" by blast})
  9.1429 -     in rtac th i THEN Simplifier.asm_full_simp_tac ss i end);
  9.1430 -
  9.1431 -fun elim_denom_tac ctxt i = REPEAT (elim_one_denom_tac ctxt i);
  9.1432 -
  9.1433 -fun sos_tac print_cert prover ctxt =
  9.1434 -  Object_Logic.full_atomize_tac THEN'
  9.1435 -  elim_denom_tac ctxt THEN'
  9.1436 -  core_sos_tac print_cert prover ctxt;
  9.1437 -
  9.1438 -end;
    10.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
    10.2 +++ b/src/HOL/Library/Sum_of_Squares.thy	Sat Jan 08 17:39:51 2011 +0100
    10.3 @@ -0,0 +1,159 @@
    10.4 +(*  Title:      HOL/Library/Sum_of_Squares.thy
    10.5 +    Author:     Amine Chaieb, University of Cambridge
    10.6 +    Author:     Philipp Meyer, TU Muenchen
    10.7 +*)
    10.8 +
    10.9 +header {* A decision method for universal multivariate real arithmetic with addition, 
   10.10 +  multiplication and ordering using semidefinite programming *}
   10.11 +
   10.12 +theory Sum_of_Squares
   10.13 +imports Complex_Main
   10.14 +uses
   10.15 +  "positivstellensatz.ML"
   10.16 +  "Sum_of_Squares/sum_of_squares.ML"
   10.17 +  "Sum_of_Squares/positivstellensatz_tools.ML"
   10.18 +  "Sum_of_Squares/sos_wrapper.ML"
   10.19 +begin
   10.20 +
   10.21 +text {*
   10.22 +  In order to use the method sos, call it with @{text "(sos
   10.23 +  remote_csdp)"} to use the remote solver.  Or install CSDP
   10.24 +  (https://projects.coin-or.org/Csdp), configure the Isabelle setting
   10.25 +  @{text CSDP_EXE}, and call it with @{text "(sos csdp)"}.  By
   10.26 +  default, sos calls @{text remote_csdp}.  This can take of the order
   10.27 +  of a minute for one sos call, because sos calls CSDP repeatedly.  If
   10.28 +  you install CSDP locally, sos calls typically takes only a few
   10.29 +  seconds.
   10.30 +  sos generates a certificate which can be used to repeat the proof
   10.31 +  without calling an external prover.
   10.32 +*}
   10.33 +
   10.34 +setup Sum_of_Squares.setup
   10.35 +setup SOS_Wrapper.setup
   10.36 +
   10.37 +text {* Tests *}
   10.38 +
   10.39 +lemma "(3::real) * x + 7 * a < 4 & 3 < 2 * x \<Longrightarrow> a < 0"
   10.40 +by (sos_cert "((R<1 + (((A<1 * R<1) * (R<2 * [1]^2)) + (((A<0 * R<1) * (R<3 * [1]^2)) + ((A<=0 * R<1) * (R<14 * [1]^2))))))")
   10.41 +
   10.42 +lemma "a1 >= 0 & a2 >= 0 \<and> (a1 * a1 + a2 * a2 = b1 * b1 + b2 * b2 + 2) \<and> (a1 * b1 + a2 * b2 = 0) --> a1 * a2 - b1 * b2 >= (0::real)"
   10.43 +by (sos_cert "(((A<0 * R<1) + (([~1/2*a1*b2 + ~1/2*a2*b1] * A=0) + (([~1/2*a1*a2 + 1/2*b1*b2] * A=1) + (((A<0 * R<1) * ((R<1/2 * [b2]^2) + (R<1/2 * [b1]^2))) + ((A<=0 * (A<=1 * R<1)) * ((R<1/2 * [b2]^2) + ((R<1/2 * [b1]^2) + ((R<1/2 * [a2]^2) + (R<1/2 * [a1]^2))))))))))")
   10.44 +
   10.45 +lemma "(3::real) * x + 7 * a < 4 & 3 < 2 * x --> a < 0"
   10.46 +by (sos_cert "((R<1 + (((A<1 * R<1) * (R<2 * [1]^2)) + (((A<0 * R<1) * (R<3 * [1]^2)) + ((A<=0 * R<1) * (R<14 * [1]^2))))))")
   10.47 +
   10.48 +lemma "(0::real) <= x & x <= 1 & 0 <= y & y <= 1  --> x^2 + y^2 < 1 |(x - 1)^2 + y^2 < 1 | x^2 + (y - 1)^2 < 1 | (x - 1)^2 + (y - 1)^2 < 1"
   10.49 +by (sos_cert "((R<1 + (((A<=3 * (A<=4 * R<1)) * (R<1 * [1]^2)) + (((A<=2 * (A<=7 * R<1)) * (R<1 * [1]^2)) + (((A<=1 * (A<=6 * R<1)) * (R<1 * [1]^2)) + ((A<=0 * (A<=5 * R<1)) * (R<1 * [1]^2)))))))")
   10.50 +
   10.51 +lemma "(0::real) <= x & 0 <= y & 0 <= z & x + y + z <= 3 --> x * y + x * z + y * z >= 3 * x * y * z"
   10.52 +by (sos_cert "(((A<0 * R<1) + (((A<0 * R<1) * (R<1/2 * [1]^2)) + (((A<=2 * R<1) * (R<1/2 * [~1*x + y]^2)) + (((A<=1 * R<1) * (R<1/2 * [~1*x + z]^2)) + (((A<=1 * (A<=2 * (A<=3 * R<1))) * (R<1/2 * [1]^2)) + (((A<=0 * R<1) * (R<1/2 * [~1*y + z]^2)) + (((A<=0 * (A<=2 * (A<=3 * R<1))) * (R<1/2 * [1]^2)) + ((A<=0 * (A<=1 * (A<=3 * R<1))) * (R<1/2 * [1]^2))))))))))")
   10.53 +
   10.54 +lemma "((x::real)^2 + y^2 + z^2 = 1) --> (x + y + z)^2 <= 3"
   10.55 +by (sos_cert "(((A<0 * R<1) + (([~3] * A=0) + (R<1 * ((R<2 * [~1/2*x + ~1/2*y + z]^2) + (R<3/2 * [~1*x + y]^2))))))")
   10.56 +
   10.57 +lemma "(w^2 + x^2 + y^2 + z^2 = 1) --> (w + x + y + z)^2 <= (4::real)"
   10.58 +by (sos_cert "(((A<0 * R<1) + (([~4] * A=0) + (R<1 * ((R<3 * [~1/3*w + ~1/3*x + ~1/3*y + z]^2) + ((R<8/3 * [~1/2*w + ~1/2*x + y]^2) + (R<2 * [~1*w + x]^2)))))))")
   10.59 +
   10.60 +lemma "(x::real) >= 1 & y >= 1 --> x * y >= x + y - 1"
   10.61 +by (sos_cert "(((A<0 * R<1) + ((A<=0 * (A<=1 * R<1)) * (R<1 * [1]^2))))")
   10.62 +
   10.63 +lemma "(x::real) > 1 & y > 1 --> x * y > x + y - 1"
   10.64 +by (sos_cert "((((A<0 * A<1) * R<1) + ((A<=0 * R<1) * (R<1 * [1]^2))))") 
   10.65 +
   10.66 +lemma "abs(x) <= 1 --> abs(64 * x^7 - 112 * x^5 + 56 * x^3 - 7 * x) <= (1::real)"
   10.67 +by (sos_cert "((((A<0 * R<1) + ((A<=1 * R<1) * (R<1 * [~8*x^3 + ~4*x^2 + 4*x + 1]^2)))) & ((((A<0 * A<1) * R<1) + ((A<=1 * (A<0 * R<1)) * (R<1 * [8*x^3 + ~4*x^2 + ~4*x + 1]^2)))))")
   10.68 +
   10.69 +(* ------------------------------------------------------------------------- *)
   10.70 +(* One component of denominator in dodecahedral example.                     *)
   10.71 +(* ------------------------------------------------------------------------- *)
   10.72 +
   10.73 +lemma "2 <= x & x <= 125841 / 50000 & 2 <= y & y <= 125841 / 50000 & 2 <= z & z <= 125841 / 50000 --> 2 * (x * z + x * y + y * z) - (x * x + y * y + z * z) >= (0::real)"
   10.74 +by (sos_cert "(((A<0 * R<1) + ((R<1 * ((R<5749028157/5000000000 * [~25000/222477*x + ~25000/222477*y + ~25000/222477*z + 1]^2) + ((R<864067/1779816 * [419113/864067*x + 419113/864067*y + z]^2) + ((R<320795/864067 * [419113/1283180*x + y]^2) + (R<1702293/5132720 * [x]^2))))) + (((A<=4 * (A<=5 * R<1)) * (R<3/2 * [1]^2)) + (((A<=3 * (A<=5 * R<1)) * (R<1/2 * [1]^2)) + (((A<=2 * (A<=4 * R<1)) * (R<1 * [1]^2)) + (((A<=2 * (A<=3 * R<1)) * (R<3/2 * [1]^2)) + (((A<=1 * (A<=5 * R<1)) * (R<1/2 * [1]^2)) + (((A<=1 * (A<=3 * R<1)) * (R<1/2 * [1]^2)) + (((A<=0 * (A<=4 * R<1)) * (R<1 * [1]^2)) + (((A<=0 * (A<=2 * R<1)) * (R<1 * [1]^2)) + ((A<=0 * (A<=1 * R<1)) * (R<3/2 * [1]^2)))))))))))))")
   10.75 +
   10.76 +(* ------------------------------------------------------------------------- *)
   10.77 +(* Over a larger but simpler interval.                                       *)
   10.78 +(* ------------------------------------------------------------------------- *)
   10.79 +
   10.80 +lemma "(2::real) <= x & x <= 4 & 2 <= y & y <= 4 & 2 <= z & z <= 4 --> 0 <= 2 * (x * z + x * y + y * z) - (x * x + y * y + z * z)"
   10.81 +by (sos_cert "((R<1 + ((R<1 * ((R<1 * [~1/6*x + ~1/6*y + ~1/6*z + 1]^2) + ((R<1/18 * [~1/2*x + ~1/2*y + z]^2) + (R<1/24 * [~1*x + y]^2)))) + (((A<0 * R<1) * (R<1/12 * [1]^2)) + (((A<=4 * (A<=5 * R<1)) * (R<1/6 * [1]^2)) + (((A<=2 * (A<=4 * R<1)) * (R<1/6 * [1]^2)) + (((A<=2 * (A<=3 * R<1)) * (R<1/6 * [1]^2)) + (((A<=0 * (A<=4 * R<1)) * (R<1/6 * [1]^2)) + (((A<=0 * (A<=2 * R<1)) * (R<1/6 * [1]^2)) + ((A<=0 * (A<=1 * R<1)) * (R<1/6 * [1]^2)))))))))))")
   10.82 +
   10.83 +(* ------------------------------------------------------------------------- *)
   10.84 +(* We can do 12. I think 12 is a sharp bound; see PP's certificate.          *)
   10.85 +(* ------------------------------------------------------------------------- *)
   10.86 +
   10.87 +lemma "2 <= (x::real) & x <= 4 & 2 <= y & y <= 4 & 2 <= z & z <= 4 --> 12 <= 2 * (x * z + x * y + y * z) - (x * x + y * y + z * z)"
   10.88 +by (sos_cert "(((A<0 * R<1) + (((A<=4 * R<1) * (R<2/3 * [1]^2)) + (((A<=4 * (A<=5 * R<1)) * (R<1 * [1]^2)) + (((A<=3 * (A<=4 * R<1)) * (R<1/3 * [1]^2)) + (((A<=2 * R<1) * (R<2/3 * [1]^2)) + (((A<=2 * (A<=5 * R<1)) * (R<1/3 * [1]^2)) + (((A<=2 * (A<=4 * R<1)) * (R<8/3 * [1]^2)) + (((A<=2 * (A<=3 * R<1)) * (R<1 * [1]^2)) + (((A<=1 * (A<=4 * R<1)) * (R<1/3 * [1]^2)) + (((A<=1 * (A<=2 * R<1)) * (R<1/3 * [1]^2)) + (((A<=0 * R<1) * (R<2/3 * [1]^2)) + (((A<=0 * (A<=5 * R<1)) * (R<1/3 * [1]^2)) + (((A<=0 * (A<=4 * R<1)) * (R<8/3 * [1]^2)) + (((A<=0 * (A<=3 * R<1)) * (R<1/3 * [1]^2)) + (((A<=0 * (A<=2 * R<1)) * (R<8/3 * [1]^2)) + ((A<=0 * (A<=1 * R<1)) * (R<1 * [1]^2))))))))))))))))))")
   10.89 +
   10.90 +(* ------------------------------------------------------------------------- *)
   10.91 +(* Inequality from sci.math (see "Leon-Sotelo, por favor").                  *)
   10.92 +(* ------------------------------------------------------------------------- *)
   10.93 +
   10.94 +lemma "0 <= (x::real) & 0 <= y & (x * y = 1) --> x + y <= x^2 + y^2"
   10.95 +by (sos_cert "(((A<0 * R<1) + (([1] * A=0) + (R<1 * ((R<1 * [~1/2*x + ~1/2*y + 1]^2) + (R<3/4 * [~1*x + y]^2))))))") 
   10.96 +
   10.97 +lemma "0 <= (x::real) & 0 <= y & (x * y = 1) --> x * y * (x + y) <= x^2 + y^2"
   10.98 +by (sos_cert "(((A<0 * R<1) + (([~1*x + ~1*y + 1] * A=0) + (R<1 * ((R<1 * [~1/2*x + ~1/2*y + 1]^2) + (R<3/4 * [~1*x + y]^2))))))") 
   10.99 +
  10.100 +lemma "0 <= (x::real) & 0 <= y --> x * y * (x + y)^2 <= (x^2 + y^2)^2"
  10.101 +by (sos_cert "(((A<0 * R<1) + (R<1 * ((R<1 * [~1/2*x^2 + y^2 + ~1/2*x*y]^2) + (R<3/4 * [~1*x^2 + x*y]^2)))))")
  10.102 +
  10.103 +lemma "(0::real) <= a & 0 <= b & 0 <= c & c * (2 * a + b)^3/ 27 <= x \<longrightarrow> c * a^2 * b <= x"
  10.104 +by (sos_cert "(((A<0 * R<1) + (((A<=3 * R<1) * (R<1 * [1]^2)) + (((A<=1 * (A<=2 * R<1)) * (R<1/27 * [~1*a + b]^2)) + ((A<=0 * (A<=2 * R<1)) * (R<8/27 * [~1*a + b]^2))))))")
  10.105 + 
  10.106 +lemma "(0::real) < x --> 0 < 1 + x + x^2"
  10.107 +by (sos_cert "((R<1 + ((R<1 * (R<1 * [x]^2)) + (((A<0 * R<1) * (R<1 * [1]^2)) + ((A<=0 * R<1) * (R<1 * [1]^2))))))")
  10.108 +
  10.109 +lemma "(0::real) <= x --> 0 < 1 + x + x^2"
  10.110 +by (sos_cert "((R<1 + ((R<1 * (R<1 * [x]^2)) + (((A<=1 * R<1) * (R<1 * [1]^2)) + ((A<=0 * R<1) * (R<1 * [1]^2))))))")
  10.111 +
  10.112 +lemma "(0::real) < 1 + x^2"
  10.113 +by (sos_cert "((R<1 + ((R<1 * (R<1 * [x]^2)) + ((A<=0 * R<1) * (R<1 * [1]^2)))))")
  10.114 +
  10.115 +lemma "(0::real) <= 1 + 2 * x + x^2"
  10.116 +by (sos_cert "(((A<0 * R<1) + (R<1 * (R<1 * [x + 1]^2))))")
  10.117 +
  10.118 +lemma "(0::real) < 1 + abs x"
  10.119 +by (sos_cert "((R<1 + (((A<=1 * R<1) * (R<1/2 * [1]^2)) + ((A<=0 * R<1) * (R<1/2 * [1]^2)))))")
  10.120 +
  10.121 +lemma "(0::real) < 1 + (1 + x)^2 * (abs x)"
  10.122 +by (sos_cert "(((R<1 + (((A<=1 * R<1) * (R<1 * [1]^2)) + ((A<=0 * R<1) * (R<1 * [x + 1]^2))))) & ((R<1 + (((A<0 * R<1) * (R<1 * [x + 1]^2)) + ((A<=0 * R<1) * (R<1 * [1]^2))))))")
  10.123 +
  10.124 +
  10.125 +
  10.126 +lemma "abs ((1::real) + x^2) = (1::real) + x^2"
  10.127 +by (sos_cert "(() & (((R<1 + ((R<1 * (R<1 * [x]^2)) + ((A<1 * R<1) * (R<1/2 * [1]^2))))) & ((R<1 + ((R<1 * (R<1 * [x]^2)) + ((A<0 * R<1) * (R<1 * [1]^2)))))))")
  10.128 +lemma "(3::real) * x + 7 * a < 4 \<and> 3 < 2 * x \<longrightarrow> a < 0"
  10.129 +by (sos_cert "((R<1 + (((A<1 * R<1) * (R<2 * [1]^2)) + (((A<0 * R<1) * (R<3 * [1]^2)) + ((A<=0 * R<1) * (R<14 * [1]^2))))))")
  10.130 +
  10.131 +lemma "(0::real) < x --> 1 < y --> y * x <= z --> x < z"
  10.132 +by (sos_cert "((((A<0 * A<1) * R<1) + (((A<=1 * R<1) * (R<1 * [1]^2)) + ((A<=0 * R<1) * (R<1 * [1]^2)))))")
  10.133 +lemma "(1::real) < x --> x^2 < y --> 1 < y"
  10.134 +by (sos_cert "((((A<0 * A<1) * R<1) + ((R<1 * ((R<1/10 * [~2*x + y + 1]^2) + (R<1/10 * [~1*x + y]^2))) + (((A<1 * R<1) * (R<1/2 * [1]^2)) + (((A<0 * R<1) * (R<1 * [x]^2)) + (((A<=0 * R<1) * ((R<1/10 * [x + 1]^2) + (R<1/10 * [x]^2))) + (((A<=0 * (A<1 * R<1)) * (R<1/5 * [1]^2)) + ((A<=0 * (A<0 * R<1)) * (R<1/5 * [1]^2)))))))))")
  10.135 +lemma "(b::real)^2 < 4 * a * c --> ~(a * x^2 + b * x + c = 0)"
  10.136 +by (sos_cert "(((A<0 * R<1) + (R<1 * (R<1 * [2*a*x + b]^2))))")
  10.137 +lemma "(b::real)^2 < 4 * a * c --> ~(a * x^2 + b * x + c = 0)"
  10.138 +by (sos_cert "(((A<0 * R<1) + (R<1 * (R<1 * [2*a*x + b]^2))))")
  10.139 +lemma "((a::real) * x^2 + b * x + c = 0) --> b^2 >= 4 * a * c"
  10.140 +by (sos_cert "(((A<0 * R<1) + (R<1 * (R<1 * [2*a*x + b]^2))))")
  10.141 +lemma "(0::real) <= b & 0 <= c & 0 <= x & 0 <= y & (x^2 = c) & (y^2 = a^2 * c + b) --> a * c <= y * x"
  10.142 +by (sos_cert "(((A<0 * (A<0 * R<1)) + (((A<=2 * (A<=3 * (A<0 * R<1))) * (R<2 * [1]^2)) + ((A<=0 * (A<=1 * R<1)) * (R<1 * [1]^2)))))")
  10.143 +lemma "abs(x - z) <= e & abs(y - z) <= e & 0 <= u & 0 <= v & (u + v = 1) --> abs((u * x + v * y) - z) <= (e::real)"
  10.144 +by (sos_cert "((((A<0 * R<1) + (((A<=3 * (A<=6 * R<1)) * (R<1 * [1]^2)) + ((A<=1 * (A<=5 * R<1)) * (R<1 * [1]^2))))) & ((((A<0 * A<1) * R<1) + (((A<=3 * (A<=5 * (A<0 * R<1))) * (R<1 * [1]^2)) + ((A<=1 * (A<=4 * (A<0 * R<1))) * (R<1 * [1]^2))))))")
  10.145 +
  10.146 +
  10.147 +(* lemma "((x::real) - y - 2 * x^4 = 0) & 0 <= x & x <= 2 & 0 <= y & y <= 3 --> y^2 - 7 * y - 12 * x + 17 >= 0" by sos *) (* Too hard?*)
  10.148 +
  10.149 +lemma "(0::real) <= x --> (1 + x + x^2)/(1 + x^2) <= 1 + x"
  10.150 +by (sos_cert "(((((A<0 * A<1) * R<1) + ((A<=0 * (A<0 * R<1)) * (R<1 * [x]^2)))) & ((R<1 + ((R<1 * (R<1 * [x]^2)) + ((A<0 * R<1) * (R<1 * [1]^2))))))")
  10.151 +
  10.152 +lemma "(0::real) <= x --> 1 - x <= 1 / (1 + x + x^2)"
  10.153 +by (sos_cert "(((R<1 + (([~4/3] * A=0) + ((R<1 * ((R<1/3 * [3/2*x + 1]^2) + (R<7/12 * [x]^2))) + ((A<=0 * R<1) * (R<1/3 * [1]^2)))))) & (((((A<0 * A<1) * R<1) + ((A<=0 * (A<0 * R<1)) * (R<1 * [x]^2)))) & ((R<1 + ((R<1 * (R<1 * [x]^2)) + (((A<0 * R<1) * (R<1 * [1]^2)) + ((A<=0 * R<1) * (R<1 * [1]^2))))))))")
  10.154 +
  10.155 +lemma "(x::real) <= 1 / 2 --> - x - 2 * x^2 <= - x / (1 - x)"
  10.156 +by (sos_cert "((((A<0 * A<1) * R<1) + ((A<=0 * (A<0 * R<1)) * (R<1 * [x]^2))))")
  10.157 +
  10.158 +lemma "4*r^2 = p^2 - 4*q & r >= (0::real) & x^2 + p*x + q = 0 --> 2*(x::real) = - p + 2*r | 2*x = -p - 2*r"
  10.159 +by (sos_cert "((((((A<0 * A<1) * R<1) + ([~4] * A=0))) & ((((A<0 * A<1) * R<1) + ([4] * A=0)))) & (((((A<0 * A<1) * R<1) + ([4] * A=0))) & ((((A<0 * A<1) * R<1) + ([~4] * A=0)))))")
  10.160 +
  10.161 +end
  10.162 +
    11.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
    11.2 +++ b/src/HOL/Library/Sum_of_Squares/etc/settings	Sat Jan 08 17:39:51 2011 +0100
    11.3 @@ -0,0 +1,3 @@
    11.4 +# -*- shell-script -*- :mode=shellscript:
    11.5 +
    11.6 +ISABELLE_SUM_OF_SQUARES="$COMPONENT"
    12.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
    12.2 +++ b/src/HOL/Library/Sum_of_Squares/neos_csdp_client	Sat Jan 08 17:39:51 2011 +0100
    12.3 @@ -0,0 +1,87 @@
    12.4 +#!/usr/bin/env python
    12.5 +import sys
    12.6 +import signal
    12.7 +import xmlrpclib
    12.8 +import time
    12.9 +import re
   12.10 +
   12.11 +# Neos server config
   12.12 +NEOS_HOST="neos.mcs.anl.gov"
   12.13 +NEOS_PORT=3332
   12.14 +
   12.15 +neos=xmlrpclib.Server("http://%s:%d" % (NEOS_HOST, NEOS_PORT))
   12.16 +
   12.17 +jobNumber = 0
   12.18 +password = ""
   12.19 +inputfile = None
   12.20 +outputfile = None
   12.21 +# interrupt handler
   12.22 +def cleanup(signal, frame):
   12.23 +  sys.stdout.write("Caught interrupt, cleaning up\n")
   12.24 +  if jobNumber != 0:
   12.25 +    neos.killJob(jobNumber, password)
   12.26 +  if inputfile != None:
   12.27 +    inputfile.close()
   12.28 +  if outputfile != None:
   12.29 +    outputfile.close()
   12.30 +  sys.exit(21)
   12.31 +
   12.32 +signal.signal(signal.SIGHUP, cleanup)
   12.33 +signal.signal(signal.SIGINT, cleanup)
   12.34 +signal.signal(signal.SIGQUIT, cleanup)
   12.35 +signal.signal(signal.SIGTERM, cleanup)
   12.36 +
   12.37 +if len(sys.argv) <> 3:
   12.38 +  sys.stderr.write("Usage: neos_csdp_client <input_filename> <output_filename>\n")
   12.39 +  sys.exit(19)
   12.40 +
   12.41 +xml_pre = "<document>\n<category>sdp</category>\n<solver>csdp</solver>\n<inputMethod>SPARSE_SDPA</inputMethod>\n<dat><![CDATA["
   12.42 +xml_post = "]]></dat>\n</document>\n"
   12.43 +xml = xml_pre
   12.44 +inputfile = open(sys.argv[1],"r")
   12.45 +buffer = 1
   12.46 +while buffer:
   12.47 +  buffer = inputfile.read()
   12.48 +  xml += buffer
   12.49 +inputfile.close()
   12.50 +xml += xml_post
   12.51 +
   12.52 +(jobNumber,password) = neos.submitJob(xml)
   12.53 +
   12.54 +if jobNumber == 0:
   12.55 +  sys.stdout.write("error submitting job: %s" % password)
   12.56 +  sys.exit(20)
   12.57 +else:
   12.58 +  sys.stdout.write("jobNumber = %d\tpassword = %s\n" % (jobNumber,password))
   12.59 +
   12.60 +offset=0
   12.61 +messages = ""
   12.62 +status="Waiting"
   12.63 +while status == "Running" or status=="Waiting":
   12.64 +  time.sleep(1)
   12.65 +  (msg,offset) = neos.getIntermediateResults(jobNumber,password,offset)
   12.66 +  messages += msg.data
   12.67 +  sys.stdout.write(msg.data)
   12.68 +  status = neos.getJobStatus(jobNumber, password)
   12.69 +
   12.70 +msg = neos.getFinalResults(jobNumber, password).data
   12.71 +sys.stdout.write("---------- Begin CSDP Output -------------\n");
   12.72 +sys.stdout.write(msg)
   12.73 +
   12.74 +# extract solution
   12.75 +result = msg.split("Solution:")
   12.76 +if len(result) > 1:
   12.77 +  solution = result[1].strip()
   12.78 +  if solution != "":
   12.79 +    outputfile = open(sys.argv[2],"w")
   12.80 +    outputfile.write(solution)
   12.81 +    outputfile.close()
   12.82 +
   12.83 +# extract return code
   12.84 +p = re.compile(r"^Error: Command exited with non-zero status (\d+)$", re.MULTILINE)
   12.85 +m = p.search(messages)
   12.86 +if m:
   12.87 +  sys.exit(int(m.group(1)))
   12.88 +else:
   12.89 +  sys.exit(0)
   12.90 +
    13.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
    13.2 +++ b/src/HOL/Library/Sum_of_Squares/positivstellensatz_tools.ML	Sat Jan 08 17:39:51 2011 +0100
    13.3 @@ -0,0 +1,156 @@
    13.4 +(*  Title:      HOL/Library/Sum_of_Squares/positivstellensatz_tools.ML
    13.5 +    Author:     Philipp Meyer, TU Muenchen
    13.6 +
    13.7 +Functions for generating a certificate from a positivstellensatz and vice
    13.8 +versa.
    13.9 +*)
   13.10 +
   13.11 +signature POSITIVSTELLENSATZ_TOOLS =
   13.12 +sig
   13.13 +  val pss_tree_to_cert : RealArith.pss_tree -> string
   13.14 +
   13.15 +  val cert_to_pss_tree : Proof.context -> string -> RealArith.pss_tree
   13.16 +end
   13.17 +
   13.18 +
   13.19 +structure PositivstellensatzTools : POSITIVSTELLENSATZ_TOOLS =
   13.20 +struct
   13.21 +
   13.22 +(*** certificate generation ***)
   13.23 +
   13.24 +fun string_of_rat r =
   13.25 +  let
   13.26 +    val (nom, den) = Rat.quotient_of_rat r
   13.27 +  in
   13.28 +    if den = 1 then string_of_int nom
   13.29 +    else string_of_int nom ^ "/" ^ string_of_int den
   13.30 +  end
   13.31 +
   13.32 +(* map polynomials to strings *)
   13.33 +
   13.34 +fun string_of_varpow x k =
   13.35 +  let
   13.36 +    val term = term_of x
   13.37 +    val name = case term of
   13.38 +      Free (n, _) => n
   13.39 +    | _ => error "Term in monomial not free variable"
   13.40 +  in
   13.41 +    if k = 1 then name else name ^ "^" ^ string_of_int k 
   13.42 +  end
   13.43 +
   13.44 +fun string_of_monomial m = 
   13.45 + if FuncUtil.Ctermfunc.is_empty m then "1" 
   13.46 + else 
   13.47 +  let 
   13.48 +   val m' = FuncUtil.dest_monomial m
   13.49 +   val vps = fold_rev (fn (x,k) => cons (string_of_varpow x k)) m' [] 
   13.50 +  in foldr1 (fn (s, t) => s ^ "*" ^ t) vps
   13.51 +  end
   13.52 +
   13.53 +fun string_of_cmonomial (m,c) =
   13.54 +  if FuncUtil.Ctermfunc.is_empty m then string_of_rat c
   13.55 +  else if c = Rat.one then string_of_monomial m
   13.56 +  else (string_of_rat c) ^ "*" ^ (string_of_monomial m);
   13.57 +
   13.58 +fun string_of_poly p = 
   13.59 + if FuncUtil.Monomialfunc.is_empty p then "0" 
   13.60 + else
   13.61 +  let 
   13.62 +   val cms = map string_of_cmonomial
   13.63 +     (sort (prod_ord FuncUtil.monomial_order (K EQUAL)) (FuncUtil.Monomialfunc.dest p))
   13.64 +  in foldr1 (fn (t1, t2) => t1 ^ " + " ^ t2) cms
   13.65 +  end;
   13.66 +
   13.67 +fun pss_to_cert (RealArith.Axiom_eq i) = "A=" ^ string_of_int i
   13.68 +  | pss_to_cert (RealArith.Axiom_le i) = "A<=" ^ string_of_int i
   13.69 +  | pss_to_cert (RealArith.Axiom_lt i) = "A<" ^ string_of_int i
   13.70 +  | pss_to_cert (RealArith.Rational_eq r) = "R=" ^ string_of_rat r
   13.71 +  | pss_to_cert (RealArith.Rational_le r) = "R<=" ^ string_of_rat r
   13.72 +  | pss_to_cert (RealArith.Rational_lt r) = "R<" ^ string_of_rat r
   13.73 +  | pss_to_cert (RealArith.Square p) = "[" ^ string_of_poly p ^ "]^2"
   13.74 +  | pss_to_cert (RealArith.Eqmul (p, pss)) = "([" ^ string_of_poly p ^ "] * " ^ pss_to_cert pss ^ ")"
   13.75 +  | pss_to_cert (RealArith.Sum (pss1, pss2)) = "(" ^ pss_to_cert pss1 ^ " + " ^ pss_to_cert pss2 ^ ")"
   13.76 +  | pss_to_cert (RealArith.Product (pss1, pss2)) = "(" ^ pss_to_cert pss1 ^ " * " ^ pss_to_cert pss2 ^ ")"
   13.77 +
   13.78 +fun pss_tree_to_cert RealArith.Trivial = "()"
   13.79 +  | pss_tree_to_cert (RealArith.Cert pss) = "(" ^ pss_to_cert pss ^ ")"
   13.80 +  | pss_tree_to_cert (RealArith.Branch (t1, t2)) = "(" ^ pss_tree_to_cert t1 ^ " & " ^ pss_tree_to_cert t2 ^ ")"
   13.81 +
   13.82 +(*** certificate parsing ***)
   13.83 +
   13.84 +(* basic parser *)
   13.85 +
   13.86 +val str = Scan.this_string
   13.87 +
   13.88 +val number = Scan.repeat1 (Scan.one Symbol.is_ascii_digit >>
   13.89 +  (fn s => ord s - ord "0")) >>
   13.90 +  foldl1 (fn (n, d) => n * 10 + d)
   13.91 +
   13.92 +val nat = number
   13.93 +val int = Scan.optional (str "~" >> K ~1) 1 -- nat >> op *;
   13.94 +val rat = int --| str "/" -- int >> Rat.rat_of_quotient
   13.95 +val rat_int = rat || int >> Rat.rat_of_int
   13.96 +
   13.97 +(* polynomial parser *)
   13.98 +
   13.99 +fun repeat_sep s f = f ::: Scan.repeat (str s |-- f)
  13.100 +
  13.101 +val parse_id = Scan.one Symbol.is_letter ::: Scan.many Symbol.is_letdig >> implode
  13.102 +
  13.103 +fun parse_varpow ctxt = parse_id -- Scan.optional (str "^" |-- nat) 1 >>
  13.104 +  (fn (x, k) => (cterm_of (ProofContext.theory_of ctxt) (Free (x, @{typ real})), k)) 
  13.105 +
  13.106 +fun parse_monomial ctxt = repeat_sep "*" (parse_varpow ctxt) >>
  13.107 +  (fn xs => fold FuncUtil.Ctermfunc.update xs FuncUtil.Ctermfunc.empty)
  13.108 +
  13.109 +fun parse_cmonomial ctxt =
  13.110 +  rat_int --| str "*" -- (parse_monomial ctxt) >> swap ||
  13.111 +  (parse_monomial ctxt) >> (fn m => (m, Rat.one)) ||
  13.112 +  rat_int >> (fn r => (FuncUtil.Ctermfunc.empty, r))
  13.113 +
  13.114 +fun parse_poly ctxt = repeat_sep "+" (parse_cmonomial ctxt) >>
  13.115 +  (fn xs => fold FuncUtil.Monomialfunc.update xs FuncUtil.Monomialfunc.empty)
  13.116 +
  13.117 +(* positivstellensatz parser *)
  13.118 +
  13.119 +val parse_axiom =
  13.120 +  (str "A=" |-- int >> RealArith.Axiom_eq) ||
  13.121 +  (str "A<=" |-- int >> RealArith.Axiom_le) ||
  13.122 +  (str "A<" |-- int >> RealArith.Axiom_lt)
  13.123 +
  13.124 +val parse_rational =
  13.125 +  (str "R=" |-- rat_int >> RealArith.Rational_eq) ||
  13.126 +  (str "R<=" |-- rat_int >> RealArith.Rational_le) ||
  13.127 +  (str "R<" |-- rat_int >> RealArith.Rational_lt)
  13.128 +
  13.129 +fun parse_cert ctxt input =
  13.130 +  let
  13.131 +    val pc = parse_cert ctxt
  13.132 +    val pp = parse_poly ctxt
  13.133 +  in
  13.134 +  (parse_axiom ||
  13.135 +   parse_rational ||
  13.136 +   str "[" |-- pp --| str "]^2" >> RealArith.Square ||
  13.137 +   str "([" |-- pp --| str "]*" -- pc --| str ")" >> RealArith.Eqmul ||
  13.138 +   str "(" |-- pc --| str "*" -- pc --| str ")" >> RealArith.Product ||
  13.139 +   str "(" |-- pc --| str "+" -- pc --| str ")" >> RealArith.Sum) input
  13.140 +  end
  13.141 +
  13.142 +fun parse_cert_tree ctxt input =
  13.143 +  let
  13.144 +    val pc = parse_cert ctxt
  13.145 +    val pt = parse_cert_tree ctxt
  13.146 +  in
  13.147 +  (str "()" >> K RealArith.Trivial ||
  13.148 +   str "(" |-- pc --| str ")" >> RealArith.Cert ||
  13.149 +   str "(" |-- pt --| str "&" -- pt --| str ")" >> RealArith.Branch) input
  13.150 +  end
  13.151 +
  13.152 +(* scanner *)
  13.153 +
  13.154 +fun cert_to_pss_tree ctxt input_str = Symbol.scanner "bad certificate" (parse_cert_tree ctxt)
  13.155 +  (filter_out Symbol.is_blank (Symbol.explode input_str))
  13.156 +
  13.157 +end
  13.158 +
  13.159 +
    14.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
    14.2 +++ b/src/HOL/Library/Sum_of_Squares/sos_wrapper.ML	Sat Jan 08 17:39:51 2011 +0100
    14.3 @@ -0,0 +1,159 @@
    14.4 +(*  Title:      HOL/Library/Sum_of_Squares/sos_wrapper.ML
    14.5 +    Author:     Philipp Meyer, TU Muenchen
    14.6 +
    14.7 +Added functionality for sums of squares, e.g. calling a remote prover.
    14.8 +*)
    14.9 +
   14.10 +signature SOS_WRAPPER =
   14.11 +sig
   14.12 +  datatype prover_result = Success | Failure | Error
   14.13 +  val setup: theory -> theory
   14.14 +  val dest_dir: string Config.T
   14.15 +  val prover_name: string Config.T
   14.16 +end
   14.17 +
   14.18 +structure SOS_Wrapper: SOS_WRAPPER =
   14.19 +struct
   14.20 +
   14.21 +datatype prover_result = Success | Failure | Error
   14.22 +
   14.23 +fun str_of_result Success = "Success"
   14.24 +  | str_of_result Failure = "Failure"
   14.25 +  | str_of_result Error = "Error"
   14.26 +
   14.27 +
   14.28 +(*** calling provers ***)
   14.29 +
   14.30 +val (dest_dir, setup_dest_dir) = Attrib.config_string "sos_dest_dir" (K "")
   14.31 +
   14.32 +fun filename dir name =
   14.33 +  let
   14.34 +    val probfile = Path.basic (name ^ serial_string ())
   14.35 +    val dir_path = Path.explode dir
   14.36 +  in
   14.37 +    if dir = "" then
   14.38 +      File.tmp_path probfile
   14.39 +    else if File.exists dir_path then
   14.40 +      Path.append dir_path probfile
   14.41 +    else error ("No such directory: " ^ dir)
   14.42 +  end
   14.43 +
   14.44 +fun run_solver ctxt name cmd find_failure input =
   14.45 +  let
   14.46 +    val _ = warning ("Calling solver: " ^ name)
   14.47 +
   14.48 +    (* create input file *)
   14.49 +    val dir = Config.get ctxt dest_dir
   14.50 +    val input_file = filename dir "sos_in"
   14.51 +    val _ = File.write input_file input
   14.52 +
   14.53 +    (* call solver *)
   14.54 +    val output_file = filename dir "sos_out"
   14.55 +    val (output, rv) =
   14.56 +      bash_output
   14.57 +       (if File.exists cmd then
   14.58 +          space_implode " "
   14.59 +            [File.shell_path cmd, File.shell_path input_file, File.shell_path output_file]
   14.60 +        else error ("Bad executable: " ^ File.platform_path cmd))
   14.61 +
   14.62 +    (* read and analyze output *)
   14.63 +    val (res, res_msg) = find_failure rv
   14.64 +    val result = if File.exists output_file then File.read output_file else ""
   14.65 +
   14.66 +    (* remove temporary files *)
   14.67 +    val _ =
   14.68 +      if dir = "" then
   14.69 +        (File.rm input_file; if File.exists output_file then File.rm output_file else ())
   14.70 +      else ()
   14.71 +
   14.72 +    val _ =
   14.73 +      if Config.get ctxt Sum_of_Squares.trace
   14.74 +      then writeln ("Solver output:\n" ^ output)
   14.75 +      else ()
   14.76 +
   14.77 +    val _ = warning (str_of_result res ^ ": " ^ res_msg)
   14.78 +  in
   14.79 +    (case res of
   14.80 +      Success => result
   14.81 +    | Failure => raise Sum_of_Squares.Failure res_msg
   14.82 +    | Error => error ("Prover failed: " ^ res_msg))
   14.83 +  end
   14.84 +
   14.85 +
   14.86 +(*** various provers ***)
   14.87 +
   14.88 +(* local csdp client *)
   14.89 +
   14.90 +fun find_csdp_failure rv =
   14.91 +  case rv of
   14.92 +    0 => (Success, "SDP solved")
   14.93 +  | 1 => (Failure, "SDP is primal infeasible")
   14.94 +  | 2 => (Failure, "SDP is dual infeasible")
   14.95 +  | 3 => (Success, "SDP solved with reduced accuracy")
   14.96 +  | 4 => (Failure, "Maximum iterations reached")
   14.97 +  | 5 => (Failure, "Stuck at edge of primal feasibility")
   14.98 +  | 6 => (Failure, "Stuck at edge of dual infeasibility")
   14.99 +  | 7 => (Failure, "Lack of progress")
  14.100 +  | 8 => (Failure, "X, Z, or O was singular")
  14.101 +  | 9 => (Failure, "Detected NaN or Inf values")
  14.102 +  | _ => (Error, "return code is " ^ string_of_int rv)
  14.103 +
  14.104 +val csdp = ("$CSDP_EXE", find_csdp_failure)
  14.105 +
  14.106 +
  14.107 +(* remote neos server *)
  14.108 +
  14.109 +fun find_neos_failure rv =
  14.110 +  case rv of
  14.111 +    20 => (Error, "error submitting job")
  14.112 +  | 21 => (Error, "interrupt")
  14.113 +  |  _ => find_csdp_failure rv
  14.114 +
  14.115 +val neos_csdp = ("$ISABELLE_SUM_OF_SQUARES/neos_csdp_client", find_neos_failure)
  14.116 +
  14.117 +
  14.118 +(* named provers *)
  14.119 +
  14.120 +fun prover "remote_csdp" = neos_csdp
  14.121 +  | prover "csdp" = csdp
  14.122 +  | prover name = error ("Unknown prover: " ^ name)
  14.123 +
  14.124 +val (prover_name, setup_prover_name) = Attrib.config_string "sos_prover_name" (K "remote_csdp")
  14.125 +
  14.126 +fun call_solver ctxt opt_name =
  14.127 +  let
  14.128 +    val name = the_default (Config.get ctxt prover_name) opt_name
  14.129 +    val (cmd, find_failure) = prover name
  14.130 +  in run_solver ctxt name (Path.explode cmd) find_failure end
  14.131 +
  14.132 +
  14.133 +(* certificate output *)
  14.134 +
  14.135 +fun output_line cert =
  14.136 +  "To repeat this proof with a certifiate use this command:\n" ^
  14.137 +    Markup.markup Markup.sendback ("by (sos_cert \"" ^ cert ^ "\")")
  14.138 +
  14.139 +val print_cert = warning o output_line o PositivstellensatzTools.pss_tree_to_cert
  14.140 +
  14.141 +
  14.142 +(* method setup *)
  14.143 +
  14.144 +fun sos_solver print method = SIMPLE_METHOD' o Sum_of_Squares.sos_tac print method
  14.145 +
  14.146 +val setup =
  14.147 +  setup_dest_dir #>
  14.148 +  setup_prover_name #>
  14.149 +  Method.setup @{binding sos}
  14.150 +    (Scan.lift (Scan.option Parse.xname)
  14.151 +      >> (fn opt_name => fn ctxt =>
  14.152 +        sos_solver print_cert
  14.153 +          (Sum_of_Squares.Prover (call_solver ctxt opt_name)) ctxt))
  14.154 +    "prove universal problems over the reals using sums of squares" #>
  14.155 +  Method.setup @{binding sos_cert}
  14.156 +    (Scan.lift Parse.string
  14.157 +      >> (fn cert => fn ctxt =>
  14.158 +        sos_solver ignore
  14.159 +          (Sum_of_Squares.Certificate (PositivstellensatzTools.cert_to_pss_tree ctxt cert)) ctxt))
  14.160 +    "prove universal problems over the reals using sums of squares with certificates"
  14.161 +
  14.162 +end
    15.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
    15.2 +++ b/src/HOL/Library/Sum_of_Squares/sum_of_squares.ML	Sat Jan 08 17:39:51 2011 +0100
    15.3 @@ -0,0 +1,1435 @@
    15.4 +(*  Title:      HOL/Library/Sum_of_Squares/sum_of_squares.ML
    15.5 +    Author:     Amine Chaieb, University of Cambridge
    15.6 +    Author:     Philipp Meyer, TU Muenchen
    15.7 +
    15.8 +A tactic for proving nonlinear inequalities.
    15.9 +*)
   15.10 +
   15.11 +signature SUM_OF_SQUARES =
   15.12 +sig
   15.13 +  datatype proof_method = Certificate of RealArith.pss_tree | Prover of string -> string
   15.14 +  val sos_tac: (RealArith.pss_tree -> unit) -> proof_method -> Proof.context -> int -> tactic
   15.15 +  val trace: bool Config.T
   15.16 +  val setup: theory -> theory
   15.17 +  exception Failure of string;
   15.18 +end
   15.19 +
   15.20 +structure Sum_of_Squares: SUM_OF_SQUARES =
   15.21 +struct
   15.22 +
   15.23 +val rat_0 = Rat.zero;
   15.24 +val rat_1 = Rat.one;
   15.25 +val rat_2 = Rat.two;
   15.26 +val rat_10 = Rat.rat_of_int 10;
   15.27 +val rat_1_2 = rat_1 // rat_2;
   15.28 +val max = Integer.max;
   15.29 +val min = Integer.min;
   15.30 +
   15.31 +val denominator_rat = Rat.quotient_of_rat #> snd #> Rat.rat_of_int;
   15.32 +val numerator_rat = Rat.quotient_of_rat #> fst #> Rat.rat_of_int;
   15.33 +fun int_of_rat a =
   15.34 +    case Rat.quotient_of_rat a of (i,1) => i | _ => error "int_of_rat: not an int";
   15.35 +fun lcm_rat x y = Rat.rat_of_int (Integer.lcm (int_of_rat x) (int_of_rat y));
   15.36 +
   15.37 +fun rat_pow r i =
   15.38 + let fun pow r i =
   15.39 +   if i = 0 then rat_1 else
   15.40 +   let val d = pow r (i div 2)
   15.41 +   in d */ d */ (if i mod 2 = 0 then rat_1 else r)
   15.42 +   end
   15.43 + in if i < 0 then pow (Rat.inv r) (~ i) else pow r i end;
   15.44 +
   15.45 +fun round_rat r =
   15.46 + let val (a,b) = Rat.quotient_of_rat (Rat.abs r)
   15.47 +     val d = a div b
   15.48 +     val s = if r </ rat_0 then (Rat.neg o Rat.rat_of_int) else Rat.rat_of_int
   15.49 +     val x2 = 2 * (a - (b * d))
   15.50 + in s (if x2 >= b then d + 1 else d) end
   15.51 +
   15.52 +val abs_rat = Rat.abs;
   15.53 +val pow2 = rat_pow rat_2;
   15.54 +val pow10 = rat_pow rat_10;
   15.55 +
   15.56 +val (trace, setup_trace) = Attrib.config_bool "sos_trace" (K false);
   15.57 +val setup = setup_trace;
   15.58 +
   15.59 +exception Sanity;
   15.60 +
   15.61 +exception Unsolvable;
   15.62 +
   15.63 +exception Failure of string;
   15.64 +
   15.65 +datatype proof_method =
   15.66 +    Certificate of RealArith.pss_tree
   15.67 +  | Prover of (string -> string)
   15.68 +
   15.69 +(* Turn a rational into a decimal string with d sig digits.                  *)
   15.70 +
   15.71 +local
   15.72 +fun normalize y =
   15.73 +  if abs_rat y </ (rat_1 // rat_10) then normalize (rat_10 */ y) - 1
   15.74 +  else if abs_rat y >=/ rat_1 then normalize (y // rat_10) + 1
   15.75 +  else 0
   15.76 + in
   15.77 +fun decimalize d x =
   15.78 +  if x =/ rat_0 then "0.0" else
   15.79 +  let
   15.80 +   val y = Rat.abs x
   15.81 +   val e = normalize y
   15.82 +   val z = pow10(~ e) */ y +/ rat_1
   15.83 +   val k = int_of_rat (round_rat(pow10 d */ z))
   15.84 +  in (if x </ rat_0 then "-0." else "0.") ^
   15.85 +     implode(tl(raw_explode(string_of_int k))) ^
   15.86 +     (if e = 0 then "" else "e"^string_of_int e)
   15.87 +  end
   15.88 +end;
   15.89 +
   15.90 +(* Iterations over numbers, and lists indexed by numbers.                    *)
   15.91 +
   15.92 +fun itern k l f a =
   15.93 +  case l of
   15.94 +    [] => a
   15.95 +  | h::t => itern (k + 1) t f (f h k a);
   15.96 +
   15.97 +fun iter (m,n) f a =
   15.98 +  if n < m then a
   15.99 +  else iter (m+1,n) f (f m a);
  15.100 +
  15.101 +(* The main types.                                                           *)
  15.102 +
  15.103 +type vector = int* Rat.rat FuncUtil.Intfunc.table;
  15.104 +
  15.105 +type matrix = (int*int)*(Rat.rat FuncUtil.Intpairfunc.table);
  15.106 +
  15.107 +fun iszero (k,r) = r =/ rat_0;
  15.108 +
  15.109 +
  15.110 +(* Vectors. Conventionally indexed 1..n.                                     *)
  15.111 +
  15.112 +fun vector_0 n = (n,FuncUtil.Intfunc.empty):vector;
  15.113 +
  15.114 +fun dim (v:vector) = fst v;
  15.115 +
  15.116 +fun vector_const c n =
  15.117 +  if c =/ rat_0 then vector_0 n
  15.118 +  else (n,fold_rev (fn k => FuncUtil.Intfunc.update (k,c)) (1 upto n) FuncUtil.Intfunc.empty) :vector;
  15.119 +
  15.120 +val vector_1 = vector_const rat_1;
  15.121 +
  15.122 +fun vector_cmul c (v:vector) =
  15.123 + let val n = dim v
  15.124 + in if c =/ rat_0 then vector_0 n
  15.125 +    else (n,FuncUtil.Intfunc.map (fn _ => fn x => c */ x) (snd v))
  15.126 + end;
  15.127 +
  15.128 +fun vector_neg (v:vector) = (fst v,FuncUtil.Intfunc.map (K Rat.neg) (snd v)) :vector;
  15.129 +
  15.130 +fun vector_add (v1:vector) (v2:vector) =
  15.131 + let val m = dim v1
  15.132 +     val n = dim v2
  15.133 + in if m <> n then error "vector_add: incompatible dimensions"
  15.134 +    else (n,FuncUtil.Intfunc.combine (curry op +/) (fn x => x =/ rat_0) (snd v1) (snd v2)) :vector
  15.135 + end;
  15.136 +
  15.137 +fun vector_sub v1 v2 = vector_add v1 (vector_neg v2);
  15.138 +
  15.139 +fun vector_dot (v1:vector) (v2:vector) =
  15.140 + let val m = dim v1
  15.141 +     val n = dim v2
  15.142 + in if m <> n then error "vector_dot: incompatible dimensions"
  15.143 +    else FuncUtil.Intfunc.fold (fn (i,x) => fn a => x +/ a)
  15.144 +        (FuncUtil.Intfunc.combine (curry op */) (fn x => x =/ rat_0) (snd v1) (snd v2)) rat_0
  15.145 + end;
  15.146 +
  15.147 +fun vector_of_list l =
  15.148 + let val n = length l
  15.149 + in (n,fold_rev2 (curry FuncUtil.Intfunc.update) (1 upto n) l FuncUtil.Intfunc.empty) :vector
  15.150 + end;
  15.151 +
  15.152 +(* Matrices; again rows and columns indexed from 1.                          *)
  15.153 +
  15.154 +fun matrix_0 (m,n) = ((m,n),FuncUtil.Intpairfunc.empty):matrix;
  15.155 +
  15.156 +fun dimensions (m:matrix) = fst m;
  15.157 +
  15.158 +fun matrix_const c (mn as (m,n)) =
  15.159 +  if m <> n then error "matrix_const: needs to be square"
  15.160 +  else if c =/ rat_0 then matrix_0 mn
  15.161 +  else (mn,fold_rev (fn k => FuncUtil.Intpairfunc.update ((k,k), c)) (1 upto n) FuncUtil.Intpairfunc.empty) :matrix;;
  15.162 +
  15.163 +val matrix_1 = matrix_const rat_1;
  15.164 +
  15.165 +fun matrix_cmul c (m:matrix) =
  15.166 + let val (i,j) = dimensions m
  15.167 + in if c =/ rat_0 then matrix_0 (i,j)
  15.168 +    else ((i,j),FuncUtil.Intpairfunc.map (fn _ => fn x => c */ x) (snd m))
  15.169 + end;
  15.170 +
  15.171 +fun matrix_neg (m:matrix) =
  15.172 +  (dimensions m, FuncUtil.Intpairfunc.map (K Rat.neg) (snd m)) :matrix;
  15.173 +
  15.174 +fun matrix_add (m1:matrix) (m2:matrix) =
  15.175 + let val d1 = dimensions m1
  15.176 +     val d2 = dimensions m2
  15.177 + in if d1 <> d2
  15.178 +     then error "matrix_add: incompatible dimensions"
  15.179 +    else (d1,FuncUtil.Intpairfunc.combine (curry op +/) (fn x => x =/ rat_0) (snd m1) (snd m2)) :matrix
  15.180 + end;;
  15.181 +
  15.182 +fun matrix_sub m1 m2 = matrix_add m1 (matrix_neg m2);
  15.183 +
  15.184 +fun row k (m:matrix) =
  15.185 + let val (i,j) = dimensions m
  15.186 + in (j,
  15.187 +   FuncUtil.Intpairfunc.fold (fn ((i,j), c) => fn a => if i = k then FuncUtil.Intfunc.update (j,c) a else a) (snd m) FuncUtil.Intfunc.empty ) : vector
  15.188 + end;
  15.189 +
  15.190 +fun column k (m:matrix) =
  15.191 +  let val (i,j) = dimensions m
  15.192 +  in (i,
  15.193 +   FuncUtil.Intpairfunc.fold (fn ((i,j), c) => fn a => if j = k then FuncUtil.Intfunc.update (i,c) a else a) (snd m)  FuncUtil.Intfunc.empty)
  15.194 +   : vector
  15.195 + end;
  15.196 +
  15.197 +fun transp (m:matrix) =
  15.198 +  let val (i,j) = dimensions m
  15.199 +  in
  15.200 +  ((j,i),FuncUtil.Intpairfunc.fold (fn ((i,j), c) => fn a => FuncUtil.Intpairfunc.update ((j,i), c) a) (snd m) FuncUtil.Intpairfunc.empty) :matrix
  15.201 + end;
  15.202 +
  15.203 +fun diagonal (v:vector) =
  15.204 + let val n = dim v
  15.205 + in ((n,n),FuncUtil.Intfunc.fold (fn (i, c) => fn a => FuncUtil.Intpairfunc.update ((i,i), c) a) (snd v) FuncUtil.Intpairfunc.empty) : matrix
  15.206 + end;
  15.207 +
  15.208 +fun matrix_of_list l =
  15.209 + let val m = length l
  15.210 + in if m = 0 then matrix_0 (0,0) else
  15.211 +   let val n = length (hd l)
  15.212 +   in ((m,n),itern 1 l (fn v => fn i => itern 1 v (fn c => fn j => FuncUtil.Intpairfunc.update ((i,j), c))) FuncUtil.Intpairfunc.empty)
  15.213 +   end
  15.214 + end;
  15.215 +
  15.216 +(* Monomials.                                                                *)
  15.217 +
  15.218 +fun monomial_eval assig m =
  15.219 +  FuncUtil.Ctermfunc.fold (fn (x, k) => fn a => a */ rat_pow (FuncUtil.Ctermfunc.apply assig x) k)
  15.220 +        m rat_1;
  15.221 +val monomial_1 = FuncUtil.Ctermfunc.empty;
  15.222 +
  15.223 +fun monomial_var x = FuncUtil.Ctermfunc.onefunc (x, 1);
  15.224 +
  15.225 +val monomial_mul =
  15.226 +  FuncUtil.Ctermfunc.combine Integer.add (K false);
  15.227 +
  15.228 +fun monomial_pow m k =
  15.229 +  if k = 0 then monomial_1
  15.230 +  else FuncUtil.Ctermfunc.map (fn _ => fn x => k * x) m;
  15.231 +
  15.232 +fun monomial_divides m1 m2 =
  15.233 +  FuncUtil.Ctermfunc.fold (fn (x, k) => fn a => FuncUtil.Ctermfunc.tryapplyd m2 x 0 >= k andalso a) m1 true;;
  15.234 +
  15.235 +fun monomial_div m1 m2 =
  15.236 + let val m = FuncUtil.Ctermfunc.combine Integer.add
  15.237 +   (fn x => x = 0) m1 (FuncUtil.Ctermfunc.map (fn _ => fn x => ~ x) m2)
  15.238 + in if FuncUtil.Ctermfunc.fold (fn (x, k) => fn a => k >= 0 andalso a) m true then m
  15.239 +  else error "monomial_div: non-divisible"
  15.240 + end;
  15.241 +
  15.242 +fun monomial_degree x m =
  15.243 +  FuncUtil.Ctermfunc.tryapplyd m x 0;;
  15.244 +
  15.245 +fun monomial_lcm m1 m2 =
  15.246 +  fold_rev (fn x => FuncUtil.Ctermfunc.update (x, max (monomial_degree x m1) (monomial_degree x m2)))
  15.247 +          (union (is_equal o FuncUtil.cterm_ord) (FuncUtil.Ctermfunc.dom m1) (FuncUtil.Ctermfunc.dom m2)) (FuncUtil.Ctermfunc.empty);
  15.248 +
  15.249 +fun monomial_multidegree m =
  15.250 + FuncUtil.Ctermfunc.fold (fn (x, k) => fn a => k + a) m 0;;
  15.251 +
  15.252 +fun monomial_variables m = FuncUtil.Ctermfunc.dom m;;
  15.253 +
  15.254 +(* Polynomials.                                                              *)
  15.255 +
  15.256 +fun eval assig p =
  15.257 +  FuncUtil.Monomialfunc.fold (fn (m, c) => fn a => a +/ c */ monomial_eval assig m) p rat_0;
  15.258 +
  15.259 +val poly_0 = FuncUtil.Monomialfunc.empty;
  15.260 +
  15.261 +fun poly_isconst p =
  15.262 +  FuncUtil.Monomialfunc.fold (fn (m, c) => fn a => FuncUtil.Ctermfunc.is_empty m andalso a) p true;
  15.263 +
  15.264 +fun poly_var x = FuncUtil.Monomialfunc.onefunc (monomial_var x,rat_1);
  15.265 +
  15.266 +fun poly_const c =
  15.267 +  if c =/ rat_0 then poly_0 else FuncUtil.Monomialfunc.onefunc(monomial_1, c);
  15.268 +
  15.269 +fun poly_cmul c p =
  15.270 +  if c =/ rat_0 then poly_0
  15.271 +  else FuncUtil.Monomialfunc.map (fn _ => fn x => c */ x) p;
  15.272 +
  15.273 +fun poly_neg p = FuncUtil.Monomialfunc.map (K Rat.neg) p;;
  15.274 +
  15.275 +fun poly_add p1 p2 =
  15.276 +  FuncUtil.Monomialfunc.combine (curry op +/) (fn x => x =/ rat_0) p1 p2;
  15.277 +
  15.278 +fun poly_sub p1 p2 = poly_add p1 (poly_neg p2);
  15.279 +
  15.280 +fun poly_cmmul (c,m) p =
  15.281 + if c =/ rat_0 then poly_0
  15.282 + else if FuncUtil.Ctermfunc.is_empty m
  15.283 +      then FuncUtil.Monomialfunc.map (fn _ => fn d => c */ d) p
  15.284 +      else FuncUtil.Monomialfunc.fold (fn (m', d) => fn a => (FuncUtil.Monomialfunc.update (monomial_mul m m', c */ d) a)) p poly_0;
  15.285 +
  15.286 +fun poly_mul p1 p2 =
  15.287 +  FuncUtil.Monomialfunc.fold (fn (m, c) => fn a => poly_add (poly_cmmul (c,m) p2) a) p1 poly_0;
  15.288 +
  15.289 +fun poly_div p1 p2 =
  15.290 + if not(poly_isconst p2)
  15.291 + then error "poly_div: non-constant" else
  15.292 + let val c = eval FuncUtil.Ctermfunc.empty p2
  15.293 + in if c =/ rat_0 then error "poly_div: division by zero"
  15.294 +    else poly_cmul (Rat.inv c) p1
  15.295 + end;
  15.296 +
  15.297 +fun poly_square p = poly_mul p p;
  15.298 +
  15.299 +fun poly_pow p k =
  15.300 + if k = 0 then poly_const rat_1
  15.301 + else if k = 1 then p
  15.302 + else let val q = poly_square(poly_pow p (k div 2)) in
  15.303 +      if k mod 2 = 1 then poly_mul p q else q end;
  15.304 +
  15.305 +fun poly_exp p1 p2 =
  15.306 +  if not(poly_isconst p2)
  15.307 +  then error "poly_exp: not a constant"
  15.308 +  else poly_pow p1 (int_of_rat (eval FuncUtil.Ctermfunc.empty p2));
  15.309 +
  15.310 +fun degree x p =
  15.311 + FuncUtil.Monomialfunc.fold (fn (m,c) => fn a => max (monomial_degree x m) a) p 0;
  15.312 +
  15.313 +fun multidegree p =
  15.314 +  FuncUtil.Monomialfunc.fold (fn (m, c) => fn a => max (monomial_multidegree m) a) p 0;
  15.315 +
  15.316 +fun poly_variables p =
  15.317 +  sort FuncUtil.cterm_ord (FuncUtil.Monomialfunc.fold_rev (fn (m, c) => union (is_equal o FuncUtil.cterm_ord) (monomial_variables m)) p []);;
  15.318 +
  15.319 +(* Order monomials for human presentation.                                   *)
  15.320 +
  15.321 +val humanorder_varpow = prod_ord FuncUtil.cterm_ord (rev_order o int_ord);
  15.322 +
  15.323 +local
  15.324 + fun ord (l1,l2) = case (l1,l2) of
  15.325 +  (_,[]) => LESS
  15.326 + | ([],_) => GREATER
  15.327 + | (h1::t1,h2::t2) =>
  15.328 +   (case humanorder_varpow (h1, h2) of
  15.329 +     LESS => LESS
  15.330 +   | EQUAL => ord (t1,t2)
  15.331 +   | GREATER => GREATER)
  15.332 +in fun humanorder_monomial m1 m2 =
  15.333 + ord (sort humanorder_varpow (FuncUtil.Ctermfunc.dest m1),
  15.334 +  sort humanorder_varpow (FuncUtil.Ctermfunc.dest m2))
  15.335 +end;
  15.336 +
  15.337 +(* Conversions to strings.                                                   *)
  15.338 +
  15.339 +fun string_of_vector min_size max_size (v:vector) =
  15.340 + let val n_raw = dim v
  15.341 + in if n_raw = 0 then "[]" else
  15.342 +  let
  15.343 +   val n = max min_size (min n_raw max_size)
  15.344 +   val xs = map (Rat.string_of_rat o (fn i => FuncUtil.Intfunc.tryapplyd (snd v) i rat_0)) (1 upto n)
  15.345 +  in "[" ^ space_implode ", " xs ^
  15.346 +  (if n_raw > max_size then ", ...]" else "]")
  15.347 +  end
  15.348 + end;
  15.349 +
  15.350 +fun string_of_matrix max_size (m:matrix) =
  15.351 + let
  15.352 +  val (i_raw,j_raw) = dimensions m
  15.353 +  val i = min max_size i_raw
  15.354 +  val j = min max_size j_raw
  15.355 +  val rstr = map (fn k => string_of_vector j j (row k m)) (1 upto i)
  15.356 + in "["^ space_implode ";\n " rstr ^
  15.357 +  (if j > max_size then "\n ...]" else "]")
  15.358 + end;
  15.359 +
  15.360 +fun string_of_term t =
  15.361 + case t of
  15.362 +   a$b => "("^(string_of_term a)^" "^(string_of_term b)^")"
  15.363 + | Abs x =>
  15.364 +    let val (xn, b) = Term.dest_abs x
  15.365 +    in "(\\"^xn^"."^(string_of_term b)^")"
  15.366 +    end
  15.367 + | Const(s,_) => s
  15.368 + | Free (s,_) => s
  15.369 + | Var((s,_),_) => s
  15.370 + | _ => error "string_of_term";
  15.371 +
  15.372 +val string_of_cterm = string_of_term o term_of;
  15.373 +
  15.374 +fun string_of_varpow x k =
  15.375 +  if k = 1 then string_of_cterm x
  15.376 +  else string_of_cterm x^"^"^string_of_int k;
  15.377 +
  15.378 +fun string_of_monomial m =
  15.379 + if FuncUtil.Ctermfunc.is_empty m then "1" else
  15.380 + let val vps = fold_rev (fn (x,k) => fn a => string_of_varpow x k :: a)
  15.381 +  (sort humanorder_varpow (FuncUtil.Ctermfunc.dest m)) []
  15.382 + in space_implode "*" vps
  15.383 + end;
  15.384 +
  15.385 +fun string_of_cmonomial (c,m) =
  15.386 + if FuncUtil.Ctermfunc.is_empty m then Rat.string_of_rat c
  15.387 + else if c =/ rat_1 then string_of_monomial m
  15.388 + else Rat.string_of_rat c ^ "*" ^ string_of_monomial m;;
  15.389 +
  15.390 +fun string_of_poly p =
  15.391 + if FuncUtil.Monomialfunc.is_empty p then "<<0>>" else
  15.392 + let
  15.393 +  val cms = sort (fn ((m1,_),(m2,_)) => humanorder_monomial m1  m2) (FuncUtil.Monomialfunc.dest p)
  15.394 +  val s = fold (fn (m,c) => fn a =>
  15.395 +             if c </ rat_0 then a ^ " - " ^ string_of_cmonomial(Rat.neg c,m)
  15.396 +             else a ^ " + " ^ string_of_cmonomial(c,m))
  15.397 +          cms ""
  15.398 +  val s1 = String.substring (s, 0, 3)
  15.399 +  val s2 = String.substring (s, 3, String.size s - 3)
  15.400 + in "<<" ^(if s1 = " + " then s2 else "-"^s2)^">>"
  15.401 + end;
  15.402 +
  15.403 +(* Conversion from HOL term.                                                 *)
  15.404 +
  15.405 +local
  15.406 + val neg_tm = @{cterm "uminus :: real => _"}
  15.407 + val add_tm = @{cterm "op + :: real => _"}
  15.408 + val sub_tm = @{cterm "op - :: real => _"}
  15.409 + val mul_tm = @{cterm "op * :: real => _"}
  15.410 + val inv_tm = @{cterm "inverse :: real => _"}
  15.411 + val div_tm = @{cterm "op / :: real => _"}
  15.412 + val pow_tm = @{cterm "op ^ :: real => _"}
  15.413 + val zero_tm = @{cterm "0:: real"}
  15.414 + val is_numeral = can (HOLogic.dest_number o term_of)
  15.415 + fun is_comb t = case t of _$_ => true | _ => false
  15.416 + fun poly_of_term tm =
  15.417 +  if tm aconvc zero_tm then poly_0
  15.418 +  else if RealArith.is_ratconst tm
  15.419 +       then poly_const(RealArith.dest_ratconst tm)
  15.420 +  else
  15.421 +  (let val (lop,r) = Thm.dest_comb tm
  15.422 +   in if lop aconvc neg_tm then poly_neg(poly_of_term r)
  15.423 +      else if lop aconvc inv_tm then
  15.424 +       let val p = poly_of_term r
  15.425 +       in if poly_isconst p
  15.426 +          then poly_const(Rat.inv (eval FuncUtil.Ctermfunc.empty p))
  15.427 +          else error "poly_of_term: inverse of non-constant polyomial"
  15.428 +       end
  15.429 +   else (let val (opr,l) = Thm.dest_comb lop
  15.430 +         in
  15.431 +          if opr aconvc pow_tm andalso is_numeral r
  15.432 +          then poly_pow (poly_of_term l) ((snd o HOLogic.dest_number o term_of) r)
  15.433 +          else if opr aconvc add_tm
  15.434 +           then poly_add (poly_of_term l) (poly_of_term r)
  15.435 +          else if opr aconvc sub_tm
  15.436 +           then poly_sub (poly_of_term l) (poly_of_term r)
  15.437 +          else if opr aconvc mul_tm
  15.438 +           then poly_mul (poly_of_term l) (poly_of_term r)
  15.439 +          else if opr aconvc div_tm
  15.440 +           then let
  15.441 +                  val p = poly_of_term l
  15.442 +                  val q = poly_of_term r
  15.443 +                in if poly_isconst q then poly_cmul (Rat.inv (eval FuncUtil.Ctermfunc.empty q)) p
  15.444 +                   else error "poly_of_term: division by non-constant polynomial"
  15.445 +                end
  15.446 +          else poly_var tm
  15.447 +
  15.448 +         end
  15.449 +         handle CTERM ("dest_comb",_) => poly_var tm)
  15.450 +   end
  15.451 +   handle CTERM ("dest_comb",_) => poly_var tm)
  15.452 +in
  15.453 +val poly_of_term = fn tm =>
  15.454 + if type_of (term_of tm) = @{typ real} then poly_of_term tm
  15.455 + else error "poly_of_term: term does not have real type"
  15.456 +end;
  15.457 +
  15.458 +(* String of vector (just a list of space-separated numbers).                *)
  15.459 +
  15.460 +fun sdpa_of_vector (v:vector) =
  15.461 + let
  15.462 +  val n = dim v
  15.463 +  val strs = map (decimalize 20 o (fn i => FuncUtil.Intfunc.tryapplyd (snd v) i rat_0)) (1 upto n)
  15.464 + in space_implode " " strs ^ "\n"
  15.465 + end;
  15.466 +
  15.467 +fun triple_int_ord ((a,b,c),(a',b',c')) =
  15.468 + prod_ord int_ord (prod_ord int_ord int_ord)
  15.469 +    ((a,(b,c)),(a',(b',c')));
  15.470 +structure Inttriplefunc = FuncFun(type key = int*int*int val ord = triple_int_ord);
  15.471 +
  15.472 +(* String for block diagonal matrix numbered k.                              *)
  15.473 +
  15.474 +fun sdpa_of_blockdiagonal k m =
  15.475 + let
  15.476 +  val pfx = string_of_int k ^" "
  15.477 +  val ents =
  15.478 +    Inttriplefunc.fold (fn ((b,i,j), c) => fn a => if i > j then a else ((b,i,j),c)::a) m []
  15.479 +  val entss = sort (triple_int_ord o pairself fst) ents
  15.480 + in  fold_rev (fn ((b,i,j),c) => fn a =>
  15.481 +     pfx ^ string_of_int b ^ " " ^ string_of_int i ^ " " ^ string_of_int j ^
  15.482 +     " " ^ decimalize 20 c ^ "\n" ^ a) entss ""
  15.483 + end;
  15.484 +
  15.485 +(* String for a matrix numbered k, in SDPA sparse format.                    *)
  15.486 +
  15.487 +fun sdpa_of_matrix k (m:matrix) =
  15.488 + let
  15.489 +  val pfx = string_of_int k ^ " 1 "
  15.490 +  val ms = FuncUtil.Intpairfunc.fold (fn ((i,j), c) => fn  a => if i > j then a else ((i,j),c)::a) (snd m) []
  15.491 +  val mss = sort ((prod_ord int_ord int_ord) o pairself fst) ms
  15.492 + in fold_rev (fn ((i,j),c) => fn a =>
  15.493 +     pfx ^ string_of_int i ^ " " ^ string_of_int j ^
  15.494 +     " " ^ decimalize 20 c ^ "\n" ^ a) mss ""
  15.495 + end;;
  15.496 +
  15.497 +(* ------------------------------------------------------------------------- *)
  15.498 +(* String in SDPA sparse format for standard SDP problem:                    *)
  15.499 +(*                                                                           *)
  15.500 +(*    X = v_1 * [M_1] + ... + v_m * [M_m] - [M_0] must be PSD                *)
  15.501 +(*    Minimize obj_1 * v_1 + ... obj_m * v_m                                 *)
  15.502 +(* ------------------------------------------------------------------------- *)
  15.503 +
  15.504 +fun sdpa_of_problem obj mats =
  15.505 + let
  15.506 +  val m = length mats - 1
  15.507 +  val (n,_) = dimensions (hd mats)
  15.508 + in
  15.509 +  string_of_int m ^ "\n" ^
  15.510 +  "1\n" ^
  15.511 +  string_of_int n ^ "\n" ^
  15.512 +  sdpa_of_vector obj ^
  15.513 +  fold_rev2 (fn k => fn m => fn a => sdpa_of_matrix (k - 1) m ^ a) (1 upto length mats) mats ""
  15.514 + end;
  15.515 +
  15.516 +fun index_char str chr pos =
  15.517 +  if pos >= String.size str then ~1
  15.518 +  else if String.sub(str,pos) = chr then pos
  15.519 +  else index_char str chr (pos + 1);
  15.520 +fun rat_of_quotient (a,b) = if b = 0 then rat_0 else Rat.rat_of_quotient (a,b);
  15.521 +fun rat_of_string s =
  15.522 + let val n = index_char s #"/" 0 in
  15.523 +  if n = ~1 then s |> Int.fromString |> the |> Rat.rat_of_int
  15.524 +  else
  15.525 +   let val SOME numer = Int.fromString(String.substring(s,0,n))
  15.526 +       val SOME den = Int.fromString (String.substring(s,n+1,String.size s - n - 1))
  15.527 +   in rat_of_quotient(numer, den)
  15.528 +   end
  15.529 + end;
  15.530 +
  15.531 +fun isspace x = (x = " ");
  15.532 +fun isnum x = member (op =) ["0","1","2","3","4","5","6","7","8","9"] x;
  15.533 +
  15.534 +(* More parser basics.                                                       *)
  15.535 +
  15.536 + val word = Scan.this_string
  15.537 + fun token s =
  15.538 +  Scan.repeat ($$ " ") |-- word s --| Scan.repeat ($$ " ")
  15.539 + val numeral = Scan.one isnum
  15.540 + val decimalint = Scan.repeat1 numeral >> (rat_of_string o implode)
  15.541 + val decimalfrac = Scan.repeat1 numeral
  15.542 +    >> (fn s => rat_of_string(implode s) // pow10 (length s))
  15.543 + val decimalsig =
  15.544 +    decimalint -- Scan.option (Scan.$$ "." |-- decimalfrac)
  15.545 +    >> (fn (h,NONE) => h | (h,SOME x) => h +/ x)
  15.546 + fun signed prs =
  15.547 +       $$ "-" |-- prs >> Rat.neg
  15.548 +    || $$ "+" |-- prs
  15.549 +    || prs;
  15.550 +
  15.551 +fun emptyin def xs = if null xs then (def,xs) else Scan.fail xs
  15.552 +
  15.553 + val exponent = ($$ "e" || $$ "E") |-- signed decimalint;
  15.554 +
  15.555 + val decimal = signed decimalsig -- (emptyin rat_0|| exponent)
  15.556 +    >> (fn (h, x) => h */ pow10 (int_of_rat x));
  15.557 +
  15.558 + fun mkparser p s =
  15.559 +  let val (x,rst) = p (raw_explode s)
  15.560 +  in if null rst then x
  15.561 +     else error "mkparser: unparsed input"
  15.562 +  end;;
  15.563 +
  15.564 +(* Parse back csdp output.                                                      *)
  15.565 +
  15.566 + fun ignore inp = ((),[])
  15.567 + fun csdpoutput inp =
  15.568 +   ((decimal -- Scan.repeat (Scan.$$ " " |-- Scan.option decimal) >>
  15.569 +    (fn (h,to) => map_filter I ((SOME h)::to))) --| ignore >> vector_of_list) inp
  15.570 + val parse_csdpoutput = mkparser csdpoutput
  15.571 +
  15.572 +(* Run prover on a problem in linear form.                       *)
  15.573 +
  15.574 +fun run_problem prover obj mats =
  15.575 +  parse_csdpoutput (prover (sdpa_of_problem obj mats))
  15.576 +
  15.577 +(* Try some apparently sensible scaling first. Note that this is purely to   *)
  15.578 +(* get a cleaner translation to floating-point, and doesn't affect any of    *)
  15.579 +(* the results, in principle. In practice it seems a lot better when there   *)
  15.580 +(* are extreme numbers in the original problem.                              *)
  15.581 +
  15.582 +  (* Version for (int*int) keys *)
  15.583 +local
  15.584 +  fun max_rat x y = if x </ y then y else x
  15.585 +  fun common_denominator fld amat acc =
  15.586 +      fld (fn (m,c) => fn a => lcm_rat (denominator_rat c) a) amat acc
  15.587 +  fun maximal_element fld amat acc =
  15.588 +    fld (fn (m,c) => fn maxa => max_rat maxa (abs_rat c)) amat acc
  15.589 +fun float_of_rat x = let val (a,b) = Rat.quotient_of_rat x
  15.590 +                     in Real.fromLargeInt a / Real.fromLargeInt b end;
  15.591 +in
  15.592 +
  15.593 +fun pi_scale_then solver (obj:vector)  mats =
  15.594 + let
  15.595 +  val cd1 = fold_rev (common_denominator FuncUtil.Intpairfunc.fold) mats (rat_1)
  15.596 +  val cd2 = common_denominator FuncUtil.Intfunc.fold (snd obj)  (rat_1)
  15.597 +  val mats' = map (FuncUtil.Intpairfunc.map (fn _ => fn x => cd1 */ x)) mats
  15.598 +  val obj' = vector_cmul cd2 obj
  15.599 +  val max1 = fold_rev (maximal_element FuncUtil.Intpairfunc.fold) mats' (rat_0)
  15.600 +  val max2 = maximal_element FuncUtil.Intfunc.fold (snd obj') (rat_0)
  15.601 +  val scal1 = pow2 (20 - trunc(Math.ln (float_of_rat max1) / Math.ln 2.0))
  15.602 +  val scal2 = pow2 (20 - trunc(Math.ln (float_of_rat max2) / Math.ln 2.0))
  15.603 +  val mats'' = map (FuncUtil.Intpairfunc.map (fn _ => fn x => x */ scal1)) mats'
  15.604 +  val obj'' = vector_cmul scal2 obj'
  15.605 + in solver obj'' mats''
  15.606 +  end
  15.607 +end;
  15.608 +
  15.609 +(* Try some apparently sensible scaling first. Note that this is purely to   *)
  15.610 +(* get a cleaner translation to floating-point, and doesn't affect any of    *)
  15.611 +(* the results, in principle. In practice it seems a lot better when there   *)
  15.612 +(* are extreme numbers in the original problem.                              *)
  15.613 +
  15.614 +  (* Version for (int*int*int) keys *)
  15.615 +local
  15.616 +  fun max_rat x y = if x </ y then y else x
  15.617 +  fun common_denominator fld amat acc =
  15.618 +      fld (fn (m,c) => fn a => lcm_rat (denominator_rat c) a) amat acc
  15.619 +  fun maximal_element fld amat acc =
  15.620 +    fld (fn (m,c) => fn maxa => max_rat maxa (abs_rat c)) amat acc
  15.621 +fun float_of_rat x = let val (a,b) = Rat.quotient_of_rat x
  15.622 +                     in Real.fromLargeInt a / Real.fromLargeInt b end;
  15.623 +fun int_of_float x = (trunc x handle Overflow => 0 | Domain => 0)
  15.624 +in
  15.625 +
  15.626 +fun tri_scale_then solver (obj:vector)  mats =
  15.627 + let
  15.628 +  val cd1 = fold_rev (common_denominator Inttriplefunc.fold) mats (rat_1)
  15.629 +  val cd2 = common_denominator FuncUtil.Intfunc.fold (snd obj)  (rat_1)
  15.630 +  val mats' = map (Inttriplefunc.map (fn _ => fn x => cd1 */ x)) mats
  15.631 +  val obj' = vector_cmul cd2 obj
  15.632 +  val max1 = fold_rev (maximal_element Inttriplefunc.fold) mats' (rat_0)
  15.633 +  val max2 = maximal_element FuncUtil.Intfunc.fold (snd obj') (rat_0)
  15.634 +  val scal1 = pow2 (20 - int_of_float(Math.ln (float_of_rat max1) / Math.ln 2.0))
  15.635 +  val scal2 = pow2 (20 - int_of_float(Math.ln (float_of_rat max2) / Math.ln 2.0))
  15.636 +  val mats'' = map (Inttriplefunc.map (fn _ => fn x => x */ scal1)) mats'
  15.637 +  val obj'' = vector_cmul scal2 obj'
  15.638 + in solver obj'' mats''
  15.639 +  end
  15.640 +end;
  15.641 +
  15.642 +(* Round a vector to "nice" rationals.                                       *)
  15.643 +
  15.644 +fun nice_rational n x = round_rat (n */ x) // n;;
  15.645 +fun nice_vector n ((d,v) : vector) =
  15.646 + (d, FuncUtil.Intfunc.fold (fn (i,c) => fn a =>
  15.647 +   let val y = nice_rational n c
  15.648 +   in if c =/ rat_0 then a
  15.649 +      else FuncUtil.Intfunc.update (i,y) a end) v FuncUtil.Intfunc.empty):vector
  15.650 +
  15.651 +fun dest_ord f x = is_equal (f x);
  15.652 +
  15.653 +(* Stuff for "equations" ((int*int*int)->num functions).                         *)
  15.654 +
  15.655 +fun tri_equation_cmul c eq =
  15.656 +  if c =/ rat_0 then Inttriplefunc.empty else Inttriplefunc.map (fn _ => fn d => c */ d) eq;
  15.657 +
  15.658 +fun tri_equation_add eq1 eq2 = Inttriplefunc.combine (curry op +/) (fn x => x =/ rat_0) eq1 eq2;
  15.659 +
  15.660 +fun tri_equation_eval assig eq =
  15.661 + let fun value v = Inttriplefunc.apply assig v
  15.662 + in Inttriplefunc.fold (fn (v, c) => fn a => a +/ value v */ c) eq rat_0
  15.663 + end;
  15.664 +
  15.665 +(* Eliminate among linear equations: return unconstrained variables and      *)
  15.666 +(* assignments for the others in terms of them. We give one pseudo-variable  *)
  15.667 +(* "one" that's used for a constant term.                                    *)
  15.668 +
  15.669 +local
  15.670 +  fun extract_first p l = case l of  (* FIXME : use find_first instead *)
  15.671 +   [] => error "extract_first"
  15.672 + | h::t => if p h then (h,t) else
  15.673 +          let val (k,s) = extract_first p t in (k,h::s) end
  15.674 +fun eliminate vars dun eqs = case vars of
  15.675 +  [] => if forall Inttriplefunc.is_empty eqs then dun
  15.676 +        else raise Unsolvable
  15.677 + | v::vs =>
  15.678 +  ((let
  15.679 +    val (eq,oeqs) = extract_first (fn e => Inttriplefunc.defined e v) eqs
  15.680 +    val a = Inttriplefunc.apply eq v
  15.681 +    val eq' = tri_equation_cmul ((Rat.neg rat_1) // a) (Inttriplefunc.delete_safe v eq)
  15.682 +    fun elim e =
  15.683 +     let val b = Inttriplefunc.tryapplyd e v rat_0
  15.684 +     in if b =/ rat_0 then e else
  15.685 +        tri_equation_add e (tri_equation_cmul (Rat.neg b // a) eq)
  15.686 +     end
  15.687 +   in eliminate vs (Inttriplefunc.update (v,eq') (Inttriplefunc.map (K elim) dun)) (map elim oeqs)
  15.688 +   end)
  15.689 +  handle Failure _ => eliminate vs dun eqs)
  15.690 +in
  15.691 +fun tri_eliminate_equations one vars eqs =
  15.692 + let
  15.693 +  val assig = eliminate vars Inttriplefunc.empty eqs
  15.694 +  val vs = Inttriplefunc.fold (fn (x, f) => fn a => remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig []
  15.695 +  in (distinct (dest_ord triple_int_ord) vs, assig)
  15.696 +  end
  15.697 +end;
  15.698 +
  15.699 +(* Eliminate all variables, in an essentially arbitrary order.               *)
  15.700 +
  15.701 +fun tri_eliminate_all_equations one =
  15.702 + let
  15.703 +  fun choose_variable eq =
  15.704 +   let val (v,_) = Inttriplefunc.choose eq
  15.705 +   in if is_equal (triple_int_ord(v,one)) then
  15.706 +      let val eq' = Inttriplefunc.delete_safe v eq
  15.707 +      in if Inttriplefunc.is_empty eq' then error "choose_variable"
  15.708 +         else fst (Inttriplefunc.choose eq')
  15.709 +      end
  15.710 +    else v
  15.711 +   end
  15.712 +  fun eliminate dun eqs = case eqs of
  15.713 +    [] => dun
  15.714 +  | eq::oeqs =>
  15.715 +    if Inttriplefunc.is_empty eq then eliminate dun oeqs else
  15.716 +    let val v = choose_variable eq
  15.717 +        val a = Inttriplefunc.apply eq v
  15.718 +        val eq' = tri_equation_cmul ((Rat.rat_of_int ~1) // a)
  15.719 +                   (Inttriplefunc.delete_safe v eq)
  15.720 +        fun elim e =
  15.721 +         let val b = Inttriplefunc.tryapplyd e v rat_0
  15.722 +         in if b =/ rat_0 then e
  15.723 +            else tri_equation_add e (tri_equation_cmul (Rat.neg b // a) eq)
  15.724 +         end
  15.725 +    in eliminate (Inttriplefunc.update(v, eq') (Inttriplefunc.map (K elim) dun))
  15.726 +                 (map elim oeqs)
  15.727 +    end
  15.728 +in fn eqs =>
  15.729 + let
  15.730 +  val assig = eliminate Inttriplefunc.empty eqs
  15.731 +  val vs = Inttriplefunc.fold (fn (x, f) => fn a => remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig []
  15.732 + in (distinct (dest_ord triple_int_ord) vs,assig)
  15.733 + end
  15.734 +end;
  15.735 +
  15.736 +(* Solve equations by assigning arbitrary numbers.                           *)
  15.737 +
  15.738 +fun tri_solve_equations one eqs =
  15.739 + let
  15.740 +  val (vars,assigs) = tri_eliminate_all_equations one eqs
  15.741 +  val vfn = fold_rev (fn v => Inttriplefunc.update(v,rat_0)) vars
  15.742 +            (Inttriplefunc.onefunc(one, Rat.rat_of_int ~1))
  15.743 +  val ass =
  15.744 +    Inttriplefunc.combine (curry op +/) (K false)
  15.745 +    (Inttriplefunc.map (K (tri_equation_eval vfn)) assigs) vfn
  15.746 + in if forall (fn e => tri_equation_eval ass e =/ rat_0) eqs
  15.747 +    then Inttriplefunc.delete_safe one ass else raise Sanity
  15.748 + end;
  15.749 +
  15.750 +(* Multiply equation-parametrized poly by regular poly and add accumulator.  *)
  15.751 +
  15.752 +fun tri_epoly_pmul p q acc =
  15.753 + FuncUtil.Monomialfunc.fold (fn (m1, c) => fn a =>
  15.754 +  FuncUtil.Monomialfunc.fold (fn (m2,e) => fn b =>
  15.755 +   let val m =  monomial_mul m1 m2
  15.756 +       val es = FuncUtil.Monomialfunc.tryapplyd b m Inttriplefunc.empty
  15.757 +   in FuncUtil.Monomialfunc.update (m,tri_equation_add (tri_equation_cmul c e) es) b
  15.758 +   end) q a) p acc ;
  15.759 +
  15.760 +(* Usual operations on equation-parametrized poly.                           *)
  15.761 +
  15.762 +fun tri_epoly_cmul c l =
  15.763 +  if c =/ rat_0 then Inttriplefunc.empty else Inttriplefunc.map (K (tri_equation_cmul c)) l;;
  15.764 +
  15.765 +val tri_epoly_neg = tri_epoly_cmul (Rat.rat_of_int ~1);
  15.766 +
  15.767 +val tri_epoly_add = Inttriplefunc.combine tri_equation_add Inttriplefunc.is_empty;
  15.768 +
  15.769 +fun tri_epoly_sub p q = tri_epoly_add p (tri_epoly_neg q);;
  15.770 +
  15.771 +(* Stuff for "equations" ((int*int)->num functions).                         *)
  15.772 +
  15.773 +fun pi_equation_cmul c eq =
  15.774 +  if c =/ rat_0 then Inttriplefunc.empty else Inttriplefunc.map (fn _ => fn d => c */ d) eq;
  15.775 +
  15.776 +fun pi_equation_add eq1 eq2 = Inttriplefunc.combine (curry op +/) (fn x => x =/ rat_0) eq1 eq2;
  15.777 +
  15.778 +fun pi_equation_eval assig eq =
  15.779 + let fun value v = Inttriplefunc.apply assig v
  15.780 + in Inttriplefunc.fold (fn (v, c) => fn a => a +/ value v */ c) eq rat_0
  15.781 + end;
  15.782 +
  15.783 +(* Eliminate among linear equations: return unconstrained variables and      *)
  15.784 +(* assignments for the others in terms of them. We give one pseudo-variable  *)
  15.785 +(* "one" that's used for a constant term.                                    *)
  15.786 +
  15.787 +local
  15.788 +fun extract_first p l = case l of
  15.789 +   [] => error "extract_first"
  15.790 + | h::t => if p h then (h,t) else
  15.791 +          let val (k,s) = extract_first p t in (k,h::s) end
  15.792 +fun eliminate vars dun eqs = case vars of
  15.793 +  [] => if forall Inttriplefunc.is_empty eqs then dun
  15.794 +        else raise Unsolvable
  15.795 + | v::vs =>
  15.796 +   let
  15.797 +    val (eq,oeqs) = extract_first (fn e => Inttriplefunc.defined e v) eqs
  15.798 +    val a = Inttriplefunc.apply eq v
  15.799 +    val eq' = pi_equation_cmul ((Rat.neg rat_1) // a) (Inttriplefunc.delete_safe v eq)
  15.800 +    fun elim e =
  15.801 +     let val b = Inttriplefunc.tryapplyd e v rat_0
  15.802 +     in if b =/ rat_0 then e else
  15.803 +        pi_equation_add e (pi_equation_cmul (Rat.neg b // a) eq)
  15.804 +     end
  15.805 +   in eliminate vs (Inttriplefunc.update (v,eq') (Inttriplefunc.map (K elim) dun)) (map elim oeqs)
  15.806 +   end
  15.807 +  handle Failure _ => eliminate vs dun eqs
  15.808 +in
  15.809 +fun pi_eliminate_equations one vars eqs =
  15.810 + let
  15.811 +  val assig = eliminate vars Inttriplefunc.empty eqs
  15.812 +  val vs = Inttriplefunc.fold (fn (x, f) => fn a => remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig []
  15.813 +  in (distinct (dest_ord triple_int_ord) vs, assig)
  15.814 +  end
  15.815 +end;
  15.816 +
  15.817 +(* Eliminate all variables, in an essentially arbitrary order.               *)
  15.818 +
  15.819 +fun pi_eliminate_all_equations one =
  15.820 + let
  15.821 +  fun choose_variable eq =
  15.822 +   let val (v,_) = Inttriplefunc.choose eq
  15.823 +   in if is_equal (triple_int_ord(v,one)) then
  15.824 +      let val eq' = Inttriplefunc.delete_safe v eq
  15.825 +      in if Inttriplefunc.is_empty eq' then error "choose_variable"
  15.826 +         else fst (Inttriplefunc.choose eq')
  15.827 +      end
  15.828 +    else v
  15.829 +   end
  15.830 +  fun eliminate dun eqs = case eqs of
  15.831 +    [] => dun
  15.832 +  | eq::oeqs =>
  15.833 +    if Inttriplefunc.is_empty eq then eliminate dun oeqs else
  15.834 +    let val v = choose_variable eq
  15.835 +        val a = Inttriplefunc.apply eq v
  15.836 +        val eq' = pi_equation_cmul ((Rat.rat_of_int ~1) // a)
  15.837 +                   (Inttriplefunc.delete_safe v eq)
  15.838 +        fun elim e =
  15.839 +         let val b = Inttriplefunc.tryapplyd e v rat_0
  15.840 +         in if b =/ rat_0 then e
  15.841 +            else pi_equation_add e (pi_equation_cmul (Rat.neg b // a) eq)
  15.842 +         end
  15.843 +    in eliminate (Inttriplefunc.update(v, eq') (Inttriplefunc.map (K elim) dun))
  15.844 +                 (map elim oeqs)
  15.845 +    end
  15.846 +in fn eqs =>
  15.847 + let
  15.848 +  val assig = eliminate Inttriplefunc.empty eqs
  15.849 +  val vs = Inttriplefunc.fold (fn (x, f) => fn a => remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig []
  15.850 + in (distinct (dest_ord triple_int_ord) vs,assig)
  15.851 + end
  15.852 +end;
  15.853 +
  15.854 +(* Solve equations by assigning arbitrary numbers.                           *)
  15.855 +
  15.856 +fun pi_solve_equations one eqs =
  15.857 + let
  15.858 +  val (vars,assigs) = pi_eliminate_all_equations one eqs
  15.859 +  val vfn = fold_rev (fn v => Inttriplefunc.update(v,rat_0)) vars
  15.860 +            (Inttriplefunc.onefunc(one, Rat.rat_of_int ~1))
  15.861 +  val ass =
  15.862 +    Inttriplefunc.combine (curry op +/) (K false)
  15.863 +    (Inttriplefunc.map (K (pi_equation_eval vfn)) assigs) vfn
  15.864 + in if forall (fn e => pi_equation_eval ass e =/ rat_0) eqs
  15.865 +    then Inttriplefunc.delete_safe one ass else raise Sanity
  15.866 + end;
  15.867 +
  15.868 +(* Multiply equation-parametrized poly by regular poly and add accumulator.  *)
  15.869 +
  15.870 +fun pi_epoly_pmul p q acc =
  15.871 + FuncUtil.Monomialfunc.fold (fn (m1, c) => fn a =>
  15.872 +  FuncUtil.Monomialfunc.fold (fn (m2,e) => fn b =>
  15.873 +   let val m =  monomial_mul m1 m2
  15.874 +       val es = FuncUtil.Monomialfunc.tryapplyd b m Inttriplefunc.empty
  15.875 +   in FuncUtil.Monomialfunc.update (m,pi_equation_add (pi_equation_cmul c e) es) b
  15.876 +   end) q a) p acc ;
  15.877 +
  15.878 +(* Usual operations on equation-parametrized poly.                           *)
  15.879 +
  15.880 +fun pi_epoly_cmul c l =
  15.881 +  if c =/ rat_0 then Inttriplefunc.empty else Inttriplefunc.map (K (pi_equation_cmul c)) l;;
  15.882 +
  15.883 +val pi_epoly_neg = pi_epoly_cmul (Rat.rat_of_int ~1);
  15.884 +
  15.885 +val pi_epoly_add = Inttriplefunc.combine pi_equation_add Inttriplefunc.is_empty;
  15.886 +
  15.887 +fun pi_epoly_sub p q = pi_epoly_add p (pi_epoly_neg q);;
  15.888 +
  15.889 +fun allpairs f l1 l2 =  fold_rev (fn x => (curry (op @)) (map (f x) l2)) l1 [];
  15.890 +
  15.891 +(* Hence produce the "relevant" monomials: those whose squares lie in the    *)
  15.892 +(* Newton polytope of the monomials in the input. (This is enough according  *)
  15.893 +(* to Reznik: "Extremal PSD forms with few terms", Duke Math. Journal,       *)
  15.894 +(* vol 45, pp. 363--374, 1978.                                               *)
  15.895 +(*                                                                           *)
  15.896 +(* These are ordered in sort of decreasing degree. In particular the         *)
  15.897 +(* constant monomial is last; this gives an order in diagonalization of the  *)
  15.898 +(* quadratic form that will tend to display constants.                       *)
  15.899 +
  15.900 +(* Diagonalize (Cholesky/LDU) the matrix corresponding to a quadratic form.  *)
  15.901 +
  15.902 +local
  15.903 +fun diagonalize n i m =
  15.904 + if FuncUtil.Intpairfunc.is_empty (snd m) then []
  15.905 + else
  15.906 +  let val a11 = FuncUtil.Intpairfunc.tryapplyd (snd m) (i,i) rat_0
  15.907 +  in if a11 </ rat_0 then raise Failure "diagonalize: not PSD"
  15.908 +    else if a11 =/ rat_0 then
  15.909 +          if FuncUtil.Intfunc.is_empty (snd (row i m)) then diagonalize n (i + 1) m
  15.910 +          else raise Failure "diagonalize: not PSD ___ "
  15.911 +    else
  15.912 +     let
  15.913 +      val v = row i m
  15.914 +      val v' = (fst v, FuncUtil.Intfunc.fold (fn (i, c) => fn a =>
  15.915 +       let val y = c // a11
  15.916 +       in if y = rat_0 then a else FuncUtil.Intfunc.update (i,y) a
  15.917 +       end)  (snd v) FuncUtil.Intfunc.empty)
  15.918 +      fun upt0 x y a = if y = rat_0 then a else FuncUtil.Intpairfunc.update (x,y) a
  15.919 +      val m' =
  15.920 +      ((n,n),
  15.921 +      iter (i+1,n) (fn j =>
  15.922 +          iter (i+1,n) (fn k =>
  15.923 +              (upt0 (j,k) (FuncUtil.Intpairfunc.tryapplyd (snd m) (j,k) rat_0 -/ FuncUtil.Intfunc.tryapplyd (snd v) j rat_0 */ FuncUtil.Intfunc.tryapplyd (snd v') k rat_0))))
  15.924 +          FuncUtil.Intpairfunc.empty)
  15.925 +     in (a11,v')::diagonalize n (i + 1) m'
  15.926 +     end
  15.927 +  end
  15.928 +in
  15.929 +fun diag m =
  15.930 + let
  15.931 +   val nn = dimensions m
  15.932 +   val n = fst nn
  15.933 + in if snd nn <> n then error "diagonalize: non-square matrix"
  15.934 +    else diagonalize n 1 m
  15.935 + end
  15.936 +end;
  15.937 +
  15.938 +fun gcd_rat a b = Rat.rat_of_int (Integer.gcd (int_of_rat a) (int_of_rat b));
  15.939 +
  15.940 +(* Adjust a diagonalization to collect rationals at the start.               *)
  15.941 +  (* FIXME : Potentially polymorphic keys, but here only: integers!! *)
  15.942 +local
  15.943 + fun upd0 x y a = if y =/ rat_0 then a else FuncUtil.Intfunc.update(x,y) a;
  15.944 + fun mapa f (d,v) =
  15.945 +  (d, FuncUtil.Intfunc.fold (fn (i,c) => fn a => upd0 i (f c) a) v FuncUtil.Intfunc.empty)
  15.946 + fun adj (c,l) =
  15.947 + let val a =
  15.948 +  FuncUtil.Intfunc.fold (fn (i,c) => fn a => lcm_rat a (denominator_rat c))
  15.949 +    (snd l) rat_1 //
  15.950 +  FuncUtil.Intfunc.fold (fn (i,c) => fn a => gcd_rat a (numerator_rat c))
  15.951 +    (snd l) rat_0
  15.952 +  in ((c // (a */ a)),mapa (fn x => a */ x) l)
  15.953 +  end
  15.954 +in
  15.955 +fun deration d = if null d then (rat_0,d) else
  15.956 + let val d' = map adj d
  15.957 +     val a = fold (lcm_rat o denominator_rat o fst) d' rat_1 //
  15.958 +          fold (gcd_rat o numerator_rat o fst) d' rat_0
  15.959 + in ((rat_1 // a),map (fn (c,l) => (a */ c,l)) d')
  15.960 + end
  15.961 +end;
  15.962 +
  15.963 +(* Enumeration of monomials with given multidegree bound.                    *)
  15.964 +
  15.965 +fun enumerate_monomials d vars =
  15.966 + if d < 0 then []
  15.967 + else if d = 0 then [FuncUtil.Ctermfunc.empty]
  15.968 + else if null vars then [monomial_1] else
  15.969 + let val alts =
  15.970 +  map_range (fn k => let val oths = enumerate_monomials (d - k) (tl vars)
  15.971 +               in map (fn ks => if k = 0 then ks else FuncUtil.Ctermfunc.update (hd vars, k) ks) oths end) (d + 1)
  15.972 + in flat alts
  15.973 + end;
  15.974 +
  15.975 +(* Enumerate products of distinct input polys with degree <= d.              *)
  15.976 +(* We ignore any constant input polynomials.                                 *)
  15.977 +(* Give the output polynomial and a record of how it was derived.            *)
  15.978 +
  15.979 +fun enumerate_products d pols =
  15.980 +if d = 0 then [(poly_const rat_1,RealArith.Rational_lt rat_1)]
  15.981 +else if d < 0 then [] else
  15.982 +case pols of
  15.983 +   [] => [(poly_const rat_1,RealArith.Rational_lt rat_1)]
  15.984 + | (p,b)::ps =>
  15.985 +    let val e = multidegree p
  15.986 +    in if e = 0 then enumerate_products d ps else
  15.987 +       enumerate_products d ps @
  15.988 +       map (fn (q,c) => (poly_mul p q,RealArith.Product(b,c)))
  15.989 +         (enumerate_products (d - e) ps)
  15.990 +    end
  15.991 +
  15.992 +(* Convert regular polynomial. Note that we treat (0,0,0) as -1.             *)
  15.993 +
  15.994 +fun epoly_of_poly p =
  15.995 +  FuncUtil.Monomialfunc.fold (fn (m,c) => fn a => FuncUtil.Monomialfunc.update (m, Inttriplefunc.onefunc ((0,0,0), Rat.neg c)) a) p FuncUtil.Monomialfunc.empty;
  15.996 +
  15.997 +(* String for block diagonal matrix numbered k.                              *)
  15.998 +
  15.999 +fun sdpa_of_blockdiagonal k m =
 15.1000 + let
 15.1001 +  val pfx = string_of_int k ^" "
 15.1002 +  val ents =
 15.1003 +    Inttriplefunc.fold
 15.1004 +      (fn ((b,i,j),c) => fn a => if i > j then a else ((b,i,j),c)::a)
 15.1005 +      m []
 15.1006 +  val entss = sort (triple_int_ord o pairself fst) ents
 15.1007 + in fold_rev (fn ((b,i,j),c) => fn a =>
 15.1008 +     pfx ^ string_of_int b ^ " " ^ string_of_int i ^ " " ^ string_of_int j ^
 15.1009 +     " " ^ decimalize 20 c ^ "\n" ^ a) entss ""
 15.1010 + end;
 15.1011 +
 15.1012 +(* SDPA for problem using block diagonal (i.e. multiple SDPs)                *)
 15.1013 +
 15.1014 +fun sdpa_of_blockproblem nblocks blocksizes obj mats =
 15.1015 + let val m = length mats - 1
 15.1016 + in
 15.1017 +  string_of_int m ^ "\n" ^
 15.1018 +  string_of_int nblocks ^ "\n" ^
 15.1019 +  (space_implode " " (map string_of_int blocksizes)) ^
 15.1020 +  "\n" ^
 15.1021 +  sdpa_of_vector obj ^
 15.1022 +  fold_rev2 (fn k => fn m => fn a => sdpa_of_blockdiagonal (k - 1) m ^ a)
 15.1023 +    (1 upto length mats) mats ""
 15.1024 + end;
 15.1025 +
 15.1026 +(* Run prover on a problem in block diagonal form.                       *)
 15.1027 +
 15.1028 +fun run_blockproblem prover nblocks blocksizes obj mats=
 15.1029 +  parse_csdpoutput (prover (sdpa_of_blockproblem nblocks blocksizes obj mats))
 15.1030 +
 15.1031 +(* 3D versions of matrix operations to consider blocks separately.           *)
 15.1032 +
 15.1033 +val bmatrix_add = Inttriplefunc.combine (curry op +/) (fn x => x =/ rat_0);
 15.1034 +fun bmatrix_cmul c bm =
 15.1035 +  if c =/ rat_0 then Inttriplefunc.empty
 15.1036 +  else Inttriplefunc.map (fn _ => fn x => c */ x) bm;
 15.1037 +
 15.1038 +val bmatrix_neg = bmatrix_cmul (Rat.rat_of_int ~1);
 15.1039 +fun bmatrix_sub m1 m2 = bmatrix_add m1 (bmatrix_neg m2);;
 15.1040 +
 15.1041 +(* Smash a block matrix into components.                                     *)
 15.1042 +
 15.1043 +fun blocks blocksizes bm =
 15.1044 + map (fn (bs,b0) =>
 15.1045 +      let val m = Inttriplefunc.fold
 15.1046 +          (fn ((b,i,j),c) => fn a => if b = b0 then FuncUtil.Intpairfunc.update ((i,j),c) a else a) bm FuncUtil.Intpairfunc.empty
 15.1047 +          val d = FuncUtil.Intpairfunc.fold (fn ((i,j),c) => fn a => max a (max i j)) m 0
 15.1048 +      in (((bs,bs),m):matrix) end)
 15.1049 + (blocksizes ~~ (1 upto length blocksizes));;
 15.1050 +
 15.1051 +(* FIXME : Get rid of this !!!*)
 15.1052 +local
 15.1053 +  fun tryfind_with msg f [] = raise Failure msg
 15.1054 +    | tryfind_with msg f (x::xs) = (f x handle Failure s => tryfind_with s f xs);
 15.1055 +in
 15.1056 +  fun tryfind f = tryfind_with "tryfind" f
 15.1057 +end
 15.1058 +
 15.1059 +(* Positiv- and Nullstellensatz. Flag "linf" forces a linear representation. *)
 15.1060 +
 15.1061 +
 15.1062 +fun real_positivnullstellensatz_general ctxt prover linf d eqs leqs pol =
 15.1063 +let
 15.1064 + val vars = fold_rev (union (op aconvc) o poly_variables)
 15.1065 +   (pol :: eqs @ map fst leqs) []
 15.1066 + val monoid = if linf then
 15.1067 +      (poly_const rat_1,RealArith.Rational_lt rat_1)::
 15.1068 +      (filter (fn (p,c) => multidegree p <= d) leqs)
 15.1069 +    else enumerate_products d leqs
 15.1070 + val nblocks = length monoid
 15.1071 + fun mk_idmultiplier k p =
 15.1072 +  let
 15.1073 +   val e = d - multidegree p
 15.1074 +   val mons = enumerate_monomials e vars
 15.1075 +   val nons = mons ~~ (1 upto length mons)
 15.1076 +  in (mons,
 15.1077 +      fold_rev (fn (m,n) => FuncUtil.Monomialfunc.update(m,Inttriplefunc.onefunc((~k,~n,n),rat_1))) nons FuncUtil.Monomialfunc.empty)
 15.1078 +  end
 15.1079 +
 15.1080 + fun mk_sqmultiplier k (p,c) =
 15.1081 +  let
 15.1082 +   val e = (d - multidegree p) div 2
 15.1083 +   val mons = enumerate_monomials e vars
 15.1084 +   val nons = mons ~~ (1 upto length mons)
 15.1085 +  in (mons,
 15.1086 +      fold_rev (fn (m1,n1) =>
 15.1087 +       fold_rev (fn (m2,n2) => fn  a =>
 15.1088 +        let val m = monomial_mul m1 m2
 15.1089 +        in if n1 > n2 then a else
 15.1090 +          let val c = if n1 = n2 then rat_1 else rat_2
 15.1091 +              val e = FuncUtil.Monomialfunc.tryapplyd a m Inttriplefunc.empty
 15.1092 +          in FuncUtil.Monomialfunc.update(m, tri_equation_add (Inttriplefunc.onefunc((k,n1,n2), c)) e) a
 15.1093 +          end
 15.1094 +        end)  nons)
 15.1095 +       nons FuncUtil.Monomialfunc.empty)
 15.1096 +  end
 15.1097 +
 15.1098 +  val (sqmonlist,sqs) = split_list (map2 mk_sqmultiplier (1 upto length monoid) monoid)
 15.1099 +  val (idmonlist,ids) =  split_list(map2 mk_idmultiplier (1 upto length eqs) eqs)
 15.1100 +  val blocksizes = map length sqmonlist
 15.1101 +  val bigsum =
 15.1102 +    fold_rev2 (fn p => fn q => fn a => tri_epoly_pmul p q a) eqs ids
 15.1103 +            (fold_rev2 (fn (p,c) => fn s => fn a => tri_epoly_pmul p s a) monoid sqs
 15.1104 +                     (epoly_of_poly(poly_neg pol)))
 15.1105 +  val eqns = FuncUtil.Monomialfunc.fold (fn (m,e) => fn a => e::a) bigsum []
 15.1106 +  val (pvs,assig) = tri_eliminate_all_equations (0,0,0) eqns
 15.1107 +  val qvars = (0,0,0)::pvs
 15.1108 +  val allassig = fold_rev (fn v => Inttriplefunc.update(v,(Inttriplefunc.onefunc(v,rat_1)))) pvs assig
 15.1109 +  fun mk_matrix v =
 15.1110 +    Inttriplefunc.fold (fn ((b,i,j), ass) => fn m =>
 15.1111 +        if b < 0 then m else
 15.1112 +         let val c = Inttriplefunc.tryapplyd ass v rat_0
 15.1113 +         in if c = rat_0 then m else
 15.1114 +            Inttriplefunc.update ((b,j,i), c) (Inttriplefunc.update ((b,i,j), c) m)
 15.1115 +         end)
 15.1116 +          allassig Inttriplefunc.empty
 15.1117 +  val diagents = Inttriplefunc.fold
 15.1118 +    (fn ((b,i,j), e) => fn a => if b > 0 andalso i = j then tri_equation_add e a else a)
 15.1119 +    allassig Inttriplefunc.empty
 15.1120 +
 15.1121 +  val mats = map mk_matrix qvars
 15.1122 +  val obj = (length pvs,
 15.1123 +            itern 1 pvs (fn v => fn i => FuncUtil.Intfunc.updatep iszero (i,Inttriplefunc.tryapplyd diagents v rat_0))
 15.1124 +                        FuncUtil.Intfunc.empty)
 15.1125 +  val raw_vec = if null pvs then vector_0 0
 15.1126 +                else tri_scale_then (run_blockproblem prover nblocks blocksizes) obj mats
 15.1127 +  fun int_element (d,v) i = FuncUtil.Intfunc.tryapplyd v i rat_0
 15.1128 +
 15.1129 +  fun find_rounding d =
 15.1130 +   let
 15.1131 +    val _ =
 15.1132 +      if Config.get ctxt trace
 15.1133 +      then writeln ("Trying rounding with limit "^Rat.string_of_rat d ^ "\n")
 15.1134 +      else ()
 15.1135 +    val vec = nice_vector d raw_vec
 15.1136 +    val blockmat = iter (1,dim vec)
 15.1137 +     (fn i => fn a => bmatrix_add (bmatrix_cmul (int_element vec i) (nth mats i)) a)
 15.1138 +     (bmatrix_neg (nth mats 0))
 15.1139 +    val allmats = blocks blocksizes blockmat
 15.1140 +   in (vec,map diag allmats)
 15.1141 +   end
 15.1142 +  val (vec,ratdias) =
 15.1143 +    if null pvs then find_rounding rat_1
 15.1144 +    else tryfind find_rounding (map Rat.rat_of_int (1 upto 31) @
 15.1145 +                                map pow2 (5 upto 66))
 15.1146 +  val newassigs =
 15.1147 +    fold_rev (fn k => Inttriplefunc.update (nth pvs (k - 1), int_element vec k))
 15.1148 +           (1 upto dim vec) (Inttriplefunc.onefunc ((0,0,0), Rat.rat_of_int ~1))
 15.1149 +  val finalassigs =
 15.1150 +    Inttriplefunc.fold (fn (v,e) => fn a => Inttriplefunc.update(v, tri_equation_eval newassigs e) a) allassig newassigs
 15.1151 +  fun poly_of_epoly p =
 15.1152 +    FuncUtil.Monomialfunc.fold (fn (v,e) => fn a => FuncUtil.Monomialfunc.updatep iszero (v,tri_equation_eval finalassigs e) a)
 15.1153 +          p FuncUtil.Monomialfunc.empty
 15.1154 +  fun  mk_sos mons =
 15.1155 +   let fun mk_sq (c,m) =
 15.1156 +    (c,fold_rev (fn k=> fn a => FuncUtil.Monomialfunc.updatep iszero (nth mons (k - 1), int_element m k) a)
 15.1157 +                 (1 upto length mons) FuncUtil.Monomialfunc.empty)
 15.1158 +   in map mk_sq
 15.1159 +   end
 15.1160 +  val sqs = map2 mk_sos sqmonlist ratdias
 15.1161 +  val cfs = map poly_of_epoly ids
 15.1162 +  val msq = filter (fn (a,b) => not (null b)) (map2 pair monoid sqs)
 15.1163 +  fun eval_sq sqs = fold_rev (fn (c,q) => poly_add (poly_cmul c (poly_mul q q))) sqs poly_0
 15.1164 +  val sanity =
 15.1165 +    fold_rev (fn ((p,c),s) => poly_add (poly_mul p (eval_sq s))) msq
 15.1166 +           (fold_rev2 (fn p => fn q => poly_add (poly_mul p q)) cfs eqs
 15.1167 +                    (poly_neg pol))
 15.1168 +
 15.1169 +in if not(FuncUtil.Monomialfunc.is_empty sanity) then raise Sanity else
 15.1170 +  (cfs,map (fn (a,b) => (snd a,b)) msq)
 15.1171 + end
 15.1172 +
 15.1173 +
 15.1174 +(* Iterative deepening.                                                      *)
 15.1175 +
 15.1176 +fun deepen f n =
 15.1177 +  (writeln ("Searching with depth limit " ^ string_of_int n);
 15.1178 +    (f n handle Failure s => (writeln ("failed with message: " ^ s); deepen f (n + 1))));
 15.1179 +
 15.1180 +
 15.1181 +(* Map back polynomials and their composites to a positivstellensatz.        *)
 15.1182 +
 15.1183 +fun cterm_of_sqterm (c,p) = RealArith.Product(RealArith.Rational_lt c,RealArith.Square p);
 15.1184 +
 15.1185 +fun cterm_of_sos (pr,sqs) = if null sqs then pr
 15.1186 +  else RealArith.Product(pr,foldr1 RealArith.Sum (map cterm_of_sqterm sqs));
 15.1187 +
 15.1188 +(* Interface to HOL.                                                         *)
 15.1189 +local
 15.1190 +  open Conv
 15.1191 +  val concl = Thm.dest_arg o cprop_of
 15.1192 +  fun simple_cterm_ord t u = Term_Ord.fast_term_ord (term_of t, term_of u) = LESS
 15.1193 +in
 15.1194 +  (* FIXME: Replace tryfind by get_first !! *)
 15.1195 +fun real_nonlinear_prover proof_method ctxt =
 15.1196 + let
 15.1197 +  val {add,mul,neg,pow,sub,main} =  Semiring_Normalizer.semiring_normalizers_ord_wrapper ctxt
 15.1198 +      (the (Semiring_Normalizer.match ctxt @{cterm "(0::real) + 1"}))
 15.1199 +     simple_cterm_ord
 15.1200 +  val (real_poly_add_conv,real_poly_mul_conv,real_poly_neg_conv,
 15.1201 +       real_poly_pow_conv,real_poly_sub_conv,real_poly_conv) = (add,mul,neg,pow,sub,main)
 15.1202 +  fun mainf cert_choice translator (eqs,les,lts) =
 15.1203 +  let
 15.1204 +   val eq0 = map (poly_of_term o Thm.dest_arg1 o concl) eqs
 15.1205 +   val le0 = map (poly_of_term o Thm.dest_arg o concl) les
 15.1206 +   val lt0 = map (poly_of_term o Thm.dest_arg o concl) lts
 15.1207 +   val eqp0 = map_index (fn (i, t) => (t,RealArith.Axiom_eq i)) eq0
 15.1208 +   val lep0 = map_index (fn (i, t) => (t,RealArith.Axiom_le i)) le0
 15.1209 +   val ltp0 = map_index (fn (i, t) => (t,RealArith.Axiom_lt i)) lt0
 15.1210 +   val (keq,eq) = List.partition (fn (p,_) => multidegree p = 0) eqp0
 15.1211 +   val (klep,lep) = List.partition (fn (p,_) => multidegree p = 0) lep0
 15.1212 +   val (kltp,ltp) = List.partition (fn (p,_) => multidegree p = 0) ltp0
 15.1213 +   fun trivial_axiom (p,ax) =
 15.1214 +    case ax of
 15.1215 +       RealArith.Axiom_eq n => if eval FuncUtil.Ctermfunc.empty p <>/ Rat.zero then nth eqs n
 15.1216 +                     else raise Failure "trivial_axiom: Not a trivial axiom"
 15.1217 +     | RealArith.Axiom_le n => if eval FuncUtil.Ctermfunc.empty p </ Rat.zero then nth les n
 15.1218 +                     else raise Failure "trivial_axiom: Not a trivial axiom"
 15.1219 +     | RealArith.Axiom_lt n => if eval FuncUtil.Ctermfunc.empty p <=/ Rat.zero then nth lts n
 15.1220 +                     else raise Failure "trivial_axiom: Not a trivial axiom"
 15.1221 +     | _ => error "trivial_axiom: Not a trivial axiom"
 15.1222 +   in
 15.1223 +  (let val th = tryfind trivial_axiom (keq @ klep @ kltp)
 15.1224 +   in
 15.1225 +    (fconv_rule (arg_conv (arg1_conv real_poly_conv) then_conv Numeral_Simprocs.field_comp_conv) th, RealArith.Trivial)
 15.1226 +   end)
 15.1227 +   handle Failure _ =>
 15.1228 +     (let val proof =
 15.1229 +       (case proof_method of Certificate certs =>
 15.1230 +         (* choose certificate *)
 15.1231 +         let
 15.1232 +           fun chose_cert [] (RealArith.Cert c) = c
 15.1233 +             | chose_cert (RealArith.Left::s) (RealArith.Branch (l, _)) = chose_cert s l
 15.1234 +             | chose_cert (RealArith.Right::s) (RealArith.Branch (_, r)) = chose_cert s r
 15.1235 +             | chose_cert _ _ = error "certificate tree in invalid form"
 15.1236 +         in
 15.1237 +           chose_cert cert_choice certs
 15.1238 +         end
 15.1239 +       | Prover prover =>
 15.1240 +         (* call prover *)
 15.1241 +         let
 15.1242 +          val pol = fold_rev poly_mul (map fst ltp) (poly_const Rat.one)
 15.1243 +          val leq = lep @ ltp
 15.1244 +          fun tryall d =
 15.1245 +           let val e = multidegree pol
 15.1246 +               val k = if e = 0 then 0 else d div e
 15.1247 +               val eq' = map fst eq
 15.1248 +           in tryfind (fn i => (d,i,real_positivnullstellensatz_general ctxt prover false d eq' leq
 15.1249 +                                 (poly_neg(poly_pow pol i))))
 15.1250 +                   (0 upto k)
 15.1251 +           end
 15.1252 +         val (d,i,(cert_ideal,cert_cone)) = deepen tryall 0
 15.1253 +         val proofs_ideal =
 15.1254 +           map2 (fn q => fn (p,ax) => RealArith.Eqmul(q,ax)) cert_ideal eq
 15.1255 +         val proofs_cone = map cterm_of_sos cert_cone
 15.1256 +         val proof_ne = if null ltp then RealArith.Rational_lt Rat.one else
 15.1257 +           let val p = foldr1 RealArith.Product (map snd ltp)
 15.1258 +           in  funpow i (fn q => RealArith.Product(p,q)) (RealArith.Rational_lt Rat.one)
 15.1259 +           end
 15.1260 +         in
 15.1261 +           foldr1 RealArith.Sum (proof_ne :: proofs_ideal @ proofs_cone)
 15.1262 +         end)
 15.1263 +     in
 15.1264 +        (translator (eqs,les,lts) proof, RealArith.Cert proof)
 15.1265 +     end)
 15.1266 +   end
 15.1267 + in mainf end
 15.1268 +end
 15.1269 +
 15.1270 +fun C f x y = f y x;
 15.1271 +  (* FIXME : This is very bad!!!*)
 15.1272 +fun subst_conv eqs t =
 15.1273 + let
 15.1274 +  val t' = fold (Thm.cabs o Thm.lhs_of) eqs t
 15.1275 + in Conv.fconv_rule (Thm.beta_conversion true) (fold (C Thm.combination) eqs (Thm.reflexive t'))
 15.1276 + end
 15.1277 +
 15.1278 +(* A wrapper that tries to substitute away variables first.                  *)
 15.1279 +
 15.1280 +local
 15.1281 + open Conv
 15.1282 +  fun simple_cterm_ord t u = Term_Ord.fast_term_ord (term_of t, term_of u) = LESS
 15.1283 + val concl = Thm.dest_arg o cprop_of
 15.1284 + val shuffle1 =
 15.1285 +   fconv_rule (rewr_conv @{lemma "(a + x == y) == (x == y - (a::real))" by (atomize (full)) (simp add: field_simps) })
 15.1286 + val shuffle2 =
 15.1287 +    fconv_rule (rewr_conv @{lemma "(x + a == y) ==  (x == y - (a::real))" by (atomize (full)) (simp add: field_simps)})
 15.1288 + fun substitutable_monomial fvs tm = case term_of tm of
 15.1289 +    Free(_,@{typ real}) => if not (member (op aconvc) fvs tm) then (Rat.one,tm)
 15.1290 +                           else raise Failure "substitutable_monomial"
 15.1291 +  | @{term "op * :: real => _"}$c$(t as Free _ ) =>
 15.1292 +     if RealArith.is_ratconst (Thm.dest_arg1 tm) andalso not (member (op aconvc) fvs (Thm.dest_arg tm))
 15.1293 +         then (RealArith.dest_ratconst (Thm.dest_arg1 tm),Thm.dest_arg tm) else raise Failure "substitutable_monomial"
 15.1294 +  | @{term "op + :: real => _"}$s$t =>
 15.1295 +       (substitutable_monomial (Thm.add_cterm_frees (Thm.dest_arg tm) fvs) (Thm.dest_arg1 tm)
 15.1296 +        handle Failure _ => substitutable_monomial (Thm.add_cterm_frees (Thm.dest_arg1 tm) fvs) (Thm.dest_arg tm))
 15.1297 +  | _ => raise Failure "substitutable_monomial"
 15.1298 +
 15.1299 +  fun isolate_variable v th =
 15.1300 +   let val w = Thm.dest_arg1 (cprop_of th)
 15.1301 +   in if v aconvc w then th
 15.1302 +      else case term_of w of
 15.1303 +           @{term "op + :: real => _"}$s$t =>
 15.1304 +              if Thm.dest_arg1 w aconvc v then shuffle2 th
 15.1305 +              else isolate_variable v (shuffle1 th)
 15.1306 +          | _ => error "isolate variable : This should not happen?"
 15.1307 +   end
 15.1308 +in
 15.1309 +
 15.1310 +fun real_nonlinear_subst_prover prover ctxt =
 15.1311 + let
 15.1312 +  val {add,mul,neg,pow,sub,main} =  Semiring_Normalizer.semiring_normalizers_ord_wrapper ctxt
 15.1313 +      (the (Semiring_Normalizer.match ctxt @{cterm "(0::real) + 1"}))
 15.1314 +     simple_cterm_ord
 15.1315 +
 15.1316 +  val (real_poly_add_conv,real_poly_mul_conv,real_poly_neg_conv,
 15.1317 +       real_poly_pow_conv,real_poly_sub_conv,real_poly_conv) = (add,mul,neg,pow,sub,main)
 15.1318 +
 15.1319 +  fun make_substitution th =
 15.1320 +   let
 15.1321 +    val (c,v) = substitutable_monomial [] (Thm.dest_arg1(concl th))
 15.1322 +    val th1 = Drule.arg_cong_rule (Thm.capply @{cterm "op * :: real => _"} (RealArith.cterm_of_rat (Rat.inv c))) (mk_meta_eq th)
 15.1323 +    val th2 = fconv_rule (binop_conv real_poly_mul_conv) th1
 15.1324 +   in fconv_rule (arg_conv real_poly_conv) (isolate_variable v th2)
 15.1325 +   end
 15.1326 +   fun oprconv cv ct =
 15.1327 +    let val g = Thm.dest_fun2 ct
 15.1328 +    in if g aconvc @{cterm "op <= :: real => _"}
 15.1329 +         orelse g aconvc @{cterm "op < :: real => _"}
 15.1330 +       then arg_conv cv ct else arg1_conv cv ct
 15.1331 +    end
 15.1332 +  fun mainf cert_choice translator =
 15.1333 +   let
 15.1334 +    fun substfirst(eqs,les,lts) =
 15.1335 +      ((let
 15.1336 +           val eth = tryfind make_substitution eqs
 15.1337 +           val modify = fconv_rule (arg_conv (oprconv(subst_conv [eth] then_conv real_poly_conv)))
 15.1338 +       in  substfirst
 15.1339 +             (filter_out (fn t => (Thm.dest_arg1 o Thm.dest_arg o cprop_of) t
 15.1340 +                                   aconvc @{cterm "0::real"}) (map modify eqs),
 15.1341 +                                   map modify les,map modify lts)
 15.1342 +       end)
 15.1343 +       handle Failure  _ => real_nonlinear_prover prover ctxt cert_choice translator (rev eqs, rev les, rev lts))
 15.1344 +    in substfirst
 15.1345 +   end
 15.1346 +
 15.1347 +
 15.1348 + in mainf
 15.1349 + end
 15.1350 +
 15.1351 +(* Overall function. *)
 15.1352 +
 15.1353 +fun real_sos prover ctxt =
 15.1354 +  RealArith.gen_prover_real_arith ctxt (real_nonlinear_subst_prover prover ctxt)
 15.1355 +end;
 15.1356 +
 15.1357 +val known_sos_constants =
 15.1358 +  [@{term "op ==>"}, @{term "Trueprop"},
 15.1359 +   @{term HOL.implies}, @{term HOL.conj}, @{term HOL.disj},
 15.1360 +   @{term "Not"}, @{term "op = :: bool => _"},
 15.1361 +   @{term "All :: (real => _) => _"}, @{term "Ex :: (real => _) => _"},
 15.1362 +   @{term "op = :: real => _"}, @{term "op < :: real => _"},
 15.1363 +   @{term "op <= :: real => _"},
 15.1364 +   @{term "op + :: real => _"}, @{term "op - :: real => _"},
 15.1365 +   @{term "op * :: real => _"}, @{term "uminus :: real => _"},
 15.1366 +   @{term "op / :: real => _"}, @{term "inverse :: real => _"},
 15.1367 +   @{term "op ^ :: real => _"}, @{term "abs :: real => _"},
 15.1368 +   @{term "min :: real => _"}, @{term "max :: real => _"},
 15.1369 +   @{term "0::real"}, @{term "1::real"}, @{term "number_of :: int => real"},
 15.1370 +   @{term "number_of :: int => nat"},
 15.1371 +   @{term "Int.Bit0"}, @{term "Int.Bit1"},
 15.1372 +   @{term "Int.Pls"}, @{term "Int.Min"}];
 15.1373 +
 15.1374 +fun check_sos kcts ct =
 15.1375 + let
 15.1376 +  val t = term_of ct
 15.1377 +  val _ = if not (null (Term.add_tfrees t [])
 15.1378 +                  andalso null (Term.add_tvars t []))
 15.1379 +          then error "SOS: not sos. Additional type varables" else ()
 15.1380 +  val fs = Term.add_frees t []
 15.1381 +  val _ = if exists (fn ((_,T)) => not (T = @{typ "real"})) fs
 15.1382 +          then error "SOS: not sos. Variables with type not real" else ()
 15.1383 +  val vs = Term.add_vars t []
 15.1384 +  val _ = if exists (fn ((_,T)) => not (T = @{typ "real"})) vs
 15.1385 +          then error "SOS: not sos. Variables with type not real" else ()
 15.1386 +  val ukcs = subtract (fn (t,p) => Const p aconv t) kcts (Term.add_consts t [])
 15.1387 +  val _ = if  null ukcs then ()
 15.1388 +              else error ("SOSO: Unknown constants in Subgoal:" ^ commas (map fst ukcs))
 15.1389 +in () end
 15.1390 +
 15.1391 +fun core_sos_tac print_cert prover = SUBPROOF (fn {concl, context, ...} =>
 15.1392 +  let
 15.1393 +    val _ = check_sos known_sos_constants concl
 15.1394 +    val (ths, certificates) = real_sos prover context (Thm.dest_arg concl)
 15.1395 +    val _ = print_cert certificates
 15.1396 +  in rtac ths 1 end)
 15.1397 +
 15.1398 +fun default_SOME f NONE v = SOME v
 15.1399 +  | default_SOME f (SOME v) _ = SOME v;
 15.1400 +
 15.1401 +fun lift_SOME f NONE a = f a
 15.1402 +  | lift_SOME f (SOME a) _ = SOME a;
 15.1403 +
 15.1404 +
 15.1405 +local
 15.1406 + val is_numeral = can (HOLogic.dest_number o term_of)
 15.1407 +in
 15.1408 +fun get_denom b ct = case term_of ct of
 15.1409 +  @{term "op / :: real => _"} $ _ $ _ =>
 15.1410 +     if is_numeral (Thm.dest_arg ct) then get_denom b (Thm.dest_arg1 ct)
 15.1411 +     else default_SOME (get_denom b) (get_denom b (Thm.dest_arg ct))   (Thm.dest_arg ct, b)
 15.1412 + | @{term "op < :: real => _"} $ _ $ _ => lift_SOME (get_denom true) (get_denom true (Thm.dest_arg ct)) (Thm.dest_arg1 ct)
 15.1413 + | @{term "op <= :: real => _"} $ _ $ _ => lift_SOME (get_denom true) (get_denom true (Thm.dest_arg ct)) (Thm.dest_arg1 ct)
 15.1414 + | _ $ _ => lift_SOME (get_denom b) (get_denom b (Thm.dest_fun ct)) (Thm.dest_arg ct)
 15.1415 + | _ => NONE
 15.1416 +end;
 15.1417 +
 15.1418 +fun elim_one_denom_tac ctxt =
 15.1419 +CSUBGOAL (fn (P,i) =>
 15.1420 + case get_denom false P of
 15.1421 +   NONE => no_tac
 15.1422 + | SOME (d,ord) =>
 15.1423 +     let
 15.1424 +      val ss = simpset_of ctxt addsimps @{thms field_simps}
 15.1425 +               addsimps [@{thm nonzero_power_divide}, @{thm power_divide}]
 15.1426 +      val th = instantiate' [] [SOME d, SOME (Thm.dest_arg P)]
 15.1427 +         (if ord then @{lemma "(d=0 --> P) & (d>0 --> P) & (d<(0::real) --> P) ==> P" by auto}
 15.1428 +          else @{lemma "(d=0 --> P) & (d ~= (0::real) --> P) ==> P" by blast})
 15.1429 +     in rtac th i THEN Simplifier.asm_full_simp_tac ss i end);
 15.1430 +
 15.1431 +fun elim_denom_tac ctxt i = REPEAT (elim_one_denom_tac ctxt i);
 15.1432 +
 15.1433 +fun sos_tac print_cert prover ctxt =
 15.1434 +  Object_Logic.full_atomize_tac THEN'
 15.1435 +  elim_denom_tac ctxt THEN'
 15.1436 +  core_sos_tac print_cert prover ctxt;
 15.1437 +
 15.1438 +end;