--- 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
+
--- /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
+
+
--- 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
--- 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;
--- 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
--- 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