sos method generates and uses proof certificates
authorPhilipp Meyer
Mon Sep 21 15:05:26 2009 +0200 (2009-09-21)
changeset 326451cc5b24f5a01
parent 32644 e4511a1b4c1b
child 32646 962b4354ed90
sos method generates and uses proof certificates
src/HOL/Library/Sum_Of_Squares.thy
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/normarith.ML
src/HOL/Library/positivstellensatz.ML
     1.1 --- a/src/HOL/Library/Sum_Of_Squares.thy	Tue Sep 22 20:25:31 2009 +0200
     1.2 +++ b/src/HOL/Library/Sum_Of_Squares.thy	Mon Sep 21 15:05:26 2009 +0200
     1.3 @@ -10,6 +10,7 @@
     1.4  uses
     1.5    ("positivstellensatz.ML")  (* duplicate use!? -- cf. Euclidian_Space.thy *)
     1.6    ("Sum_Of_Squares/sum_of_squares.ML")
     1.7 +  ("Sum_Of_Squares/positivstellensatz_tools.ML")
     1.8    ("Sum_Of_Squares/sos_wrapper.ML")
     1.9  begin
    1.10  
    1.11 @@ -22,112 +23,142 @@
    1.12    of a minute for one sos call, because sos calls CSDP repeatedly.  If
    1.13    you install CSDP locally, sos calls typically takes only a few
    1.14    seconds.
    1.15 +  sos generates a certificate which can be used to repeat the proof
    1.16 +  without calling an external prover.
    1.17  *}
    1.18  
    1.19  text {* setup sos tactic *}
    1.20  
    1.21  use "positivstellensatz.ML"
    1.22 +use "Sum_Of_Squares/positivstellensatz_tools.ML"
    1.23  use "Sum_Of_Squares/sum_of_squares.ML"
    1.24  use "Sum_Of_Squares/sos_wrapper.ML"
    1.25  
    1.26  setup SosWrapper.setup
    1.27  
    1.28 -text {* Tests -- commented since they work only when csdp is installed
    1.29 -  or take too long with remote csdps *}
    1.30 +text {* Tests *}
    1.31 +
    1.32 +lemma "(3::real) * x + 7 * a < 4 & 3 < 2 * x \<Longrightarrow> a < 0"
    1.33 +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))))))")
    1.34 +
    1.35 +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)"
    1.36 +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))))))))))")
    1.37  
    1.38 -(*
    1.39 -lemma "(3::real) * x + 7 * a < 4 & 3 < 2 * x \<Longrightarrow> a < 0" by sos
    1.40 +lemma "(3::real) * x + 7 * a < 4 & 3 < 2 * x --> a < 0"
    1.41 +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))))))")
    1.42 +
    1.43 +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"
    1.44 +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)))))))")
    1.45  
    1.46 -lemma "a1 >= 0 & a2 >= 0 \<and> (a1 * a1 + a2 * a2 = b1 * b1 + b2 * b2 + 2) \<and>
    1.47 -  (a1 * b1 + a2 * b2 = 0) --> a1 * a2 - b1 * b2 >= (0::real)" by sos
    1.48 +lemma "(0::real) <= x & 0 <= y & 0 <= z & x + y + z <= 3 --> x * y + x * z + y * z >= 3 * x * y * z"
    1.49 +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))))))))))")
    1.50 +
    1.51 +lemma "((x::real)^2 + y^2 + z^2 = 1) --> (x + y + z)^2 <= 3"
    1.52 +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))))))")
    1.53  
    1.54 -lemma "(3::real) * x + 7 * a < 4 & 3 < 2 * x --> a < 0" by sos
    1.55 +lemma "(w^2 + x^2 + y^2 + z^2 = 1) --> (w + x + y + z)^2 <= (4::real)"
    1.56 +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)))))))")
    1.57  
    1.58 -lemma "(0::real) <= x & x <= 1 & 0 <= y & y <= 1 -->
    1.59 -  x^2 + y^2 < 1 |(x - 1)^2 + y^2 < 1 | x^2 + (y - 1)^2 < 1 | (x - 1)^2 + (y - 1)^2 < 1" by sos
    1.60 +lemma "(x::real) >= 1 & y >= 1 --> x * y >= x + y - 1"
    1.61 +by (sos_cert "(((A<0 * R<1) + ((A<=0 * (A<=1 * R<1)) * (R<1 * [1]^2))))")
    1.62 +
    1.63 +lemma "(x::real) > 1 & y > 1 --> x * y > x + y - 1"
    1.64 +by (sos_cert "((((A<0 * A<1) * R<1) + ((A<=0 * R<1) * (R<1 * [1]^2))))") 
    1.65  
    1.66 -lemma "(0::real) <= x & 0 <= y & 0 <= z & x + y + z <= 3 -->
    1.67 -  x * y + x * z + y * z >= 3 * x * y * z" by sos
    1.68 +lemma "abs(x) <= 1 --> abs(64 * x^7 - 112 * x^5 + 56 * x^3 - 7 * x) <= (1::real)"
    1.69 +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)))))")
    1.70 +
    1.71 +(* ------------------------------------------------------------------------- *)
    1.72 +(* One component of denominator in dodecahedral example.                     *)
    1.73 +(* ------------------------------------------------------------------------- *)
    1.74  
    1.75 -lemma "((x::real)^2 + y^2 + z^2 = 1) --> (x + y + z)^2 <= 3" by sos
    1.76 +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)"
    1.77 +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)))))))))))))")
    1.78  
    1.79 -lemma "(w^2 + x^2 + y^2 + z^2 = 1) --> (w + x + y + z)^2 <= (4::real)" by sos
    1.80 +(* ------------------------------------------------------------------------- *)
    1.81 +(* Over a larger but simpler interval.                                       *)
    1.82 +(* ------------------------------------------------------------------------- *)
    1.83  
    1.84 -lemma "(x::real) >= 1 & y >= 1 --> x * y >= x + y - 1" by sos
    1.85 +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)"
    1.86 +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)))))))))))")
    1.87  
    1.88 -lemma "(x::real) > 1 & y > 1 --> x * y > x + y - 1" by sos; 
    1.89 +(* ------------------------------------------------------------------------- *)
    1.90 +(* We can do 12. I think 12 is a sharp bound; see PP's certificate.          *)
    1.91 +(* ------------------------------------------------------------------------- *)
    1.92 +
    1.93 +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)"
    1.94 +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))))))))))))))))))")
    1.95  
    1.96 -lemma "abs(x) <= 1 --> abs(64 * x^7 - 112 * x^5 + 56 * x^3 - 7 * x) <= (1::real)" by sos  
    1.97 +(* ------------------------------------------------------------------------- *)
    1.98 +(* Inequality from sci.math (see "Leon-Sotelo, por favor").                  *)
    1.99 +(* ------------------------------------------------------------------------- *)
   1.100  
   1.101 +lemma "0 <= (x::real) & 0 <= y & (x * y = 1) --> x + y <= x^2 + y^2"
   1.102 +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))))))") 
   1.103  
   1.104 -text {* One component of denominator in dodecahedral example. *}
   1.105 +lemma "0 <= (x::real) & 0 <= y & (x * y = 1) --> x * y * (x + y) <= x^2 + y^2"
   1.106 +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))))))") 
   1.107  
   1.108 -lemma "2 <= x & x <= 125841 / 50000 & 2 <= y & y <= 125841 / 50000 & 2 <= z &
   1.109 -  z <= 125841 / 50000 --> 2 * (x * z + x * y + y * z) - (x * x + y * y + z * z) >= (0::real)" by sos
   1.110 +lemma "0 <= (x::real) & 0 <= y --> x * y * (x + y)^2 <= (x^2 + y^2)^2"
   1.111 +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)))))")
   1.112  
   1.113 +lemma "(0::real) <= a & 0 <= b & 0 <= c & c * (2 * a + b)^3/ 27 <= x \<longrightarrow> c * a^2 * b <= x"
   1.114 +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))))))")
   1.115 + 
   1.116 +lemma "(0::real) < x --> 0 < 1 + x + x^2"
   1.117 +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))))))")
   1.118  
   1.119 -text {* Over a larger but simpler interval. *}
   1.120 +lemma "(0::real) <= x --> 0 < 1 + x + x^2"
   1.121 +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))))))")
   1.122  
   1.123 -lemma "(2::real) <= x & x <= 4 & 2 <= y & y <= 4 & 2 <= z &
   1.124 -  z <= 4 --> 0 <= 2 * (x * z + x * y + y * z) - (x * x + y * y + z * z)" by sos
   1.125 +lemma "(0::real) < 1 + x^2"
   1.126 +by (sos_cert "((R<1 + ((R<1 * (R<1 * [x]^2)) + ((A<=0 * R<1) * (R<1 * [1]^2)))))")
   1.127 +
   1.128 +lemma "(0::real) <= 1 + 2 * x + x^2"
   1.129 +by (sos_cert "(((A<0 * R<1) + (R<1 * (R<1 * [x + 1]^2))))")
   1.130  
   1.131 -text {* We can do 12. I think 12 is a sharp bound; see PP's certificate. *}
   1.132 +lemma "(0::real) < 1 + abs x"
   1.133 +by (sos_cert "((R<1 + (((A<=1 * R<1) * (R<1/2 * [1]^2)) + ((A<=0 * R<1) * (R<1/2 * [1]^2)))))")
   1.134  
   1.135 -lemma "2 <= (x::real) & x <= 4 & 2 <= y & y <= 4 & 2 <= z & z <= 4 -->
   1.136 -  12 <= 2 * (x * z + x * y + y * z) - (x * x + y * y + z * z)" by sos
   1.137 +lemma "(0::real) < 1 + (1 + x)^2 * (abs x)"
   1.138 +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))))))")
   1.139  
   1.140  
   1.141 -text {* Inequality from sci.math (see "Leon-Sotelo, por favor"). *}
   1.142  
   1.143 -lemma "0 <= (x::real) & 0 <= y & (x * y = 1) --> x + y <= x^2 + y^2" by sos 
   1.144 -
   1.145 -lemma "0 <= (x::real) & 0 <= y & (x * y = 1) --> x * y * (x + y) <= x^2 + y^2" by sos 
   1.146 -
   1.147 -lemma "0 <= (x::real) & 0 <= y --> x * y * (x + y)^2 <= (x^2 + y^2)^2" by sos
   1.148 +lemma "abs ((1::real) + x^2) = (1::real) + x^2"
   1.149 +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)))))))")
   1.150 +lemma "(3::real) * x + 7 * a < 4 \<and> 3 < 2 * x \<longrightarrow> a < 0"
   1.151 +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))))))")
   1.152  
   1.153 -lemma "(0::real) <= a & 0 <= b & 0 <= c & c * (2 * a + b)^3/ 27 <= x \<longrightarrow> c * a^2 * b <= x" by sos
   1.154 - 
   1.155 -lemma "(0::real) < x --> 0 < 1 + x + x^2" by sos
   1.156 -
   1.157 -lemma "(0::real) <= x --> 0 < 1 + x + x^2" by sos
   1.158 -
   1.159 -lemma "(0::real) < 1 + x^2" by sos
   1.160 -
   1.161 -lemma "(0::real) <= 1 + 2 * x + x^2" by sos
   1.162 -
   1.163 -lemma "(0::real) < 1 + abs x" by sos
   1.164 -
   1.165 -lemma "(0::real) < 1 + (1 + x)^2 * (abs x)" by sos
   1.166 +lemma "(0::real) < x --> 1 < y --> y * x <= z --> x < z"
   1.167 +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)))))")
   1.168 +lemma "(1::real) < x --> x^2 < y --> 1 < y"
   1.169 +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)))))))))")
   1.170 +lemma "(b::real)^2 < 4 * a * c --> ~(a * x^2 + b * x + c = 0)"
   1.171 +by (sos_cert "(((A<0 * R<1) + (R<1 * (R<1 * [2*a*x + b]^2))))")
   1.172 +lemma "(b::real)^2 < 4 * a * c --> ~(a * x^2 + b * x + c = 0)"
   1.173 +by (sos_cert "(((A<0 * R<1) + (R<1 * (R<1 * [2*a*x + b]^2))))")
   1.174 +lemma "((a::real) * x^2 + b * x + c = 0) --> b^2 >= 4 * a * c"
   1.175 +by (sos_cert "(((A<0 * R<1) + (R<1 * (R<1 * [2*a*x + b]^2))))")
   1.176 +lemma "(0::real) <= b & 0 <= c & 0 <= x & 0 <= y & (x^2 = c) & (y^2 = a^2 * c + b) --> a * c <= y * x"
   1.177 +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)))))")
   1.178 +lemma "abs(x - z) <= e & abs(y - z) <= e & 0 <= u & 0 <= v & (u + v = 1) --> abs((u * x + v * y) - z) <= (e::real)"
   1.179 +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))))))")
   1.180  
   1.181  
   1.182 -lemma "abs ((1::real) + x^2) = (1::real) + x^2" by sos
   1.183 -lemma "(3::real) * x + 7 * a < 4 \<and> 3 < 2 * x \<longrightarrow> a < 0" by sos
   1.184 -
   1.185 -lemma "(0::real) < x --> 1 < y --> y * x <= z --> x < z" by sos
   1.186 -lemma "(1::real) < x --> x^2 < y --> 1 < y" by sos
   1.187 -lemma "(b::real)^2 < 4 * a * c --> ~(a * x^2 + b * x + c = 0)" by sos
   1.188 -lemma "(b::real)^2 < 4 * a * c --> ~(a * x^2 + b * x + c = 0)" by sos
   1.189 -lemma "((a::real) * x^2 + b * x + c = 0) --> b^2 >= 4 * a * c" by sos
   1.190 -lemma "(0::real) <= b & 0 <= c & 0 <= x & 0 <= y & (x^2 = c) & (y^2 = a^2 * c + b) --> a * c <= y * x" by sos
   1.191 -lemma "abs(x - z) <= e & abs(y - z) <= e & 0 <= u & 0 <= v & (u + v = 1) -->
   1.192 -  abs((u * x + v * y) - z) <= (e::real)" by sos
   1.193 -
   1.194 -(*
   1.195 -lemma "((x::real) - y - 2 * x^4 = 0) & 0 <= x & x <= 2 & 0 <= y & y <= 3 -->
   1.196 -  y^2 - 7 * y - 12 * x + 17 >= 0" by sos  -- {* Too hard?*}
   1.197 -*)
   1.198 +(* 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?*)
   1.199  
   1.200  lemma "(0::real) <= x --> (1 + x + x^2)/(1 + x^2) <= 1 + x"
   1.201 -  by sos
   1.202 +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))))))")
   1.203  
   1.204  lemma "(0::real) <= x --> 1 - x <= 1 / (1 + x + x^2)"
   1.205 -  by sos
   1.206 +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))))))))")
   1.207  
   1.208  lemma "(x::real) <= 1 / 2 --> - x - 2 * x^2 <= - x / (1 - x)"
   1.209 -  by sos
   1.210 +by (sos_cert "((((A<0 * A<1) * R<1) + ((A<=0 * (A<0 * R<1)) * (R<1 * [x]^2))))")
   1.211  
   1.212 -lemma "4*r^2 = p^2 - 4*q & r >= (0::real) & x^2 + p*x + q = 0 -->
   1.213 -  2*(x::real) = - p + 2*r | 2*x = -p - 2*r" by sos
   1.214 -*)
   1.215 +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"
   1.216 +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)))))")
   1.217  
   1.218  end
   1.219 +
     2.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     2.2 +++ b/src/HOL/Library/Sum_Of_Squares/positivstellensatz_tools.ML	Mon Sep 21 15:05:26 2009 +0200
     2.3 @@ -0,0 +1,158 @@
     2.4 +(* Title:      positivstellensatz_tools.ML
     2.5 +   Author:     Philipp Meyer, TU Muenchen
     2.6 +
     2.7 +Functions for generating a certificate from a positivstellensatz and vice versa
     2.8 +*)
     2.9 +
    2.10 +signature POSITIVSTELLENSATZ_TOOLS =
    2.11 +sig
    2.12 +  val pss_tree_to_cert : RealArith.pss_tree -> string
    2.13 +
    2.14 +  val cert_to_pss_tree : Proof.context -> string -> RealArith.pss_tree
    2.15 +
    2.16 +end
    2.17 +
    2.18 +
    2.19 +structure PositivstellensatzTools : POSITIVSTELLENSATZ_TOOLS =
    2.20 +struct
    2.21 +
    2.22 +open RealArith FuncUtil
    2.23 +
    2.24 +(*** certificate generation ***)
    2.25 +
    2.26 +fun string_of_rat r =
    2.27 +  let
    2.28 +    val (nom, den) = Rat.quotient_of_rat r
    2.29 +  in
    2.30 +    if den = 1 then string_of_int nom
    2.31 +    else string_of_int nom ^ "/" ^ string_of_int den
    2.32 +  end
    2.33 +
    2.34 +(* map polynomials to strings *)
    2.35 +
    2.36 +fun string_of_varpow x k =
    2.37 +  let
    2.38 +    val term = term_of x
    2.39 +    val name = case term of
    2.40 +      Free (n, _) => n
    2.41 +    | _ => error "Term in monomial not free variable"
    2.42 +  in
    2.43 +    if k = 1 then name else name ^ "^" ^ string_of_int k 
    2.44 +  end
    2.45 +
    2.46 +fun string_of_monomial m = 
    2.47 + if Ctermfunc.is_undefined m then "1" 
    2.48 + else 
    2.49 +  let 
    2.50 +   val m' = dest_monomial m
    2.51 +   val vps = fold_rev (fn (x,k) => cons (string_of_varpow x k)) m' [] 
    2.52 +  in foldr1 (fn (s, t) => s ^ "*" ^ t) vps
    2.53 +  end
    2.54 +
    2.55 +fun string_of_cmonomial (m,c) =
    2.56 +  if Ctermfunc.is_undefined m then string_of_rat c
    2.57 +  else if c = Rat.one then string_of_monomial m
    2.58 +  else (string_of_rat c) ^ "*" ^ (string_of_monomial m);
    2.59 +
    2.60 +fun string_of_poly p = 
    2.61 + if Monomialfunc.is_undefined p then "0" 
    2.62 + else
    2.63 +  let 
    2.64 +   val cms = map string_of_cmonomial
    2.65 +     (sort (prod_ord monomial_order (K EQUAL)) (Monomialfunc.graph p))
    2.66 +  in foldr1 (fn (t1, t2) => t1 ^ " + " ^ t2) cms
    2.67 +  end;
    2.68 +
    2.69 +fun pss_to_cert (Axiom_eq i) = "A=" ^ string_of_int i
    2.70 +  | pss_to_cert (Axiom_le i) = "A<=" ^ string_of_int i
    2.71 +  | pss_to_cert (Axiom_lt i) = "A<" ^ string_of_int i
    2.72 +  | pss_to_cert (Rational_eq r) = "R=" ^ string_of_rat r
    2.73 +  | pss_to_cert (Rational_le r) = "R<=" ^ string_of_rat r
    2.74 +  | pss_to_cert (Rational_lt r) = "R<" ^ string_of_rat r
    2.75 +  | pss_to_cert (Square p) = "[" ^ string_of_poly p ^ "]^2"
    2.76 +  | pss_to_cert (Eqmul (p, pss)) = "([" ^ string_of_poly p ^ "] * " ^ pss_to_cert pss ^ ")"
    2.77 +  | pss_to_cert (Sum (pss1, pss2)) = "(" ^ pss_to_cert pss1 ^ " + " ^ pss_to_cert pss2 ^ ")"
    2.78 +  | pss_to_cert (Product (pss1, pss2)) = "(" ^ pss_to_cert pss1 ^ " * " ^ pss_to_cert pss2 ^ ")"
    2.79 +
    2.80 +fun pss_tree_to_cert Trivial = "()"
    2.81 +  | pss_tree_to_cert (Cert pss) = "(" ^ pss_to_cert pss ^ ")"
    2.82 +  | pss_tree_to_cert (Branch (t1, t2)) = "(" ^ pss_tree_to_cert t1 ^ " & " ^ pss_tree_to_cert t2 ^ ")"
    2.83 +
    2.84 +(*** certificate parsing ***)
    2.85 +
    2.86 +(* basic parser *)
    2.87 +
    2.88 +fun $$ k = Scan.this_string k
    2.89 +
    2.90 +val number = Scan.repeat1 (Scan.one Symbol.is_ascii_digit >>
    2.91 +  (fn s => ord s - ord "0")) >>
    2.92 +  foldl1 (fn (n, d) => n * 10 + d)
    2.93 +
    2.94 +val nat = number
    2.95 +val int = Scan.optional ($$ "~" >> K ~1) 1 -- nat >> op *;
    2.96 +val rat = int --| $$ "/" -- int >> Rat.rat_of_quotient
    2.97 +val rat_int = rat || int >> Rat.rat_of_int
    2.98 +
    2.99 +(* polynomial parser *)
   2.100 +
   2.101 +fun repeat_sep s f = f ::: Scan.repeat ($$ s |-- f)
   2.102 +
   2.103 +val parse_id = Scan.one Symbol.is_letter ::: Scan.many Symbol.is_letdig >> implode
   2.104 +
   2.105 +fun parse_varpow ctxt = parse_id -- Scan.optional ($$ "^" |-- nat) 1 >>
   2.106 +  (fn (x, k) => (cterm_of (Context.theory_of_proof ctxt) (Free (x, @{typ real})), k)) 
   2.107 +
   2.108 +fun parse_monomial ctxt = repeat_sep "*" (parse_varpow ctxt) >>
   2.109 +  foldl (uncurry Ctermfunc.update) Ctermfunc.undefined
   2.110 +
   2.111 +fun parse_cmonomial ctxt =
   2.112 +  rat_int --| $$ "*" -- (parse_monomial ctxt) >> swap ||
   2.113 +  (parse_monomial ctxt) >> (fn m => (m, Rat.one)) ||
   2.114 +  rat_int >> (fn r => (Ctermfunc.undefined, r))
   2.115 +
   2.116 +fun parse_poly ctxt = repeat_sep "+" (parse_cmonomial ctxt) >>
   2.117 +  foldl (uncurry Monomialfunc.update) Monomialfunc.undefined
   2.118 +
   2.119 +(* positivstellensatz parser *)
   2.120 +
   2.121 +val parse_axiom =
   2.122 +  ($$ "A=" |-- int >> Axiom_eq) ||
   2.123 +  ($$ "A<=" |-- int >> Axiom_le) ||
   2.124 +  ($$ "A<" |-- int >> Axiom_lt)
   2.125 +
   2.126 +val parse_rational =
   2.127 +  ($$ "R=" |-- rat_int >> Rational_eq) ||
   2.128 +  ($$ "R<=" |-- rat_int >> Rational_le) ||
   2.129 +  ($$ "R<" |-- rat_int >> Rational_lt)
   2.130 +
   2.131 +fun parse_cert ctxt input =
   2.132 +  let
   2.133 +    val pc = parse_cert ctxt
   2.134 +    val pp = parse_poly ctxt
   2.135 +  in
   2.136 +  (parse_axiom ||
   2.137 +   parse_rational ||
   2.138 +   $$ "[" |-- pp --| $$ "]^2" >> Square ||
   2.139 +   $$ "([" |-- pp --| $$ "]*" -- pc --| $$ ")" >> Eqmul ||
   2.140 +   $$ "(" |-- pc --| $$ "*" -- pc --| $$ ")" >> Product ||
   2.141 +   $$ "(" |-- pc --| $$ "+" -- pc --| $$ ")" >> Sum) input
   2.142 +  end
   2.143 +
   2.144 +fun parse_cert_tree ctxt input =
   2.145 +  let
   2.146 +    val pc = parse_cert ctxt
   2.147 +    val pt = parse_cert_tree ctxt
   2.148 +  in
   2.149 +  ($$ "()" >> K Trivial ||
   2.150 +   $$ "(" |-- pc --| $$ ")" >> Cert ||
   2.151 +   $$ "(" |-- pt --| $$ "&" -- pt --| $$ ")" >> Branch) input
   2.152 +  end
   2.153 +
   2.154 +(* scanner *)
   2.155 +
   2.156 +fun cert_to_pss_tree ctxt str = Symbol.scanner "bad certificate" (parse_cert_tree ctxt)
   2.157 +  (filter_out Symbol.is_blank (Symbol.explode str))
   2.158 +
   2.159 +end
   2.160 +
   2.161 +
     3.1 --- a/src/HOL/Library/Sum_Of_Squares/sos_wrapper.ML	Tue Sep 22 20:25:31 2009 +0200
     3.2 +++ b/src/HOL/Library/Sum_Of_Squares/sos_wrapper.ML	Mon Sep 21 15:05:26 2009 +0200
     3.3 @@ -136,13 +136,42 @@
     3.4      run_solver name (Path.explode cmd) find_failure
     3.5    end
     3.6  
     3.7 +(* certificate output *)
     3.8 +
     3.9 +fun output_line cert = "To repeat this proof with a certifiate use this command:\n" ^
    3.10 +        (Markup.markup Markup.sendback) ("by (sos_cert \"" ^ cert ^ "\")")
    3.11 +
    3.12 +fun print_certs pss_tree =
    3.13 +  let
    3.14 +    val cert = PositivstellensatzTools.pss_tree_to_cert pss_tree
    3.15 +    val str = output_line cert
    3.16 +  in
    3.17 +    Output.writeln str
    3.18 +  end
    3.19 +
    3.20  (* setup tactic *)
    3.21  
    3.22 -fun sos_solver name = (SIMPLE_METHOD' o (Sos.sos_tac (call_solver name))) 
    3.23 +fun parse_cert (ctxt, tk) =
    3.24 +  let
    3.25 +    val (str, tk') = OuterParse.string tk
    3.26 +    val cert = PositivstellensatzTools.cert_to_pss_tree (Context.proof_of ctxt) str
    3.27 +  in
    3.28 +    (cert, (ctxt, tk'))
    3.29 +  end
    3.30 +
    3.31 +fun sos_solver print method = (SIMPLE_METHOD' o (Sos.sos_tac print method)) 
    3.32  
    3.33 -val sos_method = Scan.option (Scan.lift OuterParse.xname) >> sos_solver
    3.34 +val sos_method =
    3.35 +   Scan.lift (Scan.option OuterParse.xname) >> (fn n => Sos.Prover (call_solver n)) >>
    3.36 +   sos_solver print_certs
    3.37  
    3.38 -val setup = Method.setup @{binding sos} sos_method
    3.39 -  "Prove universal problems over the reals using sums of squares"
    3.40 +val sos_cert_method =
    3.41 +  parse_cert >> Sos.Certificate >> sos_solver ignore
    3.42 +
    3.43 +val setup =
    3.44 +     Method.setup @{binding sos} sos_method
    3.45 +     "Prove universal problems over the reals using sums of squares"
    3.46 +  #> Method.setup @{binding sos_cert} sos_cert_method
    3.47 +     "Prove universal problems over the reals using sums of squares with certificates"
    3.48  
    3.49  end
     4.1 --- a/src/HOL/Library/Sum_Of_Squares/sum_of_squares.ML	Tue Sep 22 20:25:31 2009 +0200
     4.2 +++ b/src/HOL/Library/Sum_Of_Squares/sum_of_squares.ML	Mon Sep 21 15:05:26 2009 +0200
     4.3 @@ -8,7 +8,12 @@
     4.4  signature SOS =
     4.5  sig
     4.6  
     4.7 -  val sos_tac : (string -> string) -> Proof.context -> int -> Tactical.tactic
     4.8 +  datatype proof_method =
     4.9 +    Certificate of RealArith.pss_tree
    4.10 +  | Prover of (string -> string)
    4.11 +
    4.12 +  val sos_tac : (RealArith.pss_tree -> unit) ->
    4.13 +    proof_method -> Proof.context -> int -> Tactical.tactic
    4.14  
    4.15    val debugging : bool ref;
    4.16    
    4.17 @@ -18,6 +23,8 @@
    4.18  structure Sos : SOS = 
    4.19  struct
    4.20  
    4.21 +open FuncUtil;
    4.22 +
    4.23  val rat_0 = Rat.zero;
    4.24  val rat_1 = Rat.one;
    4.25  val rat_2 = Rat.two;
    4.26 @@ -59,6 +66,10 @@
    4.27  
    4.28  exception Failure of string;
    4.29  
    4.30 +datatype proof_method =
    4.31 +    Certificate of RealArith.pss_tree
    4.32 +  | Prover of (string -> string)
    4.33 +
    4.34  (* Turn a rational into a decimal string with d sig digits.                  *)
    4.35  
    4.36  local
    4.37 @@ -93,23 +104,11 @@
    4.38  
    4.39  (* The main types.                                                           *)
    4.40  
    4.41 -fun strict_ord ord (x,y) = case ord (x,y) of LESS => LESS | _ => GREATER
    4.42 -
    4.43 -structure Intpairfunc = FuncFun(type key = int*int val ord = prod_ord int_ord int_ord);
    4.44 -
    4.45  type vector = int* Rat.rat Intfunc.T;
    4.46  
    4.47  type matrix = (int*int)*(Rat.rat Intpairfunc.T);
    4.48  
    4.49 -type monomial = int Ctermfunc.T;
    4.50 -
    4.51 -val cterm_ord = (fn (s,t) => TermOrd.fast_term_ord(term_of s, term_of t))
    4.52 - fun monomial_ord (m1,m2) = list_ord (prod_ord cterm_ord int_ord) (Ctermfunc.graph m1, Ctermfunc.graph m2)
    4.53 -structure Monomialfunc = FuncFun(type key = monomial val ord = monomial_ord)
    4.54 -
    4.55 -type poly = Rat.rat Monomialfunc.T;
    4.56 -
    4.57 - fun iszero (k,r) = r =/ rat_0;
    4.58 +fun iszero (k,r) = r =/ rat_0;
    4.59  
    4.60  fun fold_rev2 f l1 l2 b =
    4.61    case (l1,l2) of
    4.62 @@ -346,10 +345,7 @@
    4.63    sort humanorder_varpow (Ctermfunc.graph m2))
    4.64  end;
    4.65  
    4.66 -fun fold1 f l =  case l of
    4.67 -   []     => error "fold1"
    4.68 - | [x]    => x
    4.69 - | (h::t) => f h (fold1 f t);
    4.70 +fun fold1 f = foldr1 (uncurry f) 
    4.71  
    4.72  (* Conversions to strings.                                                   *)
    4.73  
    4.74 @@ -404,7 +400,7 @@
    4.75   else if c =/ rat_1 then string_of_monomial m
    4.76   else Rat.string_of_rat c ^ "*" ^ string_of_monomial m;;
    4.77  
    4.78 -fun string_of_poly (p:poly) =
    4.79 +fun string_of_poly p =
    4.80   if Monomialfunc.is_undefined p then "<<0>>" else
    4.81   let 
    4.82    val cms = sort (fn ((m1,_),(m2,_)) => humanorder_monomial m1  m2) (Monomialfunc.graph p)
    4.83 @@ -481,7 +477,6 @@
    4.84   in fold1 (fn x => fn y => x ^ " " ^ y) strs ^ "\n"
    4.85   end;
    4.86  
    4.87 -fun increasing f ord (x,y) = ord (f x, f y);
    4.88  fun triple_int_ord ((a,b,c),(a',b',c')) = 
    4.89   prod_ord int_ord (prod_ord int_ord int_ord) 
    4.90      ((a,(b,c)),(a',(b',c')));
    4.91 @@ -1080,11 +1075,6 @@
    4.92    fun tryfind f = tryfind_with "tryfind" f
    4.93  end
    4.94  
    4.95 -(*
    4.96 -fun tryfind f [] = error "tryfind"
    4.97 -  | tryfind f (x::xs) = (f x handle ERROR _ => tryfind f xs);
    4.98 -*)
    4.99 -
   4.100  (* Positiv- and Nullstellensatz. Flag "linf" forces a linear representation. *)
   4.101  
   4.102   
   4.103 @@ -1210,58 +1200,14 @@
   4.104  fun deepen f n = 
   4.105    (writeln ("Searching with depth limit " ^ string_of_int n) ; (f n handle Failure s => (writeln ("failed with message: " ^ s) ; deepen f (n+1))))
   4.106  
   4.107 -(* The ordering so we can create canonical HOL polynomials.                  *)
   4.108  
   4.109 -fun dest_monomial mon = sort (increasing fst cterm_ord) (Ctermfunc.graph mon);
   4.110 -
   4.111 -fun monomial_order (m1,m2) =
   4.112 - if Ctermfunc.is_undefined m2 then LESS 
   4.113 - else if Ctermfunc.is_undefined m1 then GREATER 
   4.114 - else
   4.115 -  let val mon1 = dest_monomial m1 
   4.116 -      val mon2 = dest_monomial m2
   4.117 -      val deg1 = fold (curry op + o snd) mon1 0
   4.118 -      val deg2 = fold (curry op + o snd) mon2 0 
   4.119 -  in if deg1 < deg2 then GREATER else if deg1 > deg2 then LESS
   4.120 -     else list_ord (prod_ord cterm_ord int_ord) (mon1,mon2)
   4.121 -  end;
   4.122 -
   4.123 -fun dest_poly p =
   4.124 -  map (fn (m,c) => (c,dest_monomial m))
   4.125 -      (sort (prod_ord monomial_order (K EQUAL)) (Monomialfunc.graph p));
   4.126 -
   4.127 -(* Map back polynomials and their composites to HOL.                         *)
   4.128 +(* Map back polynomials and their composites to a positivstellensatz.        *)
   4.129  
   4.130  local
   4.131   open Thm Numeral RealArith
   4.132  in
   4.133  
   4.134 -fun cterm_of_varpow x k = if k = 1 then x else capply (capply @{cterm "op ^ :: real => _"} x) 
   4.135 -  (mk_cnumber @{ctyp nat} k)
   4.136 -
   4.137 -fun cterm_of_monomial m = 
   4.138 - if Ctermfunc.is_undefined m then @{cterm "1::real"} 
   4.139 - else 
   4.140 -  let 
   4.141 -   val m' = dest_monomial m
   4.142 -   val vps = fold_rev (fn (x,k) => cons (cterm_of_varpow x k)) m' [] 
   4.143 -  in fold1 (fn s => fn t => capply (capply @{cterm "op * :: real => _"} s) t) vps
   4.144 -  end
   4.145 -
   4.146 -fun cterm_of_cmonomial (m,c) = if Ctermfunc.is_undefined m then cterm_of_rat c
   4.147 -    else if c = Rat.one then cterm_of_monomial m
   4.148 -    else capply (capply @{cterm "op *::real => _"} (cterm_of_rat c)) (cterm_of_monomial m);
   4.149 -
   4.150 -fun cterm_of_poly p = 
   4.151 - if Monomialfunc.is_undefined p then @{cterm "0::real"} 
   4.152 - else
   4.153 -  let 
   4.154 -   val cms = map cterm_of_cmonomial
   4.155 -     (sort (prod_ord monomial_order (K EQUAL)) (Monomialfunc.graph p))
   4.156 -  in fold1 (fn t1 => fn t2 => capply(capply @{cterm "op + :: real => _"} t1) t2) cms
   4.157 -  end;
   4.158 -
   4.159 -fun cterm_of_sqterm (c,p) = Product(Rational_lt c,Square(cterm_of_poly p));
   4.160 +fun cterm_of_sqterm (c,p) = Product(Rational_lt c,Square p);
   4.161  
   4.162  fun cterm_of_sos (pr,sqs) = if null sqs then pr
   4.163    else Product(pr,fold1 (fn a => fn b => Sum(a,b)) (map cterm_of_sqterm sqs));
   4.164 @@ -1275,14 +1221,14 @@
   4.165    fun simple_cterm_ord t u = TermOrd.fast_term_ord (term_of t, term_of u) = LESS
   4.166  in
   4.167    (* FIXME: Replace tryfind by get_first !! *)
   4.168 -fun real_nonlinear_prover prover ctxt =
   4.169 +fun real_nonlinear_prover proof_method ctxt =
   4.170   let 
   4.171    val {add,mul,neg,pow,sub,main} =  Normalizer.semiring_normalizers_ord_wrapper ctxt
   4.172        (valOf (NormalizerData.match ctxt @{cterm "(0::real) + 1"})) 
   4.173       simple_cterm_ord
   4.174    val (real_poly_add_conv,real_poly_mul_conv,real_poly_neg_conv,
   4.175         real_poly_pow_conv,real_poly_sub_conv,real_poly_conv) = (add,mul,neg,pow,sub,main)
   4.176 -  fun mainf  translator (eqs,les,lts) = 
   4.177 +  fun mainf cert_choice translator (eqs,les,lts) = 
   4.178    let 
   4.179     val eq0 = map (poly_of_term o dest_arg1 o concl) eqs
   4.180     val le0 = map (poly_of_term o dest_arg o concl) les
   4.181 @@ -1303,33 +1249,49 @@
   4.182                       else raise Failure "trivial_axiom: Not a trivial axiom"
   4.183       | _ => error "trivial_axiom: Not a trivial axiom"
   4.184     in 
   4.185 -  ((let val th = tryfind trivial_axiom (keq @ klep @ kltp)
   4.186 -   in fconv_rule (arg_conv (arg1_conv real_poly_conv) then_conv field_comp_conv) th end)
   4.187 -   handle Failure _ => (
   4.188 -    let 
   4.189 -     val pol = fold_rev poly_mul (map fst ltp) (poly_const Rat.one)
   4.190 -     val leq = lep @ ltp
   4.191 -     fun tryall d =
   4.192 -      let val e = multidegree pol
   4.193 -          val k = if e = 0 then 0 else d div e
   4.194 -          val eq' = map fst eq 
   4.195 -      in tryfind (fn i => (d,i,real_positivnullstellensatz_general prover false d eq' leq
   4.196 -                            (poly_neg(poly_pow pol i))))
   4.197 -              (0 upto k)
   4.198 -      end
   4.199 -    val (d,i,(cert_ideal,cert_cone)) = deepen tryall 0
   4.200 -    val proofs_ideal =
   4.201 -      map2 (fn q => fn (p,ax) => Eqmul(cterm_of_poly q,ax)) cert_ideal eq
   4.202 -    val proofs_cone = map cterm_of_sos cert_cone
   4.203 -    val proof_ne = if null ltp then Rational_lt Rat.one else
   4.204 -      let val p = fold1 (fn s => fn t => Product(s,t)) (map snd ltp) 
   4.205 -      in  funpow i (fn q => Product(p,q)) (Rational_lt Rat.one)
   4.206 -      end
   4.207 -    val proof = fold1 (fn s => fn t => Sum(s,t))
   4.208 -                           (proof_ne :: proofs_ideal @ proofs_cone) 
   4.209 -    in writeln "Translating proof certificate to HOL";
   4.210 -       translator (eqs,les,lts) proof
   4.211 -    end))
   4.212 +  (let val th = tryfind trivial_axiom (keq @ klep @ kltp)
   4.213 +   in
   4.214 +    (fconv_rule (arg_conv (arg1_conv real_poly_conv) then_conv field_comp_conv) th, Trivial)
   4.215 +   end)
   4.216 +   handle Failure _ => 
   4.217 +     (let val proof =
   4.218 +       (case proof_method of Certificate certs =>
   4.219 +         (* choose certificate *)
   4.220 +         let
   4.221 +           fun chose_cert [] (Cert c) = c
   4.222 +             | chose_cert (Left::s) (Branch (l, _)) = chose_cert s l
   4.223 +             | chose_cert (Right::s) (Branch (_, r)) = chose_cert s r
   4.224 +             | chose_cert _ _ = error "certificate tree in invalid form"
   4.225 +         in
   4.226 +           chose_cert cert_choice certs
   4.227 +         end
   4.228 +       | Prover prover =>
   4.229 +         (* call prover *)
   4.230 +         let 
   4.231 +          val pol = fold_rev poly_mul (map fst ltp) (poly_const Rat.one)
   4.232 +          val leq = lep @ ltp
   4.233 +          fun tryall d =
   4.234 +           let val e = multidegree pol
   4.235 +               val k = if e = 0 then 0 else d div e
   4.236 +               val eq' = map fst eq 
   4.237 +           in tryfind (fn i => (d,i,real_positivnullstellensatz_general prover false d eq' leq
   4.238 +                                 (poly_neg(poly_pow pol i))))
   4.239 +                   (0 upto k)
   4.240 +           end
   4.241 +         val (d,i,(cert_ideal,cert_cone)) = deepen tryall 0
   4.242 +         val proofs_ideal =
   4.243 +           map2 (fn q => fn (p,ax) => Eqmul(q,ax)) cert_ideal eq
   4.244 +         val proofs_cone = map cterm_of_sos cert_cone
   4.245 +         val proof_ne = if null ltp then Rational_lt Rat.one else
   4.246 +           let val p = fold1 (fn s => fn t => Product(s,t)) (map snd ltp) 
   4.247 +           in  funpow i (fn q => Product(p,q)) (Rational_lt Rat.one)
   4.248 +           end
   4.249 +         in 
   4.250 +           fold1 (fn s => fn t => Sum(s,t)) (proof_ne :: proofs_ideal @ proofs_cone) 
   4.251 +         end)
   4.252 +     in
   4.253 +        (translator (eqs,les,lts) proof, Cert proof)
   4.254 +     end)
   4.255     end
   4.256   in mainf end
   4.257  end
   4.258 @@ -1396,7 +1358,7 @@
   4.259           orelse g aconvc @{cterm "op < :: real => _"} 
   4.260         then arg_conv cv ct else arg1_conv cv ct
   4.261      end
   4.262 -  fun mainf translator =
   4.263 +  fun mainf cert_choice translator =
   4.264     let 
   4.265      fun substfirst(eqs,les,lts) =
   4.266        ((let 
   4.267 @@ -1407,7 +1369,7 @@
   4.268                                     aconvc @{cterm "0::real"}) (map modify eqs),
   4.269                                     map modify les,map modify lts)
   4.270         end)
   4.271 -       handle Failure  _ => real_nonlinear_prover prover ctxt translator (rev eqs, rev les, rev lts))
   4.272 +       handle Failure  _ => real_nonlinear_prover prover ctxt cert_choice translator (rev eqs, rev les, rev lts))
   4.273      in substfirst
   4.274     end
   4.275  
   4.276 @@ -1417,7 +1379,8 @@
   4.277  
   4.278  (* Overall function. *)
   4.279  
   4.280 -fun real_sos prover ctxt t = gen_prover_real_arith ctxt (real_nonlinear_subst_prover prover ctxt) t;
   4.281 +fun real_sos prover ctxt =
   4.282 +  gen_prover_real_arith ctxt (real_nonlinear_subst_prover prover ctxt)
   4.283  end;
   4.284  
   4.285  (* A tactic *)
   4.286 @@ -1429,8 +1392,6 @@
   4.287     end
   4.288  | _ => ([],ct)
   4.289  
   4.290 -fun core_sos_conv prover ctxt t = Drule.arg_cong_rule @{cterm Trueprop} (real_sos prover ctxt (Thm.dest_arg t) RS @{thm Eq_TrueI})
   4.291 -
   4.292  val known_sos_constants = 
   4.293    [@{term "op ==>"}, @{term "Trueprop"}, 
   4.294     @{term "op -->"}, @{term "op &"}, @{term "op |"}, 
   4.295 @@ -1458,17 +1419,19 @@
   4.296    val _ = if exists (fn ((_,T)) => not (T = @{typ "real"})) fs 
   4.297            then error "SOS: not sos. Variables with type not real" else ()
   4.298    val vs = Term.add_vars t []
   4.299 -  val _ = if exists (fn ((_,T)) => not (T = @{typ "real"})) fs 
   4.300 +  val _ = if exists (fn ((_,T)) => not (T = @{typ "real"})) vs 
   4.301            then error "SOS: not sos. Variables with type not real" else ()
   4.302    val ukcs = subtract (fn (t,p) => Const p aconv t) kcts (Term.add_consts t [])
   4.303    val _ = if  null ukcs then () 
   4.304                else error ("SOSO: Unknown constants in Subgoal:" ^ commas (map fst ukcs))
   4.305  in () end
   4.306  
   4.307 -fun core_sos_tac prover ctxt = CSUBGOAL (fn (ct, i) => 
   4.308 +fun core_sos_tac print_certs prover ctxt = CSUBGOAL (fn (ct, i) => 
   4.309    let val _ = check_sos known_sos_constants ct
   4.310        val (avs, p) = strip_all ct
   4.311 -      val th = standard (fold_rev forall_intr avs (real_sos prover ctxt (Thm.dest_arg p)))
   4.312 +      val (ths, certificates) = real_sos prover ctxt (Thm.dest_arg p)
   4.313 +      val th = standard (fold_rev forall_intr avs ths)
   4.314 +      val _ = print_certs certificates
   4.315    in rtac th i end);
   4.316  
   4.317  fun default_SOME f NONE v = SOME v
   4.318 @@ -1506,7 +1469,7 @@
   4.319  
   4.320  fun elim_denom_tac ctxt i = REPEAT (elim_one_denom_tac ctxt i);
   4.321  
   4.322 -fun sos_tac prover ctxt = ObjectLogic.full_atomize_tac THEN' elim_denom_tac ctxt THEN' core_sos_tac prover ctxt
   4.323 +fun sos_tac print_certs prover ctxt = ObjectLogic.full_atomize_tac THEN' elim_denom_tac ctxt THEN' core_sos_tac print_certs prover ctxt
   4.324  
   4.325  
   4.326  end;
     5.1 --- a/src/HOL/Library/normarith.ML	Tue Sep 22 20:25:31 2009 +0200
     5.2 +++ b/src/HOL/Library/normarith.ML	Mon Sep 21 15:05:26 2009 +0200
     5.3 @@ -15,7 +15,7 @@
     5.4  structure NormArith : NORM_ARITH = 
     5.5  struct
     5.6  
     5.7 - open Conv Thm;
     5.8 + open Conv Thm FuncUtil;
     5.9   val bool_eq = op = : bool *bool -> bool
    5.10    fun dest_ratconst t = case term_of t of
    5.11     Const(@{const_name divide}, _)$a$b => Rat.rat_of_quotient(HOLogic.dest_number a |> snd, HOLogic.dest_number b |> snd)
    5.12 @@ -330,13 +330,13 @@
    5.13     val zerodests = filter
    5.14          (fn t => null (Ctermfunc.dom (vector_lincomb t))) (map snd rawdests)
    5.15  
    5.16 -  in RealArith.real_linear_prover translator
    5.17 +  in fst (RealArith.real_linear_prover translator
    5.18          (map (fn t => instantiate ([(tv_n, ctyp_of_term t)],[]) pth_zero)
    5.19              zerodests,
    5.20          map (fconv_rule (try_conv (More_Conv.top_sweep_conv (K norm_canon_conv) ctxt) then_conv
    5.21                         arg_conv (arg_conv real_poly_conv))) ges',
    5.22          map (fconv_rule (try_conv (More_Conv.top_sweep_conv (K norm_canon_conv) ctxt) then_conv 
    5.23 -                       arg_conv (arg_conv real_poly_conv))) gts)
    5.24 +                       arg_conv (arg_conv real_poly_conv))) gts))
    5.25    end
    5.26  in val real_vector_combo_prover = real_vector_combo_prover
    5.27  end;
    5.28 @@ -389,9 +389,9 @@
    5.29     val (th1,th2) = conj_pair(rawrule th)
    5.30    in th1::fconv_rule (arg_conv (arg_conv real_poly_neg_conv)) th2::acc
    5.31    end
    5.32 -in fun real_vector_prover ctxt translator (eqs,ges,gts) =
    5.33 -     real_vector_ineq_prover ctxt translator
    5.34 -         (fold_rev (splitequation ctxt) eqs ges,gts)
    5.35 +in fun real_vector_prover ctxt _ translator (eqs,ges,gts) =
    5.36 +     (real_vector_ineq_prover ctxt translator
    5.37 +         (fold_rev (splitequation ctxt) eqs ges,gts), RealArith.Trivial)
    5.38  end;
    5.39  
    5.40    fun init_conv ctxt = 
    5.41 @@ -400,7 +400,7 @@
    5.42     then_conv field_comp_conv 
    5.43     then_conv nnf_conv
    5.44  
    5.45 - fun pure ctxt = RealArith.gen_prover_real_arith ctxt (real_vector_prover ctxt);
    5.46 + fun pure ctxt = fst o RealArith.gen_prover_real_arith ctxt (real_vector_prover ctxt);
    5.47   fun norm_arith ctxt ct = 
    5.48    let 
    5.49     val ctxt' = Variable.declare_term (term_of ct) ctxt
     6.1 --- a/src/HOL/Library/positivstellensatz.ML	Tue Sep 22 20:25:31 2009 +0200
     6.2 +++ b/src/HOL/Library/positivstellensatz.ML	Mon Sep 21 15:05:26 2009 +0200
     6.3 @@ -1,10 +1,11 @@
     6.4 -(* Title:      Library/positivstellensatz
     6.5 +(* Title:      Library/Sum_Of_Squares/positivstellensatz
     6.6     Author:     Amine Chaieb, University of Cambridge
     6.7     Description: A generic arithmetic prover based on Positivstellensatz certificates --- 
     6.8      also implements Fourrier-Motzkin elimination as a special case Fourrier-Motzkin elimination.
     6.9  *)
    6.10  
    6.11  (* A functor for finite mappings based on Tables *)
    6.12 +
    6.13  signature FUNC = 
    6.14  sig
    6.15   type 'a T
    6.16 @@ -75,27 +76,54 @@
    6.17  end
    6.18  end;
    6.19  
    6.20 +(* Some standard functors and utility functions for them *)
    6.21 +
    6.22 +structure FuncUtil =
    6.23 +struct
    6.24 +
    6.25 +fun increasing f ord (x,y) = ord (f x, f y);
    6.26 +
    6.27  structure Intfunc = FuncFun(type key = int val ord = int_ord);
    6.28 +structure Ratfunc = FuncFun(type key = Rat.rat val ord = Rat.ord);
    6.29 +structure Intpairfunc = FuncFun(type key = int*int val ord = prod_ord int_ord int_ord);
    6.30  structure Symfunc = FuncFun(type key = string val ord = fast_string_ord);
    6.31  structure Termfunc = FuncFun(type key = term val ord = TermOrd.fast_term_ord);
    6.32 -structure Ctermfunc = FuncFun(type key = cterm val ord = (fn (s,t) => TermOrd.fast_term_ord(term_of s, term_of t)));
    6.33 +
    6.34 +val cterm_ord = (fn (s,t) => TermOrd.fast_term_ord(term_of s, term_of t))
    6.35 +
    6.36 +structure Ctermfunc = FuncFun(type key = cterm val ord = cterm_ord);
    6.37 +
    6.38 +type monomial = int Ctermfunc.T;
    6.39  
    6.40 -structure Ratfunc = FuncFun(type key = Rat.rat val ord = Rat.ord);
    6.41 +fun monomial_ord (m1,m2) = list_ord (prod_ord cterm_ord int_ord) (Ctermfunc.graph m1, Ctermfunc.graph m2)
    6.42 +
    6.43 +structure Monomialfunc = FuncFun(type key = monomial val ord = monomial_ord)
    6.44  
    6.45 +type poly = Rat.rat Monomialfunc.T;
    6.46 +
    6.47 +(* The ordering so we can create canonical HOL polynomials.                  *)
    6.48  
    6.49 -    (* Some useful derived rules *)
    6.50 -fun deduct_antisym_rule tha thb = 
    6.51 -    equal_intr (implies_intr (cprop_of thb) tha) 
    6.52 -     (implies_intr (cprop_of tha) thb);
    6.53 +fun dest_monomial mon = sort (increasing fst cterm_ord) (Ctermfunc.graph mon);
    6.54  
    6.55 -fun prove_hyp tha thb = 
    6.56 -  if exists (curry op aconv (concl_of tha)) (#hyps (rep_thm thb)) 
    6.57 -  then equal_elim (symmetric (deduct_antisym_rule tha thb)) tha else thb;
    6.58 +fun monomial_order (m1,m2) =
    6.59 + if Ctermfunc.is_undefined m2 then LESS 
    6.60 + else if Ctermfunc.is_undefined m1 then GREATER 
    6.61 + else
    6.62 +  let val mon1 = dest_monomial m1 
    6.63 +      val mon2 = dest_monomial m2
    6.64 +      val deg1 = fold (curry op + o snd) mon1 0
    6.65 +      val deg2 = fold (curry op + o snd) mon2 0 
    6.66 +  in if deg1 < deg2 then GREATER else if deg1 > deg2 then LESS
    6.67 +     else list_ord (prod_ord cterm_ord int_ord) (mon1,mon2)
    6.68 +  end;
    6.69  
    6.70 +end
    6.71  
    6.72 +(* positivstellensatz datatype and prover generation *)
    6.73  
    6.74  signature REAL_ARITH = 
    6.75  sig
    6.76 +  
    6.77    datatype positivstellensatz =
    6.78     Axiom_eq of int
    6.79   | Axiom_le of int
    6.80 @@ -103,34 +131,41 @@
    6.81   | Rational_eq of Rat.rat
    6.82   | Rational_le of Rat.rat
    6.83   | Rational_lt of Rat.rat
    6.84 - | Square of cterm
    6.85 - | Eqmul of cterm * positivstellensatz
    6.86 + | Square of FuncUtil.poly
    6.87 + | Eqmul of FuncUtil.poly * positivstellensatz
    6.88   | Sum of positivstellensatz * positivstellensatz
    6.89   | Product of positivstellensatz * positivstellensatz;
    6.90  
    6.91 +datatype pss_tree = Trivial | Cert of positivstellensatz | Branch of pss_tree * pss_tree
    6.92 +
    6.93 +datatype tree_choice = Left | Right
    6.94 +
    6.95 +type prover = tree_choice list -> 
    6.96 +  (thm list * thm list * thm list -> positivstellensatz -> thm) ->
    6.97 +  thm list * thm list * thm list -> thm * pss_tree
    6.98 +type cert_conv = cterm -> thm * pss_tree
    6.99 +
   6.100  val gen_gen_real_arith :
   6.101 -  Proof.context -> (Rat.rat -> Thm.cterm) * conv * conv * conv * 
   6.102 -   conv * conv * conv * conv * conv * conv * 
   6.103 -    ( (thm list * thm list * thm list -> positivstellensatz -> thm) ->
   6.104 -        thm list * thm list * thm list -> thm) -> conv
   6.105 -val real_linear_prover : 
   6.106 -  (thm list * thm list * thm list -> positivstellensatz -> thm) ->
   6.107 -   thm list * thm list * thm list -> thm
   6.108 +  Proof.context -> (Rat.rat -> Thm.cterm) * conv * conv * conv *
   6.109 +   conv * conv * conv * conv * conv * conv * prover -> cert_conv
   6.110 +val real_linear_prover : (thm list * thm list * thm list -> positivstellensatz -> thm) ->
   6.111 +  thm list * thm list * thm list -> thm * pss_tree
   6.112  
   6.113  val gen_real_arith : Proof.context ->
   6.114 -   (Rat.rat -> cterm) * conv * conv * conv * conv * conv * conv * conv *
   6.115 -   ( (thm list * thm list * thm list -> positivstellensatz -> thm) ->
   6.116 -       thm list * thm list * thm list -> thm) -> conv
   6.117 -val gen_prover_real_arith : Proof.context ->
   6.118 -   ((thm list * thm list * thm list -> positivstellensatz -> thm) ->
   6.119 -     thm list * thm list * thm list -> thm) -> conv
   6.120 -val real_arith : Proof.context -> conv
   6.121 +  (Rat.rat -> Thm.cterm) * conv * conv * conv * conv * conv * conv * conv * prover -> cert_conv
   6.122 +
   6.123 +val gen_prover_real_arith : Proof.context -> prover -> cert_conv
   6.124 +
   6.125 +val is_ratconst : Thm.cterm -> bool
   6.126 +val dest_ratconst : Thm.cterm -> Rat.rat
   6.127 +val cterm_of_rat : Rat.rat -> Thm.cterm
   6.128 +
   6.129  end
   6.130  
   6.131 -structure RealArith (* : REAL_ARITH *)=
   6.132 +structure RealArith : REAL_ARITH =
   6.133  struct
   6.134  
   6.135 - open Conv Thm;;
   6.136 + open Conv Thm FuncUtil;;
   6.137  (* ------------------------------------------------------------------------- *)
   6.138  (* Data structure for Positivstellensatz refutations.                        *)
   6.139  (* ------------------------------------------------------------------------- *)
   6.140 @@ -142,12 +177,18 @@
   6.141   | Rational_eq of Rat.rat
   6.142   | Rational_le of Rat.rat
   6.143   | Rational_lt of Rat.rat
   6.144 - | Square of cterm
   6.145 - | Eqmul of cterm * positivstellensatz
   6.146 + | Square of FuncUtil.poly
   6.147 + | Eqmul of FuncUtil.poly * positivstellensatz
   6.148   | Sum of positivstellensatz * positivstellensatz
   6.149   | Product of positivstellensatz * positivstellensatz;
   6.150           (* Theorems used in the procedure *)
   6.151  
   6.152 +datatype pss_tree = Trivial | Cert of positivstellensatz | Branch of pss_tree * pss_tree
   6.153 +datatype tree_choice = Left | Right
   6.154 +type prover = tree_choice list -> 
   6.155 +  (thm list * thm list * thm list -> positivstellensatz -> thm) ->
   6.156 +  thm list * thm list * thm list -> thm * pss_tree
   6.157 +type cert_conv = cterm -> thm * pss_tree
   6.158  
   6.159  val my_eqs = ref ([] : thm list);
   6.160  val my_les = ref ([] : thm list);
   6.161 @@ -164,6 +205,16 @@
   6.162  val my_poly_add_conv = ref no_conv;
   6.163  val my_poly_mul_conv = ref no_conv;
   6.164  
   6.165 +
   6.166 +    (* Some useful derived rules *)
   6.167 +fun deduct_antisym_rule tha thb = 
   6.168 +    equal_intr (implies_intr (cprop_of thb) tha) 
   6.169 +     (implies_intr (cprop_of tha) thb);
   6.170 +
   6.171 +fun prove_hyp tha thb = 
   6.172 +  if exists (curry op aconv (concl_of tha)) (#hyps (rep_thm thb)) 
   6.173 +  then equal_elim (symmetric (deduct_antisym_rule tha thb)) tha else thb;
   6.174 +
   6.175  fun conjunctions th = case try Conjunction.elim th of
   6.176     SOME (th1,th2) => (conjunctions th1) @ conjunctions th2
   6.177   | NONE => [th];
   6.178 @@ -292,7 +343,6 @@
   6.179   | Abs (_,_,t') => find_cterm p (Thm.dest_abs NONE t |> snd)
   6.180   | _ => raise CTERM ("find_cterm",[t]);
   6.181  
   6.182 -
   6.183      (* Some conversions-related stuff which has been forbidden entrance into Pure/conv.ML*)
   6.184  fun instantiate_cterm' ty tms = Drule.cterm_rule (Drule.instantiate' ty tms)
   6.185  fun is_comb t = case (term_of t) of _$_ => true | _ => false;
   6.186 @@ -300,6 +350,39 @@
   6.187  fun is_binop ct ct' = ct aconvc (Thm.dest_fun (Thm.dest_fun ct'))
   6.188    handle CTERM _ => false;
   6.189  
   6.190 +
   6.191 +(* Map back polynomials to HOL.                         *)
   6.192 +
   6.193 +local
   6.194 + open Thm Numeral
   6.195 +in
   6.196 +
   6.197 +fun cterm_of_varpow x k = if k = 1 then x else capply (capply @{cterm "op ^ :: real => _"} x) 
   6.198 +  (mk_cnumber @{ctyp nat} k)
   6.199 +
   6.200 +fun cterm_of_monomial m = 
   6.201 + if Ctermfunc.is_undefined m then @{cterm "1::real"} 
   6.202 + else 
   6.203 +  let 
   6.204 +   val m' = dest_monomial m
   6.205 +   val vps = fold_rev (fn (x,k) => cons (cterm_of_varpow x k)) m' [] 
   6.206 +  in foldr1 (fn (s, t) => capply (capply @{cterm "op * :: real => _"} s) t) vps
   6.207 +  end
   6.208 +
   6.209 +fun cterm_of_cmonomial (m,c) = if Ctermfunc.is_undefined m then cterm_of_rat c
   6.210 +    else if c = Rat.one then cterm_of_monomial m
   6.211 +    else capply (capply @{cterm "op *::real => _"} (cterm_of_rat c)) (cterm_of_monomial m);
   6.212 +
   6.213 +fun cterm_of_poly p = 
   6.214 + if Monomialfunc.is_undefined p then @{cterm "0::real"} 
   6.215 + else
   6.216 +  let 
   6.217 +   val cms = map cterm_of_cmonomial
   6.218 +     (sort (prod_ord monomial_order (K EQUAL)) (Monomialfunc.graph p))
   6.219 +  in foldr1 (fn (t1, t2) => capply(capply @{cterm "op + :: real => _"} t1) t2) cms
   6.220 +  end;
   6.221 +
   6.222 +end;
   6.223      (* A general real arithmetic prover *)
   6.224  
   6.225  fun gen_gen_real_arith ctxt (mk_numeric,
   6.226 @@ -368,8 +451,8 @@
   6.227        | Rational_lt x => eqT_elim(numeric_gt_conv(capply @{cterm Trueprop} 
   6.228                        (capply (capply @{cterm "op <::real => _"} @{cterm "0::real"})
   6.229                          (mk_numeric x))))
   6.230 -      | Square t => square_rule t
   6.231 -      | Eqmul(t,p) => emul_rule t (translate p)
   6.232 +      | Square pt => square_rule (cterm_of_poly pt)
   6.233 +      | Eqmul(pt,p) => emul_rule (cterm_of_poly pt) (translate p)
   6.234        | Sum(p1,p2) => add_rule (translate p1) (translate p2)
   6.235        | Product(p1,p2) => mul_rule (translate p1) (translate p2)
   6.236     in fconv_rule (first_conv [numeric_ge_conv, numeric_gt_conv, numeric_eq_conv, all_conv]) 
   6.237 @@ -394,13 +477,13 @@
   6.238         val _ = if c aconvc (concl th2) then () else error "disj_cases : conclusions not alpha convertible"
   6.239     in implies_elim (implies_elim (implies_elim (instantiate' [] (map SOME [p,q,c]) @{thm disjE}) th) (implies_intr (capply @{cterm Trueprop} p) th1)) (implies_intr (capply @{cterm Trueprop} q) th2)
   6.240     end
   6.241 - fun overall dun ths = case ths of
   6.242 + fun overall cert_choice dun ths = case ths of
   6.243    [] =>
   6.244     let 
   6.245      val (eq,ne) = List.partition (is_req o concl) dun
   6.246       val (le,nl) = List.partition (is_ge o concl) ne
   6.247       val lt = filter (is_gt o concl) nl 
   6.248 -    in prover hol_of_positivstellensatz (eq,le,lt) end
   6.249 +    in prover (rev cert_choice) hol_of_positivstellensatz (eq,le,lt) end
   6.250   | th::oths =>
   6.251     let 
   6.252      val ct = concl th 
   6.253 @@ -408,13 +491,13 @@
   6.254      if is_conj ct  then
   6.255       let 
   6.256        val (th1,th2) = conj_pair th in
   6.257 -      overall dun (th1::th2::oths) end
   6.258 +      overall cert_choice dun (th1::th2::oths) end
   6.259      else if is_disj ct then
   6.260        let 
   6.261 -       val th1 = overall dun (assume (capply @{cterm Trueprop} (dest_arg1 ct))::oths)
   6.262 -       val th2 = overall dun (assume (capply @{cterm Trueprop} (dest_arg ct))::oths)
   6.263 -      in disj_cases th th1 th2 end
   6.264 -   else overall (th::dun) oths
   6.265 +       val (th1, cert1) = overall (Left::cert_choice) dun (assume (capply @{cterm Trueprop} (dest_arg1 ct))::oths)
   6.266 +       val (th2, cert2) = overall (Right::cert_choice) dun (assume (capply @{cterm Trueprop} (dest_arg ct))::oths)
   6.267 +      in (disj_cases th th1 th2, Branch (cert1, cert2)) end
   6.268 +   else overall cert_choice (th::dun) oths
   6.269    end
   6.270    fun dest_binary b ct = if is_binop b ct then dest_binop ct 
   6.271                           else raise CTERM ("dest_binary",[b,ct])
   6.272 @@ -496,16 +579,16 @@
   6.273    val nct = capply @{cterm Trueprop} (capply @{cterm "Not"} ct)
   6.274    val th0 = (init_conv then_conv arg_conv nnf_norm_conv') nct
   6.275    val tm0 = dest_arg (rhs_of th0)
   6.276 -  val th = if tm0 aconvc @{cterm False} then equal_implies_1_rule th0 else
   6.277 +  val (th, certificates) = if tm0 aconvc @{cterm False} then (equal_implies_1_rule th0, Trivial) else
   6.278     let 
   6.279      val (evs,bod) = strip_exists tm0
   6.280      val (avs,ibod) = strip_forall bod
   6.281      val th1 = Drule.arg_cong_rule @{cterm Trueprop} (fold mk_forall avs (absremover ibod))
   6.282 -    val th2 = overall [] [specl avs (assume (rhs_of th1))]
   6.283 +    val (th2, certs) = overall [] [] [specl avs (assume (rhs_of th1))]
   6.284      val th3 = fold simple_choose evs (prove_hyp (equal_elim th1 (assume (capply @{cterm Trueprop} bod))) th2)
   6.285 -   in  Drule.implies_intr_hyps (prove_hyp (equal_elim th0 (assume nct)) th3)
   6.286 +   in (Drule.implies_intr_hyps (prove_hyp (equal_elim th0 (assume nct)) th3), certs)
   6.287     end
   6.288 -  in implies_elim (instantiate' [] [SOME ct] pth_final) th
   6.289 +  in (implies_elim (instantiate' [] [SOME ct] pth_final) th, certificates)
   6.290   end
   6.291  in f
   6.292  end;
   6.293 @@ -580,7 +663,7 @@
   6.294           val k = (Rat.neg d) */ Rat.abs c // c
   6.295           val e' = linear_cmul k e
   6.296           val t' = linear_cmul (Rat.abs c) t
   6.297 -         val p' = Eqmul(cterm_of_rat k,p)
   6.298 +         val p' = Eqmul(Monomialfunc.onefunc (Ctermfunc.undefined, k),p)
   6.299           val q' = Product(Rational_lt(Rat.abs c),q) 
   6.300          in (linear_add e' t',Sum(p',q')) 
   6.301          end 
   6.302 @@ -632,7 +715,7 @@
   6.303    val le_pols' = le_pols @ map (fn v => Ctermfunc.onefunc (v,Rat.one)) aliens
   6.304    val (_,proof) = linear_prover (eq_pols,le_pols',lt_pols)
   6.305    val le' = le @ map (fn a => instantiate' [] [SOME (dest_arg a)] @{thm real_of_nat_ge_zero}) aliens 
   6.306 - in (translator (eq,le',lt) proof) : thm
   6.307 + in ((translator (eq,le',lt) proof), Trivial)
   6.308   end
   6.309  end;
   6.310  
   6.311 @@ -698,5 +781,4 @@
   6.312      main,neg,add,mul, prover)
   6.313  end;
   6.314  
   6.315 -fun real_arith ctxt = gen_prover_real_arith ctxt real_linear_prover;
   6.316  end