# HG changeset patch
# User Philipp Meyer
# Date 1253538326 -7200
# Node ID 1cc5b24f5a01ae68c2e806206235eda6b391f1ce
# Parent  e4511a1b4c1b2014e1ff1ea84cd66a0765429b06
sos method generates and uses proof certificates

diff -r e4511a1b4c1b -r 1cc5b24f5a01 src/HOL/Library/Sum_Of_Squares.thy
--- a/src/HOL/Library/Sum_Of_Squares.thy	Tue Sep 22 20:25:31 2009 +0200
+++ b/src/HOL/Library/Sum_Of_Squares.thy	Mon Sep 21 15:05:26 2009 +0200
@@ -10,6 +10,7 @@
 uses
   ("positivstellensatz.ML")  (* duplicate use!? -- cf. Euclidian_Space.thy *)
   ("Sum_Of_Squares/sum_of_squares.ML")
+  ("Sum_Of_Squares/positivstellensatz_tools.ML")
   ("Sum_Of_Squares/sos_wrapper.ML")
 begin
 
@@ -22,112 +23,142 @@
   of a minute for one sos call, because sos calls CSDP repeatedly.  If
   you install CSDP locally, sos calls typically takes only a few
   seconds.
+  sos generates a certificate which can be used to repeat the proof
+  without calling an external prover.
 *}
 
 text {* setup sos tactic *}
 
 use "positivstellensatz.ML"
+use "Sum_Of_Squares/positivstellensatz_tools.ML"
 use "Sum_Of_Squares/sum_of_squares.ML"
 use "Sum_Of_Squares/sos_wrapper.ML"
 
 setup SosWrapper.setup
 
-text {* Tests -- commented since they work only when csdp is installed
-  or take too long with remote csdps *}
+text {* Tests *}
+
+lemma "(3::real) * x + 7 * a < 4 & 3 < 2 * x \<Longrightarrow> a < 0"
+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))))))")
+
+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)"
+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))))))))))")
 
-(*
-lemma "(3::real) * x + 7 * a < 4 & 3 < 2 * x \<Longrightarrow> a < 0" by sos
+lemma "(3::real) * x + 7 * a < 4 & 3 < 2 * x --> a < 0"
+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))))))")
+
+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"
+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)))))))")
 
-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)" by sos
+lemma "(0::real) <= x & 0 <= y & 0 <= z & x + y + z <= 3 --> x * y + x * z + y * z >= 3 * x * y * z"
+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))))))))))")
+
+lemma "((x::real)^2 + y^2 + z^2 = 1) --> (x + y + z)^2 <= 3"
+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))))))")
 
-lemma "(3::real) * x + 7 * a < 4 & 3 < 2 * x --> a < 0" by sos
+lemma "(w^2 + x^2 + y^2 + z^2 = 1) --> (w + x + y + z)^2 <= (4::real)"
+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)))))))")
 
-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" by sos
+lemma "(x::real) >= 1 & y >= 1 --> x * y >= x + y - 1"
+by (sos_cert "(((A<0 * R<1) + ((A<=0 * (A<=1 * R<1)) * (R<1 * [1]^2))))")
+
+lemma "(x::real) > 1 & y > 1 --> x * y > x + y - 1"
+by (sos_cert "((((A<0 * A<1) * R<1) + ((A<=0 * R<1) * (R<1 * [1]^2))))") 
 
-lemma "(0::real) <= x & 0 <= y & 0 <= z & x + y + z <= 3 -->
-  x * y + x * z + y * z >= 3 * x * y * z" by sos
+lemma "abs(x) <= 1 --> abs(64 * x^7 - 112 * x^5 + 56 * x^3 - 7 * x) <= (1::real)"
+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)))))")
+
+(* ------------------------------------------------------------------------- *)
+(* One component of denominator in dodecahedral example.                     *)
+(* ------------------------------------------------------------------------- *)
 
-lemma "((x::real)^2 + y^2 + z^2 = 1) --> (x + y + z)^2 <= 3" by sos
+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)"
+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)))))))))))))")
 
-lemma "(w^2 + x^2 + y^2 + z^2 = 1) --> (w + x + y + z)^2 <= (4::real)" by sos
+(* ------------------------------------------------------------------------- *)
+(* Over a larger but simpler interval.                                       *)
+(* ------------------------------------------------------------------------- *)
 
-lemma "(x::real) >= 1 & y >= 1 --> x * y >= x + y - 1" by sos
+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)"
+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)))))))))))")
 
-lemma "(x::real) > 1 & y > 1 --> x * y > x + y - 1" by sos; 
+(* ------------------------------------------------------------------------- *)
+(* We can do 12. I think 12 is a sharp bound; see PP's certificate.          *)
+(* ------------------------------------------------------------------------- *)
+
+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)"
+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))))))))))))))))))")
 
-lemma "abs(x) <= 1 --> abs(64 * x^7 - 112 * x^5 + 56 * x^3 - 7 * x) <= (1::real)" by sos  
+(* ------------------------------------------------------------------------- *)
+(* Inequality from sci.math (see "Leon-Sotelo, por favor").                  *)
+(* ------------------------------------------------------------------------- *)
 
+lemma "0 <= (x::real) & 0 <= y & (x * y = 1) --> x + y <= x^2 + y^2"
+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))))))") 
 
-text {* One component of denominator in dodecahedral example. *}
+lemma "0 <= (x::real) & 0 <= y & (x * y = 1) --> x * y * (x + y) <= x^2 + y^2"
+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))))))") 
 
-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)" by sos
+lemma "0 <= (x::real) & 0 <= y --> x * y * (x + y)^2 <= (x^2 + y^2)^2"
+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)))))")
 
+lemma "(0::real) <= a & 0 <= b & 0 <= c & c * (2 * a + b)^3/ 27 <= x \<longrightarrow> c * a^2 * b <= x"
+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))))))")
+ 
+lemma "(0::real) < x --> 0 < 1 + x + x^2"
+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))))))")
 
-text {* Over a larger but simpler interval. *}
+lemma "(0::real) <= x --> 0 < 1 + x + x^2"
+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))))))")
 
-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)" by sos
+lemma "(0::real) < 1 + x^2"
+by (sos_cert "((R<1 + ((R<1 * (R<1 * [x]^2)) + ((A<=0 * R<1) * (R<1 * [1]^2)))))")
+
+lemma "(0::real) <= 1 + 2 * x + x^2"
+by (sos_cert "(((A<0 * R<1) + (R<1 * (R<1 * [x + 1]^2))))")
 
-text {* We can do 12. I think 12 is a sharp bound; see PP's certificate. *}
+lemma "(0::real) < 1 + abs x"
+by (sos_cert "((R<1 + (((A<=1 * R<1) * (R<1/2 * [1]^2)) + ((A<=0 * R<1) * (R<1/2 * [1]^2)))))")
 
-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)" by sos
+lemma "(0::real) < 1 + (1 + x)^2 * (abs x)"
+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))))))")
 
 
-text {* Inequality from sci.math (see "Leon-Sotelo, por favor"). *}
 
