--- a/etc/components Sat Jan 08 17:30:05 2011 +0100
+++ b/etc/components Sat Jan 08 17:39:51 2011 +0100
@@ -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
--- a/src/HOL/IsaMakefile Sat Jan 08 17:30:05 2011 +0100
+++ b/src/HOL/IsaMakefile Sat Jan 08 17:39:51 2011 +0100
@@ -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 \
--- a/src/HOL/Library/Library.thy Sat Jan 08 17:30:05 2011 +0100
+++ b/src/HOL/Library/Library.thy Sat Jan 08 17:39:51 2011 +0100
@@ -56,7 +56,7 @@
Set_Algebras
SML_Quickcheck
State_Monad
- Sum_Of_Squares
+ Sum_of_Squares
Transitive_Closure_Table
Univ_Poly
While_Combinator
--- a/src/HOL/Library/Sum_Of_Squares.thy Sat Jan 08 17:30:05 2011 +0100
+++ /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 \<Longrightarrow> a < 0"
-by (sos_cert "((R<1 + (((A<1 * R<1) * (R<2 * [1]^2)) + (((A<0 * R<1) * (R<3 * [1]^2)) + ((A<=0 * R<1) * (R<14 * [1]^2))))))")
-
-lemma "a1 >= 0 & a2 >= 0 \<and> (a1 * a1 + a2 * a2 = b1 * b1 + b2 * b2 + 2) \<and> (a1 * b1 + a2 * b2 = 0) --> a1 * a2 - b1 * b2 >= (0::real)"
-by (sos_cert "(((A<0 * R<1) + (([~1/2*a1*b2 + ~1/2*a2*b1] * A=0) + (([~1/2*a1*a2 + 1/2*b1*b2] * A=1) + (((A<0 * R<1) * ((R<1/2 * [b2]^2) + (R<1/2 * [b1]^2))) + ((A<=0 * (A<=1 * R<1)) * ((R<1/2 * [b2]^2) + ((R<1/2 * [b1]^2) + ((R<1/2 * [a2]^2) + (R<1/2 * [a1]^2))))))))))")
-
-lemma "(3::real) * x + 7 * a < 4 & 3 < 2 * x --> 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 \<longrightarrow> c * a^2 * b <= x"
-by (sos_cert "(((A<0 * R<1) + (((A<=3 * R<1) * (R<1 * [1]^2)) + (((A<=1 * (A<=2 * R<1)) * (R<1/27 * [~1*a + b]^2)) + ((A<=0 * (A<=2 * R<1)) * (R<8/27 * [~1*a + b]^2))))))")
-
-lemma "(0::real) < x --> 0 < 1 + x + x^2"
-by (sos_cert "((R<1 + ((R<1 * (R<1 * [x]^2)) + (((A<0 * R<1) * (R<1 * [1]^2)) + ((A<=0 * R<1) * (R<1 * [1]^2))))))")
-
-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 \<and> 3 < 2 * x \<longrightarrow> a < 0"
-by (sos_cert "((R<1 + (((A<1 * R<1) * (R<2 * [1]^2)) + (((A<0 * R<1) * (R<3 * [1]^2)) + ((A<=0 * R<1) * (R<14 * [1]^2))))))")
-
-lemma "(0::real) < 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
-
--- a/src/HOL/Library/Sum_Of_Squares/etc/settings Sat Jan 08 17:30:05 2011 +0100
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,3 +0,0 @@
-# -*- shell-script -*- :mode=shellscript:
-
-ISABELLE_SUM_OF_SQUARES="$COMPONENT"
--- a/src/HOL/Library/Sum_Of_Squares/neos_csdp_client Sat Jan 08 17:30:05 2011 +0100
+++ /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 <input_filename> <output_filename>\n")
- sys.exit(19)
-
-xml_pre = "<document>\n<category>sdp</category>\n<solver>csdp</solver>\n<inputMethod>SPARSE_SDPA</inputMethod>\n<dat><![CDATA["
-xml_post = "]]></dat>\n</document>\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)
-
--- a/src/HOL/Library/Sum_Of_Squares/positivstellensatz_tools.ML Sat Jan 08 17:30:05 2011 +0100
+++ /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
-
-
--- a/src/HOL/Library/Sum_Of_Squares/sos_wrapper.ML Sat Jan 08 17:30:05 2011 +0100
+++ /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
--- a/src/HOL/Library/Sum_Of_Squares/sum_of_squares.ML Sat Jan 08 17:30:05 2011 +0100
+++ /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 </ rat_0 then (Rat.neg o Rat.rat_of_int) else Rat.rat_of_int
- val x2 = 2 * (a - (b * d))
- in s (if x2 >= 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 // rat_10) then normalize (rat_10 */ y) - 1
- else 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 </ rat_0 then "-0." else "0.") ^
- implode(tl(raw_explode(string_of_int k))) ^
- (if e = 0 then "" else "e"^string_of_int e)
- end
-end;
-
-(* Iterations over numbers, and lists indexed by numbers. *)
-
-fun itern k l f a =
- case l of
- [] => 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 </ rat_0 then a ^ " - " ^ string_of_cmonomial(Rat.neg c,m)
- else a ^ " + " ^ string_of_cmonomial(c,m))
- cms ""
- val s1 = String.substring (s, 0, 3)
- val s2 = String.substring (s, 3, String.size s - 3)
- in "<<" ^(if s1 = " + " then s2 else "-"^s2)^">>"
- 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 </ y then y else x
- fun common_denominator fld amat acc =
- fld (fn (m,c) => 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 </ y then y else x
- fun common_denominator fld amat acc =
- fld (fn (m,c) => 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 </ rat_0 then raise Failure "diagonalize: not PSD"
- else if a11 =/ rat_0 then
- if FuncUtil.Intfunc.is_empty (snd (row i m)) then diagonalize n (i + 1) m
- else raise Failure "diagonalize: not PSD ___ "
- else
- let
- val v = row i m
- val v' = (fst v, FuncUtil.Intfunc.fold (fn (i, c) => 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 </ Rat.zero then nth les n
- else raise Failure "trivial_axiom: Not a trivial axiom"
- | RealArith.Axiom_lt n => 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;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/src/HOL/Library/Sum_of_Squares.thy Sat Jan 08 17:39:51 2011 +0100
@@ -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 \<Longrightarrow> a < 0"
+by (sos_cert "((R<1 + (((A<1 * R<1) * (R<2 * [1]^2)) + (((A<0 * R<1) * (R<3 * [1]^2)) + ((A<=0 * R<1) * (R<14 * [1]^2))))))")
+
+lemma "a1 >= 0 & a2 >= 0 \<and> (a1 * a1 + a2 * a2 = b1 * b1 + b2 * b2 + 2) \<and> (a1 * b1 + a2 * b2 = 0) --> a1 * a2 - b1 * b2 >= (0::real)"
+by (sos_cert "(((A<0 * R<1) + (([~1/2*a1*b2 + ~1/2*a2*b1] * A=0) + (([~1/2*a1*a2 + 1/2*b1*b2] * A=1) + (((A<0 * R<1) * ((R<1/2 * [b2]^2) + (R<1/2 * [b1]^2))) + ((A<=0 * (A<=1 * R<1)) * ((R<1/2 * [b2]^2) + ((R<1/2 * [b1]^2) + ((R<1/2 * [a2]^2) + (R<1/2 * [a1]^2))))))))))")
+
+lemma "(3::real) * x + 7 * a < 4 & 3 < 2 * x --> 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 \<longrightarrow> c * a^2 * b <= x"
+by (sos_cert "(((A<0 * R<1) + (((A<=3 * R<1) * (R<1 * [1]^2)) + (((A<=1 * (A<=2 * R<1)) * (R<1/27 * [~1*a + b]^2)) + ((A<=0 * (A<=2 * R<1)) * (R<8/27 * [~1*a + b]^2))))))")
+
+lemma "(0::real) < x --> 0 < 1 + x + x^2"
+by (sos_cert "((R<1 + ((R<1 * (R<1 * [x]^2)) + (((A<0 * R<1) * (R<1 * [1]^2)) + ((A<=0 * R<1) * (R<1 * [1]^2))))))")
+
+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 \<and> 3 < 2 * x \<longrightarrow> a < 0"
+by (sos_cert "((R<1 + (((A<1 * R<1) * (R<2 * [1]^2)) + (((A<0 * R<1) * (R<3 * [1]^2)) + ((A<=0 * R<1) * (R<14 * [1]^2))))))")
+
+lemma "(0::real) < 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
+
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/src/HOL/Library/Sum_of_Squares/etc/settings Sat Jan 08 17:39:51 2011 +0100
@@ -0,0 +1,3 @@
+# -*- shell-script -*- :mode=shellscript:
+
+ISABELLE_SUM_OF_SQUARES="$COMPONENT"
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/src/HOL/Library/Sum_of_Squares/neos_csdp_client Sat Jan 08 17:39:51 2011 +0100
@@ -0,0 +1,87 @@
+#!/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 <input_filename> <output_filename>\n")
+ sys.exit(19)
+
+xml_pre = "<document>\n<category>sdp</category>\n<solver>csdp</solver>\n<inputMethod>SPARSE_SDPA</inputMethod>\n<dat><![CDATA["
+xml_post = "]]></dat>\n</document>\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)
+
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/src/HOL/Library/Sum_of_Squares/positivstellensatz_tools.ML Sat Jan 08 17:39:51 2011 +0100
@@ -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
+
+
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/src/HOL/Library/Sum_of_Squares/sos_wrapper.ML Sat Jan 08 17:39:51 2011 +0100
@@ -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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/src/HOL/Library/Sum_of_Squares/sum_of_squares.ML Sat Jan 08 17:39:51 2011 +0100
@@ -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 </ rat_0 then (Rat.neg o Rat.rat_of_int) else Rat.rat_of_int
+ val x2 = 2 * (a - (b * d))
+ in s (if x2 >= 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 // rat_10) then normalize (rat_10 */ y) - 1
+ else 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 </ rat_0 then "-0." else "0.") ^
+ implode(tl(raw_explode(string_of_int k))) ^
+ (if e = 0 then "" else "e"^string_of_int e)
+ end
+end;
+
+(* Iterations over numbers, and lists indexed by numbers. *)
+
+fun itern k l f a =
+ case l of
+ [] => 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 </ rat_0 then a ^ " - " ^ string_of_cmonomial(Rat.neg c,m)
+ else a ^ " + " ^ string_of_cmonomial(c,m))
+ cms ""
+ val s1 = String.substring (s, 0, 3)
+ val s2 = String.substring (s, 3, String.size s - 3)
+ in "<<" ^(if s1 = " + " then s2 else "-"^s2)^">>"
+ 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 </ y then y else x
+ fun common_denominator fld amat acc =
+ fld (fn (m,c) => 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 </ y then y else x
+ fun common_denominator fld amat acc =
+ fld (fn (m,c) => 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 </ rat_0 then raise Failure "diagonalize: not PSD"
+ else if a11 =/ rat_0 then
+ if FuncUtil.Intfunc.is_empty (snd (row i m)) then diagonalize n (i + 1) m
+ else raise Failure "diagonalize: not PSD ___ "
+ else
+ let
+ val v = row i m
+ val v' = (fst v, FuncUtil.Intfunc.fold (fn (i, c) => 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 </ Rat.zero then nth les n
+ else raise Failure "trivial_axiom: Not a trivial axiom"
+ | RealArith.Axiom_lt n => 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;