# HG changeset patch # User huffman # Date 1294514965 28800 # Node ID 9908cf4af3946b963e959cf5df4d640c6ea7e9e1 # Parent 655f583840d02671e0fa57ebb4bf7a0021d4cc6c# Parent fe4f0d9f9dbb512df1d709eaaec7c596e09a4618 merged diff -r 655f583840d0 -r 9908cf4af394 etc/components --- a/etc/components Sat Jan 08 11:18:09 2011 -0800 +++ b/etc/components Sat Jan 08 11:29:25 2011 -0800 @@ -14,7 +14,7 @@ src/Tools/WWW_Find src/HOL/Tools/ATP src/HOL/Mirabelle -src/HOL/Library/Sum_Of_Squares +src/HOL/Library/Sum_of_Squares src/HOL/Tools/SMT src/HOL/Tools/Predicate_Compile src/HOL/Mutabelle diff -r 655f583840d0 -r 9908cf4af394 src/HOL/Boogie/Tools/boogie_vcs.ML --- a/src/HOL/Boogie/Tools/boogie_vcs.ML Sat Jan 08 11:18:09 2011 -0800 +++ b/src/HOL/Boogie/Tools/boogie_vcs.ML Sat Jan 08 11:29:25 2011 -0800 @@ -254,7 +254,7 @@ end -(* the VC store *) +(* the VC store *) (* FIXME just one data slot (record) per program unit *) fun err_unfinished () = error "An unfinished Boogie environment is still open." @@ -278,7 +278,7 @@ type T = (serial * (term -> bool)) Ord_List.T val empty = [] val extend = I - fun merge (fs1, fs2) = fold (Ord_List.insert serial_ord) fs2 fs1 + fun merge data = Ord_List.merge serial_ord data ) fun add_assertion_filter f = diff -r 655f583840d0 -r 9908cf4af394 src/HOL/IsaMakefile --- a/src/HOL/IsaMakefile Sat Jan 08 11:18:09 2011 -0800 +++ b/src/HOL/IsaMakefile Sat Jan 08 11:29:25 2011 -0800 @@ -460,9 +460,9 @@ Library/RBT.thy Library/RBT_Impl.thy Library/README.html \ Library/Set_Algebras.thy Library/State_Monad.thy Library/Ramsey.thy \ Library/Reflection.thy Library/SML_Quickcheck.thy \ - Library/Sublist_Order.thy Library/Sum_Of_Squares.thy \ - Library/Sum_Of_Squares/sos_wrapper.ML \ - Library/Sum_Of_Squares/sum_of_squares.ML \ + Library/Sublist_Order.thy Library/Sum_of_Squares.thy \ + Library/Sum_of_Squares/sos_wrapper.ML \ + Library/Sum_of_Squares/sum_of_squares.ML \ Library/Transitive_Closure_Table.thy Library/Univ_Poly.thy \ Library/While_Combinator.thy Library/Zorn.thy \ $(SRC)/Tools/adhoc_overloading.ML Library/positivstellensatz.ML \ diff -r 655f583840d0 -r 9908cf4af394 src/HOL/Library/Eval_Witness.thy --- a/src/HOL/Library/Eval_Witness.thy Sat Jan 08 11:18:09 2011 -0800 +++ b/src/HOL/Library/Eval_Witness.thy Sat Jan 08 11:29:25 2011 -0800 @@ -45,8 +45,10 @@ instance list :: (ml_equiv) ml_equiv .. ML {* -structure Eval_Method = Proof_Data ( +structure Eval_Method = Proof_Data +( type T = unit -> bool + (* FIXME avoid user error with non-user text *) fun init _ () = error "Eval_Method" ) *} diff -r 655f583840d0 -r 9908cf4af394 src/HOL/Library/Library.thy --- a/src/HOL/Library/Library.thy Sat Jan 08 11:18:09 2011 -0800 +++ b/src/HOL/Library/Library.thy Sat Jan 08 11:29:25 2011 -0800 @@ -56,7 +56,7 @@ Set_Algebras SML_Quickcheck State_Monad - Sum_Of_Squares + Sum_of_Squares Transitive_Closure_Table Univ_Poly While_Combinator diff -r 655f583840d0 -r 9908cf4af394 src/HOL/Library/Sum_Of_Squares.thy --- a/src/HOL/Library/Sum_Of_Squares.thy Sat Jan 08 11:18:09 2011 -0800 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,159 +0,0 @@ -(* Title: HOL/Library/Sum_Of_Squares.thy - Author: Amine Chaieb, University of Cambridge - Author: Philipp Meyer, TU Muenchen -*) - -header {* A decision method for universal multivariate real arithmetic with addition, - multiplication and ordering using semidefinite programming *} - -theory Sum_Of_Squares -imports Complex_Main -uses - "positivstellensatz.ML" - "Sum_Of_Squares/sum_of_squares.ML" - "Sum_Of_Squares/positivstellensatz_tools.ML" - "Sum_Of_Squares/sos_wrapper.ML" -begin - -text {* - In order to use the method sos, call it with @{text "(sos - remote_csdp)"} to use the remote solver. Or install CSDP - (https://projects.coin-or.org/Csdp), configure the Isabelle setting - @{text CSDP_EXE}, and call it with @{text "(sos csdp)"}. By - default, sos calls @{text remote_csdp}. This can take of the order - 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. -*} - -setup Sum_Of_Squares.setup -setup SOS_Wrapper.setup - -text {* Tests *} - -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 "a1 >= 0 & a2 >= 0 \ (a1 * a1 + a2 * a2 = b1 * b1 + b2 * b2 + 2) \ (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 --> 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 "(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 "(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 "(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 "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 "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)))))))))))))") - -(* ------------------------------------------------------------------------- *) -(* Over a larger but simpler interval. *) -(* ------------------------------------------------------------------------- *) - -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)))))))))))") - -(* ------------------------------------------------------------------------- *) -(* 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))))))))))))))))))") - -(* ------------------------------------------------------------------------- *) -(* 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))))))") - -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 "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 \ 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))))))") - -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 "(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))))") - -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 "(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))))))") - - - -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 \ 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 --> 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 "((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_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_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_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_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 655f583840d0 -r 9908cf4af394 src/HOL/Library/Sum_Of_Squares/etc/settings --- a/src/HOL/Library/Sum_Of_Squares/etc/settings Sat Jan 08 11:18:09 2011 -0800 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,3 +0,0 @@ -# -*- shell-script -*- :mode=shellscript: - -ISABELLE_SUM_OF_SQUARES="$COMPONENT" diff -r 655f583840d0 -r 9908cf4af394 src/HOL/Library/Sum_Of_Squares/neos_csdp_client --- a/src/HOL/Library/Sum_Of_Squares/neos_csdp_client Sat Jan 08 11:18:09 2011 -0800 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,87 +0,0 @@ -#!/usr/bin/env python -import sys -import signal -import xmlrpclib -import time -import re - -# Neos server config -NEOS_HOST="neos.mcs.anl.gov" -NEOS_PORT=3332 - -neos=xmlrpclib.Server("http://%s:%d" % (NEOS_HOST, NEOS_PORT)) - -jobNumber = 0 -password = "" -inputfile = None -outputfile = None -# interrupt handler -def cleanup(signal, frame): - sys.stdout.write("Caught interrupt, cleaning up\n") - if jobNumber != 0: - neos.killJob(jobNumber, password) - if inputfile != None: - inputfile.close() - if outputfile != None: - outputfile.close() - sys.exit(21) - -signal.signal(signal.SIGHUP, cleanup) -signal.signal(signal.SIGINT, cleanup) -signal.signal(signal.SIGQUIT, cleanup) -signal.signal(signal.SIGTERM, cleanup) - -if len(sys.argv) <> 3: - sys.stderr.write("Usage: neos_csdp_client \n") - sys.exit(19) - -xml_pre = "\nsdp\ncsdp\nSPARSE_SDPA\n\n\n" -xml = xml_pre -inputfile = open(sys.argv[1],"r") -buffer = 1 -while buffer: - buffer = inputfile.read() - xml += buffer -inputfile.close() -xml += xml_post - -(jobNumber,password) = neos.submitJob(xml) - -if jobNumber == 0: - sys.stdout.write("error submitting job: %s" % password) - sys.exit(20) -else: - sys.stdout.write("jobNumber = %d\tpassword = %s\n" % (jobNumber,password)) - -offset=0 -messages = "" -status="Waiting" -while status == "Running" or status=="Waiting": - time.sleep(1) - (msg,offset) = neos.getIntermediateResults(jobNumber,password,offset) - messages += msg.data - sys.stdout.write(msg.data) - status = neos.getJobStatus(jobNumber, password) - -msg = neos.getFinalResults(jobNumber, password).data -sys.stdout.write("---------- Begin CSDP Output -------------\n"); -sys.stdout.write(msg) - -# extract solution -result = msg.split("Solution:") -if len(result) > 1: - solution = result[1].strip() - if solution != "": - outputfile = open(sys.argv[2],"w") - outputfile.write(solution) - outputfile.close() - -# extract return code -p = re.compile(r"^Error: Command exited with non-zero status (\d+)$", re.MULTILINE) -m = p.search(messages) -if m: - sys.exit(int(m.group(1))) -else: - sys.exit(0) - diff -r 655f583840d0 -r 9908cf4af394 src/HOL/Library/Sum_Of_Squares/positivstellensatz_tools.ML --- a/src/HOL/Library/Sum_Of_Squares/positivstellensatz_tools.ML Sat Jan 08 11:18:09 2011 -0800 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,155 +0,0 @@ -(* Title: HOL/Library/Sum_Of_Squares/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 - -(*** 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 FuncUtil.Ctermfunc.is_empty m then "1" - else - let - val m' = FuncUtil.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 FuncUtil.Ctermfunc.is_empty 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 FuncUtil.Monomialfunc.is_empty p then "0" - else - let - val cms = map string_of_cmonomial - (sort (prod_ord FuncUtil.monomial_order (K EQUAL)) (FuncUtil.Monomialfunc.dest p)) - in foldr1 (fn (t1, t2) => t1 ^ " + " ^ t2) cms - end; - -fun pss_to_cert (RealArith.Axiom_eq i) = "A=" ^ string_of_int i - | pss_to_cert (RealArith.Axiom_le i) = "A<=" ^ string_of_int i - | pss_to_cert (RealArith.Axiom_lt i) = "A<" ^ string_of_int i - | pss_to_cert (RealArith.Rational_eq r) = "R=" ^ string_of_rat r - | pss_to_cert (RealArith.Rational_le r) = "R<=" ^ string_of_rat r - | pss_to_cert (RealArith.Rational_lt r) = "R<" ^ string_of_rat r - | pss_to_cert (RealArith.Square p) = "[" ^ string_of_poly p ^ "]^2" - | pss_to_cert (RealArith.Eqmul (p, pss)) = "([" ^ string_of_poly p ^ "] * " ^ pss_to_cert pss ^ ")" - | pss_to_cert (RealArith.Sum (pss1, pss2)) = "(" ^ pss_to_cert pss1 ^ " + " ^ pss_to_cert pss2 ^ ")" - | pss_to_cert (RealArith.Product (pss1, pss2)) = "(" ^ pss_to_cert pss1 ^ " * " ^ pss_to_cert pss2 ^ ")" - -fun pss_tree_to_cert RealArith.Trivial = "()" - | pss_tree_to_cert (RealArith.Cert pss) = "(" ^ pss_to_cert pss ^ ")" - | pss_tree_to_cert (RealArith.Branch (t1, t2)) = "(" ^ pss_tree_to_cert t1 ^ " & " ^ pss_tree_to_cert t2 ^ ")" - -(*** certificate parsing ***) - -(* basic parser *) - -val str = Scan.this_string - -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 (str "~" >> K ~1) 1 -- nat >> op *; -val rat = int --| str "/" -- int >> Rat.rat_of_quotient -val rat_int = rat || int >> Rat.rat_of_int - -(* polynomial parser *) - -fun repeat_sep s f = f ::: Scan.repeat (str s |-- f) - -val parse_id = Scan.one Symbol.is_letter ::: Scan.many Symbol.is_letdig >> implode - -fun parse_varpow ctxt = parse_id -- Scan.optional (str "^" |-- nat) 1 >> - (fn (x, k) => (cterm_of (ProofContext.theory_of ctxt) (Free (x, @{typ real})), k)) - -fun parse_monomial ctxt = repeat_sep "*" (parse_varpow ctxt) >> - (fn xs => fold FuncUtil.Ctermfunc.update xs FuncUtil.Ctermfunc.empty) - -fun parse_cmonomial ctxt = - rat_int --| str "*" -- (parse_monomial ctxt) >> swap || - (parse_monomial ctxt) >> (fn m => (m, Rat.one)) || - rat_int >> (fn r => (FuncUtil.Ctermfunc.empty, r)) - -fun parse_poly ctxt = repeat_sep "+" (parse_cmonomial ctxt) >> - (fn xs => fold FuncUtil.Monomialfunc.update xs FuncUtil.Monomialfunc.empty) - -(* positivstellensatz parser *) - -val parse_axiom = - (str "A=" |-- int >> RealArith.Axiom_eq) || - (str "A<=" |-- int >> RealArith.Axiom_le) || - (str "A<" |-- int >> RealArith.Axiom_lt) - -val parse_rational = - (str "R=" |-- rat_int >> RealArith.Rational_eq) || - (str "R<=" |-- rat_int >> RealArith.Rational_le) || - (str "R<" |-- rat_int >> RealArith.Rational_lt) - -fun parse_cert ctxt input = - let - val pc = parse_cert ctxt - val pp = parse_poly ctxt - in - (parse_axiom || - parse_rational || - str "[" |-- pp --| str "]^2" >> RealArith.Square || - str "([" |-- pp --| str "]*" -- pc --| str ")" >> RealArith.Eqmul || - str "(" |-- pc --| str "*" -- pc --| str ")" >> RealArith.Product || - str "(" |-- pc --| str "+" -- pc --| str ")" >> RealArith.Sum) input - end - -fun parse_cert_tree ctxt input = - let - val pc = parse_cert ctxt - val pt = parse_cert_tree ctxt - in - (str "()" >> K RealArith.Trivial || - str "(" |-- pc --| str ")" >> RealArith.Cert || - str "(" |-- pt --| str "&" -- pt --| str ")" >> RealArith.Branch) input - end - -(* scanner *) - -fun cert_to_pss_tree ctxt input_str = Symbol.scanner "bad certificate" (parse_cert_tree ctxt) - (filter_out Symbol.is_blank (Symbol.explode input_str)) - -end - - diff -r 655f583840d0 -r 9908cf4af394 src/HOL/Library/Sum_Of_Squares/sos_wrapper.ML --- a/src/HOL/Library/Sum_Of_Squares/sos_wrapper.ML Sat Jan 08 11:18:09 2011 -0800 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,159 +0,0 @@ -(* Title: HOL/Library/Sum_Of_Squares/sos_wrapper.ML - Author: Philipp Meyer, TU Muenchen - -Added functionality for sums of squares, e.g. calling a remote prover. -*) - -signature SOS_WRAPPER = -sig - datatype prover_result = Success | Failure | Error - val setup: theory -> theory - val dest_dir: string Config.T - val prover_name: string Config.T -end - -structure SOS_Wrapper: SOS_WRAPPER = -struct - -datatype prover_result = Success | Failure | Error - -fun str_of_result Success = "Success" - | str_of_result Failure = "Failure" - | str_of_result Error = "Error" - - -(*** calling provers ***) - -val (dest_dir, setup_dest_dir) = Attrib.config_string "sos_dest_dir" (K "") - -fun filename dir name = - let - val probfile = Path.basic (name ^ serial_string ()) - val dir_path = Path.explode dir - in - if dir = "" then - File.tmp_path probfile - else if File.exists dir_path then - Path.append dir_path probfile - else error ("No such directory: " ^ dir) - end - -fun run_solver ctxt name cmd find_failure input = - let - val _ = warning ("Calling solver: " ^ name) - - (* create input file *) - val dir = Config.get ctxt dest_dir - val input_file = filename dir "sos_in" - val _ = File.write input_file input - - (* call solver *) - val output_file = filename dir "sos_out" - val (output, rv) = - bash_output - (if File.exists cmd then - space_implode " " - [File.shell_path cmd, File.shell_path input_file, File.shell_path output_file] - else error ("Bad executable: " ^ File.platform_path cmd)) - - (* read and analyze output *) - val (res, res_msg) = find_failure rv - val result = if File.exists output_file then File.read output_file else "" - - (* remove temporary files *) - val _ = - if dir = "" then - (File.rm input_file; if File.exists output_file then File.rm output_file else ()) - else () - - val _ = - if Config.get ctxt Sum_Of_Squares.trace - then writeln ("Solver output:\n" ^ output) - else () - - val _ = warning (str_of_result res ^ ": " ^ res_msg) - in - (case res of - Success => result - | Failure => raise Sum_Of_Squares.Failure res_msg - | Error => error ("Prover failed: " ^ res_msg)) - end - - -(*** various provers ***) - -(* local csdp client *) - -fun find_csdp_failure rv = - case rv of - 0 => (Success, "SDP solved") - | 1 => (Failure, "SDP is primal infeasible") - | 2 => (Failure, "SDP is dual infeasible") - | 3 => (Success, "SDP solved with reduced accuracy") - | 4 => (Failure, "Maximum iterations reached") - | 5 => (Failure, "Stuck at edge of primal feasibility") - | 6 => (Failure, "Stuck at edge of dual infeasibility") - | 7 => (Failure, "Lack of progress") - | 8 => (Failure, "X, Z, or O was singular") - | 9 => (Failure, "Detected NaN or Inf values") - | _ => (Error, "return code is " ^ string_of_int rv) - -val csdp = ("$CSDP_EXE", find_csdp_failure) - - -(* remote neos server *) - -fun find_neos_failure rv = - case rv of - 20 => (Error, "error submitting job") - | 21 => (Error, "interrupt") - | _ => find_csdp_failure rv - -val neos_csdp = ("$ISABELLE_SUM_OF_SQUARES/neos_csdp_client", find_neos_failure) - - -(* named provers *) - -fun prover "remote_csdp" = neos_csdp - | prover "csdp" = csdp - | prover name = error ("Unknown prover: " ^ name) - -val (prover_name, setup_prover_name) = Attrib.config_string "sos_prover_name" (K "remote_csdp") - -fun call_solver ctxt opt_name = - let - val name = the_default (Config.get ctxt prover_name) opt_name - val (cmd, find_failure) = prover name - in run_solver ctxt 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 ^ "\")") - -val print_cert = warning o output_line o PositivstellensatzTools.pss_tree_to_cert - - -(* method setup *) - -fun sos_solver print method = SIMPLE_METHOD' o Sum_Of_Squares.sos_tac print method - -val setup = - setup_dest_dir #> - setup_prover_name #> - Method.setup @{binding sos} - (Scan.lift (Scan.option Parse.xname) - >> (fn opt_name => fn ctxt => - sos_solver print_cert - (Sum_Of_Squares.Prover (call_solver ctxt opt_name)) ctxt)) - "prove universal problems over the reals using sums of squares" #> - Method.setup @{binding sos_cert} - (Scan.lift Parse.string - >> (fn cert => fn ctxt => - sos_solver ignore - (Sum_Of_Squares.Certificate (PositivstellensatzTools.cert_to_pss_tree ctxt cert)) ctxt)) - "prove universal problems over the reals using sums of squares with certificates" - -end diff -r 655f583840d0 -r 9908cf4af394 src/HOL/Library/Sum_Of_Squares/sum_of_squares.ML --- a/src/HOL/Library/Sum_Of_Squares/sum_of_squares.ML Sat Jan 08 11:18:09 2011 -0800 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,1435 +0,0 @@ -(* Title: HOL/Library/Sum_Of_Squares/sum_of_squares.ML - Author: Amine Chaieb, University of Cambridge - Author: Philipp Meyer, TU Muenchen - -A tactic for proving nonlinear inequalities. -*) - -signature SUM_OF_SQUARES = -sig - datatype proof_method = Certificate of RealArith.pss_tree | Prover of string -> string - val sos_tac: (RealArith.pss_tree -> unit) -> proof_method -> Proof.context -> int -> tactic - val trace: bool Config.T - val setup: theory -> theory - exception Failure of string; -end - -structure Sum_Of_Squares: SUM_OF_SQUARES = -struct - -val rat_0 = Rat.zero; -val rat_1 = Rat.one; -val rat_2 = Rat.two; -val rat_10 = Rat.rat_of_int 10; -val rat_1_2 = rat_1 // rat_2; -val max = Integer.max; -val min = Integer.min; - -val denominator_rat = Rat.quotient_of_rat #> snd #> Rat.rat_of_int; -val numerator_rat = Rat.quotient_of_rat #> fst #> Rat.rat_of_int; -fun int_of_rat a = - case Rat.quotient_of_rat a of (i,1) => i | _ => error "int_of_rat: not an int"; -fun lcm_rat x y = Rat.rat_of_int (Integer.lcm (int_of_rat x) (int_of_rat y)); - -fun rat_pow r i = - let fun pow r i = - if i = 0 then rat_1 else - let val d = pow r (i div 2) - in d */ d */ (if i mod 2 = 0 then rat_1 else r) - end - in if i < 0 then pow (Rat.inv r) (~ i) else pow r i end; - -fun round_rat r = - let val (a,b) = Rat.quotient_of_rat (Rat.abs r) - val d = a div b - val s = if r = b then d + 1 else d) end - -val abs_rat = Rat.abs; -val pow2 = rat_pow rat_2; -val pow10 = rat_pow rat_10; - -val (trace, setup_trace) = Attrib.config_bool "sos_trace" (K false); -val setup = setup_trace; - -exception Sanity; - -exception Unsolvable; - -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 -fun normalize y = - if abs_rat y =/ rat_1 then normalize (y // rat_10) + 1 - else 0 - in -fun decimalize d x = - if x =/ rat_0 then "0.0" else - let - val y = Rat.abs x - val e = normalize y - val z = pow10(~ e) */ y +/ rat_1 - val k = int_of_rat (round_rat(pow10 d */ z)) - in (if x a - | h::t => itern (k + 1) t f (f h k a); - -fun iter (m,n) f a = - if n < m then a - else iter (m+1,n) f (f m a); - -(* The main types. *) - -type vector = int* Rat.rat FuncUtil.Intfunc.table; - -type matrix = (int*int)*(Rat.rat FuncUtil.Intpairfunc.table); - -fun iszero (k,r) = r =/ rat_0; - - -(* Vectors. Conventionally indexed 1..n. *) - -fun vector_0 n = (n,FuncUtil.Intfunc.empty):vector; - -fun dim (v:vector) = fst v; - -fun vector_const c n = - if c =/ rat_0 then vector_0 n - else (n,fold_rev (fn k => FuncUtil.Intfunc.update (k,c)) (1 upto n) FuncUtil.Intfunc.empty) :vector; - -val vector_1 = vector_const rat_1; - -fun vector_cmul c (v:vector) = - let val n = dim v - in if c =/ rat_0 then vector_0 n - else (n,FuncUtil.Intfunc.map (fn _ => fn x => c */ x) (snd v)) - end; - -fun vector_neg (v:vector) = (fst v,FuncUtil.Intfunc.map (K Rat.neg) (snd v)) :vector; - -fun vector_add (v1:vector) (v2:vector) = - let val m = dim v1 - val n = dim v2 - in if m <> n then error "vector_add: incompatible dimensions" - else (n,FuncUtil.Intfunc.combine (curry op +/) (fn x => x =/ rat_0) (snd v1) (snd v2)) :vector - end; - -fun vector_sub v1 v2 = vector_add v1 (vector_neg v2); - -fun vector_dot (v1:vector) (v2:vector) = - let val m = dim v1 - val n = dim v2 - in if m <> n then error "vector_dot: incompatible dimensions" - else FuncUtil.Intfunc.fold (fn (i,x) => fn a => x +/ a) - (FuncUtil.Intfunc.combine (curry op */) (fn x => x =/ rat_0) (snd v1) (snd v2)) rat_0 - end; - -fun vector_of_list l = - let val n = length l - in (n,fold_rev2 (curry FuncUtil.Intfunc.update) (1 upto n) l FuncUtil.Intfunc.empty) :vector - end; - -(* Matrices; again rows and columns indexed from 1. *) - -fun matrix_0 (m,n) = ((m,n),FuncUtil.Intpairfunc.empty):matrix; - -fun dimensions (m:matrix) = fst m; - -fun matrix_const c (mn as (m,n)) = - if m <> n then error "matrix_const: needs to be square" - else if c =/ rat_0 then matrix_0 mn - else (mn,fold_rev (fn k => FuncUtil.Intpairfunc.update ((k,k), c)) (1 upto n) FuncUtil.Intpairfunc.empty) :matrix;; - -val matrix_1 = matrix_const rat_1; - -fun matrix_cmul c (m:matrix) = - let val (i,j) = dimensions m - in if c =/ rat_0 then matrix_0 (i,j) - else ((i,j),FuncUtil.Intpairfunc.map (fn _ => fn x => c */ x) (snd m)) - end; - -fun matrix_neg (m:matrix) = - (dimensions m, FuncUtil.Intpairfunc.map (K Rat.neg) (snd m)) :matrix; - -fun matrix_add (m1:matrix) (m2:matrix) = - let val d1 = dimensions m1 - val d2 = dimensions m2 - in if d1 <> d2 - then error "matrix_add: incompatible dimensions" - else (d1,FuncUtil.Intpairfunc.combine (curry op +/) (fn x => x =/ rat_0) (snd m1) (snd m2)) :matrix - end;; - -fun matrix_sub m1 m2 = matrix_add m1 (matrix_neg m2); - -fun row k (m:matrix) = - let val (i,j) = dimensions m - in (j, - FuncUtil.Intpairfunc.fold (fn ((i,j), c) => fn a => if i = k then FuncUtil.Intfunc.update (j,c) a else a) (snd m) FuncUtil.Intfunc.empty ) : vector - end; - -fun column k (m:matrix) = - let val (i,j) = dimensions m - in (i, - FuncUtil.Intpairfunc.fold (fn ((i,j), c) => fn a => if j = k then FuncUtil.Intfunc.update (i,c) a else a) (snd m) FuncUtil.Intfunc.empty) - : vector - end; - -fun transp (m:matrix) = - let val (i,j) = dimensions m - in - ((j,i),FuncUtil.Intpairfunc.fold (fn ((i,j), c) => fn a => FuncUtil.Intpairfunc.update ((j,i), c) a) (snd m) FuncUtil.Intpairfunc.empty) :matrix - end; - -fun diagonal (v:vector) = - let val n = dim v - in ((n,n),FuncUtil.Intfunc.fold (fn (i, c) => fn a => FuncUtil.Intpairfunc.update ((i,i), c) a) (snd v) FuncUtil.Intpairfunc.empty) : matrix - end; - -fun matrix_of_list l = - let val m = length l - in if m = 0 then matrix_0 (0,0) else - let val n = length (hd l) - in ((m,n),itern 1 l (fn v => fn i => itern 1 v (fn c => fn j => FuncUtil.Intpairfunc.update ((i,j), c))) FuncUtil.Intpairfunc.empty) - end - end; - -(* Monomials. *) - -fun monomial_eval assig m = - FuncUtil.Ctermfunc.fold (fn (x, k) => fn a => a */ rat_pow (FuncUtil.Ctermfunc.apply assig x) k) - m rat_1; -val monomial_1 = FuncUtil.Ctermfunc.empty; - -fun monomial_var x = FuncUtil.Ctermfunc.onefunc (x, 1); - -val monomial_mul = - FuncUtil.Ctermfunc.combine Integer.add (K false); - -fun monomial_pow m k = - if k = 0 then monomial_1 - else FuncUtil.Ctermfunc.map (fn _ => fn x => k * x) m; - -fun monomial_divides m1 m2 = - FuncUtil.Ctermfunc.fold (fn (x, k) => fn a => FuncUtil.Ctermfunc.tryapplyd m2 x 0 >= k andalso a) m1 true;; - -fun monomial_div m1 m2 = - let val m = FuncUtil.Ctermfunc.combine Integer.add - (fn x => x = 0) m1 (FuncUtil.Ctermfunc.map (fn _ => fn x => ~ x) m2) - in if FuncUtil.Ctermfunc.fold (fn (x, k) => fn a => k >= 0 andalso a) m true then m - else error "monomial_div: non-divisible" - end; - -fun monomial_degree x m = - FuncUtil.Ctermfunc.tryapplyd m x 0;; - -fun monomial_lcm m1 m2 = - fold_rev (fn x => FuncUtil.Ctermfunc.update (x, max (monomial_degree x m1) (monomial_degree x m2))) - (union (is_equal o FuncUtil.cterm_ord) (FuncUtil.Ctermfunc.dom m1) (FuncUtil.Ctermfunc.dom m2)) (FuncUtil.Ctermfunc.empty); - -fun monomial_multidegree m = - FuncUtil.Ctermfunc.fold (fn (x, k) => fn a => k + a) m 0;; - -fun monomial_variables m = FuncUtil.Ctermfunc.dom m;; - -(* Polynomials. *) - -fun eval assig p = - FuncUtil.Monomialfunc.fold (fn (m, c) => fn a => a +/ c */ monomial_eval assig m) p rat_0; - -val poly_0 = FuncUtil.Monomialfunc.empty; - -fun poly_isconst p = - FuncUtil.Monomialfunc.fold (fn (m, c) => fn a => FuncUtil.Ctermfunc.is_empty m andalso a) p true; - -fun poly_var x = FuncUtil.Monomialfunc.onefunc (monomial_var x,rat_1); - -fun poly_const c = - if c =/ rat_0 then poly_0 else FuncUtil.Monomialfunc.onefunc(monomial_1, c); - -fun poly_cmul c p = - if c =/ rat_0 then poly_0 - else FuncUtil.Monomialfunc.map (fn _ => fn x => c */ x) p; - -fun poly_neg p = FuncUtil.Monomialfunc.map (K Rat.neg) p;; - -fun poly_add p1 p2 = - FuncUtil.Monomialfunc.combine (curry op +/) (fn x => x =/ rat_0) p1 p2; - -fun poly_sub p1 p2 = poly_add p1 (poly_neg p2); - -fun poly_cmmul (c,m) p = - if c =/ rat_0 then poly_0 - else if FuncUtil.Ctermfunc.is_empty m - then FuncUtil.Monomialfunc.map (fn _ => fn d => c */ d) p - else FuncUtil.Monomialfunc.fold (fn (m', d) => fn a => (FuncUtil.Monomialfunc.update (monomial_mul m m', c */ d) a)) p poly_0; - -fun poly_mul p1 p2 = - FuncUtil.Monomialfunc.fold (fn (m, c) => fn a => poly_add (poly_cmmul (c,m) p2) a) p1 poly_0; - -fun poly_div p1 p2 = - if not(poly_isconst p2) - then error "poly_div: non-constant" else - let val c = eval FuncUtil.Ctermfunc.empty p2 - in if c =/ rat_0 then error "poly_div: division by zero" - else poly_cmul (Rat.inv c) p1 - end; - -fun poly_square p = poly_mul p p; - -fun poly_pow p k = - if k = 0 then poly_const rat_1 - else if k = 1 then p - else let val q = poly_square(poly_pow p (k div 2)) in - if k mod 2 = 1 then poly_mul p q else q end; - -fun poly_exp p1 p2 = - if not(poly_isconst p2) - then error "poly_exp: not a constant" - else poly_pow p1 (int_of_rat (eval FuncUtil.Ctermfunc.empty p2)); - -fun degree x p = - FuncUtil.Monomialfunc.fold (fn (m,c) => fn a => max (monomial_degree x m) a) p 0; - -fun multidegree p = - FuncUtil.Monomialfunc.fold (fn (m, c) => fn a => max (monomial_multidegree m) a) p 0; - -fun poly_variables p = - sort FuncUtil.cterm_ord (FuncUtil.Monomialfunc.fold_rev (fn (m, c) => union (is_equal o FuncUtil.cterm_ord) (monomial_variables m)) p []);; - -(* Order monomials for human presentation. *) - -val humanorder_varpow = prod_ord FuncUtil.cterm_ord (rev_order o int_ord); - -local - fun ord (l1,l2) = case (l1,l2) of - (_,[]) => LESS - | ([],_) => GREATER - | (h1::t1,h2::t2) => - (case humanorder_varpow (h1, h2) of - LESS => LESS - | EQUAL => ord (t1,t2) - | GREATER => GREATER) -in fun humanorder_monomial m1 m2 = - ord (sort humanorder_varpow (FuncUtil.Ctermfunc.dest m1), - sort humanorder_varpow (FuncUtil.Ctermfunc.dest m2)) -end; - -(* Conversions to strings. *) - -fun string_of_vector min_size max_size (v:vector) = - let val n_raw = dim v - in if n_raw = 0 then "[]" else - let - val n = max min_size (min n_raw max_size) - val xs = map (Rat.string_of_rat o (fn i => FuncUtil.Intfunc.tryapplyd (snd v) i rat_0)) (1 upto n) - in "[" ^ space_implode ", " xs ^ - (if n_raw > max_size then ", ...]" else "]") - end - end; - -fun string_of_matrix max_size (m:matrix) = - let - val (i_raw,j_raw) = dimensions m - val i = min max_size i_raw - val j = min max_size j_raw - val rstr = map (fn k => string_of_vector j j (row k m)) (1 upto i) - in "["^ space_implode ";\n " rstr ^ - (if j > max_size then "\n ...]" else "]") - end; - -fun string_of_term t = - case t of - a$b => "("^(string_of_term a)^" "^(string_of_term b)^")" - | Abs x => - let val (xn, b) = Term.dest_abs x - in "(\\"^xn^"."^(string_of_term b)^")" - end - | Const(s,_) => s - | Free (s,_) => s - | Var((s,_),_) => s - | _ => error "string_of_term"; - -val string_of_cterm = string_of_term o term_of; - -fun string_of_varpow x k = - if k = 1 then string_of_cterm x - else string_of_cterm x^"^"^string_of_int k; - -fun string_of_monomial m = - if FuncUtil.Ctermfunc.is_empty m then "1" else - let val vps = fold_rev (fn (x,k) => fn a => string_of_varpow x k :: a) - (sort humanorder_varpow (FuncUtil.Ctermfunc.dest m)) [] - in space_implode "*" vps - end; - -fun string_of_cmonomial (c,m) = - if FuncUtil.Ctermfunc.is_empty m then Rat.string_of_rat c - 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 = - if FuncUtil.Monomialfunc.is_empty p then "<<0>>" else - let - val cms = sort (fn ((m1,_),(m2,_)) => humanorder_monomial m1 m2) (FuncUtil.Monomialfunc.dest p) - val s = fold (fn (m,c) => fn a => - if c >" - end; - -(* Conversion from HOL term. *) - -local - val neg_tm = @{cterm "uminus :: real => _"} - val add_tm = @{cterm "op + :: real => _"} - val sub_tm = @{cterm "op - :: real => _"} - val mul_tm = @{cterm "op * :: real => _"} - val inv_tm = @{cterm "inverse :: real => _"} - val div_tm = @{cterm "op / :: real => _"} - val pow_tm = @{cterm "op ^ :: real => _"} - val zero_tm = @{cterm "0:: real"} - val is_numeral = can (HOLogic.dest_number o term_of) - fun is_comb t = case t of _$_ => true | _ => false - fun poly_of_term tm = - if tm aconvc zero_tm then poly_0 - else if RealArith.is_ratconst tm - then poly_const(RealArith.dest_ratconst tm) - else - (let val (lop,r) = Thm.dest_comb tm - in if lop aconvc neg_tm then poly_neg(poly_of_term r) - else if lop aconvc inv_tm then - let val p = poly_of_term r - in if poly_isconst p - then poly_const(Rat.inv (eval FuncUtil.Ctermfunc.empty p)) - else error "poly_of_term: inverse of non-constant polyomial" - end - else (let val (opr,l) = Thm.dest_comb lop - in - if opr aconvc pow_tm andalso is_numeral r - then poly_pow (poly_of_term l) ((snd o HOLogic.dest_number o term_of) r) - else if opr aconvc add_tm - then poly_add (poly_of_term l) (poly_of_term r) - else if opr aconvc sub_tm - then poly_sub (poly_of_term l) (poly_of_term r) - else if opr aconvc mul_tm - then poly_mul (poly_of_term l) (poly_of_term r) - else if opr aconvc div_tm - then let - val p = poly_of_term l - val q = poly_of_term r - in if poly_isconst q then poly_cmul (Rat.inv (eval FuncUtil.Ctermfunc.empty q)) p - else error "poly_of_term: division by non-constant polynomial" - end - else poly_var tm - - end - handle CTERM ("dest_comb",_) => poly_var tm) - end - handle CTERM ("dest_comb",_) => poly_var tm) -in -val poly_of_term = fn tm => - if type_of (term_of tm) = @{typ real} then poly_of_term tm - else error "poly_of_term: term does not have real type" -end; - -(* String of vector (just a list of space-separated numbers). *) - -fun sdpa_of_vector (v:vector) = - let - val n = dim v - val strs = map (decimalize 20 o (fn i => FuncUtil.Intfunc.tryapplyd (snd v) i rat_0)) (1 upto n) - in space_implode " " strs ^ "\n" - end; - -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'))); -structure Inttriplefunc = FuncFun(type key = int*int*int val ord = triple_int_ord); - -(* String for block diagonal matrix numbered k. *) - -fun sdpa_of_blockdiagonal k m = - let - val pfx = string_of_int k ^" " - val ents = - Inttriplefunc.fold (fn ((b,i,j), c) => fn a => if i > j then a else ((b,i,j),c)::a) m [] - val entss = sort (triple_int_ord o pairself fst) ents - in fold_rev (fn ((b,i,j),c) => fn a => - pfx ^ string_of_int b ^ " " ^ string_of_int i ^ " " ^ string_of_int j ^ - " " ^ decimalize 20 c ^ "\n" ^ a) entss "" - end; - -(* String for a matrix numbered k, in SDPA sparse format. *) - -fun sdpa_of_matrix k (m:matrix) = - let - val pfx = string_of_int k ^ " 1 " - val ms = FuncUtil.Intpairfunc.fold (fn ((i,j), c) => fn a => if i > j then a else ((i,j),c)::a) (snd m) [] - val mss = sort ((prod_ord int_ord int_ord) o pairself fst) ms - in fold_rev (fn ((i,j),c) => fn a => - pfx ^ string_of_int i ^ " " ^ string_of_int j ^ - " " ^ decimalize 20 c ^ "\n" ^ a) mss "" - end;; - -(* ------------------------------------------------------------------------- *) -(* String in SDPA sparse format for standard SDP problem: *) -(* *) -(* X = v_1 * [M_1] + ... + v_m * [M_m] - [M_0] must be PSD *) -(* Minimize obj_1 * v_1 + ... obj_m * v_m *) -(* ------------------------------------------------------------------------- *) - -fun sdpa_of_problem obj mats = - let - val m = length mats - 1 - val (n,_) = dimensions (hd mats) - in - string_of_int m ^ "\n" ^ - "1\n" ^ - string_of_int n ^ "\n" ^ - sdpa_of_vector obj ^ - fold_rev2 (fn k => fn m => fn a => sdpa_of_matrix (k - 1) m ^ a) (1 upto length mats) mats "" - end; - -fun index_char str chr pos = - if pos >= String.size str then ~1 - else if String.sub(str,pos) = chr then pos - else index_char str chr (pos + 1); -fun rat_of_quotient (a,b) = if b = 0 then rat_0 else Rat.rat_of_quotient (a,b); -fun rat_of_string s = - let val n = index_char s #"/" 0 in - if n = ~1 then s |> Int.fromString |> the |> Rat.rat_of_int - else - let val SOME numer = Int.fromString(String.substring(s,0,n)) - val SOME den = Int.fromString (String.substring(s,n+1,String.size s - n - 1)) - in rat_of_quotient(numer, den) - end - end; - -fun isspace x = (x = " "); -fun isnum x = member (op =) ["0","1","2","3","4","5","6","7","8","9"] x; - -(* More parser basics. *) - - val word = Scan.this_string - fun token s = - Scan.repeat ($$ " ") |-- word s --| Scan.repeat ($$ " ") - val numeral = Scan.one isnum - val decimalint = Scan.repeat1 numeral >> (rat_of_string o implode) - val decimalfrac = Scan.repeat1 numeral - >> (fn s => rat_of_string(implode s) // pow10 (length s)) - val decimalsig = - decimalint -- Scan.option (Scan.$$ "." |-- decimalfrac) - >> (fn (h,NONE) => h | (h,SOME x) => h +/ x) - fun signed prs = - $$ "-" |-- prs >> Rat.neg - || $$ "+" |-- prs - || prs; - -fun emptyin def xs = if null xs then (def,xs) else Scan.fail xs - - val exponent = ($$ "e" || $$ "E") |-- signed decimalint; - - val decimal = signed decimalsig -- (emptyin rat_0|| exponent) - >> (fn (h, x) => h */ pow10 (int_of_rat x)); - - fun mkparser p s = - let val (x,rst) = p (raw_explode s) - in if null rst then x - else error "mkparser: unparsed input" - end;; - -(* Parse back csdp output. *) - - fun ignore inp = ((),[]) - fun csdpoutput inp = - ((decimal -- Scan.repeat (Scan.$$ " " |-- Scan.option decimal) >> - (fn (h,to) => map_filter I ((SOME h)::to))) --| ignore >> vector_of_list) inp - val parse_csdpoutput = mkparser csdpoutput - -(* Run prover on a problem in linear form. *) - -fun run_problem prover obj mats = - parse_csdpoutput (prover (sdpa_of_problem obj mats)) - -(* Try some apparently sensible scaling first. Note that this is purely to *) -(* get a cleaner translation to floating-point, and doesn't affect any of *) -(* the results, in principle. In practice it seems a lot better when there *) -(* are extreme numbers in the original problem. *) - - (* Version for (int*int) keys *) -local - fun max_rat x y = if x fn a => lcm_rat (denominator_rat c) a) amat acc - fun maximal_element fld amat acc = - fld (fn (m,c) => fn maxa => max_rat maxa (abs_rat c)) amat acc -fun float_of_rat x = let val (a,b) = Rat.quotient_of_rat x - in Real.fromLargeInt a / Real.fromLargeInt b end; -in - -fun pi_scale_then solver (obj:vector) mats = - let - val cd1 = fold_rev (common_denominator FuncUtil.Intpairfunc.fold) mats (rat_1) - val cd2 = common_denominator FuncUtil.Intfunc.fold (snd obj) (rat_1) - val mats' = map (FuncUtil.Intpairfunc.map (fn _ => fn x => cd1 */ x)) mats - val obj' = vector_cmul cd2 obj - val max1 = fold_rev (maximal_element FuncUtil.Intpairfunc.fold) mats' (rat_0) - val max2 = maximal_element FuncUtil.Intfunc.fold (snd obj') (rat_0) - val scal1 = pow2 (20 - trunc(Math.ln (float_of_rat max1) / Math.ln 2.0)) - val scal2 = pow2 (20 - trunc(Math.ln (float_of_rat max2) / Math.ln 2.0)) - val mats'' = map (FuncUtil.Intpairfunc.map (fn _ => fn x => x */ scal1)) mats' - val obj'' = vector_cmul scal2 obj' - in solver obj'' mats'' - end -end; - -(* Try some apparently sensible scaling first. Note that this is purely to *) -(* get a cleaner translation to floating-point, and doesn't affect any of *) -(* the results, in principle. In practice it seems a lot better when there *) -(* are extreme numbers in the original problem. *) - - (* Version for (int*int*int) keys *) -local - fun max_rat x y = if x fn a => lcm_rat (denominator_rat c) a) amat acc - fun maximal_element fld amat acc = - fld (fn (m,c) => fn maxa => max_rat maxa (abs_rat c)) amat acc -fun float_of_rat x = let val (a,b) = Rat.quotient_of_rat x - in Real.fromLargeInt a / Real.fromLargeInt b end; -fun int_of_float x = (trunc x handle Overflow => 0 | Domain => 0) -in - -fun tri_scale_then solver (obj:vector) mats = - let - val cd1 = fold_rev (common_denominator Inttriplefunc.fold) mats (rat_1) - val cd2 = common_denominator FuncUtil.Intfunc.fold (snd obj) (rat_1) - val mats' = map (Inttriplefunc.map (fn _ => fn x => cd1 */ x)) mats - val obj' = vector_cmul cd2 obj - val max1 = fold_rev (maximal_element Inttriplefunc.fold) mats' (rat_0) - val max2 = maximal_element FuncUtil.Intfunc.fold (snd obj') (rat_0) - val scal1 = pow2 (20 - int_of_float(Math.ln (float_of_rat max1) / Math.ln 2.0)) - val scal2 = pow2 (20 - int_of_float(Math.ln (float_of_rat max2) / Math.ln 2.0)) - val mats'' = map (Inttriplefunc.map (fn _ => fn x => x */ scal1)) mats' - val obj'' = vector_cmul scal2 obj' - in solver obj'' mats'' - end -end; - -(* Round a vector to "nice" rationals. *) - -fun nice_rational n x = round_rat (n */ x) // n;; -fun nice_vector n ((d,v) : vector) = - (d, FuncUtil.Intfunc.fold (fn (i,c) => fn a => - let val y = nice_rational n c - in if c =/ rat_0 then a - else FuncUtil.Intfunc.update (i,y) a end) v FuncUtil.Intfunc.empty):vector - -fun dest_ord f x = is_equal (f x); - -(* Stuff for "equations" ((int*int*int)->num functions). *) - -fun tri_equation_cmul c eq = - if c =/ rat_0 then Inttriplefunc.empty else Inttriplefunc.map (fn _ => fn d => c */ d) eq; - -fun tri_equation_add eq1 eq2 = Inttriplefunc.combine (curry op +/) (fn x => x =/ rat_0) eq1 eq2; - -fun tri_equation_eval assig eq = - let fun value v = Inttriplefunc.apply assig v - in Inttriplefunc.fold (fn (v, c) => fn a => a +/ value v */ c) eq rat_0 - end; - -(* Eliminate among linear equations: return unconstrained variables and *) -(* assignments for the others in terms of them. We give one pseudo-variable *) -(* "one" that's used for a constant term. *) - -local - fun extract_first p l = case l of (* FIXME : use find_first instead *) - [] => error "extract_first" - | h::t => if p h then (h,t) else - let val (k,s) = extract_first p t in (k,h::s) end -fun eliminate vars dun eqs = case vars of - [] => if forall Inttriplefunc.is_empty eqs then dun - else raise Unsolvable - | v::vs => - ((let - val (eq,oeqs) = extract_first (fn e => Inttriplefunc.defined e v) eqs - val a = Inttriplefunc.apply eq v - val eq' = tri_equation_cmul ((Rat.neg rat_1) // a) (Inttriplefunc.delete_safe v eq) - fun elim e = - let val b = Inttriplefunc.tryapplyd e v rat_0 - in if b =/ rat_0 then e else - tri_equation_add e (tri_equation_cmul (Rat.neg b // a) eq) - end - in eliminate vs (Inttriplefunc.update (v,eq') (Inttriplefunc.map (K elim) dun)) (map elim oeqs) - end) - handle Failure _ => eliminate vs dun eqs) -in -fun tri_eliminate_equations one vars eqs = - let - val assig = eliminate vars Inttriplefunc.empty eqs - val vs = Inttriplefunc.fold (fn (x, f) => fn a => remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig [] - in (distinct (dest_ord triple_int_ord) vs, assig) - end -end; - -(* Eliminate all variables, in an essentially arbitrary order. *) - -fun tri_eliminate_all_equations one = - let - fun choose_variable eq = - let val (v,_) = Inttriplefunc.choose eq - in if is_equal (triple_int_ord(v,one)) then - let val eq' = Inttriplefunc.delete_safe v eq - in if Inttriplefunc.is_empty eq' then error "choose_variable" - else fst (Inttriplefunc.choose eq') - end - else v - end - fun eliminate dun eqs = case eqs of - [] => dun - | eq::oeqs => - if Inttriplefunc.is_empty eq then eliminate dun oeqs else - let val v = choose_variable eq - val a = Inttriplefunc.apply eq v - val eq' = tri_equation_cmul ((Rat.rat_of_int ~1) // a) - (Inttriplefunc.delete_safe v eq) - fun elim e = - let val b = Inttriplefunc.tryapplyd e v rat_0 - in if b =/ rat_0 then e - else tri_equation_add e (tri_equation_cmul (Rat.neg b // a) eq) - end - in eliminate (Inttriplefunc.update(v, eq') (Inttriplefunc.map (K elim) dun)) - (map elim oeqs) - end -in fn eqs => - let - val assig = eliminate Inttriplefunc.empty eqs - val vs = Inttriplefunc.fold (fn (x, f) => fn a => remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig [] - in (distinct (dest_ord triple_int_ord) vs,assig) - end -end; - -(* Solve equations by assigning arbitrary numbers. *) - -fun tri_solve_equations one eqs = - let - val (vars,assigs) = tri_eliminate_all_equations one eqs - val vfn = fold_rev (fn v => Inttriplefunc.update(v,rat_0)) vars - (Inttriplefunc.onefunc(one, Rat.rat_of_int ~1)) - val ass = - Inttriplefunc.combine (curry op +/) (K false) - (Inttriplefunc.map (K (tri_equation_eval vfn)) assigs) vfn - in if forall (fn e => tri_equation_eval ass e =/ rat_0) eqs - then Inttriplefunc.delete_safe one ass else raise Sanity - end; - -(* Multiply equation-parametrized poly by regular poly and add accumulator. *) - -fun tri_epoly_pmul p q acc = - FuncUtil.Monomialfunc.fold (fn (m1, c) => fn a => - FuncUtil.Monomialfunc.fold (fn (m2,e) => fn b => - let val m = monomial_mul m1 m2 - val es = FuncUtil.Monomialfunc.tryapplyd b m Inttriplefunc.empty - in FuncUtil.Monomialfunc.update (m,tri_equation_add (tri_equation_cmul c e) es) b - end) q a) p acc ; - -(* Usual operations on equation-parametrized poly. *) - -fun tri_epoly_cmul c l = - if c =/ rat_0 then Inttriplefunc.empty else Inttriplefunc.map (K (tri_equation_cmul c)) l;; - -val tri_epoly_neg = tri_epoly_cmul (Rat.rat_of_int ~1); - -val tri_epoly_add = Inttriplefunc.combine tri_equation_add Inttriplefunc.is_empty; - -fun tri_epoly_sub p q = tri_epoly_add p (tri_epoly_neg q);; - -(* Stuff for "equations" ((int*int)->num functions). *) - -fun pi_equation_cmul c eq = - if c =/ rat_0 then Inttriplefunc.empty else Inttriplefunc.map (fn _ => fn d => c */ d) eq; - -fun pi_equation_add eq1 eq2 = Inttriplefunc.combine (curry op +/) (fn x => x =/ rat_0) eq1 eq2; - -fun pi_equation_eval assig eq = - let fun value v = Inttriplefunc.apply assig v - in Inttriplefunc.fold (fn (v, c) => fn a => a +/ value v */ c) eq rat_0 - end; - -(* Eliminate among linear equations: return unconstrained variables and *) -(* assignments for the others in terms of them. We give one pseudo-variable *) -(* "one" that's used for a constant term. *) - -local -fun extract_first p l = case l of - [] => error "extract_first" - | h::t => if p h then (h,t) else - let val (k,s) = extract_first p t in (k,h::s) end -fun eliminate vars dun eqs = case vars of - [] => if forall Inttriplefunc.is_empty eqs then dun - else raise Unsolvable - | v::vs => - let - val (eq,oeqs) = extract_first (fn e => Inttriplefunc.defined e v) eqs - val a = Inttriplefunc.apply eq v - val eq' = pi_equation_cmul ((Rat.neg rat_1) // a) (Inttriplefunc.delete_safe v eq) - fun elim e = - let val b = Inttriplefunc.tryapplyd e v rat_0 - in if b =/ rat_0 then e else - pi_equation_add e (pi_equation_cmul (Rat.neg b // a) eq) - end - in eliminate vs (Inttriplefunc.update (v,eq') (Inttriplefunc.map (K elim) dun)) (map elim oeqs) - end - handle Failure _ => eliminate vs dun eqs -in -fun pi_eliminate_equations one vars eqs = - let - val assig = eliminate vars Inttriplefunc.empty eqs - val vs = Inttriplefunc.fold (fn (x, f) => fn a => remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig [] - in (distinct (dest_ord triple_int_ord) vs, assig) - end -end; - -(* Eliminate all variables, in an essentially arbitrary order. *) - -fun pi_eliminate_all_equations one = - let - fun choose_variable eq = - let val (v,_) = Inttriplefunc.choose eq - in if is_equal (triple_int_ord(v,one)) then - let val eq' = Inttriplefunc.delete_safe v eq - in if Inttriplefunc.is_empty eq' then error "choose_variable" - else fst (Inttriplefunc.choose eq') - end - else v - end - fun eliminate dun eqs = case eqs of - [] => dun - | eq::oeqs => - if Inttriplefunc.is_empty eq then eliminate dun oeqs else - let val v = choose_variable eq - val a = Inttriplefunc.apply eq v - val eq' = pi_equation_cmul ((Rat.rat_of_int ~1) // a) - (Inttriplefunc.delete_safe v eq) - fun elim e = - let val b = Inttriplefunc.tryapplyd e v rat_0 - in if b =/ rat_0 then e - else pi_equation_add e (pi_equation_cmul (Rat.neg b // a) eq) - end - in eliminate (Inttriplefunc.update(v, eq') (Inttriplefunc.map (K elim) dun)) - (map elim oeqs) - end -in fn eqs => - let - val assig = eliminate Inttriplefunc.empty eqs - val vs = Inttriplefunc.fold (fn (x, f) => fn a => remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig [] - in (distinct (dest_ord triple_int_ord) vs,assig) - end -end; - -(* Solve equations by assigning arbitrary numbers. *) - -fun pi_solve_equations one eqs = - let - val (vars,assigs) = pi_eliminate_all_equations one eqs - val vfn = fold_rev (fn v => Inttriplefunc.update(v,rat_0)) vars - (Inttriplefunc.onefunc(one, Rat.rat_of_int ~1)) - val ass = - Inttriplefunc.combine (curry op +/) (K false) - (Inttriplefunc.map (K (pi_equation_eval vfn)) assigs) vfn - in if forall (fn e => pi_equation_eval ass e =/ rat_0) eqs - then Inttriplefunc.delete_safe one ass else raise Sanity - end; - -(* Multiply equation-parametrized poly by regular poly and add accumulator. *) - -fun pi_epoly_pmul p q acc = - FuncUtil.Monomialfunc.fold (fn (m1, c) => fn a => - FuncUtil.Monomialfunc.fold (fn (m2,e) => fn b => - let val m = monomial_mul m1 m2 - val es = FuncUtil.Monomialfunc.tryapplyd b m Inttriplefunc.empty - in FuncUtil.Monomialfunc.update (m,pi_equation_add (pi_equation_cmul c e) es) b - end) q a) p acc ; - -(* Usual operations on equation-parametrized poly. *) - -fun pi_epoly_cmul c l = - if c =/ rat_0 then Inttriplefunc.empty else Inttriplefunc.map (K (pi_equation_cmul c)) l;; - -val pi_epoly_neg = pi_epoly_cmul (Rat.rat_of_int ~1); - -val pi_epoly_add = Inttriplefunc.combine pi_equation_add Inttriplefunc.is_empty; - -fun pi_epoly_sub p q = pi_epoly_add p (pi_epoly_neg q);; - -fun allpairs f l1 l2 = fold_rev (fn x => (curry (op @)) (map (f x) l2)) l1 []; - -(* Hence produce the "relevant" monomials: those whose squares lie in the *) -(* Newton polytope of the monomials in the input. (This is enough according *) -(* to Reznik: "Extremal PSD forms with few terms", Duke Math. Journal, *) -(* vol 45, pp. 363--374, 1978. *) -(* *) -(* These are ordered in sort of decreasing degree. In particular the *) -(* constant monomial is last; this gives an order in diagonalization of the *) -(* quadratic form that will tend to display constants. *) - -(* Diagonalize (Cholesky/LDU) the matrix corresponding to a quadratic form. *) - -local -fun diagonalize n i m = - if FuncUtil.Intpairfunc.is_empty (snd m) then [] - else - let val a11 = FuncUtil.Intpairfunc.tryapplyd (snd m) (i,i) rat_0 - in if a11 fn a => - let val y = c // a11 - in if y = rat_0 then a else FuncUtil.Intfunc.update (i,y) a - end) (snd v) FuncUtil.Intfunc.empty) - fun upt0 x y a = if y = rat_0 then a else FuncUtil.Intpairfunc.update (x,y) a - val m' = - ((n,n), - iter (i+1,n) (fn j => - iter (i+1,n) (fn k => - (upt0 (j,k) (FuncUtil.Intpairfunc.tryapplyd (snd m) (j,k) rat_0 -/ FuncUtil.Intfunc.tryapplyd (snd v) j rat_0 */ FuncUtil.Intfunc.tryapplyd (snd v') k rat_0)))) - FuncUtil.Intpairfunc.empty) - in (a11,v')::diagonalize n (i + 1) m' - end - end -in -fun diag m = - let - val nn = dimensions m - val n = fst nn - in if snd nn <> n then error "diagonalize: non-square matrix" - else diagonalize n 1 m - end -end; - -fun gcd_rat a b = Rat.rat_of_int (Integer.gcd (int_of_rat a) (int_of_rat b)); - -(* Adjust a diagonalization to collect rationals at the start. *) - (* FIXME : Potentially polymorphic keys, but here only: integers!! *) -local - fun upd0 x y a = if y =/ rat_0 then a else FuncUtil.Intfunc.update(x,y) a; - fun mapa f (d,v) = - (d, FuncUtil.Intfunc.fold (fn (i,c) => fn a => upd0 i (f c) a) v FuncUtil.Intfunc.empty) - fun adj (c,l) = - let val a = - FuncUtil.Intfunc.fold (fn (i,c) => fn a => lcm_rat a (denominator_rat c)) - (snd l) rat_1 // - FuncUtil.Intfunc.fold (fn (i,c) => fn a => gcd_rat a (numerator_rat c)) - (snd l) rat_0 - in ((c // (a */ a)),mapa (fn x => a */ x) l) - end -in -fun deration d = if null d then (rat_0,d) else - let val d' = map adj d - val a = fold (lcm_rat o denominator_rat o fst) d' rat_1 // - fold (gcd_rat o numerator_rat o fst) d' rat_0 - in ((rat_1 // a),map (fn (c,l) => (a */ c,l)) d') - end -end; - -(* Enumeration of monomials with given multidegree bound. *) - -fun enumerate_monomials d vars = - if d < 0 then [] - else if d = 0 then [FuncUtil.Ctermfunc.empty] - else if null vars then [monomial_1] else - let val alts = - map_range (fn k => let val oths = enumerate_monomials (d - k) (tl vars) - in map (fn ks => if k = 0 then ks else FuncUtil.Ctermfunc.update (hd vars, k) ks) oths end) (d + 1) - in flat alts - end; - -(* Enumerate products of distinct input polys with degree <= d. *) -(* We ignore any constant input polynomials. *) -(* Give the output polynomial and a record of how it was derived. *) - -fun enumerate_products d pols = -if d = 0 then [(poly_const rat_1,RealArith.Rational_lt rat_1)] -else if d < 0 then [] else -case pols of - [] => [(poly_const rat_1,RealArith.Rational_lt rat_1)] - | (p,b)::ps => - let val e = multidegree p - in if e = 0 then enumerate_products d ps else - enumerate_products d ps @ - map (fn (q,c) => (poly_mul p q,RealArith.Product(b,c))) - (enumerate_products (d - e) ps) - end - -(* Convert regular polynomial. Note that we treat (0,0,0) as -1. *) - -fun epoly_of_poly p = - FuncUtil.Monomialfunc.fold (fn (m,c) => fn a => FuncUtil.Monomialfunc.update (m, Inttriplefunc.onefunc ((0,0,0), Rat.neg c)) a) p FuncUtil.Monomialfunc.empty; - -(* String for block diagonal matrix numbered k. *) - -fun sdpa_of_blockdiagonal k m = - let - val pfx = string_of_int k ^" " - val ents = - Inttriplefunc.fold - (fn ((b,i,j),c) => fn a => if i > j then a else ((b,i,j),c)::a) - m [] - val entss = sort (triple_int_ord o pairself fst) ents - in fold_rev (fn ((b,i,j),c) => fn a => - pfx ^ string_of_int b ^ " " ^ string_of_int i ^ " " ^ string_of_int j ^ - " " ^ decimalize 20 c ^ "\n" ^ a) entss "" - end; - -(* SDPA for problem using block diagonal (i.e. multiple SDPs) *) - -fun sdpa_of_blockproblem nblocks blocksizes obj mats = - let val m = length mats - 1 - in - string_of_int m ^ "\n" ^ - string_of_int nblocks ^ "\n" ^ - (space_implode " " (map string_of_int blocksizes)) ^ - "\n" ^ - sdpa_of_vector obj ^ - fold_rev2 (fn k => fn m => fn a => sdpa_of_blockdiagonal (k - 1) m ^ a) - (1 upto length mats) mats "" - end; - -(* Run prover on a problem in block diagonal form. *) - -fun run_blockproblem prover nblocks blocksizes obj mats= - parse_csdpoutput (prover (sdpa_of_blockproblem nblocks blocksizes obj mats)) - -(* 3D versions of matrix operations to consider blocks separately. *) - -val bmatrix_add = Inttriplefunc.combine (curry op +/) (fn x => x =/ rat_0); -fun bmatrix_cmul c bm = - if c =/ rat_0 then Inttriplefunc.empty - else Inttriplefunc.map (fn _ => fn x => c */ x) bm; - -val bmatrix_neg = bmatrix_cmul (Rat.rat_of_int ~1); -fun bmatrix_sub m1 m2 = bmatrix_add m1 (bmatrix_neg m2);; - -(* Smash a block matrix into components. *) - -fun blocks blocksizes bm = - map (fn (bs,b0) => - let val m = Inttriplefunc.fold - (fn ((b,i,j),c) => fn a => if b = b0 then FuncUtil.Intpairfunc.update ((i,j),c) a else a) bm FuncUtil.Intpairfunc.empty - val d = FuncUtil.Intpairfunc.fold (fn ((i,j),c) => fn a => max a (max i j)) m 0 - in (((bs,bs),m):matrix) end) - (blocksizes ~~ (1 upto length blocksizes));; - -(* FIXME : Get rid of this !!!*) -local - fun tryfind_with msg f [] = raise Failure msg - | tryfind_with msg f (x::xs) = (f x handle Failure s => tryfind_with s f xs); -in - fun tryfind f = tryfind_with "tryfind" f -end - -(* Positiv- and Nullstellensatz. Flag "linf" forces a linear representation. *) - - -fun real_positivnullstellensatz_general ctxt prover linf d eqs leqs pol = -let - val vars = fold_rev (union (op aconvc) o poly_variables) - (pol :: eqs @ map fst leqs) [] - val monoid = if linf then - (poly_const rat_1,RealArith.Rational_lt rat_1):: - (filter (fn (p,c) => multidegree p <= d) leqs) - else enumerate_products d leqs - val nblocks = length monoid - fun mk_idmultiplier k p = - let - val e = d - multidegree p - val mons = enumerate_monomials e vars - val nons = mons ~~ (1 upto length mons) - in (mons, - fold_rev (fn (m,n) => FuncUtil.Monomialfunc.update(m,Inttriplefunc.onefunc((~k,~n,n),rat_1))) nons FuncUtil.Monomialfunc.empty) - end - - fun mk_sqmultiplier k (p,c) = - let - val e = (d - multidegree p) div 2 - val mons = enumerate_monomials e vars - val nons = mons ~~ (1 upto length mons) - in (mons, - fold_rev (fn (m1,n1) => - fold_rev (fn (m2,n2) => fn a => - let val m = monomial_mul m1 m2 - in if n1 > n2 then a else - let val c = if n1 = n2 then rat_1 else rat_2 - val e = FuncUtil.Monomialfunc.tryapplyd a m Inttriplefunc.empty - in FuncUtil.Monomialfunc.update(m, tri_equation_add (Inttriplefunc.onefunc((k,n1,n2), c)) e) a - end - end) nons) - nons FuncUtil.Monomialfunc.empty) - end - - val (sqmonlist,sqs) = split_list (map2 mk_sqmultiplier (1 upto length monoid) monoid) - val (idmonlist,ids) = split_list(map2 mk_idmultiplier (1 upto length eqs) eqs) - val blocksizes = map length sqmonlist - val bigsum = - fold_rev2 (fn p => fn q => fn a => tri_epoly_pmul p q a) eqs ids - (fold_rev2 (fn (p,c) => fn s => fn a => tri_epoly_pmul p s a) monoid sqs - (epoly_of_poly(poly_neg pol))) - val eqns = FuncUtil.Monomialfunc.fold (fn (m,e) => fn a => e::a) bigsum [] - val (pvs,assig) = tri_eliminate_all_equations (0,0,0) eqns - val qvars = (0,0,0)::pvs - val allassig = fold_rev (fn v => Inttriplefunc.update(v,(Inttriplefunc.onefunc(v,rat_1)))) pvs assig - fun mk_matrix v = - Inttriplefunc.fold (fn ((b,i,j), ass) => fn m => - if b < 0 then m else - let val c = Inttriplefunc.tryapplyd ass v rat_0 - in if c = rat_0 then m else - Inttriplefunc.update ((b,j,i), c) (Inttriplefunc.update ((b,i,j), c) m) - end) - allassig Inttriplefunc.empty - val diagents = Inttriplefunc.fold - (fn ((b,i,j), e) => fn a => if b > 0 andalso i = j then tri_equation_add e a else a) - allassig Inttriplefunc.empty - - val mats = map mk_matrix qvars - val obj = (length pvs, - itern 1 pvs (fn v => fn i => FuncUtil.Intfunc.updatep iszero (i,Inttriplefunc.tryapplyd diagents v rat_0)) - FuncUtil.Intfunc.empty) - val raw_vec = if null pvs then vector_0 0 - else tri_scale_then (run_blockproblem prover nblocks blocksizes) obj mats - fun int_element (d,v) i = FuncUtil.Intfunc.tryapplyd v i rat_0 - - fun find_rounding d = - let - val _ = - if Config.get ctxt trace - then writeln ("Trying rounding with limit "^Rat.string_of_rat d ^ "\n") - else () - val vec = nice_vector d raw_vec - val blockmat = iter (1,dim vec) - (fn i => fn a => bmatrix_add (bmatrix_cmul (int_element vec i) (nth mats i)) a) - (bmatrix_neg (nth mats 0)) - val allmats = blocks blocksizes blockmat - in (vec,map diag allmats) - end - val (vec,ratdias) = - if null pvs then find_rounding rat_1 - else tryfind find_rounding (map Rat.rat_of_int (1 upto 31) @ - map pow2 (5 upto 66)) - val newassigs = - fold_rev (fn k => Inttriplefunc.update (nth pvs (k - 1), int_element vec k)) - (1 upto dim vec) (Inttriplefunc.onefunc ((0,0,0), Rat.rat_of_int ~1)) - val finalassigs = - Inttriplefunc.fold (fn (v,e) => fn a => Inttriplefunc.update(v, tri_equation_eval newassigs e) a) allassig newassigs - fun poly_of_epoly p = - FuncUtil.Monomialfunc.fold (fn (v,e) => fn a => FuncUtil.Monomialfunc.updatep iszero (v,tri_equation_eval finalassigs e) a) - p FuncUtil.Monomialfunc.empty - fun mk_sos mons = - let fun mk_sq (c,m) = - (c,fold_rev (fn k=> fn a => FuncUtil.Monomialfunc.updatep iszero (nth mons (k - 1), int_element m k) a) - (1 upto length mons) FuncUtil.Monomialfunc.empty) - in map mk_sq - end - val sqs = map2 mk_sos sqmonlist ratdias - val cfs = map poly_of_epoly ids - val msq = filter (fn (a,b) => not (null b)) (map2 pair monoid sqs) - fun eval_sq sqs = fold_rev (fn (c,q) => poly_add (poly_cmul c (poly_mul q q))) sqs poly_0 - val sanity = - fold_rev (fn ((p,c),s) => poly_add (poly_mul p (eval_sq s))) msq - (fold_rev2 (fn p => fn q => poly_add (poly_mul p q)) cfs eqs - (poly_neg pol)) - -in if not(FuncUtil.Monomialfunc.is_empty sanity) then raise Sanity else - (cfs,map (fn (a,b) => (snd a,b)) msq) - end - - -(* Iterative deepening. *) - -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)))); - - -(* Map back polynomials and their composites to a positivstellensatz. *) - -fun cterm_of_sqterm (c,p) = RealArith.Product(RealArith.Rational_lt c,RealArith.Square p); - -fun cterm_of_sos (pr,sqs) = if null sqs then pr - else RealArith.Product(pr,foldr1 RealArith.Sum (map cterm_of_sqterm sqs)); - -(* Interface to HOL. *) -local - open Conv - val concl = Thm.dest_arg o cprop_of - fun simple_cterm_ord t u = Term_Ord.fast_term_ord (term_of t, term_of u) = LESS -in - (* FIXME: Replace tryfind by get_first !! *) -fun real_nonlinear_prover proof_method ctxt = - let - val {add,mul,neg,pow,sub,main} = Semiring_Normalizer.semiring_normalizers_ord_wrapper ctxt - (the (Semiring_Normalizer.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 cert_choice translator (eqs,les,lts) = - let - val eq0 = map (poly_of_term o Thm.dest_arg1 o concl) eqs - val le0 = map (poly_of_term o Thm.dest_arg o concl) les - val lt0 = map (poly_of_term o Thm.dest_arg o concl) lts - val eqp0 = map_index (fn (i, t) => (t,RealArith.Axiom_eq i)) eq0 - val lep0 = map_index (fn (i, t) => (t,RealArith.Axiom_le i)) le0 - val ltp0 = map_index (fn (i, t) => (t,RealArith.Axiom_lt i)) lt0 - val (keq,eq) = List.partition (fn (p,_) => multidegree p = 0) eqp0 - val (klep,lep) = List.partition (fn (p,_) => multidegree p = 0) lep0 - val (kltp,ltp) = List.partition (fn (p,_) => multidegree p = 0) ltp0 - fun trivial_axiom (p,ax) = - case ax of - RealArith.Axiom_eq n => if eval FuncUtil.Ctermfunc.empty p <>/ Rat.zero then nth eqs n - else raise Failure "trivial_axiom: Not a trivial axiom" - | RealArith.Axiom_le n => if eval FuncUtil.Ctermfunc.empty p if eval FuncUtil.Ctermfunc.empty p <=/ Rat.zero then nth lts n - 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 Numeral_Simprocs.field_comp_conv) th, RealArith.Trivial) - end) - handle Failure _ => - (let val proof = - (case proof_method of Certificate certs => - (* choose certificate *) - let - fun chose_cert [] (RealArith.Cert c) = c - | chose_cert (RealArith.Left::s) (RealArith.Branch (l, _)) = chose_cert s l - | chose_cert (RealArith.Right::s) (RealArith.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 ctxt 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) => RealArith.Eqmul(q,ax)) cert_ideal eq - val proofs_cone = map cterm_of_sos cert_cone - val proof_ne = if null ltp then RealArith.Rational_lt Rat.one else - let val p = foldr1 RealArith.Product (map snd ltp) - in funpow i (fn q => RealArith.Product(p,q)) (RealArith.Rational_lt Rat.one) - end - in - foldr1 RealArith.Sum (proof_ne :: proofs_ideal @ proofs_cone) - end) - in - (translator (eqs,les,lts) proof, RealArith.Cert proof) - end) - end - in mainf end -end - -fun C f x y = f y x; - (* FIXME : This is very bad!!!*) -fun subst_conv eqs t = - let - val t' = fold (Thm.cabs o Thm.lhs_of) eqs t - in Conv.fconv_rule (Thm.beta_conversion true) (fold (C Thm.combination) eqs (Thm.reflexive t')) - end - -(* A wrapper that tries to substitute away variables first. *) - -local - open Conv - fun simple_cterm_ord t u = Term_Ord.fast_term_ord (term_of t, term_of u) = LESS - val concl = Thm.dest_arg o cprop_of - val shuffle1 = - fconv_rule (rewr_conv @{lemma "(a + x == y) == (x == y - (a::real))" by (atomize (full)) (simp add: field_simps) }) - val shuffle2 = - fconv_rule (rewr_conv @{lemma "(x + a == y) == (x == y - (a::real))" by (atomize (full)) (simp add: field_simps)}) - fun substitutable_monomial fvs tm = case term_of tm of - Free(_,@{typ real}) => if not (member (op aconvc) fvs tm) then (Rat.one,tm) - else raise Failure "substitutable_monomial" - | @{term "op * :: real => _"}$c$(t as Free _ ) => - if RealArith.is_ratconst (Thm.dest_arg1 tm) andalso not (member (op aconvc) fvs (Thm.dest_arg tm)) - then (RealArith.dest_ratconst (Thm.dest_arg1 tm),Thm.dest_arg tm) else raise Failure "substitutable_monomial" - | @{term "op + :: real => _"}$s$t => - (substitutable_monomial (Thm.add_cterm_frees (Thm.dest_arg tm) fvs) (Thm.dest_arg1 tm) - handle Failure _ => substitutable_monomial (Thm.add_cterm_frees (Thm.dest_arg1 tm) fvs) (Thm.dest_arg tm)) - | _ => raise Failure "substitutable_monomial" - - fun isolate_variable v th = - let val w = Thm.dest_arg1 (cprop_of th) - in if v aconvc w then th - else case term_of w of - @{term "op + :: real => _"}$s$t => - if Thm.dest_arg1 w aconvc v then shuffle2 th - else isolate_variable v (shuffle1 th) - | _ => error "isolate variable : This should not happen?" - end -in - -fun real_nonlinear_subst_prover prover ctxt = - let - val {add,mul,neg,pow,sub,main} = Semiring_Normalizer.semiring_normalizers_ord_wrapper ctxt - (the (Semiring_Normalizer.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 make_substitution th = - let - val (c,v) = substitutable_monomial [] (Thm.dest_arg1(concl th)) - val th1 = Drule.arg_cong_rule (Thm.capply @{cterm "op * :: real => _"} (RealArith.cterm_of_rat (Rat.inv c))) (mk_meta_eq th) - val th2 = fconv_rule (binop_conv real_poly_mul_conv) th1 - in fconv_rule (arg_conv real_poly_conv) (isolate_variable v th2) - end - fun oprconv cv ct = - let val g = Thm.dest_fun2 ct - in if g aconvc @{cterm "op <= :: real => _"} - orelse g aconvc @{cterm "op < :: real => _"} - then arg_conv cv ct else arg1_conv cv ct - end - fun mainf cert_choice translator = - let - fun substfirst(eqs,les,lts) = - ((let - val eth = tryfind make_substitution eqs - val modify = fconv_rule (arg_conv (oprconv(subst_conv [eth] then_conv real_poly_conv))) - in substfirst - (filter_out (fn t => (Thm.dest_arg1 o Thm.dest_arg o cprop_of) t - aconvc @{cterm "0::real"}) (map modify eqs), - map modify les,map modify lts) - end) - handle Failure _ => real_nonlinear_prover prover ctxt cert_choice translator (rev eqs, rev les, rev lts)) - in substfirst - end - - - in mainf - end - -(* Overall function. *) - -fun real_sos prover ctxt = - RealArith.gen_prover_real_arith ctxt (real_nonlinear_subst_prover prover ctxt) -end; - -val known_sos_constants = - [@{term "op ==>"}, @{term "Trueprop"}, - @{term HOL.implies}, @{term HOL.conj}, @{term HOL.disj}, - @{term "Not"}, @{term "op = :: bool => _"}, - @{term "All :: (real => _) => _"}, @{term "Ex :: (real => _) => _"}, - @{term "op = :: real => _"}, @{term "op < :: real => _"}, - @{term "op <= :: real => _"}, - @{term "op + :: real => _"}, @{term "op - :: real => _"}, - @{term "op * :: real => _"}, @{term "uminus :: real => _"}, - @{term "op / :: real => _"}, @{term "inverse :: real => _"}, - @{term "op ^ :: real => _"}, @{term "abs :: real => _"}, - @{term "min :: real => _"}, @{term "max :: real => _"}, - @{term "0::real"}, @{term "1::real"}, @{term "number_of :: int => real"}, - @{term "number_of :: int => nat"}, - @{term "Int.Bit0"}, @{term "Int.Bit1"}, - @{term "Int.Pls"}, @{term "Int.Min"}]; - -fun check_sos kcts ct = - let - val t = term_of ct - val _ = if not (null (Term.add_tfrees t []) - andalso null (Term.add_tvars t [])) - then error "SOS: not sos. Additional type varables" else () - val fs = Term.add_frees t [] - 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"})) 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 print_cert prover = SUBPROOF (fn {concl, context, ...} => - let - val _ = check_sos known_sos_constants concl - val (ths, certificates) = real_sos prover context (Thm.dest_arg concl) - val _ = print_cert certificates - in rtac ths 1 end) - -fun default_SOME f NONE v = SOME v - | default_SOME f (SOME v) _ = SOME v; - -fun lift_SOME f NONE a = f a - | lift_SOME f (SOME a) _ = SOME a; - - -local - val is_numeral = can (HOLogic.dest_number o term_of) -in -fun get_denom b ct = case term_of ct of - @{term "op / :: real => _"} $ _ $ _ => - if is_numeral (Thm.dest_arg ct) then get_denom b (Thm.dest_arg1 ct) - else default_SOME (get_denom b) (get_denom b (Thm.dest_arg ct)) (Thm.dest_arg ct, b) - | @{term "op < :: real => _"} $ _ $ _ => lift_SOME (get_denom true) (get_denom true (Thm.dest_arg ct)) (Thm.dest_arg1 ct) - | @{term "op <= :: real => _"} $ _ $ _ => lift_SOME (get_denom true) (get_denom true (Thm.dest_arg ct)) (Thm.dest_arg1 ct) - | _ $ _ => lift_SOME (get_denom b) (get_denom b (Thm.dest_fun ct)) (Thm.dest_arg ct) - | _ => NONE -end; - -fun elim_one_denom_tac ctxt = -CSUBGOAL (fn (P,i) => - case get_denom false P of - NONE => no_tac - | SOME (d,ord) => - let - val ss = simpset_of ctxt addsimps @{thms field_simps} - addsimps [@{thm nonzero_power_divide}, @{thm power_divide}] - val th = instantiate' [] [SOME d, SOME (Thm.dest_arg P)] - (if ord then @{lemma "(d=0 --> P) & (d>0 --> P) & (d<(0::real) --> P) ==> P" by auto} - else @{lemma "(d=0 --> P) & (d ~= (0::real) --> P) ==> P" by blast}) - in rtac th i THEN Simplifier.asm_full_simp_tac ss i end); - -fun elim_denom_tac ctxt i = REPEAT (elim_one_denom_tac ctxt i); - -fun sos_tac print_cert prover ctxt = - Object_Logic.full_atomize_tac THEN' - elim_denom_tac ctxt THEN' - core_sos_tac print_cert prover ctxt; - -end; diff -r 655f583840d0 -r 9908cf4af394 src/HOL/Library/Sum_of_Squares.thy --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/HOL/Library/Sum_of_Squares.thy Sat Jan 08 11:29:25 2011 -0800 @@ -0,0 +1,159 @@ +(* Title: HOL/Library/Sum_of_Squares.thy + Author: Amine Chaieb, University of Cambridge + Author: Philipp Meyer, TU Muenchen +*) + +header {* A decision method for universal multivariate real arithmetic with addition, + multiplication and ordering using semidefinite programming *} + +theory Sum_of_Squares +imports Complex_Main +uses + "positivstellensatz.ML" + "Sum_of_Squares/sum_of_squares.ML" + "Sum_of_Squares/positivstellensatz_tools.ML" + "Sum_of_Squares/sos_wrapper.ML" +begin + +text {* + In order to use the method sos, call it with @{text "(sos + remote_csdp)"} to use the remote solver. Or install CSDP + (https://projects.coin-or.org/Csdp), configure the Isabelle setting + @{text CSDP_EXE}, and call it with @{text "(sos csdp)"}. By + default, sos calls @{text remote_csdp}. This can take of the order + 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. +*} + +setup Sum_of_Squares.setup +setup SOS_Wrapper.setup + +text {* Tests *} + +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 "a1 >= 0 & a2 >= 0 \ (a1 * a1 + a2 * a2 = b1 * b1 + b2 * b2 + 2) \ (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 --> 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 "(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 "(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 "(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 "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 "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)))))))))))))") + +(* ------------------------------------------------------------------------- *) +(* Over a larger but simpler interval. *) +(* ------------------------------------------------------------------------- *) + +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)))))))))))") + +(* ------------------------------------------------------------------------- *) +(* 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))))))))))))))))))") + +(* ------------------------------------------------------------------------- *) +(* 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))))))") + +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 "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 \ 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))))))") + +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 "(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))))") + +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 "(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))))))") + + + +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 \ 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 --> 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 "((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_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_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_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_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 655f583840d0 -r 9908cf4af394 src/HOL/Library/Sum_of_Squares/etc/settings --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/HOL/Library/Sum_of_Squares/etc/settings Sat Jan 08 11:29:25 2011 -0800 @@ -0,0 +1,5 @@ +# -*- shell-script -*- :mode=shellscript: + +ISABELLE_SUM_OF_SQUARES="$COMPONENT" +NEOS_SERVER="http://neos-server.org:3332" + diff -r 655f583840d0 -r 9908cf4af394 src/HOL/Library/Sum_of_Squares/neos_csdp_client --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/HOL/Library/Sum_of_Squares/neos_csdp_client Sat Jan 08 11:29:25 2011 -0800 @@ -0,0 +1,85 @@ +#!/usr/bin/env python +import sys +import signal +import xmlrpclib +import time +import re +import os + +# Neos server config +neos = xmlrpclib.Server(os.getenv("NEOS_SERVER")) + +jobNumber = 0 +password = "" +inputfile = None +outputfile = None +# interrupt handler +def cleanup(signal, frame): + sys.stdout.write("Caught interrupt, cleaning up\n") + if jobNumber != 0: + neos.killJob(jobNumber, password) + if inputfile != None: + inputfile.close() + if outputfile != None: + outputfile.close() + sys.exit(21) + +signal.signal(signal.SIGHUP, cleanup) +signal.signal(signal.SIGINT, cleanup) +signal.signal(signal.SIGQUIT, cleanup) +signal.signal(signal.SIGTERM, cleanup) + +if len(sys.argv) <> 3: + sys.stderr.write("Usage: neos_csdp_client \n") + sys.exit(19) + +xml_pre = "\nsdp\ncsdp\nSPARSE_SDPA\n\n\n" +xml = xml_pre +inputfile = open(sys.argv[1],"r") +buffer = 1 +while buffer: + buffer = inputfile.read() + xml += buffer +inputfile.close() +xml += xml_post + +(jobNumber,password) = neos.submitJob(xml) + +if jobNumber == 0: + sys.stdout.write("error submitting job: %s" % password) + sys.exit(20) +else: + sys.stdout.write("jobNumber = %d\tpassword = %s\n" % (jobNumber,password)) + +offset=0 +messages = "" +status="Waiting" +while status == "Running" or status=="Waiting": + time.sleep(1) + (msg,offset) = neos.getIntermediateResults(jobNumber,password,offset) + messages += msg.data + sys.stdout.write(msg.data) + status = neos.getJobStatus(jobNumber, password) + +msg = neos.getFinalResults(jobNumber, password).data +sys.stdout.write("---------- Begin CSDP Output -------------\n"); +sys.stdout.write(msg) + +# extract solution +result = msg.split("Solution:") +if len(result) > 1: + solution = result[1].strip() + if solution != "": + outputfile = open(sys.argv[2],"w") + outputfile.write(solution) + outputfile.close() + +# extract return code +p = re.compile(r"^Error: Command exited with non-zero status (\d+)$", re.MULTILINE) +m = p.search(messages) +if m: + sys.exit(int(m.group(1))) +else: + sys.exit(0) + diff -r 655f583840d0 -r 9908cf4af394 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 Sat Jan 08 11:29:25 2011 -0800 @@ -0,0 +1,156 @@ +(* Title: HOL/Library/Sum_of_Squares/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 + +(*** 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 FuncUtil.Ctermfunc.is_empty m then "1" + else + let + val m' = FuncUtil.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 FuncUtil.Ctermfunc.is_empty 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 FuncUtil.Monomialfunc.is_empty p then "0" + else + let + val cms = map string_of_cmonomial + (sort (prod_ord FuncUtil.monomial_order (K EQUAL)) (FuncUtil.Monomialfunc.dest p)) + in foldr1 (fn (t1, t2) => t1 ^ " + " ^ t2) cms + end; + +fun pss_to_cert (RealArith.Axiom_eq i) = "A=" ^ string_of_int i + | pss_to_cert (RealArith.Axiom_le i) = "A<=" ^ string_of_int i + | pss_to_cert (RealArith.Axiom_lt i) = "A<" ^ string_of_int i + | pss_to_cert (RealArith.Rational_eq r) = "R=" ^ string_of_rat r + | pss_to_cert (RealArith.Rational_le r) = "R<=" ^ string_of_rat r + | pss_to_cert (RealArith.Rational_lt r) = "R<" ^ string_of_rat r + | pss_to_cert (RealArith.Square p) = "[" ^ string_of_poly p ^ "]^2" + | pss_to_cert (RealArith.Eqmul (p, pss)) = "([" ^ string_of_poly p ^ "] * " ^ pss_to_cert pss ^ ")" + | pss_to_cert (RealArith.Sum (pss1, pss2)) = "(" ^ pss_to_cert pss1 ^ " + " ^ pss_to_cert pss2 ^ ")" + | pss_to_cert (RealArith.Product (pss1, pss2)) = "(" ^ pss_to_cert pss1 ^ " * " ^ pss_to_cert pss2 ^ ")" + +fun pss_tree_to_cert RealArith.Trivial = "()" + | pss_tree_to_cert (RealArith.Cert pss) = "(" ^ pss_to_cert pss ^ ")" + | pss_tree_to_cert (RealArith.Branch (t1, t2)) = "(" ^ pss_tree_to_cert t1 ^ " & " ^ pss_tree_to_cert t2 ^ ")" + +(*** certificate parsing ***) + +(* basic parser *) + +val str = Scan.this_string + +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 (str "~" >> K ~1) 1 -- nat >> op *; +val rat = int --| str "/" -- int >> Rat.rat_of_quotient +val rat_int = rat || int >> Rat.rat_of_int + +(* polynomial parser *) + +fun repeat_sep s f = f ::: Scan.repeat (str s |-- f) + +val parse_id = Scan.one Symbol.is_letter ::: Scan.many Symbol.is_letdig >> implode + +fun parse_varpow ctxt = parse_id -- Scan.optional (str "^" |-- nat) 1 >> + (fn (x, k) => (cterm_of (ProofContext.theory_of ctxt) (Free (x, @{typ real})), k)) + +fun parse_monomial ctxt = repeat_sep "*" (parse_varpow ctxt) >> + (fn xs => fold FuncUtil.Ctermfunc.update xs FuncUtil.Ctermfunc.empty) + +fun parse_cmonomial ctxt = + rat_int --| str "*" -- (parse_monomial ctxt) >> swap || + (parse_monomial ctxt) >> (fn m => (m, Rat.one)) || + rat_int >> (fn r => (FuncUtil.Ctermfunc.empty, r)) + +fun parse_poly ctxt = repeat_sep "+" (parse_cmonomial ctxt) >> + (fn xs => fold FuncUtil.Monomialfunc.update xs FuncUtil.Monomialfunc.empty) + +(* positivstellensatz parser *) + +val parse_axiom = + (str "A=" |-- int >> RealArith.Axiom_eq) || + (str "A<=" |-- int >> RealArith.Axiom_le) || + (str "A<" |-- int >> RealArith.Axiom_lt) + +val parse_rational = + (str "R=" |-- rat_int >> RealArith.Rational_eq) || + (str "R<=" |-- rat_int >> RealArith.Rational_le) || + (str "R<" |-- rat_int >> RealArith.Rational_lt) + +fun parse_cert ctxt input = + let + val pc = parse_cert ctxt + val pp = parse_poly ctxt + in + (parse_axiom || + parse_rational || + str "[" |-- pp --| str "]^2" >> RealArith.Square || + str "([" |-- pp --| str "]*" -- pc --| str ")" >> RealArith.Eqmul || + str "(" |-- pc --| str "*" -- pc --| str ")" >> RealArith.Product || + str "(" |-- pc --| str "+" -- pc --| str ")" >> RealArith.Sum) input + end + +fun parse_cert_tree ctxt input = + let + val pc = parse_cert ctxt + val pt = parse_cert_tree ctxt + in + (str "()" >> K RealArith.Trivial || + str "(" |-- pc --| str ")" >> RealArith.Cert || + str "(" |-- pt --| str "&" -- pt --| str ")" >> RealArith.Branch) input + end + +(* scanner *) + +fun cert_to_pss_tree ctxt input_str = Symbol.scanner "bad certificate" (parse_cert_tree ctxt) + (filter_out Symbol.is_blank (Symbol.explode input_str)) + +end + + diff -r 655f583840d0 -r 9908cf4af394 src/HOL/Library/Sum_of_Squares/sos_wrapper.ML --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/HOL/Library/Sum_of_Squares/sos_wrapper.ML Sat Jan 08 11:29:25 2011 -0800 @@ -0,0 +1,159 @@ +(* Title: HOL/Library/Sum_of_Squares/sos_wrapper.ML + Author: Philipp Meyer, TU Muenchen + +Added functionality for sums of squares, e.g. calling a remote prover. +*) + +signature SOS_WRAPPER = +sig + datatype prover_result = Success | Failure | Error + val setup: theory -> theory + val dest_dir: string Config.T + val prover_name: string Config.T +end + +structure SOS_Wrapper: SOS_WRAPPER = +struct + +datatype prover_result = Success | Failure | Error + +fun str_of_result Success = "Success" + | str_of_result Failure = "Failure" + | str_of_result Error = "Error" + + +(*** calling provers ***) + +val (dest_dir, setup_dest_dir) = Attrib.config_string "sos_dest_dir" (K "") + +fun filename dir name = + let + val probfile = Path.basic (name ^ serial_string ()) + val dir_path = Path.explode dir + in + if dir = "" then + File.tmp_path probfile + else if File.exists dir_path then + Path.append dir_path probfile + else error ("No such directory: " ^ dir) + end + +fun run_solver ctxt name cmd find_failure input = + let + val _ = warning ("Calling solver: " ^ name) + + (* create input file *) + val dir = Config.get ctxt dest_dir + val input_file = filename dir "sos_in" + val _ = File.write input_file input + + (* call solver *) + val output_file = filename dir "sos_out" + val (output, rv) = + bash_output + (if File.exists cmd then + space_implode " " + [File.shell_path cmd, File.shell_path input_file, File.shell_path output_file] + else error ("Bad executable: " ^ File.platform_path cmd)) + + (* read and analyze output *) + val (res, res_msg) = find_failure rv + val result = if File.exists output_file then File.read output_file else "" + + (* remove temporary files *) + val _ = + if dir = "" then + (File.rm input_file; if File.exists output_file then File.rm output_file else ()) + else () + + val _ = + if Config.get ctxt Sum_of_Squares.trace + then writeln ("Solver output:\n" ^ output) + else () + + val _ = warning (str_of_result res ^ ": " ^ res_msg) + in + (case res of + Success => result + | Failure => raise Sum_of_Squares.Failure res_msg + | Error => error ("Prover failed: " ^ res_msg)) + end + + +(*** various provers ***) + +(* local csdp client *) + +fun find_csdp_failure rv = + case rv of + 0 => (Success, "SDP solved") + | 1 => (Failure, "SDP is primal infeasible") + | 2 => (Failure, "SDP is dual infeasible") + | 3 => (Success, "SDP solved with reduced accuracy") + | 4 => (Failure, "Maximum iterations reached") + | 5 => (Failure, "Stuck at edge of primal feasibility") + | 6 => (Failure, "Stuck at edge of dual infeasibility") + | 7 => (Failure, "Lack of progress") + | 8 => (Failure, "X, Z, or O was singular") + | 9 => (Failure, "Detected NaN or Inf values") + | _ => (Error, "return code is " ^ string_of_int rv) + +val csdp = ("$CSDP_EXE", find_csdp_failure) + + +(* remote neos server *) + +fun find_neos_failure rv = + case rv of + 20 => (Error, "error submitting job") + | 21 => (Error, "interrupt") + | _ => find_csdp_failure rv + +val neos_csdp = ("$ISABELLE_SUM_OF_SQUARES/neos_csdp_client", find_neos_failure) + + +(* named provers *) + +fun prover "remote_csdp" = neos_csdp + | prover "csdp" = csdp + | prover name = error ("Unknown prover: " ^ name) + +val (prover_name, setup_prover_name) = Attrib.config_string "sos_prover_name" (K "remote_csdp") + +fun call_solver ctxt opt_name = + let + val name = the_default (Config.get ctxt prover_name) opt_name + val (cmd, find_failure) = prover name + in run_solver ctxt 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 ^ "\")") + +val print_cert = warning o output_line o PositivstellensatzTools.pss_tree_to_cert + + +(* method setup *) + +fun sos_solver print method = SIMPLE_METHOD' o Sum_of_Squares.sos_tac print method + +val setup = + setup_dest_dir #> + setup_prover_name #> + Method.setup @{binding sos} + (Scan.lift (Scan.option Parse.xname) + >> (fn opt_name => fn ctxt => + sos_solver print_cert + (Sum_of_Squares.Prover (call_solver ctxt opt_name)) ctxt)) + "prove universal problems over the reals using sums of squares" #> + Method.setup @{binding sos_cert} + (Scan.lift Parse.string + >> (fn cert => fn ctxt => + sos_solver ignore + (Sum_of_Squares.Certificate (PositivstellensatzTools.cert_to_pss_tree ctxt cert)) ctxt)) + "prove universal problems over the reals using sums of squares with certificates" + +end diff -r 655f583840d0 -r 9908cf4af394 src/HOL/Library/Sum_of_Squares/sum_of_squares.ML --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/HOL/Library/Sum_of_Squares/sum_of_squares.ML Sat Jan 08 11:29:25 2011 -0800 @@ -0,0 +1,1435 @@ +(* Title: HOL/Library/Sum_of_Squares/sum_of_squares.ML + Author: Amine Chaieb, University of Cambridge + Author: Philipp Meyer, TU Muenchen + +A tactic for proving nonlinear inequalities. +*) + +signature SUM_OF_SQUARES = +sig + datatype proof_method = Certificate of RealArith.pss_tree | Prover of string -> string + val sos_tac: (RealArith.pss_tree -> unit) -> proof_method -> Proof.context -> int -> tactic + val trace: bool Config.T + val setup: theory -> theory + exception Failure of string; +end + +structure Sum_of_Squares: SUM_OF_SQUARES = +struct + +val rat_0 = Rat.zero; +val rat_1 = Rat.one; +val rat_2 = Rat.two; +val rat_10 = Rat.rat_of_int 10; +val rat_1_2 = rat_1 // rat_2; +val max = Integer.max; +val min = Integer.min; + +val denominator_rat = Rat.quotient_of_rat #> snd #> Rat.rat_of_int; +val numerator_rat = Rat.quotient_of_rat #> fst #> Rat.rat_of_int; +fun int_of_rat a = + case Rat.quotient_of_rat a of (i,1) => i | _ => error "int_of_rat: not an int"; +fun lcm_rat x y = Rat.rat_of_int (Integer.lcm (int_of_rat x) (int_of_rat y)); + +fun rat_pow r i = + let fun pow r i = + if i = 0 then rat_1 else + let val d = pow r (i div 2) + in d */ d */ (if i mod 2 = 0 then rat_1 else r) + end + in if i < 0 then pow (Rat.inv r) (~ i) else pow r i end; + +fun round_rat r = + let val (a,b) = Rat.quotient_of_rat (Rat.abs r) + val d = a div b + val s = if r = b then d + 1 else d) end + +val abs_rat = Rat.abs; +val pow2 = rat_pow rat_2; +val pow10 = rat_pow rat_10; + +val (trace, setup_trace) = Attrib.config_bool "sos_trace" (K false); +val setup = setup_trace; + +exception Sanity; + +exception Unsolvable; + +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 +fun normalize y = + if abs_rat y =/ rat_1 then normalize (y // rat_10) + 1 + else 0 + in +fun decimalize d x = + if x =/ rat_0 then "0.0" else + let + val y = Rat.abs x + val e = normalize y + val z = pow10(~ e) */ y +/ rat_1 + val k = int_of_rat (round_rat(pow10 d */ z)) + in (if x a + | h::t => itern (k + 1) t f (f h k a); + +fun iter (m,n) f a = + if n < m then a + else iter (m+1,n) f (f m a); + +(* The main types. *) + +type vector = int* Rat.rat FuncUtil.Intfunc.table; + +type matrix = (int*int)*(Rat.rat FuncUtil.Intpairfunc.table); + +fun iszero (k,r) = r =/ rat_0; + + +(* Vectors. Conventionally indexed 1..n. *) + +fun vector_0 n = (n,FuncUtil.Intfunc.empty):vector; + +fun dim (v:vector) = fst v; + +fun vector_const c n = + if c =/ rat_0 then vector_0 n + else (n,fold_rev (fn k => FuncUtil.Intfunc.update (k,c)) (1 upto n) FuncUtil.Intfunc.empty) :vector; + +val vector_1 = vector_const rat_1; + +fun vector_cmul c (v:vector) = + let val n = dim v + in if c =/ rat_0 then vector_0 n + else (n,FuncUtil.Intfunc.map (fn _ => fn x => c */ x) (snd v)) + end; + +fun vector_neg (v:vector) = (fst v,FuncUtil.Intfunc.map (K Rat.neg) (snd v)) :vector; + +fun vector_add (v1:vector) (v2:vector) = + let val m = dim v1 + val n = dim v2 + in if m <> n then error "vector_add: incompatible dimensions" + else (n,FuncUtil.Intfunc.combine (curry op +/) (fn x => x =/ rat_0) (snd v1) (snd v2)) :vector + end; + +fun vector_sub v1 v2 = vector_add v1 (vector_neg v2); + +fun vector_dot (v1:vector) (v2:vector) = + let val m = dim v1 + val n = dim v2 + in if m <> n then error "vector_dot: incompatible dimensions" + else FuncUtil.Intfunc.fold (fn (i,x) => fn a => x +/ a) + (FuncUtil.Intfunc.combine (curry op */) (fn x => x =/ rat_0) (snd v1) (snd v2)) rat_0 + end; + +fun vector_of_list l = + let val n = length l + in (n,fold_rev2 (curry FuncUtil.Intfunc.update) (1 upto n) l FuncUtil.Intfunc.empty) :vector + end; + +(* Matrices; again rows and columns indexed from 1. *) + +fun matrix_0 (m,n) = ((m,n),FuncUtil.Intpairfunc.empty):matrix; + +fun dimensions (m:matrix) = fst m; + +fun matrix_const c (mn as (m,n)) = + if m <> n then error "matrix_const: needs to be square" + else if c =/ rat_0 then matrix_0 mn + else (mn,fold_rev (fn k => FuncUtil.Intpairfunc.update ((k,k), c)) (1 upto n) FuncUtil.Intpairfunc.empty) :matrix;; + +val matrix_1 = matrix_const rat_1; + +fun matrix_cmul c (m:matrix) = + let val (i,j) = dimensions m + in if c =/ rat_0 then matrix_0 (i,j) + else ((i,j),FuncUtil.Intpairfunc.map (fn _ => fn x => c */ x) (snd m)) + end; + +fun matrix_neg (m:matrix) = + (dimensions m, FuncUtil.Intpairfunc.map (K Rat.neg) (snd m)) :matrix; + +fun matrix_add (m1:matrix) (m2:matrix) = + let val d1 = dimensions m1 + val d2 = dimensions m2 + in if d1 <> d2 + then error "matrix_add: incompatible dimensions" + else (d1,FuncUtil.Intpairfunc.combine (curry op +/) (fn x => x =/ rat_0) (snd m1) (snd m2)) :matrix + end;; + +fun matrix_sub m1 m2 = matrix_add m1 (matrix_neg m2); + +fun row k (m:matrix) = + let val (i,j) = dimensions m + in (j, + FuncUtil.Intpairfunc.fold (fn ((i,j), c) => fn a => if i = k then FuncUtil.Intfunc.update (j,c) a else a) (snd m) FuncUtil.Intfunc.empty ) : vector + end; + +fun column k (m:matrix) = + let val (i,j) = dimensions m + in (i, + FuncUtil.Intpairfunc.fold (fn ((i,j), c) => fn a => if j = k then FuncUtil.Intfunc.update (i,c) a else a) (snd m) FuncUtil.Intfunc.empty) + : vector + end; + +fun transp (m:matrix) = + let val (i,j) = dimensions m + in + ((j,i),FuncUtil.Intpairfunc.fold (fn ((i,j), c) => fn a => FuncUtil.Intpairfunc.update ((j,i), c) a) (snd m) FuncUtil.Intpairfunc.empty) :matrix + end; + +fun diagonal (v:vector) = + let val n = dim v + in ((n,n),FuncUtil.Intfunc.fold (fn (i, c) => fn a => FuncUtil.Intpairfunc.update ((i,i), c) a) (snd v) FuncUtil.Intpairfunc.empty) : matrix + end; + +fun matrix_of_list l = + let val m = length l + in if m = 0 then matrix_0 (0,0) else + let val n = length (hd l) + in ((m,n),itern 1 l (fn v => fn i => itern 1 v (fn c => fn j => FuncUtil.Intpairfunc.update ((i,j), c))) FuncUtil.Intpairfunc.empty) + end + end; + +(* Monomials. *) + +fun monomial_eval assig m = + FuncUtil.Ctermfunc.fold (fn (x, k) => fn a => a */ rat_pow (FuncUtil.Ctermfunc.apply assig x) k) + m rat_1; +val monomial_1 = FuncUtil.Ctermfunc.empty; + +fun monomial_var x = FuncUtil.Ctermfunc.onefunc (x, 1); + +val monomial_mul = + FuncUtil.Ctermfunc.combine Integer.add (K false); + +fun monomial_pow m k = + if k = 0 then monomial_1 + else FuncUtil.Ctermfunc.map (fn _ => fn x => k * x) m; + +fun monomial_divides m1 m2 = + FuncUtil.Ctermfunc.fold (fn (x, k) => fn a => FuncUtil.Ctermfunc.tryapplyd m2 x 0 >= k andalso a) m1 true;; + +fun monomial_div m1 m2 = + let val m = FuncUtil.Ctermfunc.combine Integer.add + (fn x => x = 0) m1 (FuncUtil.Ctermfunc.map (fn _ => fn x => ~ x) m2) + in if FuncUtil.Ctermfunc.fold (fn (x, k) => fn a => k >= 0 andalso a) m true then m + else error "monomial_div: non-divisible" + end; + +fun monomial_degree x m = + FuncUtil.Ctermfunc.tryapplyd m x 0;; + +fun monomial_lcm m1 m2 = + fold_rev (fn x => FuncUtil.Ctermfunc.update (x, max (monomial_degree x m1) (monomial_degree x m2))) + (union (is_equal o FuncUtil.cterm_ord) (FuncUtil.Ctermfunc.dom m1) (FuncUtil.Ctermfunc.dom m2)) (FuncUtil.Ctermfunc.empty); + +fun monomial_multidegree m = + FuncUtil.Ctermfunc.fold (fn (x, k) => fn a => k + a) m 0;; + +fun monomial_variables m = FuncUtil.Ctermfunc.dom m;; + +(* Polynomials. *) + +fun eval assig p = + FuncUtil.Monomialfunc.fold (fn (m, c) => fn a => a +/ c */ monomial_eval assig m) p rat_0; + +val poly_0 = FuncUtil.Monomialfunc.empty; + +fun poly_isconst p = + FuncUtil.Monomialfunc.fold (fn (m, c) => fn a => FuncUtil.Ctermfunc.is_empty m andalso a) p true; + +fun poly_var x = FuncUtil.Monomialfunc.onefunc (monomial_var x,rat_1); + +fun poly_const c = + if c =/ rat_0 then poly_0 else FuncUtil.Monomialfunc.onefunc(monomial_1, c); + +fun poly_cmul c p = + if c =/ rat_0 then poly_0 + else FuncUtil.Monomialfunc.map (fn _ => fn x => c */ x) p; + +fun poly_neg p = FuncUtil.Monomialfunc.map (K Rat.neg) p;; + +fun poly_add p1 p2 = + FuncUtil.Monomialfunc.combine (curry op +/) (fn x => x =/ rat_0) p1 p2; + +fun poly_sub p1 p2 = poly_add p1 (poly_neg p2); + +fun poly_cmmul (c,m) p = + if c =/ rat_0 then poly_0 + else if FuncUtil.Ctermfunc.is_empty m + then FuncUtil.Monomialfunc.map (fn _ => fn d => c */ d) p + else FuncUtil.Monomialfunc.fold (fn (m', d) => fn a => (FuncUtil.Monomialfunc.update (monomial_mul m m', c */ d) a)) p poly_0; + +fun poly_mul p1 p2 = + FuncUtil.Monomialfunc.fold (fn (m, c) => fn a => poly_add (poly_cmmul (c,m) p2) a) p1 poly_0; + +fun poly_div p1 p2 = + if not(poly_isconst p2) + then error "poly_div: non-constant" else + let val c = eval FuncUtil.Ctermfunc.empty p2 + in if c =/ rat_0 then error "poly_div: division by zero" + else poly_cmul (Rat.inv c) p1 + end; + +fun poly_square p = poly_mul p p; + +fun poly_pow p k = + if k = 0 then poly_const rat_1 + else if k = 1 then p + else let val q = poly_square(poly_pow p (k div 2)) in + if k mod 2 = 1 then poly_mul p q else q end; + +fun poly_exp p1 p2 = + if not(poly_isconst p2) + then error "poly_exp: not a constant" + else poly_pow p1 (int_of_rat (eval FuncUtil.Ctermfunc.empty p2)); + +fun degree x p = + FuncUtil.Monomialfunc.fold (fn (m,c) => fn a => max (monomial_degree x m) a) p 0; + +fun multidegree p = + FuncUtil.Monomialfunc.fold (fn (m, c) => fn a => max (monomial_multidegree m) a) p 0; + +fun poly_variables p = + sort FuncUtil.cterm_ord (FuncUtil.Monomialfunc.fold_rev (fn (m, c) => union (is_equal o FuncUtil.cterm_ord) (monomial_variables m)) p []);; + +(* Order monomials for human presentation. *) + +val humanorder_varpow = prod_ord FuncUtil.cterm_ord (rev_order o int_ord); + +local + fun ord (l1,l2) = case (l1,l2) of + (_,[]) => LESS + | ([],_) => GREATER + | (h1::t1,h2::t2) => + (case humanorder_varpow (h1, h2) of + LESS => LESS + | EQUAL => ord (t1,t2) + | GREATER => GREATER) +in fun humanorder_monomial m1 m2 = + ord (sort humanorder_varpow (FuncUtil.Ctermfunc.dest m1), + sort humanorder_varpow (FuncUtil.Ctermfunc.dest m2)) +end; + +(* Conversions to strings. *) + +fun string_of_vector min_size max_size (v:vector) = + let val n_raw = dim v + in if n_raw = 0 then "[]" else + let + val n = max min_size (min n_raw max_size) + val xs = map (Rat.string_of_rat o (fn i => FuncUtil.Intfunc.tryapplyd (snd v) i rat_0)) (1 upto n) + in "[" ^ space_implode ", " xs ^ + (if n_raw > max_size then ", ...]" else "]") + end + end; + +fun string_of_matrix max_size (m:matrix) = + let + val (i_raw,j_raw) = dimensions m + val i = min max_size i_raw + val j = min max_size j_raw + val rstr = map (fn k => string_of_vector j j (row k m)) (1 upto i) + in "["^ space_implode ";\n " rstr ^ + (if j > max_size then "\n ...]" else "]") + end; + +fun string_of_term t = + case t of + a$b => "("^(string_of_term a)^" "^(string_of_term b)^")" + | Abs x => + let val (xn, b) = Term.dest_abs x + in "(\\"^xn^"."^(string_of_term b)^")" + end + | Const(s,_) => s + | Free (s,_) => s + | Var((s,_),_) => s + | _ => error "string_of_term"; + +val string_of_cterm = string_of_term o term_of; + +fun string_of_varpow x k = + if k = 1 then string_of_cterm x + else string_of_cterm x^"^"^string_of_int k; + +fun string_of_monomial m = + if FuncUtil.Ctermfunc.is_empty m then "1" else + let val vps = fold_rev (fn (x,k) => fn a => string_of_varpow x k :: a) + (sort humanorder_varpow (FuncUtil.Ctermfunc.dest m)) [] + in space_implode "*" vps + end; + +fun string_of_cmonomial (c,m) = + if FuncUtil.Ctermfunc.is_empty m then Rat.string_of_rat c + 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 = + if FuncUtil.Monomialfunc.is_empty p then "<<0>>" else + let + val cms = sort (fn ((m1,_),(m2,_)) => humanorder_monomial m1 m2) (FuncUtil.Monomialfunc.dest p) + val s = fold (fn (m,c) => fn a => + if c >" + end; + +(* Conversion from HOL term. *) + +local + val neg_tm = @{cterm "uminus :: real => _"} + val add_tm = @{cterm "op + :: real => _"} + val sub_tm = @{cterm "op - :: real => _"} + val mul_tm = @{cterm "op * :: real => _"} + val inv_tm = @{cterm "inverse :: real => _"} + val div_tm = @{cterm "op / :: real => _"} + val pow_tm = @{cterm "op ^ :: real => _"} + val zero_tm = @{cterm "0:: real"} + val is_numeral = can (HOLogic.dest_number o term_of) + fun is_comb t = case t of _$_ => true | _ => false + fun poly_of_term tm = + if tm aconvc zero_tm then poly_0 + else if RealArith.is_ratconst tm + then poly_const(RealArith.dest_ratconst tm) + else + (let val (lop,r) = Thm.dest_comb tm + in if lop aconvc neg_tm then poly_neg(poly_of_term r) + else if lop aconvc inv_tm then + let val p = poly_of_term r + in if poly_isconst p + then poly_const(Rat.inv (eval FuncUtil.Ctermfunc.empty p)) + else error "poly_of_term: inverse of non-constant polyomial" + end + else (let val (opr,l) = Thm.dest_comb lop + in + if opr aconvc pow_tm andalso is_numeral r + then poly_pow (poly_of_term l) ((snd o HOLogic.dest_number o term_of) r) + else if opr aconvc add_tm + then poly_add (poly_of_term l) (poly_of_term r) + else if opr aconvc sub_tm + then poly_sub (poly_of_term l) (poly_of_term r) + else if opr aconvc mul_tm + then poly_mul (poly_of_term l) (poly_of_term r) + else if opr aconvc div_tm + then let + val p = poly_of_term l + val q = poly_of_term r + in if poly_isconst q then poly_cmul (Rat.inv (eval FuncUtil.Ctermfunc.empty q)) p + else error "poly_of_term: division by non-constant polynomial" + end + else poly_var tm + + end + handle CTERM ("dest_comb",_) => poly_var tm) + end + handle CTERM ("dest_comb",_) => poly_var tm) +in +val poly_of_term = fn tm => + if type_of (term_of tm) = @{typ real} then poly_of_term tm + else error "poly_of_term: term does not have real type" +end; + +(* String of vector (just a list of space-separated numbers). *) + +fun sdpa_of_vector (v:vector) = + let + val n = dim v + val strs = map (decimalize 20 o (fn i => FuncUtil.Intfunc.tryapplyd (snd v) i rat_0)) (1 upto n) + in space_implode " " strs ^ "\n" + end; + +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'))); +structure Inttriplefunc = FuncFun(type key = int*int*int val ord = triple_int_ord); + +(* String for block diagonal matrix numbered k. *) + +fun sdpa_of_blockdiagonal k m = + let + val pfx = string_of_int k ^" " + val ents = + Inttriplefunc.fold (fn ((b,i,j), c) => fn a => if i > j then a else ((b,i,j),c)::a) m [] + val entss = sort (triple_int_ord o pairself fst) ents + in fold_rev (fn ((b,i,j),c) => fn a => + pfx ^ string_of_int b ^ " " ^ string_of_int i ^ " " ^ string_of_int j ^ + " " ^ decimalize 20 c ^ "\n" ^ a) entss "" + end; + +(* String for a matrix numbered k, in SDPA sparse format. *) + +fun sdpa_of_matrix k (m:matrix) = + let + val pfx = string_of_int k ^ " 1 " + val ms = FuncUtil.Intpairfunc.fold (fn ((i,j), c) => fn a => if i > j then a else ((i,j),c)::a) (snd m) [] + val mss = sort ((prod_ord int_ord int_ord) o pairself fst) ms + in fold_rev (fn ((i,j),c) => fn a => + pfx ^ string_of_int i ^ " " ^ string_of_int j ^ + " " ^ decimalize 20 c ^ "\n" ^ a) mss "" + end;; + +(* ------------------------------------------------------------------------- *) +(* String in SDPA sparse format for standard SDP problem: *) +(* *) +(* X = v_1 * [M_1] + ... + v_m * [M_m] - [M_0] must be PSD *) +(* Minimize obj_1 * v_1 + ... obj_m * v_m *) +(* ------------------------------------------------------------------------- *) + +fun sdpa_of_problem obj mats = + let + val m = length mats - 1 + val (n,_) = dimensions (hd mats) + in + string_of_int m ^ "\n" ^ + "1\n" ^ + string_of_int n ^ "\n" ^ + sdpa_of_vector obj ^ + fold_rev2 (fn k => fn m => fn a => sdpa_of_matrix (k - 1) m ^ a) (1 upto length mats) mats "" + end; + +fun index_char str chr pos = + if pos >= String.size str then ~1 + else if String.sub(str,pos) = chr then pos + else index_char str chr (pos + 1); +fun rat_of_quotient (a,b) = if b = 0 then rat_0 else Rat.rat_of_quotient (a,b); +fun rat_of_string s = + let val n = index_char s #"/" 0 in + if n = ~1 then s |> Int.fromString |> the |> Rat.rat_of_int + else + let val SOME numer = Int.fromString(String.substring(s,0,n)) + val SOME den = Int.fromString (String.substring(s,n+1,String.size s - n - 1)) + in rat_of_quotient(numer, den) + end + end; + +fun isspace x = (x = " "); +fun isnum x = member (op =) ["0","1","2","3","4","5","6","7","8","9"] x; + +(* More parser basics. *) + + val word = Scan.this_string + fun token s = + Scan.repeat ($$ " ") |-- word s --| Scan.repeat ($$ " ") + val numeral = Scan.one isnum + val decimalint = Scan.repeat1 numeral >> (rat_of_string o implode) + val decimalfrac = Scan.repeat1 numeral + >> (fn s => rat_of_string(implode s) // pow10 (length s)) + val decimalsig = + decimalint -- Scan.option (Scan.$$ "." |-- decimalfrac) + >> (fn (h,NONE) => h | (h,SOME x) => h +/ x) + fun signed prs = + $$ "-" |-- prs >> Rat.neg + || $$ "+" |-- prs + || prs; + +fun emptyin def xs = if null xs then (def,xs) else Scan.fail xs + + val exponent = ($$ "e" || $$ "E") |-- signed decimalint; + + val decimal = signed decimalsig -- (emptyin rat_0|| exponent) + >> (fn (h, x) => h */ pow10 (int_of_rat x)); + + fun mkparser p s = + let val (x,rst) = p (raw_explode s) + in if null rst then x + else error "mkparser: unparsed input" + end;; + +(* Parse back csdp output. *) + + fun ignore inp = ((),[]) + fun csdpoutput inp = + ((decimal -- Scan.repeat (Scan.$$ " " |-- Scan.option decimal) >> + (fn (h,to) => map_filter I ((SOME h)::to))) --| ignore >> vector_of_list) inp + val parse_csdpoutput = mkparser csdpoutput + +(* Run prover on a problem in linear form. *) + +fun run_problem prover obj mats = + parse_csdpoutput (prover (sdpa_of_problem obj mats)) + +(* Try some apparently sensible scaling first. Note that this is purely to *) +(* get a cleaner translation to floating-point, and doesn't affect any of *) +(* the results, in principle. In practice it seems a lot better when there *) +(* are extreme numbers in the original problem. *) + + (* Version for (int*int) keys *) +local + fun max_rat x y = if x fn a => lcm_rat (denominator_rat c) a) amat acc + fun maximal_element fld amat acc = + fld (fn (m,c) => fn maxa => max_rat maxa (abs_rat c)) amat acc +fun float_of_rat x = let val (a,b) = Rat.quotient_of_rat x + in Real.fromLargeInt a / Real.fromLargeInt b end; +in + +fun pi_scale_then solver (obj:vector) mats = + let + val cd1 = fold_rev (common_denominator FuncUtil.Intpairfunc.fold) mats (rat_1) + val cd2 = common_denominator FuncUtil.Intfunc.fold (snd obj) (rat_1) + val mats' = map (FuncUtil.Intpairfunc.map (fn _ => fn x => cd1 */ x)) mats + val obj' = vector_cmul cd2 obj + val max1 = fold_rev (maximal_element FuncUtil.Intpairfunc.fold) mats' (rat_0) + val max2 = maximal_element FuncUtil.Intfunc.fold (snd obj') (rat_0) + val scal1 = pow2 (20 - trunc(Math.ln (float_of_rat max1) / Math.ln 2.0)) + val scal2 = pow2 (20 - trunc(Math.ln (float_of_rat max2) / Math.ln 2.0)) + val mats'' = map (FuncUtil.Intpairfunc.map (fn _ => fn x => x */ scal1)) mats' + val obj'' = vector_cmul scal2 obj' + in solver obj'' mats'' + end +end; + +(* Try some apparently sensible scaling first. Note that this is purely to *) +(* get a cleaner translation to floating-point, and doesn't affect any of *) +(* the results, in principle. In practice it seems a lot better when there *) +(* are extreme numbers in the original problem. *) + + (* Version for (int*int*int) keys *) +local + fun max_rat x y = if x fn a => lcm_rat (denominator_rat c) a) amat acc + fun maximal_element fld amat acc = + fld (fn (m,c) => fn maxa => max_rat maxa (abs_rat c)) amat acc +fun float_of_rat x = let val (a,b) = Rat.quotient_of_rat x + in Real.fromLargeInt a / Real.fromLargeInt b end; +fun int_of_float x = (trunc x handle Overflow => 0 | Domain => 0) +in + +fun tri_scale_then solver (obj:vector) mats = + let + val cd1 = fold_rev (common_denominator Inttriplefunc.fold) mats (rat_1) + val cd2 = common_denominator FuncUtil.Intfunc.fold (snd obj) (rat_1) + val mats' = map (Inttriplefunc.map (fn _ => fn x => cd1 */ x)) mats + val obj' = vector_cmul cd2 obj + val max1 = fold_rev (maximal_element Inttriplefunc.fold) mats' (rat_0) + val max2 = maximal_element FuncUtil.Intfunc.fold (snd obj') (rat_0) + val scal1 = pow2 (20 - int_of_float(Math.ln (float_of_rat max1) / Math.ln 2.0)) + val scal2 = pow2 (20 - int_of_float(Math.ln (float_of_rat max2) / Math.ln 2.0)) + val mats'' = map (Inttriplefunc.map (fn _ => fn x => x */ scal1)) mats' + val obj'' = vector_cmul scal2 obj' + in solver obj'' mats'' + end +end; + +(* Round a vector to "nice" rationals. *) + +fun nice_rational n x = round_rat (n */ x) // n;; +fun nice_vector n ((d,v) : vector) = + (d, FuncUtil.Intfunc.fold (fn (i,c) => fn a => + let val y = nice_rational n c + in if c =/ rat_0 then a + else FuncUtil.Intfunc.update (i,y) a end) v FuncUtil.Intfunc.empty):vector + +fun dest_ord f x = is_equal (f x); + +(* Stuff for "equations" ((int*int*int)->num functions). *) + +fun tri_equation_cmul c eq = + if c =/ rat_0 then Inttriplefunc.empty else Inttriplefunc.map (fn _ => fn d => c */ d) eq; + +fun tri_equation_add eq1 eq2 = Inttriplefunc.combine (curry op +/) (fn x => x =/ rat_0) eq1 eq2; + +fun tri_equation_eval assig eq = + let fun value v = Inttriplefunc.apply assig v + in Inttriplefunc.fold (fn (v, c) => fn a => a +/ value v */ c) eq rat_0 + end; + +(* Eliminate among linear equations: return unconstrained variables and *) +(* assignments for the others in terms of them. We give one pseudo-variable *) +(* "one" that's used for a constant term. *) + +local + fun extract_first p l = case l of (* FIXME : use find_first instead *) + [] => error "extract_first" + | h::t => if p h then (h,t) else + let val (k,s) = extract_first p t in (k,h::s) end +fun eliminate vars dun eqs = case vars of + [] => if forall Inttriplefunc.is_empty eqs then dun + else raise Unsolvable + | v::vs => + ((let + val (eq,oeqs) = extract_first (fn e => Inttriplefunc.defined e v) eqs + val a = Inttriplefunc.apply eq v + val eq' = tri_equation_cmul ((Rat.neg rat_1) // a) (Inttriplefunc.delete_safe v eq) + fun elim e = + let val b = Inttriplefunc.tryapplyd e v rat_0 + in if b =/ rat_0 then e else + tri_equation_add e (tri_equation_cmul (Rat.neg b // a) eq) + end + in eliminate vs (Inttriplefunc.update (v,eq') (Inttriplefunc.map (K elim) dun)) (map elim oeqs) + end) + handle Failure _ => eliminate vs dun eqs) +in +fun tri_eliminate_equations one vars eqs = + let + val assig = eliminate vars Inttriplefunc.empty eqs + val vs = Inttriplefunc.fold (fn (x, f) => fn a => remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig [] + in (distinct (dest_ord triple_int_ord) vs, assig) + end +end; + +(* Eliminate all variables, in an essentially arbitrary order. *) + +fun tri_eliminate_all_equations one = + let + fun choose_variable eq = + let val (v,_) = Inttriplefunc.choose eq + in if is_equal (triple_int_ord(v,one)) then + let val eq' = Inttriplefunc.delete_safe v eq + in if Inttriplefunc.is_empty eq' then error "choose_variable" + else fst (Inttriplefunc.choose eq') + end + else v + end + fun eliminate dun eqs = case eqs of + [] => dun + | eq::oeqs => + if Inttriplefunc.is_empty eq then eliminate dun oeqs else + let val v = choose_variable eq + val a = Inttriplefunc.apply eq v + val eq' = tri_equation_cmul ((Rat.rat_of_int ~1) // a) + (Inttriplefunc.delete_safe v eq) + fun elim e = + let val b = Inttriplefunc.tryapplyd e v rat_0 + in if b =/ rat_0 then e + else tri_equation_add e (tri_equation_cmul (Rat.neg b // a) eq) + end + in eliminate (Inttriplefunc.update(v, eq') (Inttriplefunc.map (K elim) dun)) + (map elim oeqs) + end +in fn eqs => + let + val assig = eliminate Inttriplefunc.empty eqs + val vs = Inttriplefunc.fold (fn (x, f) => fn a => remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig [] + in (distinct (dest_ord triple_int_ord) vs,assig) + end +end; + +(* Solve equations by assigning arbitrary numbers. *) + +fun tri_solve_equations one eqs = + let + val (vars,assigs) = tri_eliminate_all_equations one eqs + val vfn = fold_rev (fn v => Inttriplefunc.update(v,rat_0)) vars + (Inttriplefunc.onefunc(one, Rat.rat_of_int ~1)) + val ass = + Inttriplefunc.combine (curry op +/) (K false) + (Inttriplefunc.map (K (tri_equation_eval vfn)) assigs) vfn + in if forall (fn e => tri_equation_eval ass e =/ rat_0) eqs + then Inttriplefunc.delete_safe one ass else raise Sanity + end; + +(* Multiply equation-parametrized poly by regular poly and add accumulator. *) + +fun tri_epoly_pmul p q acc = + FuncUtil.Monomialfunc.fold (fn (m1, c) => fn a => + FuncUtil.Monomialfunc.fold (fn (m2,e) => fn b => + let val m = monomial_mul m1 m2 + val es = FuncUtil.Monomialfunc.tryapplyd b m Inttriplefunc.empty + in FuncUtil.Monomialfunc.update (m,tri_equation_add (tri_equation_cmul c e) es) b + end) q a) p acc ; + +(* Usual operations on equation-parametrized poly. *) + +fun tri_epoly_cmul c l = + if c =/ rat_0 then Inttriplefunc.empty else Inttriplefunc.map (K (tri_equation_cmul c)) l;; + +val tri_epoly_neg = tri_epoly_cmul (Rat.rat_of_int ~1); + +val tri_epoly_add = Inttriplefunc.combine tri_equation_add Inttriplefunc.is_empty; + +fun tri_epoly_sub p q = tri_epoly_add p (tri_epoly_neg q);; + +(* Stuff for "equations" ((int*int)->num functions). *) + +fun pi_equation_cmul c eq = + if c =/ rat_0 then Inttriplefunc.empty else Inttriplefunc.map (fn _ => fn d => c */ d) eq; + +fun pi_equation_add eq1 eq2 = Inttriplefunc.combine (curry op +/) (fn x => x =/ rat_0) eq1 eq2; + +fun pi_equation_eval assig eq = + let fun value v = Inttriplefunc.apply assig v + in Inttriplefunc.fold (fn (v, c) => fn a => a +/ value v */ c) eq rat_0 + end; + +(* Eliminate among linear equations: return unconstrained variables and *) +(* assignments for the others in terms of them. We give one pseudo-variable *) +(* "one" that's used for a constant term. *) + +local +fun extract_first p l = case l of + [] => error "extract_first" + | h::t => if p h then (h,t) else + let val (k,s) = extract_first p t in (k,h::s) end +fun eliminate vars dun eqs = case vars of + [] => if forall Inttriplefunc.is_empty eqs then dun + else raise Unsolvable + | v::vs => + let + val (eq,oeqs) = extract_first (fn e => Inttriplefunc.defined e v) eqs + val a = Inttriplefunc.apply eq v + val eq' = pi_equation_cmul ((Rat.neg rat_1) // a) (Inttriplefunc.delete_safe v eq) + fun elim e = + let val b = Inttriplefunc.tryapplyd e v rat_0 + in if b =/ rat_0 then e else + pi_equation_add e (pi_equation_cmul (Rat.neg b // a) eq) + end + in eliminate vs (Inttriplefunc.update (v,eq') (Inttriplefunc.map (K elim) dun)) (map elim oeqs) + end + handle Failure _ => eliminate vs dun eqs +in +fun pi_eliminate_equations one vars eqs = + let + val assig = eliminate vars Inttriplefunc.empty eqs + val vs = Inttriplefunc.fold (fn (x, f) => fn a => remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig [] + in (distinct (dest_ord triple_int_ord) vs, assig) + end +end; + +(* Eliminate all variables, in an essentially arbitrary order. *) + +fun pi_eliminate_all_equations one = + let + fun choose_variable eq = + let val (v,_) = Inttriplefunc.choose eq + in if is_equal (triple_int_ord(v,one)) then + let val eq' = Inttriplefunc.delete_safe v eq + in if Inttriplefunc.is_empty eq' then error "choose_variable" + else fst (Inttriplefunc.choose eq') + end + else v + end + fun eliminate dun eqs = case eqs of + [] => dun + | eq::oeqs => + if Inttriplefunc.is_empty eq then eliminate dun oeqs else + let val v = choose_variable eq + val a = Inttriplefunc.apply eq v + val eq' = pi_equation_cmul ((Rat.rat_of_int ~1) // a) + (Inttriplefunc.delete_safe v eq) + fun elim e = + let val b = Inttriplefunc.tryapplyd e v rat_0 + in if b =/ rat_0 then e + else pi_equation_add e (pi_equation_cmul (Rat.neg b // a) eq) + end + in eliminate (Inttriplefunc.update(v, eq') (Inttriplefunc.map (K elim) dun)) + (map elim oeqs) + end +in fn eqs => + let + val assig = eliminate Inttriplefunc.empty eqs + val vs = Inttriplefunc.fold (fn (x, f) => fn a => remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig [] + in (distinct (dest_ord triple_int_ord) vs,assig) + end +end; + +(* Solve equations by assigning arbitrary numbers. *) + +fun pi_solve_equations one eqs = + let + val (vars,assigs) = pi_eliminate_all_equations one eqs + val vfn = fold_rev (fn v => Inttriplefunc.update(v,rat_0)) vars + (Inttriplefunc.onefunc(one, Rat.rat_of_int ~1)) + val ass = + Inttriplefunc.combine (curry op +/) (K false) + (Inttriplefunc.map (K (pi_equation_eval vfn)) assigs) vfn + in if forall (fn e => pi_equation_eval ass e =/ rat_0) eqs + then Inttriplefunc.delete_safe one ass else raise Sanity + end; + +(* Multiply equation-parametrized poly by regular poly and add accumulator. *) + +fun pi_epoly_pmul p q acc = + FuncUtil.Monomialfunc.fold (fn (m1, c) => fn a => + FuncUtil.Monomialfunc.fold (fn (m2,e) => fn b => + let val m = monomial_mul m1 m2 + val es = FuncUtil.Monomialfunc.tryapplyd b m Inttriplefunc.empty + in FuncUtil.Monomialfunc.update (m,pi_equation_add (pi_equation_cmul c e) es) b + end) q a) p acc ; + +(* Usual operations on equation-parametrized poly. *) + +fun pi_epoly_cmul c l = + if c =/ rat_0 then Inttriplefunc.empty else Inttriplefunc.map (K (pi_equation_cmul c)) l;; + +val pi_epoly_neg = pi_epoly_cmul (Rat.rat_of_int ~1); + +val pi_epoly_add = Inttriplefunc.combine pi_equation_add Inttriplefunc.is_empty; + +fun pi_epoly_sub p q = pi_epoly_add p (pi_epoly_neg q);; + +fun allpairs f l1 l2 = fold_rev (fn x => (curry (op @)) (map (f x) l2)) l1 []; + +(* Hence produce the "relevant" monomials: those whose squares lie in the *) +(* Newton polytope of the monomials in the input. (This is enough according *) +(* to Reznik: "Extremal PSD forms with few terms", Duke Math. Journal, *) +(* vol 45, pp. 363--374, 1978. *) +(* *) +(* These are ordered in sort of decreasing degree. In particular the *) +(* constant monomial is last; this gives an order in diagonalization of the *) +(* quadratic form that will tend to display constants. *) + +(* Diagonalize (Cholesky/LDU) the matrix corresponding to a quadratic form. *) + +local +fun diagonalize n i m = + if FuncUtil.Intpairfunc.is_empty (snd m) then [] + else + let val a11 = FuncUtil.Intpairfunc.tryapplyd (snd m) (i,i) rat_0 + in if a11 fn a => + let val y = c // a11 + in if y = rat_0 then a else FuncUtil.Intfunc.update (i,y) a + end) (snd v) FuncUtil.Intfunc.empty) + fun upt0 x y a = if y = rat_0 then a else FuncUtil.Intpairfunc.update (x,y) a + val m' = + ((n,n), + iter (i+1,n) (fn j => + iter (i+1,n) (fn k => + (upt0 (j,k) (FuncUtil.Intpairfunc.tryapplyd (snd m) (j,k) rat_0 -/ FuncUtil.Intfunc.tryapplyd (snd v) j rat_0 */ FuncUtil.Intfunc.tryapplyd (snd v') k rat_0)))) + FuncUtil.Intpairfunc.empty) + in (a11,v')::diagonalize n (i + 1) m' + end + end +in +fun diag m = + let + val nn = dimensions m + val n = fst nn + in if snd nn <> n then error "diagonalize: non-square matrix" + else diagonalize n 1 m + end +end; + +fun gcd_rat a b = Rat.rat_of_int (Integer.gcd (int_of_rat a) (int_of_rat b)); + +(* Adjust a diagonalization to collect rationals at the start. *) + (* FIXME : Potentially polymorphic keys, but here only: integers!! *) +local + fun upd0 x y a = if y =/ rat_0 then a else FuncUtil.Intfunc.update(x,y) a; + fun mapa f (d,v) = + (d, FuncUtil.Intfunc.fold (fn (i,c) => fn a => upd0 i (f c) a) v FuncUtil.Intfunc.empty) + fun adj (c,l) = + let val a = + FuncUtil.Intfunc.fold (fn (i,c) => fn a => lcm_rat a (denominator_rat c)) + (snd l) rat_1 // + FuncUtil.Intfunc.fold (fn (i,c) => fn a => gcd_rat a (numerator_rat c)) + (snd l) rat_0 + in ((c // (a */ a)),mapa (fn x => a */ x) l) + end +in +fun deration d = if null d then (rat_0,d) else + let val d' = map adj d + val a = fold (lcm_rat o denominator_rat o fst) d' rat_1 // + fold (gcd_rat o numerator_rat o fst) d' rat_0 + in ((rat_1 // a),map (fn (c,l) => (a */ c,l)) d') + end +end; + +(* Enumeration of monomials with given multidegree bound. *) + +fun enumerate_monomials d vars = + if d < 0 then [] + else if d = 0 then [FuncUtil.Ctermfunc.empty] + else if null vars then [monomial_1] else + let val alts = + map_range (fn k => let val oths = enumerate_monomials (d - k) (tl vars) + in map (fn ks => if k = 0 then ks else FuncUtil.Ctermfunc.update (hd vars, k) ks) oths end) (d + 1) + in flat alts + end; + +(* Enumerate products of distinct input polys with degree <= d. *) +(* We ignore any constant input polynomials. *) +(* Give the output polynomial and a record of how it was derived. *) + +fun enumerate_products d pols = +if d = 0 then [(poly_const rat_1,RealArith.Rational_lt rat_1)] +else if d < 0 then [] else +case pols of + [] => [(poly_const rat_1,RealArith.Rational_lt rat_1)] + | (p,b)::ps => + let val e = multidegree p + in if e = 0 then enumerate_products d ps else + enumerate_products d ps @ + map (fn (q,c) => (poly_mul p q,RealArith.Product(b,c))) + (enumerate_products (d - e) ps) + end + +(* Convert regular polynomial. Note that we treat (0,0,0) as -1. *) + +fun epoly_of_poly p = + FuncUtil.Monomialfunc.fold (fn (m,c) => fn a => FuncUtil.Monomialfunc.update (m, Inttriplefunc.onefunc ((0,0,0), Rat.neg c)) a) p FuncUtil.Monomialfunc.empty; + +(* String for block diagonal matrix numbered k. *) + +fun sdpa_of_blockdiagonal k m = + let + val pfx = string_of_int k ^" " + val ents = + Inttriplefunc.fold + (fn ((b,i,j),c) => fn a => if i > j then a else ((b,i,j),c)::a) + m [] + val entss = sort (triple_int_ord o pairself fst) ents + in fold_rev (fn ((b,i,j),c) => fn a => + pfx ^ string_of_int b ^ " " ^ string_of_int i ^ " " ^ string_of_int j ^ + " " ^ decimalize 20 c ^ "\n" ^ a) entss "" + end; + +(* SDPA for problem using block diagonal (i.e. multiple SDPs) *) + +fun sdpa_of_blockproblem nblocks blocksizes obj mats = + let val m = length mats - 1 + in + string_of_int m ^ "\n" ^ + string_of_int nblocks ^ "\n" ^ + (space_implode " " (map string_of_int blocksizes)) ^ + "\n" ^ + sdpa_of_vector obj ^ + fold_rev2 (fn k => fn m => fn a => sdpa_of_blockdiagonal (k - 1) m ^ a) + (1 upto length mats) mats "" + end; + +(* Run prover on a problem in block diagonal form. *) + +fun run_blockproblem prover nblocks blocksizes obj mats= + parse_csdpoutput (prover (sdpa_of_blockproblem nblocks blocksizes obj mats)) + +(* 3D versions of matrix operations to consider blocks separately. *) + +val bmatrix_add = Inttriplefunc.combine (curry op +/) (fn x => x =/ rat_0); +fun bmatrix_cmul c bm = + if c =/ rat_0 then Inttriplefunc.empty + else Inttriplefunc.map (fn _ => fn x => c */ x) bm; + +val bmatrix_neg = bmatrix_cmul (Rat.rat_of_int ~1); +fun bmatrix_sub m1 m2 = bmatrix_add m1 (bmatrix_neg m2);; + +(* Smash a block matrix into components. *) + +fun blocks blocksizes bm = + map (fn (bs,b0) => + let val m = Inttriplefunc.fold + (fn ((b,i,j),c) => fn a => if b = b0 then FuncUtil.Intpairfunc.update ((i,j),c) a else a) bm FuncUtil.Intpairfunc.empty + val d = FuncUtil.Intpairfunc.fold (fn ((i,j),c) => fn a => max a (max i j)) m 0 + in (((bs,bs),m):matrix) end) + (blocksizes ~~ (1 upto length blocksizes));; + +(* FIXME : Get rid of this !!!*) +local + fun tryfind_with msg f [] = raise Failure msg + | tryfind_with msg f (x::xs) = (f x handle Failure s => tryfind_with s f xs); +in + fun tryfind f = tryfind_with "tryfind" f +end + +(* Positiv- and Nullstellensatz. Flag "linf" forces a linear representation. *) + + +fun real_positivnullstellensatz_general ctxt prover linf d eqs leqs pol = +let + val vars = fold_rev (union (op aconvc) o poly_variables) + (pol :: eqs @ map fst leqs) [] + val monoid = if linf then + (poly_const rat_1,RealArith.Rational_lt rat_1):: + (filter (fn (p,c) => multidegree p <= d) leqs) + else enumerate_products d leqs + val nblocks = length monoid + fun mk_idmultiplier k p = + let + val e = d - multidegree p + val mons = enumerate_monomials e vars + val nons = mons ~~ (1 upto length mons) + in (mons, + fold_rev (fn (m,n) => FuncUtil.Monomialfunc.update(m,Inttriplefunc.onefunc((~k,~n,n),rat_1))) nons FuncUtil.Monomialfunc.empty) + end + + fun mk_sqmultiplier k (p,c) = + let + val e = (d - multidegree p) div 2 + val mons = enumerate_monomials e vars + val nons = mons ~~ (1 upto length mons) + in (mons, + fold_rev (fn (m1,n1) => + fold_rev (fn (m2,n2) => fn a => + let val m = monomial_mul m1 m2 + in if n1 > n2 then a else + let val c = if n1 = n2 then rat_1 else rat_2 + val e = FuncUtil.Monomialfunc.tryapplyd a m Inttriplefunc.empty + in FuncUtil.Monomialfunc.update(m, tri_equation_add (Inttriplefunc.onefunc((k,n1,n2), c)) e) a + end + end) nons) + nons FuncUtil.Monomialfunc.empty) + end + + val (sqmonlist,sqs) = split_list (map2 mk_sqmultiplier (1 upto length monoid) monoid) + val (idmonlist,ids) = split_list(map2 mk_idmultiplier (1 upto length eqs) eqs) + val blocksizes = map length sqmonlist + val bigsum = + fold_rev2 (fn p => fn q => fn a => tri_epoly_pmul p q a) eqs ids + (fold_rev2 (fn (p,c) => fn s => fn a => tri_epoly_pmul p s a) monoid sqs + (epoly_of_poly(poly_neg pol))) + val eqns = FuncUtil.Monomialfunc.fold (fn (m,e) => fn a => e::a) bigsum [] + val (pvs,assig) = tri_eliminate_all_equations (0,0,0) eqns + val qvars = (0,0,0)::pvs + val allassig = fold_rev (fn v => Inttriplefunc.update(v,(Inttriplefunc.onefunc(v,rat_1)))) pvs assig + fun mk_matrix v = + Inttriplefunc.fold (fn ((b,i,j), ass) => fn m => + if b < 0 then m else + let val c = Inttriplefunc.tryapplyd ass v rat_0 + in if c = rat_0 then m else + Inttriplefunc.update ((b,j,i), c) (Inttriplefunc.update ((b,i,j), c) m) + end) + allassig Inttriplefunc.empty + val diagents = Inttriplefunc.fold + (fn ((b,i,j), e) => fn a => if b > 0 andalso i = j then tri_equation_add e a else a) + allassig Inttriplefunc.empty + + val mats = map mk_matrix qvars + val obj = (length pvs, + itern 1 pvs (fn v => fn i => FuncUtil.Intfunc.updatep iszero (i,Inttriplefunc.tryapplyd diagents v rat_0)) + FuncUtil.Intfunc.empty) + val raw_vec = if null pvs then vector_0 0 + else tri_scale_then (run_blockproblem prover nblocks blocksizes) obj mats + fun int_element (d,v) i = FuncUtil.Intfunc.tryapplyd v i rat_0 + + fun find_rounding d = + let + val _ = + if Config.get ctxt trace + then writeln ("Trying rounding with limit "^Rat.string_of_rat d ^ "\n") + else () + val vec = nice_vector d raw_vec + val blockmat = iter (1,dim vec) + (fn i => fn a => bmatrix_add (bmatrix_cmul (int_element vec i) (nth mats i)) a) + (bmatrix_neg (nth mats 0)) + val allmats = blocks blocksizes blockmat + in (vec,map diag allmats) + end + val (vec,ratdias) = + if null pvs then find_rounding rat_1 + else tryfind find_rounding (map Rat.rat_of_int (1 upto 31) @ + map pow2 (5 upto 66)) + val newassigs = + fold_rev (fn k => Inttriplefunc.update (nth pvs (k - 1), int_element vec k)) + (1 upto dim vec) (Inttriplefunc.onefunc ((0,0,0), Rat.rat_of_int ~1)) + val finalassigs = + Inttriplefunc.fold (fn (v,e) => fn a => Inttriplefunc.update(v, tri_equation_eval newassigs e) a) allassig newassigs + fun poly_of_epoly p = + FuncUtil.Monomialfunc.fold (fn (v,e) => fn a => FuncUtil.Monomialfunc.updatep iszero (v,tri_equation_eval finalassigs e) a) + p FuncUtil.Monomialfunc.empty + fun mk_sos mons = + let fun mk_sq (c,m) = + (c,fold_rev (fn k=> fn a => FuncUtil.Monomialfunc.updatep iszero (nth mons (k - 1), int_element m k) a) + (1 upto length mons) FuncUtil.Monomialfunc.empty) + in map mk_sq + end + val sqs = map2 mk_sos sqmonlist ratdias + val cfs = map poly_of_epoly ids + val msq = filter (fn (a,b) => not (null b)) (map2 pair monoid sqs) + fun eval_sq sqs = fold_rev (fn (c,q) => poly_add (poly_cmul c (poly_mul q q))) sqs poly_0 + val sanity = + fold_rev (fn ((p,c),s) => poly_add (poly_mul p (eval_sq s))) msq + (fold_rev2 (fn p => fn q => poly_add (poly_mul p q)) cfs eqs + (poly_neg pol)) + +in if not(FuncUtil.Monomialfunc.is_empty sanity) then raise Sanity else + (cfs,map (fn (a,b) => (snd a,b)) msq) + end + + +(* Iterative deepening. *) + +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)))); + + +(* Map back polynomials and their composites to a positivstellensatz. *) + +fun cterm_of_sqterm (c,p) = RealArith.Product(RealArith.Rational_lt c,RealArith.Square p); + +fun cterm_of_sos (pr,sqs) = if null sqs then pr + else RealArith.Product(pr,foldr1 RealArith.Sum (map cterm_of_sqterm sqs)); + +(* Interface to HOL. *) +local + open Conv + val concl = Thm.dest_arg o cprop_of + fun simple_cterm_ord t u = Term_Ord.fast_term_ord (term_of t, term_of u) = LESS +in + (* FIXME: Replace tryfind by get_first !! *) +fun real_nonlinear_prover proof_method ctxt = + let + val {add,mul,neg,pow,sub,main} = Semiring_Normalizer.semiring_normalizers_ord_wrapper ctxt + (the (Semiring_Normalizer.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 cert_choice translator (eqs,les,lts) = + let + val eq0 = map (poly_of_term o Thm.dest_arg1 o concl) eqs + val le0 = map (poly_of_term o Thm.dest_arg o concl) les + val lt0 = map (poly_of_term o Thm.dest_arg o concl) lts + val eqp0 = map_index (fn (i, t) => (t,RealArith.Axiom_eq i)) eq0 + val lep0 = map_index (fn (i, t) => (t,RealArith.Axiom_le i)) le0 + val ltp0 = map_index (fn (i, t) => (t,RealArith.Axiom_lt i)) lt0 + val (keq,eq) = List.partition (fn (p,_) => multidegree p = 0) eqp0 + val (klep,lep) = List.partition (fn (p,_) => multidegree p = 0) lep0 + val (kltp,ltp) = List.partition (fn (p,_) => multidegree p = 0) ltp0 + fun trivial_axiom (p,ax) = + case ax of + RealArith.Axiom_eq n => if eval FuncUtil.Ctermfunc.empty p <>/ Rat.zero then nth eqs n + else raise Failure "trivial_axiom: Not a trivial axiom" + | RealArith.Axiom_le n => if eval FuncUtil.Ctermfunc.empty p if eval FuncUtil.Ctermfunc.empty p <=/ Rat.zero then nth lts n + 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 Numeral_Simprocs.field_comp_conv) th, RealArith.Trivial) + end) + handle Failure _ => + (let val proof = + (case proof_method of Certificate certs => + (* choose certificate *) + let + fun chose_cert [] (RealArith.Cert c) = c + | chose_cert (RealArith.Left::s) (RealArith.Branch (l, _)) = chose_cert s l + | chose_cert (RealArith.Right::s) (RealArith.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 ctxt 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) => RealArith.Eqmul(q,ax)) cert_ideal eq + val proofs_cone = map cterm_of_sos cert_cone + val proof_ne = if null ltp then RealArith.Rational_lt Rat.one else + let val p = foldr1 RealArith.Product (map snd ltp) + in funpow i (fn q => RealArith.Product(p,q)) (RealArith.Rational_lt Rat.one) + end + in + foldr1 RealArith.Sum (proof_ne :: proofs_ideal @ proofs_cone) + end) + in + (translator (eqs,les,lts) proof, RealArith.Cert proof) + end) + end + in mainf end +end + +fun C f x y = f y x; + (* FIXME : This is very bad!!!*) +fun subst_conv eqs t = + let + val t' = fold (Thm.cabs o Thm.lhs_of) eqs t + in Conv.fconv_rule (Thm.beta_conversion true) (fold (C Thm.combination) eqs (Thm.reflexive t')) + end + +(* A wrapper that tries to substitute away variables first. *) + +local + open Conv + fun simple_cterm_ord t u = Term_Ord.fast_term_ord (term_of t, term_of u) = LESS + val concl = Thm.dest_arg o cprop_of + val shuffle1 = + fconv_rule (rewr_conv @{lemma "(a + x == y) == (x == y - (a::real))" by (atomize (full)) (simp add: field_simps) }) + val shuffle2 = + fconv_rule (rewr_conv @{lemma "(x + a == y) == (x == y - (a::real))" by (atomize (full)) (simp add: field_simps)}) + fun substitutable_monomial fvs tm = case term_of tm of + Free(_,@{typ real}) => if not (member (op aconvc) fvs tm) then (Rat.one,tm) + else raise Failure "substitutable_monomial" + | @{term "op * :: real => _"}$c$(t as Free _ ) => + if RealArith.is_ratconst (Thm.dest_arg1 tm) andalso not (member (op aconvc) fvs (Thm.dest_arg tm)) + then (RealArith.dest_ratconst (Thm.dest_arg1 tm),Thm.dest_arg tm) else raise Failure "substitutable_monomial" + | @{term "op + :: real => _"}$s$t => + (substitutable_monomial (Thm.add_cterm_frees (Thm.dest_arg tm) fvs) (Thm.dest_arg1 tm) + handle Failure _ => substitutable_monomial (Thm.add_cterm_frees (Thm.dest_arg1 tm) fvs) (Thm.dest_arg tm)) + | _ => raise Failure "substitutable_monomial" + + fun isolate_variable v th = + let val w = Thm.dest_arg1 (cprop_of th) + in if v aconvc w then th + else case term_of w of + @{term "op + :: real => _"}$s$t => + if Thm.dest_arg1 w aconvc v then shuffle2 th + else isolate_variable v (shuffle1 th) + | _ => error "isolate variable : This should not happen?" + end +in + +fun real_nonlinear_subst_prover prover ctxt = + let + val {add,mul,neg,pow,sub,main} = Semiring_Normalizer.semiring_normalizers_ord_wrapper ctxt + (the (Semiring_Normalizer.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 make_substitution th = + let + val (c,v) = substitutable_monomial [] (Thm.dest_arg1(concl th)) + val th1 = Drule.arg_cong_rule (Thm.capply @{cterm "op * :: real => _"} (RealArith.cterm_of_rat (Rat.inv c))) (mk_meta_eq th) + val th2 = fconv_rule (binop_conv real_poly_mul_conv) th1 + in fconv_rule (arg_conv real_poly_conv) (isolate_variable v th2) + end + fun oprconv cv ct = + let val g = Thm.dest_fun2 ct + in if g aconvc @{cterm "op <= :: real => _"} + orelse g aconvc @{cterm "op < :: real => _"} + then arg_conv cv ct else arg1_conv cv ct + end + fun mainf cert_choice translator = + let + fun substfirst(eqs,les,lts) = + ((let + val eth = tryfind make_substitution eqs + val modify = fconv_rule (arg_conv (oprconv(subst_conv [eth] then_conv real_poly_conv))) + in substfirst + (filter_out (fn t => (Thm.dest_arg1 o Thm.dest_arg o cprop_of) t + aconvc @{cterm "0::real"}) (map modify eqs), + map modify les,map modify lts) + end) + handle Failure _ => real_nonlinear_prover prover ctxt cert_choice translator (rev eqs, rev les, rev lts)) + in substfirst + end + + + in mainf + end + +(* Overall function. *) + +fun real_sos prover ctxt = + RealArith.gen_prover_real_arith ctxt (real_nonlinear_subst_prover prover ctxt) +end; + +val known_sos_constants = + [@{term "op ==>"}, @{term "Trueprop"}, + @{term HOL.implies}, @{term HOL.conj}, @{term HOL.disj}, + @{term "Not"}, @{term "op = :: bool => _"}, + @{term "All :: (real => _) => _"}, @{term "Ex :: (real => _) => _"}, + @{term "op = :: real => _"}, @{term "op < :: real => _"}, + @{term "op <= :: real => _"}, + @{term "op + :: real => _"}, @{term "op - :: real => _"}, + @{term "op * :: real => _"}, @{term "uminus :: real => _"}, + @{term "op / :: real => _"}, @{term "inverse :: real => _"}, + @{term "op ^ :: real => _"}, @{term "abs :: real => _"}, + @{term "min :: real => _"}, @{term "max :: real => _"}, + @{term "0::real"}, @{term "1::real"}, @{term "number_of :: int => real"}, + @{term "number_of :: int => nat"}, + @{term "Int.Bit0"}, @{term "Int.Bit1"}, + @{term "Int.Pls"}, @{term "Int.Min"}]; + +fun check_sos kcts ct = + let + val t = term_of ct + val _ = if not (null (Term.add_tfrees t []) + andalso null (Term.add_tvars t [])) + then error "SOS: not sos. Additional type varables" else () + val fs = Term.add_frees t [] + 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"})) 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 print_cert prover = SUBPROOF (fn {concl, context, ...} => + let + val _ = check_sos known_sos_constants concl + val (ths, certificates) = real_sos prover context (Thm.dest_arg concl) + val _ = print_cert certificates + in rtac ths 1 end) + +fun default_SOME f NONE v = SOME v + | default_SOME f (SOME v) _ = SOME v; + +fun lift_SOME f NONE a = f a + | lift_SOME f (SOME a) _ = SOME a; + + +local + val is_numeral = can (HOLogic.dest_number o term_of) +in +fun get_denom b ct = case term_of ct of + @{term "op / :: real => _"} $ _ $ _ => + if is_numeral (Thm.dest_arg ct) then get_denom b (Thm.dest_arg1 ct) + else default_SOME (get_denom b) (get_denom b (Thm.dest_arg ct)) (Thm.dest_arg ct, b) + | @{term "op < :: real => _"} $ _ $ _ => lift_SOME (get_denom true) (get_denom true (Thm.dest_arg ct)) (Thm.dest_arg1 ct) + | @{term "op <= :: real => _"} $ _ $ _ => lift_SOME (get_denom true) (get_denom true (Thm.dest_arg ct)) (Thm.dest_arg1 ct) + | _ $ _ => lift_SOME (get_denom b) (get_denom b (Thm.dest_fun ct)) (Thm.dest_arg ct) + | _ => NONE +end; + +fun elim_one_denom_tac ctxt = +CSUBGOAL (fn (P,i) => + case get_denom false P of + NONE => no_tac + | SOME (d,ord) => + let + val ss = simpset_of ctxt addsimps @{thms field_simps} + addsimps [@{thm nonzero_power_divide}, @{thm power_divide}] + val th = instantiate' [] [SOME d, SOME (Thm.dest_arg P)] + (if ord then @{lemma "(d=0 --> P) & (d>0 --> P) & (d<(0::real) --> P) ==> P" by auto} + else @{lemma "(d=0 --> P) & (d ~= (0::real) --> P) ==> P" by blast}) + in rtac th i THEN Simplifier.asm_full_simp_tac ss i end); + +fun elim_denom_tac ctxt i = REPEAT (elim_one_denom_tac ctxt i); + +fun sos_tac print_cert prover ctxt = + Object_Logic.full_atomize_tac THEN' + elim_denom_tac ctxt THEN' + core_sos_tac print_cert prover ctxt; + +end; diff -r 655f583840d0 -r 9908cf4af394 src/HOL/Mirabelle/Tools/mirabelle.ML --- a/src/HOL/Mirabelle/Tools/mirabelle.ML Sat Jan 08 11:18:09 2011 -0800 +++ b/src/HOL/Mirabelle/Tools/mirabelle.ML Sat Jan 08 11:29:25 2011 -0800 @@ -66,7 +66,7 @@ type T = (int * run_action * done_action) list val empty = [] val extend = I - fun merge data = Library.merge (K true) data (* FIXME ?!? *) + fun merge data = Library.merge (K true) data (* FIXME exponential blowup because of (K true) *) ) diff -r 655f583840d0 -r 9908cf4af394 src/HOL/Statespace/state_fun.ML --- a/src/HOL/Statespace/state_fun.ML Sat Jan 08 11:18:09 2011 -0800 +++ b/src/HOL/Statespace/state_fun.ML Sat Jan 08 11:29:25 2011 -0800 @@ -101,7 +101,7 @@ val empty = (empty_ss, empty_ss, false); val extend = I; fun merge ((ss1, ex_ss1, b1), (ss2, ex_ss2, b2)) = - (merge_ss (ss1, ss2), merge_ss (ex_ss1, ex_ss2), b1 orelse b2); + (merge_ss (ss1, ss2), merge_ss (ex_ss1, ex_ss2), b1 orelse b2 (* FIXME odd merge *)); ); val init_state_fun_data = diff -r 655f583840d0 -r 9908cf4af394 src/HOL/Statespace/state_space.ML --- a/src/HOL/Statespace/state_space.ML Sat Jan 08 11:18:09 2011 -0800 +++ b/src/HOL/Statespace/state_space.ML Sat Jan 08 11:29:25 2011 -0800 @@ -115,7 +115,7 @@ {declinfo=declinfo2, distinctthm=distinctthm2, silent=silent2}) : T = {declinfo = Termtab.merge (K true) (declinfo1, declinfo2), distinctthm = Symtab.merge (K true) (distinctthm1, distinctthm2), - silent = silent1 andalso silent2} + silent = silent1 andalso silent2 (* FIXME odd merge *)} ); fun make_namespace_data declinfo distinctthm silent = diff -r 655f583840d0 -r 9908cf4af394 src/HOL/Tools/Function/partial_function.ML --- a/src/HOL/Tools/Function/partial_function.ML Sat Jan 08 11:18:09 2011 -0800 +++ b/src/HOL/Tools/Function/partial_function.ML Sat Jan 08 11:29:25 2011 -0800 @@ -27,7 +27,7 @@ type T = ((term * term) * thm) Symtab.table; val empty = Symtab.empty; val extend = I; - fun merge (a, b) = Symtab.merge (K true) (a, b); + fun merge data = Symtab.merge (K true) data; ) fun init fixp mono fixp_eq phi = diff -r 655f583840d0 -r 9908cf4af394 src/HOL/Tools/Function/scnp_solve.ML --- a/src/HOL/Tools/Function/scnp_solve.ML Sat Jan 08 11:18:09 2011 -0800 +++ b/src/HOL/Tools/Function/scnp_solve.ML Sat Jan 08 11:29:25 2011 -0800 @@ -51,11 +51,11 @@ (** Propositional formulas **) -val Not = PropLogic.Not and And = PropLogic.And and Or = PropLogic.Or -val BoolVar = PropLogic.BoolVar +val Not = Prop_Logic.Not and And = Prop_Logic.And and Or = Prop_Logic.Or +val BoolVar = Prop_Logic.BoolVar fun Implies (p, q) = Or (Not p, q) fun Equiv (p, q) = And (Implies (p, q), Implies (q, p)) -val all = PropLogic.all +val all = Prop_Logic.all (* finite indexed quantifiers: @@ -64,7 +64,7 @@ 0<=i fn y => 2 * y + (if assign (TAG (i, j) x) then 1 else 0)) (bits - 1 downto 0) 0) diff -r 655f583840d0 -r 9908cf4af394 src/HOL/Tools/Nitpick/nitpick_hol.ML --- a/src/HOL/Tools/Nitpick/nitpick_hol.ML Sat Jan 08 11:18:09 2011 -0800 +++ b/src/HOL/Tools/Nitpick/nitpick_hol.ML Sat Jan 08 11:29:25 2011 -0800 @@ -275,7 +275,8 @@ datatype boxability = InConstr | InSel | InExpr | InPair | InFunLHS | InFunRHS1 | InFunRHS2 -structure Data = Generic_Data( +structure Data = Generic_Data +( type T = {frac_types: (string * (string * string) list) list, codatatypes: (string * (string * styp list)) list} val empty = {frac_types = [], codatatypes = []} @@ -283,7 +284,8 @@ fun merge ({frac_types = fs1, codatatypes = cs1}, {frac_types = fs2, codatatypes = cs2}) : T = {frac_types = AList.merge (op =) (K true) (fs1, fs2), - codatatypes = AList.merge (op =) (K true) (cs1, cs2)}) + codatatypes = AList.merge (op =) (K true) (cs1, cs2)} +) val name_sep = "$" val numeral_prefix = nitpick_prefix ^ "num" ^ name_sep diff -r 655f583840d0 -r 9908cf4af394 src/HOL/Tools/Nitpick/nitpick_isar.ML --- a/src/HOL/Tools/Nitpick/nitpick_isar.ML Sat Jan 08 11:18:09 2011 -0800 +++ b/src/HOL/Tools/Nitpick/nitpick_isar.ML Sat Jan 08 11:29:25 2011 -0800 @@ -130,11 +130,13 @@ else [(name, value)] -structure Data = Theory_Data( +structure Data = Theory_Data +( type T = raw_param list val empty = map (apsnd single) default_default_params val extend = I - fun merge (x, y) = AList.merge (op =) (K true) (x, y)) + fun merge data = AList.merge (op =) (K true) data +) val set_default_raw_param = Data.map o fold (AList.update (op =)) o normalize_raw_param diff -r 655f583840d0 -r 9908cf4af394 src/HOL/Tools/Nitpick/nitpick_model.ML --- a/src/HOL/Tools/Nitpick/nitpick_model.ML Sat Jan 08 11:18:09 2011 -0800 +++ b/src/HOL/Tools/Nitpick/nitpick_model.ML Sat Jan 08 11:29:25 2011 -0800 @@ -63,11 +63,13 @@ type term_postprocessor = Proof.context -> string -> (typ -> term list) -> typ -> term -> term -structure Data = Generic_Data( +structure Data = Generic_Data +( type T = (typ * term_postprocessor) list val empty = [] val extend = I - fun merge (x, y) = AList.merge (op =) (K true) (x, y)) + fun merge data = AList.merge (op =) (K true) data +) fun xsym s s' () = if print_mode_active Symbol.xsymbolsN then s else s' diff -r 655f583840d0 -r 9908cf4af394 src/HOL/Tools/Nitpick/nitpick_mono.ML --- a/src/HOL/Tools/Nitpick/nitpick_mono.ML Sat Jan 08 11:18:09 2011 -0800 +++ b/src/HOL/Tools/Nitpick/nitpick_mono.ML Sat Jan 08 11:29:25 2011 -0800 @@ -20,8 +20,6 @@ open Nitpick_Util open Nitpick_HOL -structure PL = PropLogic - datatype sign = Plus | Minus type var = int @@ -474,36 +472,36 @@ val bools_from_annotation = AList.lookup (op =) bool_table #> the val annotation_from_bools = AList.find (op =) bool_table #> the_single -fun prop_for_bool b = if b then PL.True else PL.False +fun prop_for_bool b = if b then Prop_Logic.True else Prop_Logic.False fun prop_for_bool_var_equality (v1, v2) = - PL.SAnd (PL.SOr (PL.BoolVar v1, PL.SNot (PL.BoolVar v2)), - PL.SOr (PL.SNot (PL.BoolVar v1), PL.BoolVar v2)) + Prop_Logic.SAnd (Prop_Logic.SOr (Prop_Logic.BoolVar v1, Prop_Logic.SNot (Prop_Logic.BoolVar v2)), + Prop_Logic.SOr (Prop_Logic.SNot (Prop_Logic.BoolVar v1), Prop_Logic.BoolVar v2)) fun prop_for_assign (x, a) = let val (b1, b2) = bools_from_annotation a in - PL.SAnd (PL.BoolVar (fst_var x) |> not b1 ? PL.SNot, - PL.BoolVar (snd_var x) |> not b2 ? PL.SNot) + Prop_Logic.SAnd (Prop_Logic.BoolVar (fst_var x) |> not b1 ? Prop_Logic.SNot, + Prop_Logic.BoolVar (snd_var x) |> not b2 ? Prop_Logic.SNot) end fun prop_for_assign_literal (x, (Plus, a)) = prop_for_assign (x, a) - | prop_for_assign_literal (x, (Minus, a)) = PL.SNot (prop_for_assign (x, a)) + | prop_for_assign_literal (x, (Minus, a)) = Prop_Logic.SNot (prop_for_assign (x, a)) fun prop_for_atom_assign (A a', a) = prop_for_bool (a = a') | prop_for_atom_assign (V x, a) = prop_for_assign_literal (x, (Plus, a)) fun prop_for_atom_equality (aa1, A a2) = prop_for_atom_assign (aa1, a2) | prop_for_atom_equality (A a1, aa2) = prop_for_atom_assign (aa2, a1) | prop_for_atom_equality (V x1, V x2) = - PL.SAnd (prop_for_bool_var_equality (pairself fst_var (x1, x2)), + Prop_Logic.SAnd (prop_for_bool_var_equality (pairself fst_var (x1, x2)), prop_for_bool_var_equality (pairself snd_var (x1, x2))) -val prop_for_assign_clause = PL.exists o map prop_for_assign_literal +val prop_for_assign_clause = Prop_Logic.exists o map prop_for_assign_literal fun prop_for_exists_var_assign_literal xs a = - PL.exists (map (fn x => prop_for_assign_literal (x, (Plus, a))) xs) + Prop_Logic.exists (map (fn x => prop_for_assign_literal (x, (Plus, a))) xs) fun prop_for_comp (aa1, aa2, Eq, []) = - PL.SAnd (prop_for_comp (aa1, aa2, Leq, []), + Prop_Logic.SAnd (prop_for_comp (aa1, aa2, Leq, []), prop_for_comp (aa2, aa1, Leq, [])) | prop_for_comp (aa1, aa2, Neq, []) = - PL.SNot (prop_for_comp (aa1, aa2, Eq, [])) + Prop_Logic.SNot (prop_for_comp (aa1, aa2, Eq, [])) | prop_for_comp (aa1, aa2, Leq, []) = - PL.SOr (prop_for_atom_equality (aa1, aa2), prop_for_atom_assign (aa2, Gen)) + Prop_Logic.SOr (prop_for_atom_equality (aa1, aa2), prop_for_atom_assign (aa2, Gen)) | prop_for_comp (aa1, aa2, cmp, xs) = - PL.SOr (prop_for_exists_var_assign_literal xs Gen, + Prop_Logic.SOr (prop_for_exists_var_assign_literal xs Gen, prop_for_comp (aa1, aa2, cmp, [])) fun extract_assigns max_var assigns asgs = @@ -542,11 +540,11 @@ SOME (extract_assigns max_var assigns asgs |> tap print_solution) val _ = print_problem comps clauses val prop = - PL.all (map prop_for_comp comps @ map prop_for_assign_clause clauses) + Prop_Logic.all (map prop_for_comp comps @ map prop_for_assign_clause clauses) in - if PL.eval (K false) prop then + if Prop_Logic.eval (K false) prop then do_assigns (K (SOME false)) - else if PL.eval (K true) prop then + else if Prop_Logic.eval (K true) prop then do_assigns (K (SOME true)) else let diff -r 655f583840d0 -r 9908cf4af394 src/HOL/Tools/Predicate_Compile/code_prolog.ML --- a/src/HOL/Tools/Predicate_Compile/code_prolog.ML Sat Jan 08 11:18:09 2011 -0800 +++ b/src/HOL/Tools/Predicate_Compile/code_prolog.ML Sat Jan 08 11:29:25 2011 -0800 @@ -81,7 +81,7 @@ fun merge_global_limit (NONE, NONE) = NONE | merge_global_limit (NONE, SOME n) = SOME n | merge_global_limit (SOME n, NONE) = SOME n - | merge_global_limit (SOME n, SOME m) = SOME (Int.max (n, m)) + | merge_global_limit (SOME n, SOME m) = SOME (Int.max (n, m)) (* FIXME odd merge *) structure Options = Theory_Data ( @@ -96,7 +96,7 @@ {ensure_groundness = ensure_groundness2, limit_globally = limit_globally2, limited_types = limited_types2, limited_predicates = limited_predicates2, replacing = replacing2, manual_reorder = manual_reorder2}) = - {ensure_groundness = ensure_groundness1 orelse ensure_groundness2, + {ensure_groundness = ensure_groundness1 orelse ensure_groundness2 (* FIXME odd merge *), limit_globally = merge_global_limit (limit_globally1, limit_globally2), limited_types = AList.merge (op =) (K true) (limited_types1, limited_types2), limited_predicates = AList.merge (op =) (K true) (limited_predicates1, limited_predicates2), @@ -122,9 +122,7 @@ type T = system_configuration val empty = {timeout = seconds 10.0, prolog_system = SWI_PROLOG} val extend = I; - fun merge ({timeout = timeout1, prolog_system = prolog_system1}, - {timeout = timeout2, prolog_system = prolog_system2}) = - {timeout = timeout1, prolog_system = prolog_system1} + fun merge (a, _) = a ) (* general string functions *) diff -r 655f583840d0 -r 9908cf4af394 src/HOL/Tools/Predicate_Compile/predicate_compile_core.ML --- a/src/HOL/Tools/Predicate_Compile/predicate_compile_core.ML Sat Jan 08 11:18:09 2011 -0800 +++ b/src/HOL/Tools/Predicate_Compile/predicate_compile_core.ML Sat Jan 08 11:29:25 2011 -0800 @@ -1625,44 +1625,60 @@ (* transformation for code generation *) -structure Pred_Result = Proof_Data ( +(* FIXME just one data slot (record) per program unit *) + +structure Pred_Result = Proof_Data +( type T = unit -> term Predicate.pred + (* FIXME avoid user error with non-user text *) fun init _ () = error "Pred_Result" ); val put_pred_result = Pred_Result.put; -structure Pred_Random_Result = Proof_Data ( +structure Pred_Random_Result = Proof_Data +( type T = unit -> int * int -> term Predicate.pred * (int * int) + (* FIXME avoid user error with non-user text *) fun init _ () = error "Pred_Random_Result" ); val put_pred_random_result = Pred_Random_Result.put; -structure Dseq_Result = Proof_Data ( +structure Dseq_Result = Proof_Data +( type T = unit -> term DSequence.dseq + (* FIXME avoid user error with non-user text *) fun init _ () = error "Dseq_Result" ); val put_dseq_result = Dseq_Result.put; -structure Dseq_Random_Result = Proof_Data ( +structure Dseq_Random_Result = Proof_Data +( type T = unit -> int -> int -> int * int -> term DSequence.dseq * (int * int) + (* FIXME avoid user error with non-user text *) fun init _ () = error "Dseq_Random_Result" ); val put_dseq_random_result = Dseq_Random_Result.put; -structure New_Dseq_Result = Proof_Data ( +structure New_Dseq_Result = Proof_Data +( type T = unit -> int -> term Lazy_Sequence.lazy_sequence + (* FIXME avoid user error with non-user text *) fun init _ () = error "New_Dseq_Random_Result" ); val put_new_dseq_result = New_Dseq_Result.put; -structure Lseq_Random_Result = Proof_Data ( +structure Lseq_Random_Result = Proof_Data +( type T = unit -> int -> int -> int * int -> int -> term Lazy_Sequence.lazy_sequence + (* FIXME avoid user error with non-user text *) fun init _ () = error "Lseq_Random_Result" ); val put_lseq_random_result = Lseq_Random_Result.put; -structure Lseq_Random_Stats_Result = Proof_Data ( +structure Lseq_Random_Stats_Result = Proof_Data +( type T = unit -> int -> int -> int * int -> int -> (term * int) Lazy_Sequence.lazy_sequence + (* FIXME avoid user error with non-user text *) fun init _ () = error "Lseq_Random_Stats_Result" ); val put_lseq_random_stats_result = Lseq_Random_Stats_Result.put; diff -r 655f583840d0 -r 9908cf4af394 src/HOL/Tools/Predicate_Compile/predicate_compile_quickcheck.ML --- a/src/HOL/Tools/Predicate_Compile/predicate_compile_quickcheck.ML Sat Jan 08 11:18:09 2011 -0800 +++ b/src/HOL/Tools/Predicate_Compile/predicate_compile_quickcheck.ML Sat Jan 08 11:29:25 2011 -0800 @@ -35,26 +35,36 @@ open Predicate_Compile_Aux; -structure Pred_Result = Proof_Data ( +(* FIXME just one data slot (record) per program unit *) + +structure Pred_Result = Proof_Data +( type T = unit -> int -> int -> int -> int * int -> term list Predicate.pred + (* FIXME avoid user error with non-user text *) fun init _ () = error "Pred_Result" ); val put_pred_result = Pred_Result.put; -structure Dseq_Result = Proof_Data ( +structure Dseq_Result = Proof_Data +( type T = unit -> int -> int -> int * int -> term list DSequence.dseq * (int * int) + (* FIXME avoid user error with non-user text *) fun init _ () = error "Dseq_Result" ); val put_dseq_result = Dseq_Result.put; -structure Lseq_Result = Proof_Data ( +structure Lseq_Result = Proof_Data +( type T = unit -> int -> int -> int * int -> int -> term list Lazy_Sequence.lazy_sequence + (* FIXME avoid user error with non-user text *) fun init _ () = error "Lseq_Result" ); val put_lseq_result = Lseq_Result.put; -structure New_Dseq_Result = Proof_Data ( +structure New_Dseq_Result = Proof_Data +( type T = unit -> int -> term list Lazy_Sequence.lazy_sequence + (* FIXME avoid user error with non-user text *) fun init _ () = error "New_Dseq_Random_Result" ); val put_new_dseq_result = New_Dseq_Result.put; diff -r 655f583840d0 -r 9908cf4af394 src/HOL/Tools/Qelim/cooper.ML --- a/src/HOL/Tools/Qelim/cooper.ML Sat Jan 08 11:18:09 2011 -0800 +++ b/src/HOL/Tools/Qelim/cooper.ML Sat Jan 08 11:29:25 2011 -0800 @@ -51,7 +51,7 @@ ( type T = simpset * term list; val empty = (HOL_ss, allowed_consts); - val extend = I; + val extend = I; fun merge ((ss1, ts1), (ss2, ts2)) = (merge_ss (ss1, ss2), Library.merge (op aconv) (ts1, ts2)); ); diff -r 655f583840d0 -r 9908cf4af394 src/HOL/Tools/Quotient/quotient_info.ML --- a/src/HOL/Tools/Quotient/quotient_info.ML Sat Jan 08 11:18:09 2011 -0800 +++ b/src/HOL/Tools/Quotient/quotient_info.ML Sat Jan 08 11:29:25 2011 -0800 @@ -55,6 +55,8 @@ (** data containers **) +(* FIXME just one data slot (record) per program unit *) + (* info about map- and rel-functions for a type *) type maps_info = {mapfun: string, relmap: string} diff -r 655f583840d0 -r 9908cf4af394 src/HOL/Tools/SMT/smt_builtin.ML --- a/src/HOL/Tools/SMT/smt_builtin.ML Sat Jan 08 11:18:09 2011 -0800 +++ b/src/HOL/Tools/SMT/smt_builtin.ML Sat Jan 08 11:29:25 2011 -0800 @@ -70,7 +70,7 @@ (Ord_List.insert typ_ord (perhaps (try Logic.varifyT_global) T, f)) fun merge_ttab ttabp = - SMT_Utils.dict_merge (uncurry (Ord_List.union typ_ord) o swap) ttabp + SMT_Utils.dict_merge (Ord_List.merge typ_ord) ttabp fun lookup_ttab ctxt ttab T = let fun match (U, _) = Sign.typ_instance (ProofContext.theory_of ctxt) (T, U) diff -r 655f583840d0 -r 9908cf4af394 src/HOL/Tools/SMT/smt_config.ML --- a/src/HOL/Tools/SMT/smt_config.ML Sat Jan 08 11:18:09 2011 -0800 +++ b/src/HOL/Tools/SMT/smt_config.ML Sat Jan 08 11:29:25 2011 -0800 @@ -68,7 +68,7 @@ type T = (solver_info * string list) Symtab.table * string option val empty = (Symtab.empty, NONE) val extend = I - fun merge ((ss1, s), (ss2, _)) = (Symtab.merge (K true) (ss1, ss2), s) + fun merge ((ss1, s), (ss2, _)) = (Symtab.merge (K true) (ss1, ss2), s (* FIXME merge options!? *)) ) fun set_solver_options (name, options) = @@ -227,7 +227,7 @@ type T = Cache_IO.cache option val empty = NONE val extend = I - fun merge (s, _) = s + fun merge (s, _) = s (* FIXME merge options!? *) ) val get_certificates_path = diff -r 655f583840d0 -r 9908cf4af394 src/HOL/Tools/SMT/smt_solver.ML --- a/src/HOL/Tools/SMT/smt_solver.ML Sat Jan 08 11:18:09 2011 -0800 +++ b/src/HOL/Tools/SMT/smt_solver.ML Sat Jan 08 11:29:25 2011 -0800 @@ -176,6 +176,7 @@ val extend = I fun merge data = Symtab.merge (K true) data handle Symtab.DUP name => error ("Duplicate SMT solver: " ^ quote name) + (* FIXME never happens because of (K true) *) ) local diff -r 655f583840d0 -r 9908cf4af394 src/HOL/Tools/SMT/smtlib_interface.ML --- a/src/HOL/Tools/SMT/smtlib_interface.ML Sat Jan 08 11:18:09 2011 -0800 +++ b/src/HOL/Tools/SMT/smtlib_interface.ML Sat Jan 08 11:29:25 2011 -0800 @@ -81,7 +81,7 @@ type T = (int * (term list -> string option)) list val empty = [] val extend = I - fun merge (bs1, bs2) = Ord_List.union fst_int_ord bs2 bs1 + fun merge data = Ord_List.merge fst_int_ord data ) fun add_logic pf = Logics.map (Ord_List.insert fst_int_ord pf) diff -r 655f583840d0 -r 9908cf4af394 src/HOL/Tools/SMT/z3_interface.ML --- a/src/HOL/Tools/SMT/z3_interface.ML Sat Jan 08 11:18:09 2011 -0800 +++ b/src/HOL/Tools/SMT/z3_interface.ML Sat Jan 08 11:29:25 2011 -0800 @@ -119,7 +119,7 @@ type T = (int * mk_builtins) list val empty = [] val extend = I - fun merge (bs1, bs2) = Ord_List.union fst_int_ord bs2 bs1 + fun merge data = Ord_List.merge fst_int_ord data ) fun add_mk_builtins mk = diff -r 655f583840d0 -r 9908cf4af394 src/HOL/Tools/Sledgehammer/sledgehammer_isar.ML --- a/src/HOL/Tools/Sledgehammer/sledgehammer_isar.ML Sat Jan 08 11:18:09 2011 -0800 +++ b/src/HOL/Tools/Sledgehammer/sledgehammer_isar.ML Sat Jan 08 11:29:25 2011 -0800 @@ -133,11 +133,13 @@ | _ => value) | NONE => (name, value) -structure Data = Theory_Data( +structure Data = Theory_Data +( type T = raw_param list val empty = default_default_params |> map (apsnd single) val extend = I - fun merge p : T = AList.merge (op =) (K true) p) + fun merge data : T = AList.merge (op =) (K true) data +) fun remotify_prover_if_available_and_not_already_remote ctxt name = if String.isPrefix remote_prefix name then diff -r 655f583840d0 -r 9908cf4af394 src/HOL/Tools/code_evaluation.ML --- a/src/HOL/Tools/code_evaluation.ML Sat Jan 08 11:18:09 2011 -0800 +++ b/src/HOL/Tools/code_evaluation.ML Sat Jan 08 11:29:25 2011 -0800 @@ -156,8 +156,10 @@ (** evaluation **) -structure Evaluation = Proof_Data ( +structure Evaluation = Proof_Data +( type T = unit -> term + (* FIXME avoid user error with non-user text *) fun init _ () = error "Evaluation" ); val put_term = Evaluation.put; diff -r 655f583840d0 -r 9908cf4af394 src/HOL/Tools/inductive_codegen.ML --- a/src/HOL/Tools/inductive_codegen.ML Sat Jan 08 11:18:09 2011 -0800 +++ b/src/HOL/Tools/inductive_codegen.ML Sat Jan 08 11:29:25 2011 -0800 @@ -31,8 +31,9 @@ val empty = {intros = Symtab.empty, graph = Graph.empty, eqns = Symtab.empty}; val extend = I; - fun merge ({intros=intros1, graph=graph1, eqns=eqns1}, - {intros=intros2, graph=graph2, eqns=eqns2}) : T = + fun merge + ({intros = intros1, graph = graph1, eqns = eqns1}, + {intros = intros2, graph = graph2, eqns = eqns2}) : T = {intros = merge_rules (intros1, intros2), graph = Graph.merge (K true) (graph1, graph2), eqns = merge_rules (eqns1, eqns2)}; diff -r 655f583840d0 -r 9908cf4af394 src/HOL/Tools/inductive_set.ML --- a/src/HOL/Tools/inductive_set.ML Sat Jan 08 11:18:09 2011 -0800 +++ b/src/HOL/Tools/inductive_set.ML Sat Jan 08 11:29:25 2011 -0800 @@ -175,8 +175,8 @@ set_arities = set_arities2, pred_arities = pred_arities2}) : T = {to_set_simps = Thm.merge_thms (to_set_simps1, to_set_simps2), to_pred_simps = Thm.merge_thms (to_pred_simps1, to_pred_simps2), - set_arities = Symtab.merge_list op = (set_arities1, set_arities2), - pred_arities = Symtab.merge_list op = (pred_arities1, pred_arities2)}; + set_arities = Symtab.merge_list (op =) (set_arities1, set_arities2), + pred_arities = Symtab.merge_list (op =) (pred_arities1, pred_arities2)}; ); fun name_type_of (Free p) = SOME p diff -r 655f583840d0 -r 9908cf4af394 src/HOL/Tools/prop_logic.ML --- a/src/HOL/Tools/prop_logic.ML Sat Jan 08 11:18:09 2011 -0800 +++ b/src/HOL/Tools/prop_logic.ML Sat Jan 08 11:29:25 2011 -0800 @@ -42,7 +42,7 @@ val term_of_prop_formula: prop_formula -> term end; -structure PropLogic : PROP_LOGIC = +structure Prop_Logic : PROP_LOGIC = struct (* ------------------------------------------------------------------------- *) diff -r 655f583840d0 -r 9908cf4af394 src/HOL/Tools/quickcheck_generators.ML --- a/src/HOL/Tools/quickcheck_generators.ML Sat Jan 08 11:18:09 2011 -0800 +++ b/src/HOL/Tools/quickcheck_generators.ML Sat Jan 08 11:29:25 2011 -0800 @@ -307,14 +307,20 @@ (** building and compiling generator expressions **) -structure Counterexample = Proof_Data ( +(* FIXME just one data slot (record) per program unit *) + +structure Counterexample = Proof_Data +( type T = unit -> int -> int * int -> term list option * (int * int) + (* FIXME avoid user error with non-user text *) fun init _ () = error "Counterexample" ); val put_counterexample = Counterexample.put; -structure Counterexample_Report = Proof_Data ( +structure Counterexample_Report = Proof_Data +( type T = unit -> int -> seed -> (term list option * (bool list * bool)) * seed + (* FIXME avoid user error with non-user text *) fun init _ () = error "Counterexample_Report" ); val put_counterexample_report = Counterexample_Report.put; diff -r 655f583840d0 -r 9908cf4af394 src/HOL/Tools/refute.ML --- a/src/HOL/Tools/refute.ML Sat Jan 08 11:18:09 2011 -0800 +++ b/src/HOL/Tools/refute.ML Sat Jan 08 11:29:25 2011 -0800 @@ -77,7 +77,7 @@ structure Refute : REFUTE = struct -open PropLogic; +open Prop_Logic; (* We use 'REFUTE' only for internal error conditions that should *) (* never occur in the first place (i.e. errors caused by bugs in our *) @@ -221,7 +221,7 @@ {interpreters = in2, printers = pr2, parameters = pa2}) : T = {interpreters = AList.merge (op =) (K true) (in1, in2), printers = AList.merge (op =) (K true) (pr1, pr2), - parameters = Symtab.merge (op=) (pa1, pa2)}; + parameters = Symtab.merge (op =) (pa1, pa2)}; ); val get_data = Data.get o ProofContext.theory_of; @@ -1118,8 +1118,8 @@ (* make 'u' either true or false, and make all axioms true, and *) (* add the well-formedness side condition *) val fm_u = (if negate then toFalse else toTrue) (hd intrs) - val fm_ax = PropLogic.all (map toTrue (tl intrs)) - val fm = PropLogic.all [#wellformed args, fm_ax, fm_u] + val fm_ax = Prop_Logic.all (map toTrue (tl intrs)) + val fm = Prop_Logic.all [#wellformed args, fm_ax, fm_u] val _ = (if satsolver = "dpll" orelse satsolver = "enumerate" then warning ("Using SAT solver " ^ quote satsolver ^ @@ -1394,22 +1394,22 @@ (case i1 of Leaf xs => (case i2 of - Leaf ys => PropLogic.dot_product (xs, ys) (* defined and equal *) + Leaf ys => Prop_Logic.dot_product (xs, ys) (* defined and equal *) | Node _ => raise REFUTE ("make_equality", "second interpretation is higher")) | Node xs => (case i2 of Leaf _ => raise REFUTE ("make_equality", "first interpretation is higher") - | Node ys => PropLogic.all (map equal (xs ~~ ys)))) + | Node ys => Prop_Logic.all (map equal (xs ~~ ys)))) (* interpretation * interpretation -> prop_formula *) fun not_equal (i1, i2) = (case i1 of Leaf xs => (case i2 of (* defined and not equal *) - Leaf ys => PropLogic.all ((PropLogic.exists xs) - :: (PropLogic.exists ys) + Leaf ys => Prop_Logic.all ((Prop_Logic.exists xs) + :: (Prop_Logic.exists ys) :: (map (fn (x,y) => SOr (SNot x, SNot y)) (xs ~~ ys))) | Node _ => raise REFUTE ("make_equality", "second interpretation is higher")) @@ -1417,7 +1417,7 @@ (case i2 of Leaf _ => raise REFUTE ("make_equality", "first interpretation is higher") - | Node ys => PropLogic.exists (map not_equal (xs ~~ ys)))) + | Node ys => Prop_Logic.exists (map not_equal (xs ~~ ys)))) in (* a value may be undefined; therefore 'not_equal' is not just the *) (* negation of 'equal' *) @@ -1443,15 +1443,15 @@ Leaf xs => (case i2 of (* defined and equal, or both undefined *) - Leaf ys => SOr (PropLogic.dot_product (xs, ys), - SAnd (PropLogic.all (map SNot xs), PropLogic.all (map SNot ys))) + Leaf ys => SOr (Prop_Logic.dot_product (xs, ys), + SAnd (Prop_Logic.all (map SNot xs), Prop_Logic.all (map SNot ys))) | Node _ => raise REFUTE ("make_def_equality", "second interpretation is higher")) | Node xs => (case i2 of Leaf _ => raise REFUTE ("make_def_equality", "first interpretation is higher") - | Node ys => PropLogic.all (map equal (xs ~~ ys)))) + | Node ys => Prop_Logic.all (map equal (xs ~~ ys)))) (* interpretation *) val eq = equal (i1, i2) in @@ -1499,7 +1499,7 @@ (* interpretation -> prop_formula list *) fun interpretation_to_prop_formula_list (Leaf xs) = xs | interpretation_to_prop_formula_list (Node trees) = - map PropLogic.all (pick_all + map Prop_Logic.all (pick_all (map interpretation_to_prop_formula_list trees)) in case i1 of @@ -1571,7 +1571,7 @@ val intr = Leaf fms (* prop_formula list -> prop_formula *) fun one_of_two_false [] = True - | one_of_two_false (x::xs) = SAnd (PropLogic.all (map (fn x' => + | one_of_two_false (x::xs) = SAnd (Prop_Logic.all (map (fn x' => SOr (SNot x, SNot x')) xs), one_of_two_false xs) (* prop_formula *) val wf = one_of_two_false fms @@ -1672,8 +1672,8 @@ Node xs => (* 3-valued logic *) let - val fmTrue = PropLogic.all (map toTrue xs) - val fmFalse = PropLogic.exists (map toFalse xs) + val fmTrue = Prop_Logic.all (map toTrue xs) + val fmFalse = Prop_Logic.exists (map toFalse xs) in SOME (Leaf [fmTrue, fmFalse], m, a) end @@ -1701,8 +1701,8 @@ let val (i1, m1, a1) = interpret ctxt model args t1 val (i2, m2, a2) = interpret ctxt m1 a1 t2 - val fmTrue = PropLogic.SOr (toFalse i1, toTrue i2) - val fmFalse = PropLogic.SAnd (toTrue i1, toFalse i2) + val fmTrue = Prop_Logic.SOr (toFalse i1, toTrue i2) + val fmFalse = Prop_Logic.SAnd (toTrue i1, toFalse i2) in SOME (Leaf [fmTrue, fmFalse], m2, a2) end @@ -1735,8 +1735,8 @@ Node xs => (* 3-valued logic *) let - val fmTrue = PropLogic.all (map toTrue xs) - val fmFalse = PropLogic.exists (map toFalse xs) + val fmTrue = Prop_Logic.all (map toTrue xs) + val fmFalse = Prop_Logic.exists (map toFalse xs) in SOME (Leaf [fmTrue, fmFalse], m, a) end @@ -1754,8 +1754,8 @@ Node xs => (* 3-valued logic *) let - val fmTrue = PropLogic.exists (map toTrue xs) - val fmFalse = PropLogic.all (map toFalse xs) + val fmTrue = Prop_Logic.exists (map toTrue xs) + val fmFalse = Prop_Logic.all (map toFalse xs) in SOME (Leaf [fmTrue, fmFalse], m, a) end @@ -1781,8 +1781,8 @@ let val (i1, m1, a1) = interpret ctxt model args t1 val (i2, m2, a2) = interpret ctxt m1 a1 t2 - val fmTrue = PropLogic.SAnd (toTrue i1, toTrue i2) - val fmFalse = PropLogic.SOr (toFalse i1, toFalse i2) + val fmTrue = Prop_Logic.SAnd (toTrue i1, toTrue i2) + val fmFalse = Prop_Logic.SOr (toFalse i1, toFalse i2) in SOME (Leaf [fmTrue, fmFalse], m2, a2) end @@ -1798,8 +1798,8 @@ let val (i1, m1, a1) = interpret ctxt model args t1 val (i2, m2, a2) = interpret ctxt m1 a1 t2 - val fmTrue = PropLogic.SOr (toTrue i1, toTrue i2) - val fmFalse = PropLogic.SAnd (toFalse i1, toFalse i2) + val fmTrue = Prop_Logic.SOr (toTrue i1, toTrue i2) + val fmFalse = Prop_Logic.SAnd (toFalse i1, toFalse i2) in SOME (Leaf [fmTrue, fmFalse], m2, a2) end @@ -1815,8 +1815,8 @@ let val (i1, m1, a1) = interpret ctxt model args t1 val (i2, m2, a2) = interpret ctxt m1 a1 t2 - val fmTrue = PropLogic.SOr (toFalse i1, toTrue i2) - val fmFalse = PropLogic.SAnd (toTrue i1, toFalse i2) + val fmTrue = Prop_Logic.SOr (toFalse i1, toTrue i2) + val fmFalse = Prop_Logic.SAnd (toTrue i1, toFalse i2) in SOME (Leaf [fmTrue, fmFalse], m2, a2) end @@ -1890,7 +1890,7 @@ val intr = Leaf fms (* prop_formula list -> prop_formula *) fun one_of_two_false [] = True - | one_of_two_false (x::xs) = SAnd (PropLogic.all (map (fn x' => + | one_of_two_false (x::xs) = SAnd (Prop_Logic.all (map (fn x' => SOr (SNot x, SNot x')) xs), one_of_two_false xs) (* prop_formula *) val wf = one_of_two_false fms @@ -2961,7 +2961,7 @@ strip_leading_quote x ^ string_of_int i (* interpretation -> int *) fun index_from_interpretation (Leaf xs) = - find_index (PropLogic.eval assignment) xs + find_index (Prop_Logic.eval assignment) xs | index_from_interpretation _ = raise REFUTE ("stlc_printer", "interpretation for ground type is not a leaf") @@ -3045,7 +3045,7 @@ (* the index of the element in the datatype *) val element = (case intr of - Leaf xs => find_index (PropLogic.eval assignment) xs + Leaf xs => find_index (Prop_Logic.eval assignment) xs | Node _ => raise REFUTE ("IDT_printer", "interpretation is not a leaf")) in diff -r 655f583840d0 -r 9908cf4af394 src/HOL/Tools/sat_funcs.ML --- a/src/HOL/Tools/sat_funcs.ML Sat Jan 08 11:18:09 2011 -0800 +++ b/src/HOL/Tools/sat_funcs.ML Sat Jan 08 11:29:25 2011 -0800 @@ -266,13 +266,13 @@ (* a 'prop_formula' (just for tracing) *) (* ------------------------------------------------------------------------- *) -fun string_of_prop_formula PropLogic.True = "True" - | string_of_prop_formula PropLogic.False = "False" - | string_of_prop_formula (PropLogic.BoolVar i) = "x" ^ string_of_int i - | string_of_prop_formula (PropLogic.Not fm) = "~" ^ string_of_prop_formula fm - | string_of_prop_formula (PropLogic.Or (fm1, fm2)) = +fun string_of_prop_formula Prop_Logic.True = "True" + | string_of_prop_formula Prop_Logic.False = "False" + | string_of_prop_formula (Prop_Logic.BoolVar i) = "x" ^ string_of_int i + | string_of_prop_formula (Prop_Logic.Not fm) = "~" ^ string_of_prop_formula fm + | string_of_prop_formula (Prop_Logic.Or (fm1, fm2)) = "(" ^ string_of_prop_formula fm1 ^ " v " ^ string_of_prop_formula fm2 ^ ")" - | string_of_prop_formula (PropLogic.And (fm1, fm2)) = + | string_of_prop_formula (Prop_Logic.And (fm1, fm2)) = "(" ^ string_of_prop_formula fm1 ^ " & " ^ string_of_prop_formula fm2 ^ ")"; (* ------------------------------------------------------------------------- *) @@ -313,15 +313,15 @@ tracing ("Sorted non-trivial clauses:\n" ^ cat_lines (map (Syntax.string_of_term ctxt o Thm.term_of) sorted_clauses)) else () - (* translate clauses from HOL terms to PropLogic.prop_formula *) + (* translate clauses from HOL terms to Prop_Logic.prop_formula *) val (fms, atom_table) = - fold_map (PropLogic.prop_formula_of_term o HOLogic.dest_Trueprop o Thm.term_of) + fold_map (Prop_Logic.prop_formula_of_term o HOLogic.dest_Trueprop o Thm.term_of) sorted_clauses Termtab.empty val _ = if !trace_sat then tracing ("Invoking SAT solver on clauses:\n" ^ cat_lines (map string_of_prop_formula fms)) else () - val fm = PropLogic.all fms + val fm = Prop_Logic.all fms (* unit -> Thm.thm *) fun make_quick_and_dirty_thm () = let diff -r 655f583840d0 -r 9908cf4af394 src/HOL/Tools/sat_solver.ML --- a/src/HOL/Tools/sat_solver.ML Sat Jan 08 11:18:09 2011 -0800 +++ b/src/HOL/Tools/sat_solver.ML Sat Jan 08 11:29:25 2011 -0800 @@ -14,16 +14,16 @@ datatype result = SATISFIABLE of assignment | UNSATISFIABLE of proof option | UNKNOWN - type solver = PropLogic.prop_formula -> result + type solver = Prop_Logic.prop_formula -> result (* auxiliary functions to create external SAT solvers *) - val write_dimacs_cnf_file : Path.T -> PropLogic.prop_formula -> unit - val write_dimacs_sat_file : Path.T -> PropLogic.prop_formula -> unit - val read_std_result_file : Path.T -> string * string * string -> result - val make_external_solver : string -> (PropLogic.prop_formula -> unit) -> - (unit -> result) -> solver + val write_dimacs_cnf_file : Path.T -> Prop_Logic.prop_formula -> unit + val write_dimacs_sat_file : Path.T -> Prop_Logic.prop_formula -> unit + val read_std_result_file : Path.T -> string * string * string -> result + val make_external_solver : string -> (Prop_Logic.prop_formula -> unit) -> + (unit -> result) -> solver - val read_dimacs_cnf_file : Path.T -> PropLogic.prop_formula + val read_dimacs_cnf_file : Path.T -> Prop_Logic.prop_formula (* generic solver interface *) val solvers : (string * solver) list Unsynchronized.ref @@ -34,7 +34,7 @@ structure SatSolver : SAT_SOLVER = struct - open PropLogic; + open Prop_Logic; (* ------------------------------------------------------------------------- *) (* should be raised by an external SAT solver to indicate that the solver is *) @@ -279,8 +279,6 @@ (* make_external_solver: call 'writefn', execute 'cmd', call 'readfn' *) (* ------------------------------------------------------------------------- *) - (* string -> (PropLogic.prop_formula -> unit) -> (unit -> result) -> solver *) - fun make_external_solver cmd writefn readfn fm = (writefn fm; bash cmd; readfn ()); @@ -289,8 +287,6 @@ (* a SAT problem given in DIMACS CNF format *) (* ------------------------------------------------------------------------- *) - (* Path.T -> PropLogic.prop_formula *) - fun read_dimacs_cnf_file path = let (* string list -> string list *) @@ -323,27 +319,24 @@ | (0::tl) => xs1 :: clauses tl | _ => raise Fail "SatSolver.clauses" end - (* int -> PropLogic.prop_formula *) fun literal_from_int i = (i<>0 orelse error "variable index in DIMACS CNF file is 0"; if i>0 then - PropLogic.BoolVar i + Prop_Logic.BoolVar i else - PropLogic.Not (PropLogic.BoolVar (~i))) - (* PropLogic.prop_formula list -> PropLogic.prop_formula *) + Prop_Logic.Not (Prop_Logic.BoolVar (~i))) fun disjunction [] = error "empty clause in DIMACS CNF file" | disjunction (x::xs) = (case xs of [] => x - | _ => PropLogic.Or (x, disjunction xs)) - (* PropLogic.prop_formula list -> PropLogic.prop_formula *) + | _ => Prop_Logic.Or (x, disjunction xs)) fun conjunction [] = error "no clause in DIMACS CNF file" | conjunction (x::xs) = (case xs of [] => x - | _ => PropLogic.And (x, conjunction xs)) + | _ => Prop_Logic.And (x, conjunction xs)) in (conjunction o (map disjunction) @@ -405,7 +398,7 @@ fun enum_solver fm = let (* int list *) - val indices = PropLogic.indices fm + val indices = Prop_Logic.indices fm (* int list -> int list -> int list option *) (* binary increment: list 'xs' of current bits, list 'ys' of all bits (lower bits first) *) fun next_list _ ([]:int list) = @@ -424,7 +417,7 @@ member (op =) xs i (* int list -> SatSolver.result *) fun solver_loop xs = - if PropLogic.eval (assignment_from_list xs) fm then + if Prop_Logic.eval (assignment_from_list xs) fm then SatSolver.SATISFIABLE (SOME o (assignment_from_list xs)) else (case next_list xs indices of @@ -446,16 +439,16 @@ let local - open PropLogic + open Prop_Logic in fun dpll_solver fm = let - (* We could use 'PropLogic.defcnf fm' instead of 'PropLogic.nnf fm' *) + (* We could use 'Prop_Logic.defcnf fm' instead of 'Prop_Logic.nnf fm' *) (* but that sometimes leads to worse performance due to the *) (* introduction of additional variables. *) - val fm' = PropLogic.nnf fm + val fm' = Prop_Logic.nnf fm (* int list *) - val indices = PropLogic.indices fm' + val indices = Prop_Logic.indices fm' (* int list -> int -> prop_formula *) fun partial_var_eval [] i = BoolVar i | partial_var_eval (x::xs) i = if x=i then True else if x=(~i) then False else partial_var_eval xs i @@ -587,7 +580,7 @@ fun readfn () = SatSolver.read_std_result_file outpath ("SAT", "", "UNSAT") val _ = if File.exists inpath then warning ("overwriting existing file " ^ quote (Path.implode inpath)) else () val _ = if File.exists outpath then warning ("overwriting existing file " ^ quote (Path.implode outpath)) else () - val cnf = PropLogic.defcnf fm + val cnf = Prop_Logic.defcnf fm val result = SatSolver.make_external_solver cmd writefn readfn cnf val _ = try File.rm inpath val _ = try File.rm outpath @@ -599,16 +592,16 @@ handle IO.Io _ => raise INVALID_PROOF "Could not read file \"result.prf\"" (* representation of clauses as ordered lists of literals (with duplicates removed) *) (* prop_formula -> int list *) - fun clause_to_lit_list (PropLogic.Or (fm1, fm2)) = + fun clause_to_lit_list (Prop_Logic.Or (fm1, fm2)) = Ord_List.union int_ord (clause_to_lit_list fm1) (clause_to_lit_list fm2) - | clause_to_lit_list (PropLogic.BoolVar i) = + | clause_to_lit_list (Prop_Logic.BoolVar i) = [i] - | clause_to_lit_list (PropLogic.Not (PropLogic.BoolVar i)) = + | clause_to_lit_list (Prop_Logic.Not (Prop_Logic.BoolVar i)) = [~i] | clause_to_lit_list _ = raise INVALID_PROOF "Error: invalid clause in CNF formula." (* prop_formula -> int *) - fun cnf_number_of_clauses (PropLogic.And (fm1, fm2)) = + fun cnf_number_of_clauses (Prop_Logic.And (fm1, fm2)) = cnf_number_of_clauses fm1 + cnf_number_of_clauses fm2 | cnf_number_of_clauses _ = 1 @@ -617,7 +610,7 @@ val clauses = Array.array (number_of_clauses, []) (* initialize the 'clauses' array *) (* prop_formula * int -> int *) - fun init_array (PropLogic.And (fm1, fm2), n) = + fun init_array (Prop_Logic.And (fm1, fm2), n) = init_array (fm2, init_array (fm1, n)) | init_array (fm, n) = (Array.update (clauses, n, clause_to_lit_list fm); n+1) @@ -768,7 +761,7 @@ val inpath = File.tmp_path (Path.explode ("isabelle" ^ serial_str ^ ".cnf")) val outpath = File.tmp_path (Path.explode ("result" ^ serial_str)) val cmd = getenv "MINISAT_HOME" ^ "/minisat " ^ File.shell_path inpath ^ " -r " ^ File.shell_path outpath ^ " > /dev/null" - fun writefn fm = SatSolver.write_dimacs_cnf_file inpath (PropLogic.defcnf fm) + fun writefn fm = SatSolver.write_dimacs_cnf_file inpath (Prop_Logic.defcnf fm) fun readfn () = SatSolver.read_std_result_file outpath ("SAT", "", "UNSAT") val _ = if File.exists inpath then warning ("overwriting existing file " ^ quote (Path.implode inpath)) else () val _ = if File.exists outpath then warning ("overwriting existing file " ^ quote (Path.implode outpath)) else () @@ -806,9 +799,9 @@ (* string list *) val proof_lines = ((split_lines o File.read) (Path.explode "resolve_trace")) handle IO.Io _ => raise INVALID_PROOF "Could not read file \"resolve_trace\"" - (* PropLogic.prop_formula -> int *) - fun cnf_number_of_clauses (PropLogic.And (fm1, fm2)) = cnf_number_of_clauses fm1 + cnf_number_of_clauses fm2 - | cnf_number_of_clauses _ = 1 + fun cnf_number_of_clauses (Prop_Logic.And (fm1, fm2)) = + cnf_number_of_clauses fm1 + cnf_number_of_clauses fm2 + | cnf_number_of_clauses _ = 1 (* string -> int *) fun int_from_string s = ( case Int.fromString s of @@ -848,7 +841,7 @@ val _ = if !clause_offset = ~1 then clause_offset := (case Inttab.max_key (!clause_table) of SOME id => id - | NONE => cnf_number_of_clauses (PropLogic.defcnf fm) - 1 (* the first clause ID is 0, not 1 *)) + | NONE => cnf_number_of_clauses (Prop_Logic.defcnf fm) - 1 (* the first clause ID is 0, not 1 *)) else () val vid = int_from_string id @@ -927,7 +920,7 @@ val inpath = File.tmp_path (Path.explode ("isabelle" ^ serial_str ^ ".cnf")) val outpath = File.tmp_path (Path.explode ("result" ^ serial_str)) val cmd = getenv "ZCHAFF_HOME" ^ "/zchaff " ^ File.shell_path inpath ^ " > " ^ File.shell_path outpath - fun writefn fm = SatSolver.write_dimacs_cnf_file inpath (PropLogic.defcnf fm) + fun writefn fm = SatSolver.write_dimacs_cnf_file inpath (Prop_Logic.defcnf fm) fun readfn () = SatSolver.read_std_result_file outpath ("Instance Satisfiable", "", "Instance Unsatisfiable") val _ = if File.exists inpath then warning ("overwriting existing file " ^ quote (Path.implode inpath)) else () val _ = if File.exists outpath then warning ("overwriting existing file " ^ quote (Path.implode outpath)) else () @@ -954,7 +947,7 @@ val outpath = File.tmp_path (Path.explode ("result" ^ serial_str)) val exec = getenv "BERKMIN_EXE" val cmd = getenv "BERKMIN_HOME" ^ "/" ^ (if exec = "" then "BerkMin561" else exec) ^ " " ^ File.shell_path inpath ^ " > " ^ File.shell_path outpath - fun writefn fm = SatSolver.write_dimacs_cnf_file inpath (PropLogic.defcnf fm) + fun writefn fm = SatSolver.write_dimacs_cnf_file inpath (Prop_Logic.defcnf fm) fun readfn () = SatSolver.read_std_result_file outpath ("Satisfiable !!", "solution =", "UNSATISFIABLE !!") val _ = if File.exists inpath then warning ("overwriting existing file " ^ quote (Path.implode inpath)) else () val _ = if File.exists outpath then warning ("overwriting existing file " ^ quote (Path.implode outpath)) else () @@ -980,7 +973,7 @@ val inpath = File.tmp_path (Path.explode ("isabelle" ^ serial_str ^ ".cnf")) val outpath = File.tmp_path (Path.explode ("result" ^ serial_str)) val cmd = getenv "JERUSAT_HOME" ^ "/Jerusat1.3 " ^ File.shell_path inpath ^ " > " ^ File.shell_path outpath - fun writefn fm = SatSolver.write_dimacs_cnf_file inpath (PropLogic.defcnf fm) + fun writefn fm = SatSolver.write_dimacs_cnf_file inpath (Prop_Logic.defcnf fm) fun readfn () = SatSolver.read_std_result_file outpath ("s SATISFIABLE", "v ", "s UNSATISFIABLE") val _ = if File.exists inpath then warning ("overwriting existing file " ^ quote (Path.implode inpath)) else () val _ = if File.exists outpath then warning ("overwriting existing file " ^ quote (Path.implode outpath)) else () diff -r 655f583840d0 -r 9908cf4af394 src/HOL/Tools/smallvalue_generators.ML --- a/src/HOL/Tools/smallvalue_generators.ML Sat Jan 08 11:18:09 2011 -0800 +++ b/src/HOL/Tools/smallvalue_generators.ML Sat Jan 08 11:29:25 2011 -0800 @@ -210,8 +210,10 @@ (** building and compiling generator expressions **) -structure Counterexample = Proof_Data ( +structure Counterexample = Proof_Data +( type T = unit -> int -> term list option + (* FIXME avoid user error with non-user text *) fun init _ () = error "Counterexample" ); val put_counterexample = Counterexample.put; diff -r 655f583840d0 -r 9908cf4af394 src/HOL/Tools/type_lifting.ML --- a/src/HOL/Tools/type_lifting.ML Sat Jan 08 11:18:09 2011 -0800 +++ b/src/HOL/Tools/type_lifting.ML Sat Jan 08 11:29:25 2011 -0800 @@ -27,11 +27,12 @@ type entry = { mapper: term, variances: (sort * (bool * bool)) list, comp: thm, id: thm }; -structure Data = Generic_Data( +structure Data = Generic_Data +( type T = entry list Symtab.table val empty = Symtab.empty - fun merge (xy : T * T) = Symtab.merge (K true) xy val extend = I + fun merge data = Symtab.merge (K true) data ); val entries = Data.get o Context.Proof; diff -r 655f583840d0 -r 9908cf4af394 src/HOL/ex/SAT_Examples.thy --- a/src/HOL/ex/SAT_Examples.thy Sat Jan 08 11:18:09 2011 -0800 +++ b/src/HOL/ex/SAT_Examples.thy Sat Jan 08 11:29:25 2011 -0800 @@ -533,12 +533,10 @@ fun benchmark dimacsfile = let val prop_fm = SatSolver.read_dimacs_cnf_file (Path.explode dimacsfile) - fun and_to_list (PropLogic.And (fm1, fm2)) acc = - and_to_list fm2 (fm1 :: acc) - | and_to_list fm acc = - rev (fm :: acc) + fun and_to_list (Prop_Logic.And (fm1, fm2)) acc = and_to_list fm2 (fm1 :: acc) + | and_to_list fm acc = rev (fm :: acc) val clauses = and_to_list prop_fm [] - val terms = map (HOLogic.mk_Trueprop o PropLogic.term_of_prop_formula) + val terms = map (HOLogic.mk_Trueprop o Prop_Logic.term_of_prop_formula) clauses val cterms = map (Thm.cterm_of @{theory}) terms val timer = start_timing () diff -r 655f583840d0 -r 9908cf4af394 src/Pure/General/ord_list.ML --- a/src/Pure/General/ord_list.ML Sat Jan 08 11:18:09 2011 -0800 +++ b/src/Pure/General/ord_list.ML Sat Jan 08 11:29:25 2011 -0800 @@ -14,6 +14,7 @@ val remove: ('b * 'a -> order) -> 'b -> 'a T -> 'a T val subset: ('b * 'a -> order) -> 'b T * 'a T -> bool val union: ('a * 'a -> order) -> 'a T -> 'a T -> 'a T + val merge: ('a * 'a -> order) -> 'a T * 'a T -> 'a T val inter: ('b * 'a -> order) -> 'b T -> 'a T -> 'a T val subtract: ('b * 'a -> order) -> 'b T -> 'a T -> 'a T end; @@ -85,6 +86,8 @@ | GREATER => y :: unio lst1 ys); in if pointer_eq (list1, list2) then list1 else handle_same (unio list1) list2 end; +fun merge ord (list1, list2) = union ord list2 list1; + (*intersection: filter second list for elements present in first list*) fun inter ord list1 list2 = let diff -r 655f583840d0 -r 9908cf4af394 src/Pure/Isar/code.ML --- a/src/Pure/Isar/code.ML Sat Jan 08 11:18:09 2011 -0800 +++ b/src/Pure/Isar/code.ML Sat Jan 08 11:29:25 2011 -0800 @@ -263,7 +263,7 @@ type T = spec * (data * theory_ref) option Synchronized.var; val empty = (make_spec (false, (((Symtab.empty, Symtab.empty), Symtab.empty), (Symtab.empty, (Symtab.empty, Symtab.empty)))), empty_dataref ()); - val extend = I + val extend = I (* FIXME empty_dataref!?! *) fun merge ((spec1, _), (spec2, _)) = (merge_spec (spec1, spec2), empty_dataref ()); ); @@ -1238,7 +1238,7 @@ (** infrastructure **) -(* c.f. src/HOL/Tools/recfun_codegen.ML *) +(* cf. src/HOL/Tools/recfun_codegen.ML *) structure Code_Target_Attr = Theory_Data ( diff -r 655f583840d0 -r 9908cf4af394 src/Pure/Isar/spec_rules.ML --- a/src/Pure/Isar/spec_rules.ML Sat Jan 08 11:18:09 2011 -0800 +++ b/src/Pure/Isar/spec_rules.ML Sat Jan 08 11:29:25 2011 -0800 @@ -37,7 +37,7 @@ eq_list Thm.eq_thm_prop (ths1, ths2)) (#1 o #2); val extend = I; - fun merge data = Item_Net.merge data; + val merge = Item_Net.merge; ); val get = Item_Net.content o Rules.get o Context.Proof; diff -r 655f583840d0 -r 9908cf4af394 src/Pure/Isar/specification.ML --- a/src/Pure/Isar/specification.ML Sat Jan 08 11:18:09 2011 -0800 +++ b/src/Pure/Isar/specification.ML Sat Jan 08 11:29:25 2011 -0800 @@ -366,7 +366,7 @@ type T = ((bool -> Proof.state -> Proof.state) * stamp) list; val empty = []; val extend = I; - fun merge hooks : T = Library.merge (eq_snd op =) hooks; + fun merge data : T = Library.merge (eq_snd op =) data; ); fun gen_theorem schematic prep_att prep_stmt diff -r 655f583840d0 -r 9908cf4af394 src/Pure/context_position.ML --- a/src/Pure/context_position.ML Sat Jan 08 11:18:09 2011 -0800 +++ b/src/Pure/context_position.ML Sat Jan 08 11:29:25 2011 -0800 @@ -10,6 +10,7 @@ val set_visible: bool -> Proof.context -> Proof.context val restore_visible: Proof.context -> Proof.context -> Proof.context val if_visible: Proof.context -> ('a -> unit) -> 'a -> unit + val if_visible_proof: Context.generic -> ('a -> unit) -> 'a -> unit val reported_text: Proof.context -> Position.T -> Markup.T -> string -> string val report_text: Proof.context -> Position.T -> Markup.T -> string -> unit val report: Proof.context -> Position.T -> Markup.T -> unit @@ -25,6 +26,9 @@ fun if_visible ctxt f x = if is_visible ctxt then f x else (); +fun if_visible_proof (Context.Proof ctxt) f x = if_visible ctxt f x + | if_visible_proof _ _ _ = (); + fun reported_text ctxt pos markup txt = if is_visible ctxt then Position.reported_text pos markup txt else ""; diff -r 655f583840d0 -r 9908cf4af394 src/Pure/raw_simplifier.ML --- a/src/Pure/raw_simplifier.ML Sat Jan 08 11:18:09 2011 -0800 +++ b/src/Pure/raw_simplifier.ML Sat Jan 08 11:29:25 2011 -0800 @@ -295,7 +295,7 @@ fun if_visible (Simpset ({context, ...}, _)) f x = (case context of - SOME ctxt => if Context_Position.is_visible ctxt then f x else () + SOME ctxt => Context_Position.if_visible ctxt f x | NONE => ()); local diff -r 655f583840d0 -r 9908cf4af394 src/Tools/Code/code_runtime.ML --- a/src/Tools/Code/code_runtime.ML Sat Jan 08 11:18:09 2011 -0800 +++ b/src/Tools/Code/code_runtime.ML Sat Jan 08 11:29:25 2011 -0800 @@ -138,6 +138,7 @@ structure Truth_Result = Proof_Data ( type T = unit -> truth + (* FIXME avoid user error with non-user text *) fun init _ () = error "Truth_Result" ); val put_truth = Truth_Result.put; @@ -369,8 +370,8 @@ ( type T = string list val empty = [] + val extend = I fun merge data : T = Library.merge (op =) data - val extend = I ); fun notify_val (string, value) = diff -r 655f583840d0 -r 9908cf4af394 src/Tools/adhoc_overloading.ML --- a/src/Tools/adhoc_overloading.ML Sat Jan 08 11:18:09 2011 -0800 +++ b/src/Tools/adhoc_overloading.ML Sat Jan 08 11:29:25 2011 -0800 @@ -53,10 +53,11 @@ if ext_name1 = ext_name2 then ext_name1 else duplicate_variant_err int_name ext_name1; - fun merge ({internalize=int1, externalize=ext1}, - {internalize=int2, externalize=ext2}) = - {internalize=Symtab.join (K (Library.merge (op =))) (int1, int2), - externalize=Symtab.join merge_ext (ext1, ext2)}; + fun merge + ({internalize = int1, externalize = ext1}, + {internalize = int2, externalize = ext2}) : T = + {internalize = Symtab.join (K (Library.merge (op =))) (int1, int2), + externalize = Symtab.join merge_ext (ext1, ext2)}; ); fun map_tables f g = diff -r 655f583840d0 -r 9908cf4af394 src/Tools/nbe.ML --- a/src/Tools/nbe.ML Sat Jan 08 11:18:09 2011 -0800 +++ b/src/Tools/nbe.ML Sat Jan 08 11:29:25 2011 -0800 @@ -229,8 +229,10 @@ (* nbe specific syntax and sandbox communication *) -structure Univs = Proof_Data ( +structure Univs = Proof_Data +( type T = unit -> Univ list -> Univ list + (* FIXME avoid user error with non-user text *) fun init _ () = error "Univs" ); val put_result = Univs.put; diff -r 655f583840d0 -r 9908cf4af394 src/Tools/quickcheck.ML --- a/src/Tools/quickcheck.ML Sat Jan 08 11:18:09 2011 -0800 +++ b/src/Tools/quickcheck.ML Sat Jan 08 11:29:25 2011 -0800 @@ -93,9 +93,10 @@ fun map_test_params' f (Test_Params {default_type, expect}) = make_test_params (f (default_type, expect)); fun merge_test_params - (Test_Params {default_type = default_type1, expect = expect1}, - Test_Params {default_type = default_type2, expect = expect2}) = - make_test_params (merge (op =) (default_type1, default_type2), merge_expectation (expect1, expect2)); + (Test_Params {default_type = default_type1, expect = expect1}, + Test_Params {default_type = default_type2, expect = expect2}) = + make_test_params + (merge (op =) (default_type1, default_type2), merge_expectation (expect1, expect2)); structure Data = Generic_Data (