-lemma "0 <= (x::real) & 0 <= y & (x * y = 1) --> x + y <= x^2 + y^2" by sos 
-
-lemma "0 <= (x::real) & 0 <= y & (x * y = 1) --> x * y * (x + y) <= x^2 + y^2" by sos 
-
-lemma "0 <= (x::real) & 0 <= y --> x * y * (x + y)^2 <= (x^2 + y^2)^2" by sos
+lemma "abs ((1::real) + x^2) = (1::real) + x^2"
+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)))))))")
+lemma "(3::real) * x + 7 * a < 4 \<and> 3 < 2 * x \<longrightarrow> a < 0"
+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))))))")
 
-lemma "(0::real) <= a & 0 <= b & 0 <= c & c * (2 * a + b)^3/ 27 <= x \<longrightarrow> c * a^2 * b <= x" by sos
- 
-lemma "(0::real) < x --> 0 < 1 + x + x^2" by sos
-
-lemma "(0::real) <= x --> 0 < 1 + x + x^2" by sos
-
-lemma "(0::real) < 1 + x^2" by sos
-
-lemma "(0::real) <= 1 + 2 * x + x^2" by sos
-
-lemma "(0::real) < 1 + abs x" by sos
-
-lemma "(0::real) < 1 + (1 + x)^2 * (abs x)" by sos
+lemma "(0::real) < x --> 1 < y --> y * x <= z --> x < z"
+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)))))")
+lemma "(1::real) < x --> x^2 < y --> 1 < y"
+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)))))))))")
+lemma "(b::real)^2 < 4 * a * c --> ~(a * x^2 + b * x + c = 0)"
+by (sos_cert "(((A<0 * R<1) + (R<1 * (R<1 * [2*a*x + b]^2))))")
+lemma "(b::real)^2 < 4 * a * c --> ~(a * x^2 + b * x + c = 0)"
+by (sos_cert "(((A<0 * R<1) + (R<1 * (R<1 * [2*a*x + b]^2))))")
+lemma "((a::real) * x^2 + b * x + c = 0) --> b^2 >= 4 * a * c"
+by (sos_cert "(((A<0 * R<1) + (R<1 * (R<1 * [2*a*x + b]^2))))")
+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_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)))))")
+lemma "abs(x - z) <= e & abs(y - z) <= e & 0 <= u & 0 <= v & (u + v = 1) --> abs((u * x + v * y) - z) <= (e::real)"
+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))))))")
 
 
-lemma "abs ((1::real) + x^2) = (1::real) + x^2" by sos
-lemma "(3::real) * x + 7 * a < 4 \<and> 3 < 2 * x \<longrightarrow> a < 0" by sos
-
-lemma "(0::real) < x --> 1 < y --> y * x <= z --> x < z" by sos
-lemma "(1::real) < x --> x^2 < y --> 1 < y" by sos
-lemma "(b::real)^2 < 4 * a * c --> ~(a * x^2 + b * x + c = 0)" by sos
-lemma "(b::real)^2 < 4 * a * c --> ~(a * x^2 + b * x + c = 0)" by sos
-lemma "((a::real) * x^2 + b * x + c = 0) --> b^2 >= 4 * a * c" by sos
-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
-lemma "abs(x - z) <= e & abs(y - z) <= e & 0 <= u & 0 <= v & (u + v = 1) -->
-  abs((u * x + v * y) - z) <= (e::real)" by sos
-
-(*
-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?*}
-*)
+(* 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?*)
 
 lemma "(0::real) <= x --> (1 + x + x^2)/(1 + x^2) <= 1 + x"
-  by sos
+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))))))")
 
 lemma "(0::real) <= x --> 1 - x <= 1 / (1 + x + x^2)"
-  by sos
+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))))))))")
 
 lemma "(x::real) <= 1 / 2 --> - x - 2 * x^2 <= - x / (1 - x)"
-  by sos
+by (sos_cert "((((A<0 * A<1) * R<1) + ((A<=0 * (A<0 * R<1)) * (R<1 * [x]^2))))")
 
-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" by sos
-*)
+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"
+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)))))")
 
 end
+
diff -r e4511a1b4c1b -r 1cc5b24f5a01 src/HOL/Library/Sum_Of_Squares/positivstellensatz_tools.ML
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/HOL/Library/Sum_Of_Squares/positivstellensatz_tools.ML	Mon Sep 21 15:05:26 2009 +0200
@@ -0,0 +1,158 @@
+(* Title:      positivstellensatz_tools.ML
+   Author:     Philipp Meyer, TU Muenchen
+
+Functions for generating a certificate from a positivstellensatz and vice versa
+*)
+
+signature POSITIVSTELLENSATZ_TOOLS =
+sig
+  val pss_tree_to_cert : RealArith.pss_tree -> string
+
+  val cert_to_pss_tree : Proof.context -> string -> RealArith.pss_tree
+
+end
+
+
+structure PositivstellensatzTools : POSITIVSTELLENSATZ_TOOLS =
+struct
+
+open RealArith FuncUtil
+
+(*** certificate generation ***)
+
+fun string_of_rat r =
+  let
+    val (nom, den) = Rat.quotient_of_rat r
+  in
+    if den = 1 then string_of_int nom
+    else string_of_int nom ^ "/" ^ string_of_int den
+  end
+
+(* map polynomials to strings *)
+
+fun string_of_varpow x k =
+  let
+    val term = term_of x
+    val name = case term of
+      Free (n, _) => n
+    | _ => error "Term in monomial not free variable"
+  in
+    if k = 1 then name else name ^ "^" ^ string_of_int k 
+  end
+
+fun string_of_monomial m = 
+ if Ctermfunc.is_undefined m then "1" 
+ else 
+  let 
+   val m' = dest_monomial m
+   val vps = fold_rev (fn (x,k) => cons (string_of_varpow x k)) m' [] 
+  in foldr1 (fn (s, t) => s ^ "*" ^ t) vps
+  end
+
+fun string_of_cmonomial (m,c) =
+  if Ctermfunc.is_undefined m then string_of_rat c
+  else if c = Rat.one then string_of_monomial m
+  else (string_of_rat c) ^ "*" ^ (string_of_monomial m);
+
+fun string_of_poly p = 
+ if Monomialfunc.is_undefined p then "0" 
+ else
+  let 
+   val cms = map string_of_cmonomial
+     (sort (prod_ord monomial_order (K EQUAL)) (Monomialfunc.graph p))
+  in foldr1 (fn (t1, t2) => t1 ^ " + " ^ t2) cms
+  end;
+
+fun pss_to_cert (Axiom_eq i) = "A=" ^ string_of_int i
+  | pss_to_cert (Axiom_le i) = "A<=" ^ string_of_int i
+  | pss_to_cert (Axiom_lt i) = "A<" ^ string_of_int i
+  | pss_to_cert (Rational_eq r) = "R=" ^ string_of_rat r
+  | pss_to_cert (Rational_le r) = "R<=" ^ string_of_rat r
+  | pss_to_cert (Rational_lt r) = "R<" ^ string_of_rat r
+  | pss_to_cert (Square p) = "[" ^ string_of_poly p ^ "]^2"
+  | pss_to_cert (Eqmul (p, pss)) = "([" ^ string_of_poly p ^ "] * " ^ pss_to_cert pss ^ ")"
+  | pss_to_cert (Sum (pss1, pss2)) = "(" ^ pss_to_cert pss1 ^ " + " ^ pss_to_cert pss2 ^ ")"
+  | pss_to_cert (Product (pss1, pss2)) = "(" ^ pss_to_cert pss1 ^ " * " ^ pss_to_cert pss2 ^ ")"
+
+fun pss_tree_to_cert Trivial = "()"
+  | pss_tree_to_cert (Cert pss) = "(" ^ pss_to_cert pss ^ ")"
+  | pss_tree_to_cert (Branch (t1, t2)) = "(" ^ pss_tree_to_cert t1 ^ " & " ^ pss_tree_to_cert t2 ^ ")"
+
+(*** certificate parsing ***)
+
+(* basic parser *)
+
+fun $$ k = Scan.this_string k
+
+val number = Scan.repeat1 (Scan.one Symbol.is_ascii_digit >>
+  (fn s => ord s - ord "0")) >>
+  foldl1 (fn (n, d) => n * 10 + d)
+
+val nat = number
+val int = Scan.optional ($$ "~" >> K ~1) 1 -- nat >> op *;
+val rat = int --| $$ "/" -- int >> Rat.rat_of_quotient
+val rat_int = rat || int >> Rat.rat_of_int
+
+(* polynomial parser *)
+
+fun repeat_sep s f = f ::: Scan.repeat ($$ s |-- f)
+
+val parse_id = Scan.one Symbol.is_letter ::: Scan.many Symbol.is_letdig >> implode
+
+fun parse_varpow ctxt = parse_id -- Scan.optional ($$ "^" |-- nat) 1 >>
+  (fn (x, k) => (cterm_of (Context.theory_of_proof ctxt) (Free (x, @{typ real})), k)) 
+
+fun parse_monomial ctxt = repeat_sep "*" (parse_varpow ctxt) >>
+  foldl (uncurry Ctermfunc.update) Ctermfunc.undefined
+
+fun parse_cmonomial ctxt =
+  rat_int --| $$ "*" -- (parse_monomial ctxt) >> swap ||
+  (parse_monomial ctxt) >> (fn m => (m, Rat.one)) ||
+  rat_int >> (fn r => (Ctermfunc.undefined, r))
+
+fun parse_poly ctxt = repeat_sep "+" (parse_cmonomial ctxt) >>
+  foldl (uncurry Monomialfunc.update) Monomialfunc.undefined
+
+(* positivstellensatz parser *)
+
+val parse_axiom =
+  ($$ "A=" |-- int >> Axiom_eq) ||
+  ($$ "A<=" |-- int >> Axiom_le) ||
+  ($$ "A<" |-- int >> Axiom_lt)
+
+val parse_rational =
+  ($$ "R=" |-- rat_int >> Rational_eq) ||
+  ($$ "R<=" |-- rat_int >> Rational_le) ||
+  ($$ "R<" |-- rat_int >> Rational_lt)
+
+fun parse_cert ctxt input =
+  let
+    val pc = parse_cert ctxt
+    val pp = parse_poly ctxt
+  in
+  (parse_axiom ||
+   parse_rational ||
+   $$ "[" |-- pp --| $$ "]^2" >> Square ||
+   $$ "([" |-- pp --| $$ "]*" -- pc --| $$ ")" >> Eqmul ||
+   $$ "(" |-- pc --| $$ "*" -- pc --| $$ ")" >> Product ||
+   $$ "(" |-- pc --| $$ "+" -- pc --| $$ ")" >> Sum) input
+  end
+
+fun parse_cert_tree ctxt input =
+  let
+    val pc = parse_cert ctxt
+    val pt = parse_cert_tree ctxt
+  in
+  ($$ "()" >> K Trivial ||
+   $$ "(" |-- pc --| $$ ")" >> Cert ||
+   $$ "(" |-- pt --| $$ "&" -- pt --| $$ ")" >> Branch) input
+  end
+
+(* scanner *)
+
+fun cert_to_pss_tree ctxt str = Symbol.scanner "bad certificate" (parse_cert_tree ctxt)
+  (filter_out Symbol.is_blank (Symbol.explode str))
+
+end
+
+
diff -r e4511a1b4c1b -r 1cc5b24f5a01 src/HOL/Library/Sum_Of_Squares/sos_wrapper.ML
--- a/src/HOL/Library/Sum_Of_Squares/sos_wrapper.ML	Tue Sep 22 20:25:31 2009 +0200
+++ b/src/HOL/Library/Sum_Of_Squares/sos_wrapper.ML	Mon Sep 21 15:05:26 2009 +0200
@@ -136,13 +136,42 @@
     run_solver name (Path.explode cmd) find_failure
   end
 
+(* certificate output *)
+
+fun output_line cert = "To repeat this proof with a certifiate use this command:\n" ^
+        (Markup.markup Markup.sendback) ("by (sos_cert \"" ^ cert ^ "\")")
+
+fun print_certs pss_tree =
+  let
+    val cert = PositivstellensatzTools.pss_tree_to_cert pss_tree
+    val str = output_line cert
+  in
+    Output.writeln str
+  end
+
 (* setup tactic *)
 
-fun sos_solver name = (SIMPLE_METHOD' o (Sos.sos_tac (call_solver name))) 
+fun parse_cert (ctxt, tk) =
+  let
+    val (str, tk') = OuterParse.string tk
+    val cert = PositivstellensatzTools.cert_to_pss_tree (Context.proof_of ctxt) str
+  in
+    (cert, (ctxt, tk'))
+  end
+
+fun sos_solver print method = (SIMPLE_METHOD' o (Sos.sos_tac print method)) 
 
-val sos_method = Scan.option (Scan.lift OuterParse.xname) >> sos_solver
+val sos_method =
+   Scan.lift (Scan.option OuterParse.xname) >> (fn n => Sos.Prover (call_solver n)) >>
+   sos_solver print_certs
 
-val setup = Method.setup @{binding sos} sos_method
-  "Prove universal problems over the reals using sums of squares"
+val sos_cert_method =
+  parse_cert >> Sos.Certificate >> sos_solver ignore
+
+val setup =
+     Method.setup @{binding sos} sos_method
+     "Prove universal problems over the reals using sums of squares"
+  #> Method.setup @{binding sos_cert} sos_cert_method
+     "Prove universal problems over the reals using sums of squares with certificates"
 
 end
diff -r e4511a1b4c1b -r 1cc5b24f5a01 src/HOL/Library/Sum_Of_Squares/sum_of_squares.ML
--- a/src/HOL/Library/Sum_Of_Squares/sum_of_squares.ML	Tue Sep 22 20:25:31 2009 +0200
+++ b/src/HOL/Library/Sum_Of_Squares/sum_of_squares.ML	Mon Sep 21 15:05:26 2009 +0200
@@ -8,7 +8,12 @@
 signature SOS =
 sig
 
-  val sos_tac : (string -> string) -> Proof.context -> int -> Tactical.tactic
+  datatype proof_method =
+    Certificate of RealArith.pss_tree
+  | Prover of (string -> string)
+
+  val sos_tac : (RealArith.pss_tree -> unit) ->
+    proof_method -> Proof.context -> int -> Tactical.tactic
 
   val debugging : bool ref;
   
@@ -18,6 +23,8 @@
 structure Sos : SOS = 
 struct
 
+open FuncUtil;
+
 val rat_0 = Rat.zero;
 val rat_1 = Rat.one;
 val rat_2 = Rat.two;
@@ -59,6 +66,10 @@
 
 exception Failure of string;
 
+datatype proof_method =
+    Certificate of RealArith.pss_tree
+  | Prover of (string -> string)
+
 (* Turn a rational into a decimal string with d sig digits.                  *)
 
 local
@@ -93,23 +104,11 @@
 
 (* The main types.                                                           *)
 
-fun strict_ord ord (x,y) = case ord (x,y) of LESS => LESS | _ => GREATER
-
-structure Intpairfunc = FuncFun(type key = int*int val ord = prod_ord int_ord int_ord);
-
 type vector = int* Rat.rat Intfunc.T;
 
 type matrix = (int*int)*(Rat.rat Intpairfunc.T);
 
-type monomial = int Ctermfunc.T;
-
-val cterm_ord = (fn (s,t) => TermOrd.fast_term_ord(term_of s, term_of t))
- fun monomial_ord (m1,m2) = list_ord (prod_ord cterm_ord int_ord) (Ctermfunc.graph m1, Ctermfunc.graph m2)
-structure Monomialfunc = FuncFun(type key = monomial val ord = monomial_ord)
-
-type poly = Rat.rat Monomialfunc.T;
-
- fun iszero (k,r) = r =/ rat_0;
+fun iszero (k,r) = r =/ rat_0;
 
 fun fold_rev2 f l1 l2 b =
   case (l1,l2) of
@@ -346,10 +345,7 @@
   sort humanorder_varpow (Ctermfunc.graph m2))
 end;
 
-fun fold1 f l =  case l of
-   []     => error "fold1"
- | [x]    => x
- | (h::t) => f h (fold1 f t);
+fun fold1 f = foldr1 (uncurry f) 
 
 (* Conversions to strings.                                                   *)
 
@@ -404,7 +400,7 @@
  else if c =/ rat_1 then string_of_monomial m
  else Rat.string_of_rat c ^ "*" ^ string_of_monomial m;;
 
-fun string_of_poly (p:poly) =
+fun string_of_poly p =
  if Monomialfunc.is_undefined p then "<<0>>" else
  let 
   val cms = sort (fn ((m1,_),(m2,_)) => humanorder_monomial m1  m2) (Monomialfunc.graph p)
@@ -481,7 +477,6 @@
  in fold1 (fn x => fn y => x ^ " " ^ y) strs ^ "\n"
  end;
 
-fun increasing f ord (x,y) = ord (f x, f y);
 fun triple_int_ord ((a,b,c),(a',b',c')) = 
  prod_ord int_ord (prod_ord int_ord int_ord) 
     ((a,(b,c)),(a',(b',c')));
@@ -1080,11 +1075,6 @@
   fun tryfind f = tryfind_with "tryfind" f
 end
 
-(*
-fun tryfind f [] = error "tryfind"
-  | tryfind f (x::xs) = (f x handle ERROR _ => tryfind f xs);
-*)
-
 (* Positiv- and Nullstellensatz. Flag "linf" forces a linear representation. *)
 
  
@@ -1210,58 +1200,14 @@
 fun deepen f n = 
   (writeln ("Searching with depth limit " ^ string_of_int n) ; (f n handle Failure s => (writeln ("failed with message: " ^ s) ; deepen f (n+1))))
 
-(* The ordering so we can create canonical HOL polynomials.                  *)
 
-fun dest_monomial mon = sort (increasing fst cterm_ord) (Ctermfunc.graph mon);
-
-fun monomial_order (m1,m2) =
- if Ctermfunc.is_undefined m2 then LESS 
- else if Ctermfunc.is_undefined m1 then GREATER 
- else
-  let val mon1 = dest_monomial m1 
-      val mon2 = dest_monomial m2
-      val deg1 = fold (curry op + o snd) mon1 0
-      val deg2 = fold (curry op + o snd) mon2 0 
-  in if deg1 < deg2 then GREATER else if deg1 > deg2 then LESS
-     else list_ord (prod_ord cterm_ord int_ord) (mon1,mon2)
-  end;
-
-fun dest_poly p =
-  map (fn (m,c) => (c,dest_monomial m))
-      (sort (prod_ord monomial_order (K EQUAL)) (Monomialfunc.graph p));
-
-(* Map back polynomials and their composites to HOL.                         *)
+(* Map back polynomials and their composites to a positivstellensatz.        *)
 
 local
  open Thm Numeral RealArith
 in
 
-fun cterm_of_varpow x k = if k = 1 then x else capply (capply @{cterm "op ^ :: real => _"} x) 
-  (mk_cnumber @{ctyp nat} k)
-
-fun cterm_of_monomial m = 
- if Ctermfunc.is_undefined m then @{cterm "1::real"} 
- else 
-  let 
-   val m' = dest_monomial m
-   val vps = fold_rev (fn (x,k) => cons (cterm_of_varpow x k)) m' [] 
-  in fold1 (fn s => fn t => capply (capply @{cterm "op * :: real => _"} s) t) vps
-  end
-
-fun cterm_of_cmonomial (m,c) = if Ctermfunc.is_undefined m then cterm_of_rat c
-    else if c = Rat.one then cterm_of_monomial m
-    else capply (capply @{cterm "op *::real => _"} (cterm_of_rat c)) (cterm_of_monomial m);
-
-fun cterm_of_poly p = 
- if Monomialfunc.is_undefined p then @{cterm "0::real"} 
- else
-  let 
-   val cms = map cterm_of_cmonomial
-     (sort (prod_ord monomial_order (K EQUAL)) (Monomialfunc.graph p))
-  in fold1 (fn t1 => fn t2 => capply(capply @{cterm "op + :: real => _"} t1) t2) cms
-  end;
-
-fun cterm_of_sqterm (c,p) = Product(Rational_lt c,Square(cterm_of_poly p));
+fun cterm_of_sqterm (c,p) = Product(Rational_lt c,Square p);
 
 fun cterm_of_sos (pr,sqs) = if null sqs then pr
   else Product(pr,fold1 (fn a => fn b => Sum(a,b)) (map cterm_of_sqterm sqs));
@@ -1275,14 +1221,14 @@
   fun simple_cterm_ord t u = TermOrd.fast_term_ord (term_of t, term_of u) = LESS
 in
   (* FIXME: Replace tryfind by get_first !! *)
-fun real_nonlinear_prover prover ctxt =
+fun real_nonlinear_prover proof_method ctxt =
  let 
   val {add,mul,neg,pow,sub,main} =  Normalizer.semiring_normalizers_ord_wrapper ctxt
       (valOf (NormalizerData.match ctxt @{cterm "(0::real) + 1"})) 
      simple_cterm_ord
   val (real_poly_add_conv,real_poly_mul_conv,real_poly_neg_conv,
        real_poly_pow_conv,real_poly_sub_conv,real_poly_conv) = (add,mul,neg,pow,sub,main)
-  fun mainf  translator (eqs,les,lts) = 
+  fun mainf cert_choice translator (eqs,les,lts) = 
   let 
    val eq0 = map (poly_of_term o dest_arg1 o concl) eqs
    val le0 = map (poly_of_term o dest_arg o concl) les
@@ -1303,33 +1249,49 @@
                      else raise Failure "trivial_axiom: Not a trivial axiom"
      | _ => error "trivial_axiom: Not a trivial axiom"
    in 
-  ((let val th = tryfind trivial_axiom (keq @ klep @ kltp)
-   in fconv_rule (arg_conv (arg1_conv real_poly_conv) then_conv field_comp_conv) th end)
-   handle Failure _ => (
-    let 
-     val pol = fold_rev poly_mul (map fst ltp) (poly_const Rat.one)
-     val leq = lep @ ltp
-     fun tryall d =
-      let val e = multidegree pol
-          val k = if e = 0 then 0 else d div e
-          val eq' = map fst eq 
-      in tryfind (fn i => (d,i,real_positivnullstellensatz_general prover false d eq' leq
-                            (poly_neg(poly_pow pol i))))
-              (0 upto k)
-      end
-    val (d,i,(cert_ideal,cert_cone)) = deepen tryall 0
-    val proofs_ideal =
-      map2 (fn q => fn (p,ax) => Eqmul(cterm_of_poly q,ax)) cert_ideal eq
-    val proofs_cone = map cterm_of_sos cert_cone
-    val proof_ne = if null ltp then Rational_lt Rat.one else
-      let val p = fold1 (fn s => fn t => Product(s,t)) (map snd ltp) 
-      in  funpow i (fn q => Product(p,q)) (Rational_lt Rat.one)
-      end
-    val proof = fold1 (fn s => fn t => Sum(s,t))
-                           (proof_ne :: proofs_ideal @ proofs_cone) 
-    in writeln "Translating proof certificate to HOL";
-       translator (eqs,les,lts) proof
-    end))
+  (let val th = tryfind trivial_axiom (keq @ klep @ kltp)
+   in
+    (fconv_rule (arg_conv (arg1_conv real_poly_conv) then_conv field_comp_conv) th, Trivial)
+   end)
+   handle Failure _ => 
+     (let val proof =
+       (case proof_method of Certificate certs =>
+         (* choose certificate *)
+         let
+           fun chose_cert [] (Cert c) = c
+             | chose_cert (Left::s) (Branch (l, _)) = chose_cert s l
+             | chose_cert (Right::s) (Branch (_, r)) = chose_cert s r
+             | chose_cert _ _ = error "certificate tree in invalid form"
+         in
+           chose_cert cert_choice certs
+         end
+       | Prover prover =>
+         (* call prover *)
+         let 
+          val pol = fold_rev poly_mul (map fst ltp) (poly_const Rat.one)
+          val leq = lep @ ltp
+          fun tryall d =
+           let val e = multidegree pol
+               val k = if e = 0 then 0 else d div e
+               val eq' = map fst eq 
+           in tryfind (fn i => (d,i,real_positivnullstellensatz_general prover false d eq' leq
+                                 (poly_neg(poly_pow pol i))))
+                   (0 upto k)
+           end
+         val (d,i,(cert_ideal,cert_cone)) = deepen tryall 0
+         val proofs_ideal =
+           map2 (fn q => fn (p,ax) => Eqmul(q,ax)) cert_ideal eq
+         val proofs_cone = map cterm_of_sos cert_cone
+         val proof_ne = if null ltp then Rational_lt Rat.one else
+           let val p = fold1 (fn s => fn t => Product(s,t)) (map snd ltp) 
+           in  funpow i (fn q => Product(p,q)) (Rational_lt Rat.one)
+           end
+         in 
+           fold1 (fn s => fn t => Sum(s,t)) (proof_ne :: proofs_ideal @ proofs_cone) 
+         end)
+     in
+        (translator (eqs,les,lts) proof, Cert proof)
+     end)
    end
  in mainf end
 end
@@ -1396,7 +1358,7 @@
          orelse g aconvc @{cterm "op < :: real => _"} 
        then arg_conv cv ct else arg1_conv cv ct
     end
-  fun mainf translator =
+  fun mainf cert_choice translator =
    let 
     fun substfirst(eqs,les,lts) =
       ((let 
@@ -1407,7 +1369,7 @@
                                    aconvc @{cterm "0::real"}) (map modify eqs),
                                    map modify les,map modify lts)
        end)
-       handle Failure  _ => real_nonlinear_prover prover ctxt translator (rev eqs, rev les, rev lts))
+       handle Failure  _ => real_nonlinear_prover prover ctxt cert_choice translator (rev eqs, rev les, rev lts))
     in substfirst
    end
 
@@ -1417,7 +1379,8 @@
 
 (* Overall function. *)
 
-fun real_sos prover ctxt t = gen_prover_real_arith ctxt (real_nonlinear_subst_prover prover ctxt) t;
+fun real_sos prover ctxt =
+  gen_prover_real_arith ctxt (real_nonlinear_subst_prover prover ctxt)
 end;
 
 (* A tactic *)
@@ -1429,8 +1392,6 @@
    end
 | _ => ([],ct)
 
-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})
-
 val known_sos_constants = 
   [@{term "op ==>"}, @{term "Trueprop"}, 
    @{term "op -->"}, @{term "op &"}, @{term "op |"}, 
@@ -1458,17 +1419,19 @@
   val _ = if exists (fn ((_,T)) => not (T = @{typ "real"})) fs 
           then error "SOS: not sos. Variables with type not real" else ()
   val vs = Term.add_vars t []
-  val _ = if exists (fn ((_,T)) => not (T = @{typ "real"})) fs 
+  val _ = if exists (fn ((_,T)) => not (T = @{typ "real"})) vs 
           then error "SOS: not sos. Variables with type not real" else ()
   val ukcs = subtract (fn (t,p) => Const p aconv t) kcts (Term.add_consts t [])
   val _ = if  null ukcs then () 
               else error ("SOSO: Unknown constants in Subgoal:" ^ commas (map fst ukcs))
 in () end
 
-fun core_sos_tac prover ctxt = CSUBGOAL (fn (ct, i) => 
+fun core_sos_tac print_certs prover ctxt = CSUBGOAL (fn (ct, i) => 
   let val _ = check_sos known_sos_constants ct
       val (avs, p) = strip_all ct
-      val th = standard (fold_rev forall_intr avs (real_sos prover ctxt (Thm.dest_arg p)))
+      val (ths, certificates) = real_sos prover ctxt (Thm.dest_arg p)
+      val th = standard (fold_rev forall_intr avs ths)
+      val _ = print_certs certificates
   in rtac th i end);
 
 fun default_SOME f NONE v = SOME v
@@ -1506,7 +1469,7 @@
 
 fun elim_denom_tac ctxt i = REPEAT (elim_one_denom_tac ctxt i);
 
-fun sos_tac prover ctxt = ObjectLogic.full_atomize_tac THEN' elim_denom_tac ctxt THEN' core_sos_tac prover ctxt
+fun sos_tac print_certs prover ctxt = ObjectLogic.full_atomize_tac THEN' elim_denom_tac ctxt THEN' core_sos_tac print_certs prover ctxt
 
 
 end;
diff -r e4511a1b4c1b -r 1cc5b24f5a01 src/HOL/Library/normarith.ML
--- a/src/HOL/Library/normarith.ML	Tue Sep 22 20:25:31 2009 +0200
+++ b/src/HOL/Library/normarith.ML	Mon Sep 21 15:05:26 2009 +0200
@@ -15,7 +15,7 @@
 structure NormArith : NORM_ARITH = 
 struct
 
- open Conv Thm;
+ open Conv Thm FuncUtil;
  val bool_eq = op = : bool *bool -> bool
   fun dest_ratconst t = case term_of t of
    Const(@{const_name divide}, _)$a$b => Rat.rat_of_quotient(HOLogic.dest_number a |> snd, HOLogic.dest_number b |> snd)
@@ -330,13 +330,13 @@
    val zerodests = filter
         (fn t => null (Ctermfunc.dom (vector_lincomb t))) (map snd rawdests)
 
-  in RealArith.real_linear_prover translator
+  in fst (RealArith.real_linear_prover translator
         (map (fn t => instantiate ([(tv_n, ctyp_of_term t)],[]) pth_zero)
             zerodests,
         map (fconv_rule (try_conv (More_Conv.top_sweep_conv (K norm_canon_conv) ctxt) then_conv
                        arg_conv (arg_conv real_poly_conv))) ges',
         map (fconv_rule (try_conv (More_Conv.top_sweep_conv (K norm_canon_conv) ctxt) then_conv 
-                       arg_conv (arg_conv real_poly_conv))) gts)
+                       arg_conv (arg_conv real_poly_conv))) gts))
   end
 in val real_vector_combo_prover = real_vector_combo_prover
 end;
@@ -389,9 +389,9 @@
    val (th1,th2) = conj_pair(rawrule th)
   in th1::fconv_rule (arg_conv (arg_conv real_poly_neg_conv)) th2::acc
   end
-in fun real_vector_prover ctxt translator (eqs,ges,gts) =
-     real_vector_ineq_prover ctxt translator
-         (fold_rev (splitequation ctxt) eqs ges,gts)
+in fun real_vector_prover ctxt _ translator (eqs,ges,gts) =
+     (real_vector_ineq_prover ctxt translator
+         (fold_rev (splitequation ctxt) eqs ges,gts), RealArith.Trivial)
 end;
 
   fun init_conv ctxt = 
@@ -400,7 +400,7 @@
    then_conv field_comp_conv 
    then_conv nnf_conv
 
- fun pure ctxt = RealArith.gen_prover_real_arith ctxt (real_vector_prover ctxt);
+ fun pure ctxt = fst o RealArith.gen_prover_real_arith ctxt (real_vector_prover ctxt);
  fun norm_arith ctxt ct = 
   let 
    val ctxt' = Variable.declare_term (term_of ct) ctxt
diff -r e4511a1b4c1b -r 1cc5b24f5a01 src/HOL/Library/positivstellensatz.ML
--- a/src/HOL/Library/positivstellensatz.ML	Tue Sep 22 20:25:31 2009 +0200
+++ b/src/HOL/Library/positivstellensatz.ML	Mon Sep 21 15:05:26 2009 +0200
@@ -1,10 +1,11 @@
-(* Title:      Library/positivstellensatz
+(* Title:      Library/Sum_Of_Squares/positivstellensatz
    Author:     Amine Chaieb, University of Cambridge
    Description: A generic arithmetic prover based on Positivstellensatz certificates --- 
     also implements Fourrier-Motzkin elimination as a special case Fourrier-Motzkin elimination.
 *)
 
 (* A functor for finite mappings based on Tables *)
+
 signature FUNC = 
 sig
  type 'a T
@@ -75,27 +76,54 @@
 end
 end;
 
+(* Some standard functors and utility functions for them *)
+
+structure FuncUtil =
+struct
+
+fun increasing f ord (x,y) = ord (f x, f y);
+
 structure Intfunc = FuncFun(type key = int val ord = int_ord);
+structure Ratfunc = FuncFun(type key = Rat.rat val ord = Rat.ord);
+structure Intpairfunc = FuncFun(type key = int*int val ord = prod_ord int_ord int_ord);
 structure Symfunc = FuncFun(type key = string val ord = fast_string_ord);
 structure Termfunc = FuncFun(type key = term val ord = TermOrd.fast_term_ord);
-structure Ctermfunc = FuncFun(type key = cterm val ord = (fn (s,t) => TermOrd.fast_term_ord(term_of s, term_of t)));
+
+val cterm_ord = (fn (s,t) => TermOrd.fast_term_ord(term_of s, term_of t))
+
+structure Ctermfunc = FuncFun(type key = cterm val ord = cterm_ord);
+
+type monomial = int Ctermfunc.T;
 
-structure Ratfunc = FuncFun(type key = Rat.rat val ord = Rat.ord);
+fun monomial_ord (m1,m2) = list_ord (prod_ord cterm_ord int_ord) (Ctermfunc.graph m1, Ctermfunc.graph m2)
+
+structure Monomialfunc = FuncFun(type key = monomial val ord = monomial_ord)
 
+type poly = Rat.rat Monomialfunc.T;
+
+(* The ordering so we can create canonical HOL polynomials.                  *)
 
-    (* Some useful derived rules *)
-fun deduct_antisym_rule tha thb = 
-    equal_intr (implies_intr (cprop_of thb) tha) 
-     (implies_intr (cprop_of tha) thb);
+fun dest_monomial mon = sort (increasing fst cterm_ord) (Ctermfunc.graph mon);
 
-fun prove_hyp tha thb = 
-  if exists (curry op aconv (concl_of tha)) (#hyps (rep_thm thb)) 
-  then equal_elim (symmetric (deduct_antisym_rule tha thb)) tha else thb;
+fun monomial_order (m1,m2) =
+ if Ctermfunc.is_undefined m2 then LESS 
+ else if Ctermfunc.is_undefined m1 then GREATER 
+ else
+  let val mon1 = dest_monomial m1 
+      val mon2 = dest_monomial m2
+      val deg1 = fold (curry op + o snd) mon1 0
+      val deg2 = fold (curry op + o snd) mon2 0 
+  in if deg1 < deg2 then GREATER else if deg1 > deg2 then LESS
+     else list_ord (prod_ord cterm_ord int_ord) (mon1,mon2)
+  end;
 
+end
 
+(* positivstellensatz datatype and prover generation *)
 
 signature REAL_ARITH = 
 sig
+  
   datatype positivstellensatz =
    Axiom_eq of int
  | Axiom_le of int
@@ -103,34 +131,41 @@
  | Rational_eq of Rat.rat
  | Rational_le of Rat.rat
  | Rational_lt of Rat.rat
- | Square of cterm
- | Eqmul of cterm * positivstellensatz
+ | Square of FuncUtil.poly
+ | Eqmul of FuncUtil.poly * positivstellensatz
  | Sum of positivstellensatz * positivstellensatz
  | Product of positivstellensatz * positivstellensatz;
 
+datatype pss_tree = Trivial | Cert of positivstellensatz | Branch of pss_tree * pss_tree
+
+datatype tree_choice = Left | Right
+
+type prover = tree_choice list -> 
+  (thm list * thm list * thm list -> positivstellensatz -> thm) ->
+  thm list * thm list * thm list -> thm * pss_tree
+type cert_conv = cterm -> thm * pss_tree
+
 val gen_gen_real_arith :
-  Proof.context -> (Rat.rat -> Thm.cterm) * conv * conv * conv * 
-   conv * conv * conv * conv * conv * conv * 
-    ( (thm list * thm list * thm list -> positivstellensatz -> thm) ->
-        thm list * thm list * thm list -> thm) -> conv
-val real_linear_prover : 
-  (thm list * thm list * thm list -> positivstellensatz -> thm) ->
-   thm list * thm list * thm list -> thm
+  Proof.context -> (Rat.rat -> Thm.cterm) * conv * conv * conv *
+   conv * conv * conv * conv * conv * conv * prover -> cert_conv
+val real_linear_prover : (thm list * thm list * thm list -> positivstellensatz -> thm) ->
+  thm list * thm list * thm list -> thm * pss_tree
 
 val gen_real_arith : Proof.context ->
-   (Rat.rat -> cterm) * conv * conv * conv * conv * conv * conv * conv *
-   ( (thm list * thm list * thm list -> positivstellensatz -> thm) ->
-       thm list * thm list * thm list -> thm) -> conv
-val gen_prover_real_arith : Proof.context ->
-   ((thm list * thm list * thm list -> positivstellensatz -> thm) ->
-     thm list * thm list * thm list -> thm) -> conv
-val real_arith : Proof.context -> conv
+  (Rat.rat -> Thm.cterm) * conv * conv * conv * conv * conv * conv * conv * prover -> cert_conv
+
+val gen_prover_real_arith : Proof.context -> prover -> cert_conv
+
+val is_ratconst : Thm.cterm -> bool
+val dest_ratconst : Thm.cterm -> Rat.rat
+val cterm_of_rat : Rat.rat -> Thm.cterm
+
 end
 
-structure RealArith (* : REAL_ARITH *)=
+structure RealArith : REAL_ARITH =
 struct
 
- open Conv Thm;;
+ open Conv Thm FuncUtil;;
 (* ------------------------------------------------------------------------- *)
 (* Data structure for Positivstellensatz refutations.                        *)
 (* ------------------------------------------------------------------------- *)
@@ -142,12 +177,18 @@
  | Rational_eq of Rat.rat
  | Rational_le of Rat.rat
  | Rational_lt of Rat.rat
- | Square of cterm
- | Eqmul of cterm * positivstellensatz
+ | Square of FuncUtil.poly
+ | Eqmul of FuncUtil.poly * positivstellensatz
  | Sum of positivstellensatz * positivstellensatz
  | Product of positivstellensatz * positivstellensatz;
          (* Theorems used in the procedure *)
 
+datatype pss_tree = Trivial | Cert of positivstellensatz | Branch of pss_tree * pss_tree
+datatype tree_choice = Left | Right
+type prover = tree_choice list -> 
+  (thm list * thm list * thm list -> positivstellensatz -> thm) ->
+  thm list * thm list * thm list -> thm * pss_tree
+type cert_conv = cterm -> thm * pss_tree
 
 val my_eqs = ref ([] : thm list);
 val my_les = ref ([] : thm list);
@@ -164,6 +205,16 @@
 val my_poly_add_conv = ref no_conv;
 val my_poly_mul_conv = ref no_conv;
 
+
+    (* Some useful derived rules *)
+fun deduct_antisym_rule tha thb = 
+    equal_intr (implies_intr (cprop_of thb) tha) 
+     (implies_intr (cprop_of tha) thb);
+
+fun prove_hyp tha thb = 
+  if exists (curry op aconv (concl_of tha)) (#hyps (rep_thm thb)) 
+  then equal_elim (symmetric (deduct_antisym_rule tha thb)) tha else thb;
+
 fun conjunctions th = case try Conjunction.elim th of
    SOME (th1,th2) => (conjunctions th1) @ conjunctions th2
  | NONE => [th];
@@ -292,7 +343,6 @@
  | Abs (_,_,t') => find_cterm p (Thm.dest_abs NONE t |> snd)
  | _ => raise CTERM ("find_cterm",[t]);
 
-
     (* Some conversions-related stuff which has been forbidden entrance into Pure/conv.ML*)
 fun instantiate_cterm' ty tms = Drule.cterm_rule (Drule.instantiate' ty tms)
 fun is_comb t = case (term_of t) of _$_ => true | _ => false;
@@ -300,6 +350,39 @@
 fun is_binop ct ct' = ct aconvc (Thm.dest_fun (Thm.dest_fun ct'))
   handle CTERM _ => false;
 
+
+(* Map back polynomials to HOL.                         *)
+
+local
+ open Thm Numeral
+in
+
+fun cterm_of_varpow x k = if k = 1 then x else capply (capply @{cterm "op ^ :: real => _"} x) 
+  (mk_cnumber @{ctyp nat} k)
+
+fun cterm_of_monomial m = 
+ if Ctermfunc.is_undefined m then @{cterm "1::real"} 
+ else 
+  let 
+   val m' = dest_monomial m
+   val vps = fold_rev (fn (x,k) => cons (cterm_of_varpow x k)) m' [] 
+  in foldr1 (fn (s, t) => capply (capply @{cterm "op * :: real => _"} s) t) vps
+  end
+
+fun cterm_of_cmonomial (m,c) = if Ctermfunc.is_undefined m then cterm_of_rat c
+    else if c = Rat.one then cterm_of_monomial m
+    else capply (capply @{cterm "op *::real => _"} (cterm_of_rat c)) (cterm_of_monomial m);
+
+fun cterm_of_poly p = 
+ if Monomialfunc.is_undefined p then @{cterm "0::real"} 
+ else
+  let 
+   val cms = map cterm_of_cmonomial
+     (sort (prod_ord monomial_order (K EQUAL)) (Monomialfunc.graph p))
+  in foldr1 (fn (t1, t2) => capply(capply @{cterm "op + :: real => _"} t1) t2) cms
+  end;
+
+end;
     (* A general real arithmetic prover *)
 
 fun gen_gen_real_arith ctxt (mk_numeric,
@@ -368,8 +451,8 @@
       | Rational_lt x => eqT_elim(numeric_gt_conv(capply @{cterm Trueprop} 
                       (capply (capply @{cterm "op <::real => _"} @{cterm "0::real"})
                         (mk_numeric x))))
-      | Square t => square_rule t
-      | Eqmul(t,p) => emul_rule t (translate p)
+      | Square pt => square_rule (cterm_of_poly pt)
+      | Eqmul(pt,p) => emul_rule (cterm_of_poly pt) (translate p)
       | Sum(p1,p2) => add_rule (translate p1) (translate p2)
       | Product(p1,p2) => mul_rule (translate p1) (translate p2)
    in fconv_rule (first_conv [numeric_ge_conv, numeric_gt_conv, numeric_eq_conv, all_conv]) 
@@ -394,13 +477,13 @@
        val _ = if c aconvc (concl th2) then () else error "disj_cases : conclusions not alpha convertible"
    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)
    end
- fun overall dun ths = case ths of
+ fun overall cert_choice dun ths = case ths of
   [] =>
    let 
     val (eq,ne) = List.partition (is_req o concl) dun
      val (le,nl) = List.partition (is_ge o concl) ne
      val lt = filter (is_gt o concl) nl 
-    in prover hol_of_positivstellensatz (eq,le,lt) end
+    in prover (rev cert_choice) hol_of_positivstellensatz (eq,le,lt) end
  | th::oths =>
    let 
     val ct = concl th 
@@ -408,13 +491,13 @@
     if is_conj ct  then
      let 
       val (th1,th2) = conj_pair th in
-      overall dun (th1::th2::oths) end
+      overall cert_choice dun (th1::th2::oths) end
     else if is_disj ct then
       let 
-       val th1 = overall dun (assume (capply @{cterm Trueprop} (dest_arg1 ct))::oths)
-       val th2 = overall dun (assume (capply @{cterm Trueprop} (dest_arg ct))::oths)
-      in disj_cases th th1 th2 end
-   else overall (th::dun) oths
+       val (th1, cert1) = overall (Left::cert_choice) dun (assume (capply @{cterm Trueprop} (dest_arg1 ct))::oths)
+       val (th2, cert2) = overall (Right::cert_choice) dun (assume (capply @{cterm Trueprop} (dest_arg ct))::oths)
+      in (disj_cases th th1 th2, Branch (cert1, cert2)) end
+   else overall cert_choice (th::dun) oths
   end
   fun dest_binary b ct = if is_binop b ct then dest_binop ct 
                          else raise CTERM ("dest_binary",[b,ct])
@@ -496,16 +579,16 @@
   val nct = capply @{cterm Trueprop} (capply @{cterm "Not"} ct)
   val th0 = (init_conv then_conv arg_conv nnf_norm_conv') nct
   val tm0 = dest_arg (rhs_of th0)
-  val th = if tm0 aconvc @{cterm False} then equal_implies_1_rule th0 else
+  val (th, certificates) = if tm0 aconvc @{cterm False} then (equal_implies_1_rule th0, Trivial) else
    let 
     val (evs,bod) = strip_exists tm0
     val (avs,ibod) = strip_forall bod
     val th1 = Drule.arg_cong_rule @{cterm Trueprop} (fold mk_forall avs (absremover ibod))
-    val th2 = overall [] [specl avs (assume (rhs_of th1))]
+    val (th2, certs) = overall [] [] [specl avs (assume (rhs_of th1))]
     val th3 = fold simple_choose evs (prove_hyp (equal_elim th1 (assume (capply @{cterm Trueprop} bod))) th2)
-   in  Drule.implies_intr_hyps (prove_hyp (equal_elim th0 (assume nct)) th3)
+   in (Drule.implies_intr_hyps (prove_hyp (equal_elim th0 (assume nct)) th3), certs)
    end
-  in implies_elim (instantiate' [] [SOME ct] pth_final) th
+  in (implies_elim (instantiate' [] [SOME ct] pth_final) th, certificates)
  end
 in f
 end;
@@ -580,7 +663,7 @@
          val k = (Rat.neg d) */ Rat.abs c // c
          val e' = linear_cmul k e
          val t' = linear_cmul (Rat.abs c) t
-         val p' = Eqmul(cterm_of_rat k,p)
+         val p' = Eqmul(Monomialfunc.onefunc (Ctermfunc.undefined, k),p)
          val q' = Product(Rational_lt(Rat.abs c),q) 
         in (linear_add e' t',Sum(p',q')) 
         end 
@@ -632,7 +715,7 @@
   val le_pols' = le_pols @ map (fn v => Ctermfunc.onefunc (v,Rat.one)) aliens
   val (_,proof) = linear_prover (eq_pols,le_pols',lt_pols)
   val le' = le @ map (fn a => instantiate' [] [SOME (dest_arg a)] @{thm real_of_nat_ge_zero}) aliens 
- in (translator (eq,le',lt) proof) : thm
+ in ((translator (eq,le',lt) proof), Trivial)
  end
 end;
 
@@ -698,5 +781,4 @@
     main,neg,add,mul, prover)
 end;
 
-fun real_arith ctxt = gen_prover_real_arith ctxt real_linear_prover;
 end