misc changes to SOS by Philipp Meyer:
authorwenzelm
Thu Aug 06 19:51:59 2009 +0200 (2009-08-06)
changeset 32332bc5cec7b2be6
parent 32331 e60684ecaf3d
child 32333 d4cb904cc63c
misc changes to SOS by Philipp Meyer:
CSDP_EXE as central setting;
separate component src/HOL/Library/Sum_Of_Squares;
misc tuning and rearrangement of neos_csdp_client;
more robust treatment of shell paths;
debugging depends on local flag;
removed unused parts;
etc/components
etc/settings
lib/scripts/neos/NeosCSDPClient.py
lib/scripts/neos/config.py
src/HOL/IsaMakefile
src/HOL/Library/Sum_Of_Squares.thy
src/HOL/Library/Sum_Of_Squares/etc/settings
src/HOL/Library/Sum_Of_Squares/neos_csdp_client
src/HOL/Library/Sum_Of_Squares/sos_wrapper.ML
src/HOL/Library/Sum_Of_Squares/sum_of_squares.ML
src/HOL/Library/sos_wrapper.ML
src/HOL/Library/sum_of_squares.ML
     1.1 --- a/etc/components	Wed Aug 05 17:10:10 2009 +0200
     1.2 +++ b/etc/components	Thu Aug 06 19:51:59 2009 +0200
     1.3 @@ -12,4 +12,5 @@
     1.4  src/Sequents
     1.5  #misc components
     1.6  src/HOL/Tools/ATP_Manager
     1.7 +src/HOL/Library/Sum_Of_Squares
     1.8  
     2.1 --- a/etc/settings	Wed Aug 05 17:10:10 2009 +0200
     2.2 +++ b/etc/settings	Thu Aug 06 19:51:59 2009 +0200
     2.3 @@ -222,7 +222,6 @@
     2.4  #JEDIT_JAVA_OPTIONS="-server -Xms128m -Xmx512m"
     2.5  JEDIT_OPTIONS="-reuseview -noserver -nobackground"
     2.6  
     2.7 -
     2.8  ###
     2.9  ### External reasoning tools
    2.10  ###
    2.11 @@ -274,6 +273,9 @@
    2.12  # Jerusat 1.3 (SAT Solver, cf. Isabelle/src/HOL/Tools/sat_solver.ML)
    2.13  #JERUSAT_HOME=/usr/local/bin
    2.14  
    2.15 +# CSDP (SDP Solver, cf. Isabelle/src/HOL/Library/Sum_of_Squares/sos_wrapper.ML)
    2.16 +#CSDP_EXE=csdp
    2.17 +
    2.18  # For configuring HOL/Matrix/cplex
    2.19  # LP_SOLVER is the default solver. It can be changed during runtime via Cplex.set_solver.
    2.20  # First option: use the commercial cplex solver
     3.1 --- a/lib/scripts/neos/NeosCSDPClient.py	Wed Aug 05 17:10:10 2009 +0200
     3.2 +++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
     3.3 @@ -1,56 +0,0 @@
     3.4 -#!/usr/bin/env python
     3.5 -import sys
     3.6 -import xmlrpclib
     3.7 -import time
     3.8 -import re
     3.9 -
    3.10 -from config import Variables
    3.11 -
    3.12 -if len(sys.argv) < 3 or len(sys.argv) > 3:
    3.13 -  sys.stderr.write("Usage: NeosCSDPClient <input_filename> <output_filename>\n")
    3.14 -  sys.exit(1)
    3.15 -
    3.16 -neos=xmlrpclib.Server("http://%s:%d" % (Variables.NEOS_HOST, Variables.NEOS_PORT))
    3.17 -
    3.18 -xmlfile = open(sys.argv[1],"r")
    3.19 -xml_pre = "<document>\n<category>sdp</category>\n<solver>csdp</solver>\n<inputMethod>SPARSE_SDPA</inputMethod>\n<dat><![CDATA["
    3.20 -xml_post = "]]></dat>\n</document>\n"
    3.21 -xml = xml_pre
    3.22 -buffer = 1
    3.23 -while buffer:
    3.24 -  buffer = xmlfile.read()
    3.25 -  xml += buffer
    3.26 -xmlfile.close()
    3.27 -xml += xml_post
    3.28 -
    3.29 -(jobNumber,password) = neos.submitJob(xml)
    3.30 -
    3.31 -if jobNumber == 0:
    3.32 -  sys.stdout.write("error submitting job: %s" % password)
    3.33 -  sys.exit(-1)
    3.34 -else:
    3.35 -  sys.stdout.write("jobNumber = %d\tpassword = %s\n" % (jobNumber,password))
    3.36 -
    3.37 -offset=0
    3.38 -status="Waiting"
    3.39 -while status == "Running" or status=="Waiting":
    3.40 -  time.sleep(1)
    3.41 -  (msg,offset) = neos.getIntermediateResults(jobNumber,password,offset)
    3.42 -  sys.stdout.write(msg.data)
    3.43 -  status = neos.getJobStatus(jobNumber, password)
    3.44 -
    3.45 -msg = neos.getFinalResults(jobNumber, password).data
    3.46 -result = msg.split("Solution:")
    3.47 -
    3.48 -sys.stdout.write(result[0])
    3.49 -if len(result) > 1:
    3.50 -  plain_msg = result[1].strip()
    3.51 -  if plain_msg != "":
    3.52 -    output = open(sys.argv[2],"w")
    3.53 -    output.write(plain_msg)
    3.54 -    output.close()
    3.55 -    sys.exit(0)
    3.56 -
    3.57 -sys.exit(2)
    3.58 -
    3.59 -
     4.1 --- a/lib/scripts/neos/config.py	Wed Aug 05 17:10:10 2009 +0200
     4.2 +++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
     4.3 @@ -1,3 +0,0 @@
     4.4 -class Variables:
     4.5 -  NEOS_HOST="neos.mcs.anl.gov"
     4.6 -  NEOS_PORT=3332
     5.1 --- a/src/HOL/IsaMakefile	Wed Aug 05 17:10:10 2009 +0200
     5.2 +++ b/src/HOL/IsaMakefile	Thu Aug 06 19:51:59 2009 +0200
     5.3 @@ -316,49 +316,45 @@
     5.4  
     5.5  HOL-Library: HOL $(LOG)/HOL-Library.gz
     5.6  
     5.7 -
     5.8  $(LOG)/HOL-Library.gz: $(OUT)/HOL Library/SetsAndFunctions.thy		\
     5.9 -  Library/Abstract_Rat.thy \
    5.10 -  Library/BigO.thy Library/ContNotDenum.thy Library/Efficient_Nat.thy	\
    5.11 -  Library/Euclidean_Space.thy Library/Sum_Of_Squares.thy Library/positivstellensatz.ML	\
    5.12 -  Library/Fset.thy  Library/Convex_Euclidean_Space.thy \
    5.13 -  Library/sum_of_squares.ML Library/Glbs.thy Library/normarith.ML \
    5.14 -  Library/Executable_Set.thy Library/Infinite_Set.thy			\
    5.15 -  Library/FuncSet.thy Library/Permutations.thy Library/Determinants.thy\
    5.16 -  Library/Bit.thy Library/Topology_Euclidean_Space.thy \
    5.17 -  Library/Finite_Cartesian_Product.thy \
    5.18 -  Library/FrechetDeriv.thy Library/Fraction_Field.thy\
    5.19 -  Library/Fundamental_Theorem_Algebra.thy \
    5.20 -  Library/Inner_Product.thy Library/Kleene_Algebra.thy Library/Lattice_Syntax.thy \
    5.21 -  Library/Legacy_GCD.thy \
    5.22 -  Library/Library.thy Library/List_Prefix.thy Library/List_Set.thy Library/State_Monad.thy	\
    5.23 -  Library/Nat_Int_Bij.thy Library/Multiset.thy Library/Permutation.thy	\
    5.24 -  Library/Primes.thy Library/Pocklington.thy Library/Quotient.thy	\
    5.25 -  Library/Quicksort.thy Library/Nat_Infinity.thy Library/Word.thy	\
    5.26 -  Library/README.html Library/Continuity.thy Library/Order_Relation.thy \
    5.27 -  Library/Nested_Environment.thy Library/Ramsey.thy Library/Zorn.thy	\
    5.28 -  Library/Library/ROOT.ML Library/Library/document/root.tex		\
    5.29 -  Library/Library/document/root.bib Library/While_Combinator.thy	\
    5.30 -  Library/Product_ord.thy Library/Char_nat.thy Library/Char_ord.thy	\
    5.31 -  Library/Option_ord.thy Library/Sublist_Order.thy			\
    5.32 -  Library/List_lexord.thy Library/Commutative_Ring.thy			\
    5.33 -  Library/comm_ring.ML Library/Coinductive_List.thy			\
    5.34 -  Library/AssocList.thy	Library/Formal_Power_Series.thy	\
    5.35 -  Library/Binomial.thy Library/Eval_Witness.thy				\
    5.36 -  Library/Code_Char.thy				\
    5.37 +  Library/Abstract_Rat.thy Library/BigO.thy Library/ContNotDenum.thy	\
    5.38 +  Library/Efficient_Nat.thy Library/Euclidean_Space.thy			\
    5.39 +  Library/Sum_Of_Squares.thy Library/Sum_Of_Squares/sos_wrapper.ML	\
    5.40 +  Library/Sum_Of_Squares/sum_of_squares.ML Library/Fset.thy		\
    5.41 +  Library/Convex_Euclidean_Space.thy Library/Glbs.thy			\
    5.42 +  Library/normarith.ML Library/Executable_Set.thy			\
    5.43 +  Library/Infinite_Set.thy Library/FuncSet.thy				\
    5.44 +  Library/Permutations.thy Library/Determinants.thy Library/Bit.thy	\
    5.45 +  Library/Topology_Euclidean_Space.thy					\
    5.46 +  Library/Finite_Cartesian_Product.thy Library/FrechetDeriv.thy		\
    5.47 +  Library/Fraction_Field.thy Library/Fundamental_Theorem_Algebra.thy	\
    5.48 +  Library/Inner_Product.thy Library/Kleene_Algebra.thy			\
    5.49 +  Library/Lattice_Syntax.thy Library/Legacy_GCD.thy			\
    5.50 +  Library/Library.thy Library/List_Prefix.thy Library/List_Set.thy	\
    5.51 +  Library/State_Monad.thy Library/Nat_Int_Bij.thy Library/Multiset.thy	\
    5.52 +  Library/Permutation.thy Library/Primes.thy Library/Pocklington.thy	\
    5.53 +  Library/Quotient.thy Library/Quicksort.thy Library/Nat_Infinity.thy	\
    5.54 +  Library/Word.thy Library/README.html Library/Continuity.thy		\
    5.55 +  Library/Order_Relation.thy Library/Nested_Environment.thy		\
    5.56 +  Library/Ramsey.thy Library/Zorn.thy Library/Library/ROOT.ML		\
    5.57 +  Library/Library/document/root.tex Library/Library/document/root.bib	\
    5.58 +  Library/While_Combinator.thy Library/Product_ord.thy			\
    5.59 +  Library/Char_nat.thy Library/Char_ord.thy Library/Option_ord.thy	\
    5.60 +  Library/Sublist_Order.thy Library/List_lexord.thy			\
    5.61 +  Library/Commutative_Ring.thy Library/comm_ring.ML			\
    5.62 +  Library/Coinductive_List.thy Library/AssocList.thy			\
    5.63 +  Library/Formal_Power_Series.thy Library/Binomial.thy			\
    5.64 +  Library/Eval_Witness.thy Library/Code_Char.thy			\
    5.65    Library/Code_Char_chr.thy Library/Code_Integer.thy			\
    5.66    Library/Mapping.thy Library/Numeral_Type.thy Library/Reflection.thy	\
    5.67    Library/Boolean_Algebra.thy Library/Countable.thy			\
    5.68    Library/Diagonalize.thy Library/RBT.thy Library/Univ_Poly.thy		\
    5.69 -  Library/Poly_Deriv.thy \
    5.70 -  Library/Polynomial.thy \
    5.71 -  Library/Preorder.thy \
    5.72 -  Library/Product_plus.thy \
    5.73 -  Library/Product_Vector.thy \
    5.74 -  Library/Tree.thy \
    5.75 -  Library/Enum.thy Library/Float.thy $(SRC)/Tools/float.ML $(SRC)/HOL/Tools/float_arith.ML \
    5.76 -  Library/reify_data.ML Library/reflection.ML \
    5.77 -  Library/LaTeXsugar.thy Library/OptionalSugar.thy
    5.78 +  Library/Poly_Deriv.thy Library/Polynomial.thy Library/Preorder.thy	\
    5.79 +  Library/Product_plus.thy Library/Product_Vector.thy Library/Tree.thy	\
    5.80 +  Library/Enum.thy Library/Float.thy $(SRC)/Tools/float.ML		\
    5.81 +  $(SRC)/HOL/Tools/float_arith.ML Library/positivstellensatz.ML		\
    5.82 +  Library/reify_data.ML Library/reflection.ML Library/LaTeXsugar.thy	\
    5.83 +  Library/OptionalSugar.thy
    5.84  	@cd Library; $(ISABELLE_TOOL) usedir $(OUT)/HOL Library
    5.85  
    5.86  
     6.1 --- a/src/HOL/Library/Sum_Of_Squares.thy	Wed Aug 05 17:10:10 2009 +0200
     6.2 +++ b/src/HOL/Library/Sum_Of_Squares.thy	Thu Aug 06 19:51:59 2009 +0200
     6.3 @@ -2,9 +2,9 @@
     6.4     Author:     Amine Chaieb, University of Cambridge
     6.5  
     6.6  In order to use the method sos, call it with (sos remote_csdp) to use the remote solver
     6.7 -or install CSDP (https://projects.coin-or.org/Csdp/), put the executable csdp on your path,
     6.8 -and call it with (sos csdp). By default, sos calls remote_csdp. This can take of the order
     6.9 -of a minute for one sos call, because sos calls CSDP repeatedly.
    6.10 +or install CSDP (https://projects.coin-or.org/Csdp/), set the Isabelle environment
    6.11 +variable CSDP_EXE and call it with (sos csdp). By default, sos calls remote_csdp.
    6.12 +This can take of the order of a minute for one sos call, because sos calls CSDP repeatedly.
    6.13  If you install CSDP locally, sos calls typically takes only a few seconds.
    6.14  
    6.15  *)
    6.16 @@ -13,11 +13,19 @@
    6.17            multiplication and ordering using semidefinite programming*}
    6.18  
    6.19  theory Sum_Of_Squares
    6.20 -  imports Complex_Main (* "~~/src/HOL/Decision_Procs/Dense_Linear_Order" *)
    6.21 -  uses "positivstellensatz.ML" "sum_of_squares.ML" "sos_wrapper.ML"
    6.22 -  begin
    6.23 +imports Complex_Main (* "~~/src/HOL/Decision_Procs/Dense_Linear_Order" *)
    6.24 +uses
    6.25 +  ("positivstellensatz.ML")
    6.26 +  ("Sum_Of_Squares/sum_of_squares.ML")
    6.27 +  ("Sum_Of_Squares/sos_wrapper.ML")
    6.28 +begin
    6.29  
    6.30  (* setup sos tactic *)
    6.31 +
    6.32 +use "positivstellensatz.ML"
    6.33 +use "Sum_Of_Squares/sum_of_squares.ML"
    6.34 +use "Sum_Of_Squares/sos_wrapper.ML"
    6.35 +
    6.36  setup SosWrapper.setup
    6.37  
    6.38  text{* Tests -- commented since they work only when csdp is installed  or take too long with remote csdps *}
     7.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     7.2 +++ b/src/HOL/Library/Sum_Of_Squares/etc/settings	Thu Aug 06 19:51:59 2009 +0200
     7.3 @@ -0,0 +1,1 @@
     7.4 +ISABELLE_SUM_OF_SQUARES="$COMPONENT"
     8.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     8.2 +++ b/src/HOL/Library/Sum_Of_Squares/neos_csdp_client	Thu Aug 06 19:51:59 2009 +0200
     8.3 @@ -0,0 +1,65 @@
     8.4 +#!/usr/bin/env python
     8.5 +import sys
     8.6 +import xmlrpclib
     8.7 +import time
     8.8 +import re
     8.9 +
    8.10 +# Neos server config
    8.11 +NEOS_HOST="neos.mcs.anl.gov"
    8.12 +NEOS_PORT=3332
    8.13 +
    8.14 +if len(sys.argv) < 3 or len(sys.argv) > 3:
    8.15 +  sys.stderr.write("Usage: NeosCSDPClient <input_filename> <output_filename>\n")
    8.16 +  sys.exit(1)
    8.17 +
    8.18 +neos=xmlrpclib.Server("http://%s:%d" % (NEOS_HOST, NEOS_PORT))
    8.19 +
    8.20 +inputfile = open(sys.argv[1],"r")
    8.21 +xml_pre = "<document>\n<category>sdp</category>\n<solver>csdp</solver>\n<inputMethod>SPARSE_SDPA</inputMethod>\n<dat><![CDATA["
    8.22 +xml_post = "]]></dat>\n</document>\n"
    8.23 +xml = xml_pre
    8.24 +buffer = 1
    8.25 +while buffer:
    8.26 +  buffer = inputfile.read()
    8.27 +  xml += buffer
    8.28 +inputfile.close()
    8.29 +xml += xml_post
    8.30 +
    8.31 +(jobNumber,password) = neos.submitJob(xml)
    8.32 +
    8.33 +if jobNumber == 0:
    8.34 +  sys.stdout.write("error submitting job: %s" % password)
    8.35 +  sys.exit(20)
    8.36 +else:
    8.37 +  sys.stdout.write("jobNumber = %d\tpassword = %s\n" % (jobNumber,password))
    8.38 +
    8.39 +offset=0
    8.40 +messages = ""
    8.41 +status="Waiting"
    8.42 +while status == "Running" or status=="Waiting":
    8.43 +  time.sleep(1)
    8.44 +  (msg,offset) = neos.getIntermediateResults(jobNumber,password,offset)
    8.45 +  messages += msg.data
    8.46 +  sys.stdout.write(msg.data)
    8.47 +  status = neos.getJobStatus(jobNumber, password)
    8.48 +
    8.49 +msg = neos.getFinalResults(jobNumber, password).data
    8.50 +sys.stdout.write("---------- Begin CSDP Output -------------\n");
    8.51 +sys.stdout.write(msg)
    8.52 +
    8.53 +# extract solution
    8.54 +result = msg.split("Solution:")
    8.55 +if len(result) > 1:
    8.56 +  output = open(sys.argv[2],"w")
    8.57 +  output.write(result[1].strip())
    8.58 +  output.close()
    8.59 +
    8.60 +# extract return code
    8.61 +p = re.compile(r"^Error: Command exited with non-zero status (\d+)$", re.MULTILINE)
    8.62 +m = p.search(messages)
    8.63 +if m:
    8.64 +  sys.exit(int(m.group(1)))
    8.65 +else:
    8.66 +  sys.exit(0)
    8.67 +
    8.68 +
     9.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     9.2 +++ b/src/HOL/Library/Sum_Of_Squares/sos_wrapper.ML	Thu Aug 06 19:51:59 2009 +0200
     9.3 @@ -0,0 +1,140 @@
     9.4 +(* Title:      sos_wrapper.ML
     9.5 +   Author:     Philipp Meyer, TU Muenchen
     9.6 +
     9.7 +Added functionality for sums of squares, e.g. calling a remote prover
     9.8 +*)
     9.9 +
    9.10 +signature SOS_WRAPPER =
    9.11 +sig
    9.12 +
    9.13 +  datatype prover_result = Success | Failure | Error
    9.14 +
    9.15 +  val setup: theory -> theory
    9.16 +  val destdir: string ref
    9.17 +end
    9.18 +
    9.19 +structure SosWrapper : SOS_WRAPPER=
    9.20 +struct
    9.21 +
    9.22 +datatype prover_result = Success | Failure | Error
    9.23 +fun str_of_result Success = "Success"
    9.24 +  | str_of_result Failure = "Failure"
    9.25 +  | str_of_result Error = "Error"
    9.26 +
    9.27 +(*** output control ***)
    9.28 +
    9.29 +fun debug s = if ! Sos.debugging then Output.writeln s else ()
    9.30 +val write = Output.priority
    9.31 +
    9.32 +(*** calling provers ***)
    9.33 +
    9.34 +val destdir = ref ""
    9.35 +
    9.36 +fun filename dir name =
    9.37 +  let
    9.38 +    val probfile = Path.basic (name ^ serial_string ())
    9.39 +    val dir_path = Path.explode dir
    9.40 +  in
    9.41 +    if dir = "" then
    9.42 +      File.tmp_path probfile
    9.43 +    else
    9.44 +      if File.exists dir_path then
    9.45 +        Path.append dir_path probfile
    9.46 +      else
    9.47 +        error ("No such directory: " ^ dir)
    9.48 +  end
    9.49 +
    9.50 +fun run_solver name cmd find_failure input =
    9.51 +  let
    9.52 +    val _ = write ("Calling solver: " ^ name)
    9.53 +
    9.54 +    (* create input file *)
    9.55 +    val dir = ! destdir
    9.56 +    val input_file = filename dir "sos_in" 
    9.57 +    val _ = File.write input_file input
    9.58 +
    9.59 +    (* call solver *)
    9.60 +    val output_file = filename dir "sos_out"
    9.61 +    val (output, rv) = system_out (
    9.62 +      if File.exists cmd then space_implode " "
    9.63 +        [File.shell_path cmd, File.platform_path input_file, File.platform_path output_file]
    9.64 +      else error ("Bad executable: " ^ File.shell_path cmd))
    9.65 + 
    9.66 +    (* read and analysize output *)
    9.67 +    val (res, res_msg) = find_failure rv
    9.68 +    val result = if File.exists output_file then File.read output_file else ""
    9.69 +
    9.70 +    (* remove temporary files *)
    9.71 +    val _ = if dir = "" then
    9.72 +        (File.rm input_file ; if File.exists output_file then File.rm output_file else ())
    9.73 +        else ()
    9.74 +
    9.75 +    val _ = debug ("Solver output:\n" ^ output)
    9.76 +
    9.77 +    val _ = write (str_of_result res ^ ": " ^ res_msg)
    9.78 +  in
    9.79 +    case res of
    9.80 +      Success => result
    9.81 +    | Failure => raise Sos.Failure res_msg
    9.82 +    | Error => error ("Prover failed: " ^ res_msg)
    9.83 +  end
    9.84 +
    9.85 +(*** various provers ***)
    9.86 +
    9.87 +(* local csdp client *)
    9.88 +
    9.89 +fun find_csdp_failure rv =
    9.90 +  case rv of
    9.91 +    0 => (Success, "SDP solved")
    9.92 +  | 1 => (Failure, "SDP is primal infeasible")
    9.93 +  | 2 => (Failure, "SDP is dual infeasible")
    9.94 +  | 3 => (Success, "SDP solved with reduced accuracy")
    9.95 +  | 4 => (Failure, "Maximum iterations reached")
    9.96 +  | 5 => (Failure, "Stuck at edge of primal feasibility")
    9.97 +  | 6 => (Failure, "Stuck at edge of dual infeasibility")
    9.98 +  | 7 => (Failure, "Lack of progress")
    9.99 +  | 8 => (Failure, "X, Z, or O was singular")
   9.100 +  | 9 => (Failure, "Detected NaN or Inf values")
   9.101 +  | _ => (Error, "return code is " ^ string_of_int rv)
   9.102 +
   9.103 +val csdp = ("$CSDP_EXE", find_csdp_failure)
   9.104 +
   9.105 +(* remote neos server *)
   9.106 +
   9.107 +fun find_neos_failure rv =
   9.108 +  case rv of
   9.109 +    20 => (Error, "error submitting job")
   9.110 +  | 21 => (Error, "no solution")
   9.111 +  |  _ => find_csdp_failure rv
   9.112 +
   9.113 +val neos_csdp = ("$ISABELLE_SUM_OF_SQUARES/neos_csdp_client", find_neos_failure)
   9.114 +
   9.115 +(* save provers in table *)
   9.116 +
   9.117 +val provers =
   9.118 +     Symtab.make [("remote_csdp", neos_csdp),("csdp", csdp)]
   9.119 +
   9.120 +fun get_prover name =
   9.121 +  case Symtab.lookup provers name of
   9.122 +    SOME prover => prover
   9.123 +  | NONE => error ("unknown prover: " ^ name)
   9.124 +
   9.125 +fun call_solver name =
   9.126 +  let
   9.127 +    val (cmd, find_failure) = get_prover name
   9.128 +  in
   9.129 +    run_solver name (Path.explode cmd) find_failure
   9.130 +  end
   9.131 +
   9.132 +(* setup tactic *)
   9.133 +
   9.134 +val def_solver = "remote_csdp"
   9.135 +
   9.136 +fun sos_solver name = (SIMPLE_METHOD' o (Sos.sos_tac (call_solver name))) 
   9.137 +
   9.138 +val sos_method = Scan.optional (Scan.lift OuterParse.xname) def_solver >> sos_solver
   9.139 +
   9.140 +val setup = Method.setup @{binding sos} sos_method
   9.141 +  "Prove universal problems over the reals using sums of squares"
   9.142 +
   9.143 +end
    10.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
    10.2 +++ b/src/HOL/Library/Sum_Of_Squares/sum_of_squares.ML	Thu Aug 06 19:51:59 2009 +0200
    10.3 @@ -0,0 +1,1512 @@
    10.4 +(* Title:      sum_of_squares.ML
    10.5 +   Authors:    Amine Chaieb, University of Cambridge
    10.6 +               Philipp Meyer, TU Muenchen
    10.7 +
    10.8 +A tactic for proving nonlinear inequalities
    10.9 +*)
   10.10 +
   10.11 +signature SOS =
   10.12 +sig
   10.13 +
   10.14 +  val sos_tac : (string -> string) -> Proof.context -> int -> Tactical.tactic
   10.15 +
   10.16 +  val debugging : bool ref;
   10.17 +  
   10.18 +  exception Failure of string;
   10.19 +end
   10.20 +
   10.21 +structure Sos : SOS = 
   10.22 +struct
   10.23 +
   10.24 +val rat_0 = Rat.zero;
   10.25 +val rat_1 = Rat.one;
   10.26 +val rat_2 = Rat.two;
   10.27 +val rat_10 = Rat.rat_of_int 10;
   10.28 +val rat_1_2 = rat_1 // rat_2;
   10.29 +val max = curry IntInf.max;
   10.30 +val min = curry IntInf.min;
   10.31 +
   10.32 +val denominator_rat = Rat.quotient_of_rat #> snd #> Rat.rat_of_int;
   10.33 +val numerator_rat = Rat.quotient_of_rat #> fst #> Rat.rat_of_int;
   10.34 +fun int_of_rat a = 
   10.35 +    case Rat.quotient_of_rat a of (i,1) => i | _ => error "int_of_rat: not an int";
   10.36 +fun lcm_rat x y = Rat.rat_of_int (Integer.lcm (int_of_rat x) (int_of_rat y));
   10.37 +
   10.38 +fun rat_pow r i = 
   10.39 + let fun pow r i = 
   10.40 +   if i = 0 then rat_1 else 
   10.41 +   let val d = pow r (i div 2)
   10.42 +   in d */ d */ (if i mod 2 = 0 then rat_1 else r)
   10.43 +   end
   10.44 + in if i < 0 then pow (Rat.inv r) (~ i) else pow r i end;
   10.45 +
   10.46 +fun round_rat r = 
   10.47 + let val (a,b) = Rat.quotient_of_rat (Rat.abs r)
   10.48 +     val d = a div b
   10.49 +     val s = if r </ rat_0 then (Rat.neg o Rat.rat_of_int) else Rat.rat_of_int
   10.50 +     val x2 = 2 * (a - (b * d))
   10.51 + in s (if x2 >= b then d + 1 else d) end
   10.52 +
   10.53 +val abs_rat = Rat.abs;
   10.54 +val pow2 = rat_pow rat_2;
   10.55 +val pow10 = rat_pow rat_10;
   10.56 +
   10.57 +val debugging = ref false;
   10.58 +
   10.59 +exception Sanity;
   10.60 +
   10.61 +exception Unsolvable;
   10.62 +
   10.63 +exception Failure of string;
   10.64 +
   10.65 +(* Turn a rational into a decimal string with d sig digits.                  *)
   10.66 +
   10.67 +local
   10.68 +fun normalize y =
   10.69 +  if abs_rat y </ (rat_1 // rat_10) then normalize (rat_10 */ y) - 1
   10.70 +  else if abs_rat y >=/ rat_1 then normalize (y // rat_10) + 1
   10.71 +  else 0 
   10.72 + in
   10.73 +fun decimalize d x =
   10.74 +  if x =/ rat_0 then "0.0" else
   10.75 +  let 
   10.76 +   val y = Rat.abs x
   10.77 +   val e = normalize y
   10.78 +   val z = pow10(~ e) */ y +/ rat_1
   10.79 +   val k = int_of_rat (round_rat(pow10 d */ z))
   10.80 +  in (if x </ rat_0 then "-0." else "0.") ^
   10.81 +     implode(tl(explode(string_of_int k))) ^
   10.82 +     (if e = 0 then "" else "e"^string_of_int e)
   10.83 +  end
   10.84 +end;
   10.85 +
   10.86 +(* Iterations over numbers, and lists indexed by numbers.                    *)
   10.87 +
   10.88 +fun itern k l f a =
   10.89 +  case l of
   10.90 +    [] => a
   10.91 +  | h::t => itern (k + 1) t f (f h k a);
   10.92 +
   10.93 +fun iter (m,n) f a =
   10.94 +  if n < m then a
   10.95 +  else iter (m+1,n) f (f m a);
   10.96 +
   10.97 +(* The main types.                                                           *)
   10.98 +
   10.99 +fun strict_ord ord (x,y) = case ord (x,y) of LESS => LESS | _ => GREATER
  10.100 +
  10.101 +structure Intpairfunc = FuncFun(type key = int*int val ord = prod_ord int_ord int_ord);
  10.102 +
  10.103 +type vector = int* Rat.rat Intfunc.T;
  10.104 +
  10.105 +type matrix = (int*int)*(Rat.rat Intpairfunc.T);
  10.106 +
  10.107 +type monomial = int Ctermfunc.T;
  10.108 +
  10.109 +val cterm_ord = (fn (s,t) => TermOrd.fast_term_ord(term_of s, term_of t))
  10.110 + fun monomial_ord (m1,m2) = list_ord (prod_ord cterm_ord int_ord) (Ctermfunc.graph m1, Ctermfunc.graph m2)
  10.111 +structure Monomialfunc = FuncFun(type key = monomial val ord = monomial_ord)
  10.112 +
  10.113 +type poly = Rat.rat Monomialfunc.T;
  10.114 +
  10.115 + fun iszero (k,r) = r =/ rat_0;
  10.116 +
  10.117 +fun fold_rev2 f l1 l2 b =
  10.118 +  case (l1,l2) of
  10.119 +    ([],[]) => b
  10.120 +  | (h1::t1,h2::t2) => f h1 h2 (fold_rev2 f t1 t2 b)
  10.121 +  | _ => error "fold_rev2";
  10.122 + 
  10.123 +(* Vectors. Conventionally indexed 1..n.                                     *)
  10.124 +
  10.125 +fun vector_0 n = (n,Intfunc.undefined):vector;
  10.126 +
  10.127 +fun dim (v:vector) = fst v;
  10.128 +
  10.129 +fun vector_const c n =
  10.130 +  if c =/ rat_0 then vector_0 n
  10.131 +  else (n,fold_rev (fn k => Intfunc.update (k,c)) (1 upto n) Intfunc.undefined) :vector;
  10.132 +
  10.133 +val vector_1 = vector_const rat_1;
  10.134 +
  10.135 +fun vector_cmul c (v:vector) =
  10.136 + let val n = dim v 
  10.137 + in if c =/ rat_0 then vector_0 n
  10.138 +    else (n,Intfunc.mapf (fn x => c */ x) (snd v))
  10.139 + end;
  10.140 +
  10.141 +fun vector_neg (v:vector) = (fst v,Intfunc.mapf Rat.neg (snd v)) :vector;
  10.142 +
  10.143 +fun vector_add (v1:vector) (v2:vector) =
  10.144 + let val m = dim v1  
  10.145 +     val n = dim v2 
  10.146 + in if m <> n then error "vector_add: incompatible dimensions"
  10.147 +    else (n,Intfunc.combine (curry op +/) (fn x => x =/ rat_0) (snd v1) (snd v2)) :vector 
  10.148 + end;
  10.149 +
  10.150 +fun vector_sub v1 v2 = vector_add v1 (vector_neg v2);
  10.151 +
  10.152 +fun vector_dot (v1:vector) (v2:vector) =
  10.153 + let val m = dim v1 
  10.154 +     val n = dim v2 
  10.155 + in if m <> n then error "vector_dot: incompatible dimensions" 
  10.156 +    else Intfunc.fold (fn (i,x) => fn a => x +/ a) 
  10.157 +        (Intfunc.combine (curry op */) (fn x => x =/ rat_0) (snd v1) (snd v2)) rat_0
  10.158 + end;
  10.159 +
  10.160 +fun vector_of_list l =
  10.161 + let val n = length l 
  10.162 + in (n,fold_rev2 (curry Intfunc.update) (1 upto n) l Intfunc.undefined) :vector
  10.163 + end;
  10.164 +
  10.165 +(* Matrices; again rows and columns indexed from 1.                          *)
  10.166 +
  10.167 +fun matrix_0 (m,n) = ((m,n),Intpairfunc.undefined):matrix;
  10.168 +
  10.169 +fun dimensions (m:matrix) = fst m;
  10.170 +
  10.171 +fun matrix_const c (mn as (m,n)) =
  10.172 +  if m <> n then error "matrix_const: needs to be square"
  10.173 +  else if c =/ rat_0 then matrix_0 mn
  10.174 +  else (mn,fold_rev (fn k => Intpairfunc.update ((k,k), c)) (1 upto n) Intpairfunc.undefined) :matrix;;
  10.175 +
  10.176 +val matrix_1 = matrix_const rat_1;
  10.177 +
  10.178 +fun matrix_cmul c (m:matrix) =
  10.179 + let val (i,j) = dimensions m 
  10.180 + in if c =/ rat_0 then matrix_0 (i,j)
  10.181 +    else ((i,j),Intpairfunc.mapf (fn x => c */ x) (snd m))
  10.182 + end;
  10.183 +
  10.184 +fun matrix_neg (m:matrix) = 
  10.185 +  (dimensions m, Intpairfunc.mapf Rat.neg (snd m)) :matrix;
  10.186 +
  10.187 +fun matrix_add (m1:matrix) (m2:matrix) =
  10.188 + let val d1 = dimensions m1 
  10.189 +     val d2 = dimensions m2 
  10.190 + in if d1 <> d2 
  10.191 +     then error "matrix_add: incompatible dimensions"
  10.192 +    else (d1,Intpairfunc.combine (curry op +/) (fn x => x =/ rat_0) (snd m1) (snd m2)) :matrix
  10.193 + end;;
  10.194 +
  10.195 +fun matrix_sub m1 m2 = matrix_add m1 (matrix_neg m2);
  10.196 +
  10.197 +fun row k (m:matrix) =
  10.198 + let val (i,j) = dimensions m 
  10.199 + in (j,
  10.200 +   Intpairfunc.fold (fn ((i,j), c) => fn a => if i = k then Intfunc.update (j,c) a else a) (snd m) Intfunc.undefined ) : vector
  10.201 + end;
  10.202 +
  10.203 +fun column k (m:matrix) =
  10.204 +  let val (i,j) = dimensions m 
  10.205 +  in (i,
  10.206 +   Intpairfunc.fold (fn ((i,j), c) => fn a => if j = k then Intfunc.update (i,c) a else a) (snd m)  Intfunc.undefined)
  10.207 +   : vector
  10.208 + end;
  10.209 +
  10.210 +fun transp (m:matrix) =
  10.211 +  let val (i,j) = dimensions m 
  10.212 +  in
  10.213 +  ((j,i),Intpairfunc.fold (fn ((i,j), c) => fn a => Intpairfunc.update ((j,i), c) a) (snd m) Intpairfunc.undefined) :matrix
  10.214 + end;
  10.215 +
  10.216 +fun diagonal (v:vector) =
  10.217 + let val n = dim v 
  10.218 + in ((n,n),Intfunc.fold (fn (i, c) => fn a => Intpairfunc.update ((i,i), c) a) (snd v) Intpairfunc.undefined) : matrix
  10.219 + end;
  10.220 +
  10.221 +fun matrix_of_list l =
  10.222 + let val m = length l 
  10.223 + in if m = 0 then matrix_0 (0,0) else
  10.224 +   let val n = length (hd l) 
  10.225 +   in ((m,n),itern 1 l (fn v => fn i => itern 1 v (fn c => fn j => Intpairfunc.update ((i,j), c))) Intpairfunc.undefined)
  10.226 +   end
  10.227 + end;
  10.228 +
  10.229 +(* Monomials.                                                                *)
  10.230 +
  10.231 +fun monomial_eval assig (m:monomial) =
  10.232 +  Ctermfunc.fold (fn (x, k) => fn a => a */ rat_pow (Ctermfunc.apply assig x) k)
  10.233 +        m rat_1;
  10.234 +val monomial_1 = (Ctermfunc.undefined:monomial);
  10.235 +
  10.236 +fun monomial_var x = Ctermfunc.onefunc (x, 1) :monomial;
  10.237 +
  10.238 +val (monomial_mul:monomial->monomial->monomial) =
  10.239 +  Ctermfunc.combine (curry op +) (K false);
  10.240 +
  10.241 +fun monomial_pow (m:monomial) k =
  10.242 +  if k = 0 then monomial_1
  10.243 +  else Ctermfunc.mapf (fn x => k * x) m;
  10.244 +
  10.245 +fun monomial_divides (m1:monomial) (m2:monomial) =
  10.246 +  Ctermfunc.fold (fn (x, k) => fn a => Ctermfunc.tryapplyd m2 x 0 >= k andalso a) m1 true;;
  10.247 +
  10.248 +fun monomial_div (m1:monomial) (m2:monomial) =
  10.249 + let val m = Ctermfunc.combine (curry op +) 
  10.250 +   (fn x => x = 0) m1 (Ctermfunc.mapf (fn x => ~ x) m2) 
  10.251 + in if Ctermfunc.fold (fn (x, k) => fn a => k >= 0 andalso a) m true then m
  10.252 +  else error "monomial_div: non-divisible"
  10.253 + end;
  10.254 +
  10.255 +fun monomial_degree x (m:monomial) = 
  10.256 +  Ctermfunc.tryapplyd m x 0;;
  10.257 +
  10.258 +fun monomial_lcm (m1:monomial) (m2:monomial) =
  10.259 +  fold_rev (fn x => Ctermfunc.update (x, max (monomial_degree x m1) (monomial_degree x m2)))
  10.260 +          (gen_union (is_equal o  cterm_ord) (Ctermfunc.dom m1, Ctermfunc.dom m2)) (Ctermfunc.undefined :monomial);
  10.261 +
  10.262 +fun monomial_multidegree (m:monomial) = 
  10.263 + Ctermfunc.fold (fn (x, k) => fn a => k + a) m 0;;
  10.264 +
  10.265 +fun monomial_variables m = Ctermfunc.dom m;;
  10.266 +
  10.267 +(* Polynomials.                                                              *)
  10.268 +
  10.269 +fun eval assig (p:poly) =
  10.270 +  Monomialfunc.fold (fn (m, c) => fn a => a +/ c */ monomial_eval assig m) p rat_0;
  10.271 +
  10.272 +val poly_0 = (Monomialfunc.undefined:poly);
  10.273 +
  10.274 +fun poly_isconst (p:poly) = 
  10.275 +  Monomialfunc.fold (fn (m, c) => fn a => Ctermfunc.is_undefined m andalso a) p true;
  10.276 +
  10.277 +fun poly_var x = Monomialfunc.onefunc (monomial_var x,rat_1) :poly;
  10.278 +
  10.279 +fun poly_const c =
  10.280 +  if c =/ rat_0 then poly_0 else Monomialfunc.onefunc(monomial_1, c);
  10.281 +
  10.282 +fun poly_cmul c (p:poly) =
  10.283 +  if c =/ rat_0 then poly_0
  10.284 +  else Monomialfunc.mapf (fn x => c */ x) p;
  10.285 +
  10.286 +fun poly_neg (p:poly) = (Monomialfunc.mapf Rat.neg p :poly);;
  10.287 +
  10.288 +fun poly_add (p1:poly) (p2:poly) =
  10.289 +  (Monomialfunc.combine (curry op +/) (fn x => x =/ rat_0) p1 p2 :poly);
  10.290 +
  10.291 +fun poly_sub p1 p2 = poly_add p1 (poly_neg p2);
  10.292 +
  10.293 +fun poly_cmmul (c,m) (p:poly) =
  10.294 + if c =/ rat_0 then poly_0
  10.295 + else if Ctermfunc.is_undefined m 
  10.296 +      then Monomialfunc.mapf (fn d => c */ d) p
  10.297 +      else Monomialfunc.fold (fn (m', d) => fn a => (Monomialfunc.update (monomial_mul m m', c */ d) a)) p poly_0;
  10.298 +
  10.299 +fun poly_mul (p1:poly) (p2:poly) =
  10.300 +  Monomialfunc.fold (fn (m, c) => fn a => poly_add (poly_cmmul (c,m) p2) a) p1 poly_0;
  10.301 +
  10.302 +fun poly_div (p1:poly) (p2:poly) =
  10.303 + if not(poly_isconst p2) 
  10.304 + then error "poly_div: non-constant" else
  10.305 + let val c = eval Ctermfunc.undefined p2 
  10.306 + in if c =/ rat_0 then error "poly_div: division by zero"
  10.307 +    else poly_cmul (Rat.inv c) p1
  10.308 + end;
  10.309 +
  10.310 +fun poly_square p = poly_mul p p;
  10.311 +
  10.312 +fun poly_pow p k =
  10.313 + if k = 0 then poly_const rat_1
  10.314 + else if k = 1 then p
  10.315 + else let val q = poly_square(poly_pow p (k div 2)) in
  10.316 +      if k mod 2 = 1 then poly_mul p q else q end;
  10.317 +
  10.318 +fun poly_exp p1 p2 =
  10.319 +  if not(poly_isconst p2) 
  10.320 +  then error "poly_exp: not a constant" 
  10.321 +  else poly_pow p1 (int_of_rat (eval Ctermfunc.undefined p2));
  10.322 +
  10.323 +fun degree x (p:poly) = 
  10.324 + Monomialfunc.fold (fn (m,c) => fn a => max (monomial_degree x m) a) p 0;
  10.325 +
  10.326 +fun multidegree (p:poly) =
  10.327 +  Monomialfunc.fold (fn (m, c) => fn a => max (monomial_multidegree m) a) p 0;
  10.328 +
  10.329 +fun poly_variables (p:poly) =
  10.330 +  sort cterm_ord (Monomialfunc.fold_rev (fn (m, c) => curry (gen_union (is_equal o  cterm_ord)) (monomial_variables m)) p []);;
  10.331 +
  10.332 +(* Order monomials for human presentation.                                   *)
  10.333 +
  10.334 +fun cterm_ord (t,t') = TermOrd.fast_term_ord (term_of t, term_of t');
  10.335 +
  10.336 +val humanorder_varpow = prod_ord cterm_ord (rev_order o int_ord);
  10.337 +
  10.338 +local
  10.339 + fun ord (l1,l2) = case (l1,l2) of
  10.340 +  (_,[]) => LESS 
  10.341 + | ([],_) => GREATER
  10.342 + | (h1::t1,h2::t2) => 
  10.343 +   (case humanorder_varpow (h1, h2) of 
  10.344 +     LESS => LESS
  10.345 +   | EQUAL => ord (t1,t2)
  10.346 +   | GREATER => GREATER)
  10.347 +in fun humanorder_monomial m1 m2 = 
  10.348 + ord (sort humanorder_varpow (Ctermfunc.graph m1),
  10.349 +  sort humanorder_varpow (Ctermfunc.graph m2))
  10.350 +end;
  10.351 +
  10.352 +fun fold1 f l =  case l of
  10.353 +   []     => error "fold1"
  10.354 + | [x]    => x
  10.355 + | (h::t) => f h (fold1 f t);
  10.356 +
  10.357 +(* Conversions to strings.                                                   *)
  10.358 +
  10.359 +fun string_of_vector min_size max_size (v:vector) =
  10.360 + let val n_raw = dim v 
  10.361 + in if n_raw = 0 then "[]" else
  10.362 +  let 
  10.363 +   val n = max min_size (min n_raw max_size) 
  10.364 +   val xs = map (Rat.string_of_rat o (fn i => Intfunc.tryapplyd (snd v) i rat_0)) (1 upto n) 
  10.365 +  in "[" ^ fold1 (fn s => fn t => s ^ ", " ^ t) xs ^
  10.366 +  (if n_raw > max_size then ", ...]" else "]")
  10.367 +  end
  10.368 + end;
  10.369 +
  10.370 +fun string_of_matrix max_size (m:matrix) =
  10.371 + let 
  10.372 +  val (i_raw,j_raw) = dimensions m
  10.373 +  val i = min max_size i_raw 
  10.374 +  val j = min max_size j_raw
  10.375 +  val rstr = map (fn k => string_of_vector j j (row k m)) (1 upto i) 
  10.376 + in "["^ fold1 (fn s => fn t => s^";\n "^t) rstr ^
  10.377 +  (if j > max_size then "\n ...]" else "]")
  10.378 + end;
  10.379 +
  10.380 +fun string_of_term t = 
  10.381 + case t of
  10.382 +   a$b => "("^(string_of_term a)^" "^(string_of_term b)^")"
  10.383 + | Abs x => 
  10.384 +    let val (xn, b) = Term.dest_abs x
  10.385 +    in "(\\"^xn^"."^(string_of_term b)^")"
  10.386 +    end
  10.387 + | Const(s,_) => s
  10.388 + | Free (s,_) => s
  10.389 + | Var((s,_),_) => s
  10.390 + | _ => error "string_of_term";
  10.391 +
  10.392 +val string_of_cterm = string_of_term o term_of;
  10.393 +
  10.394 +fun string_of_varpow x k =
  10.395 +  if k = 1 then string_of_cterm x 
  10.396 +  else string_of_cterm x^"^"^string_of_int k;
  10.397 +
  10.398 +fun string_of_monomial m =
  10.399 + if Ctermfunc.is_undefined m then "1" else
  10.400 + let val vps = fold_rev (fn (x,k) => fn a => string_of_varpow x k :: a)
  10.401 +  (sort humanorder_varpow (Ctermfunc.graph m)) [] 
  10.402 + in fold1 (fn s => fn t => s^"*"^t) vps
  10.403 + end;
  10.404 +
  10.405 +fun string_of_cmonomial (c,m) =
  10.406 + if Ctermfunc.is_undefined m then Rat.string_of_rat c
  10.407 + else if c =/ rat_1 then string_of_monomial m
  10.408 + else Rat.string_of_rat c ^ "*" ^ string_of_monomial m;;
  10.409 +
  10.410 +fun string_of_poly (p:poly) =
  10.411 + if Monomialfunc.is_undefined p then "<<0>>" else
  10.412 + let 
  10.413 +  val cms = sort (fn ((m1,_),(m2,_)) => humanorder_monomial m1  m2) (Monomialfunc.graph p)
  10.414 +  val s = fold (fn (m,c) => fn a =>
  10.415 +             if c </ rat_0 then a ^ " - " ^ string_of_cmonomial(Rat.neg c,m)
  10.416 +             else a ^ " + " ^ string_of_cmonomial(c,m))
  10.417 +          cms ""
  10.418 +  val s1 = String.substring (s, 0, 3)
  10.419 +  val s2 = String.substring (s, 3, String.size s - 3) 
  10.420 + in "<<" ^(if s1 = " + " then s2 else "-"^s2)^">>"
  10.421 + end;
  10.422 +
  10.423 +(* Conversion from HOL term.                                                 *)
  10.424 +
  10.425 +local
  10.426 + val neg_tm = @{cterm "uminus :: real => _"}
  10.427 + val add_tm = @{cterm "op + :: real => _"}
  10.428 + val sub_tm = @{cterm "op - :: real => _"}
  10.429 + val mul_tm = @{cterm "op * :: real => _"}
  10.430 + val inv_tm = @{cterm "inverse :: real => _"}
  10.431 + val div_tm = @{cterm "op / :: real => _"}
  10.432 + val pow_tm = @{cterm "op ^ :: real => _"}
  10.433 + val zero_tm = @{cterm "0:: real"}
  10.434 + val is_numeral = can (HOLogic.dest_number o term_of)
  10.435 + fun is_comb t = case t of _$_ => true | _ => false
  10.436 + fun poly_of_term tm =
  10.437 +  if tm aconvc zero_tm then poly_0
  10.438 +  else if RealArith.is_ratconst tm 
  10.439 +       then poly_const(RealArith.dest_ratconst tm)
  10.440 +  else 
  10.441 +  (let val (lop,r) = Thm.dest_comb tm
  10.442 +   in if lop aconvc neg_tm then poly_neg(poly_of_term r)
  10.443 +      else if lop aconvc inv_tm then
  10.444 +       let val p = poly_of_term r 
  10.445 +       in if poly_isconst p 
  10.446 +          then poly_const(Rat.inv (eval Ctermfunc.undefined p))
  10.447 +          else error "poly_of_term: inverse of non-constant polyomial"
  10.448 +       end
  10.449 +   else (let val (opr,l) = Thm.dest_comb lop
  10.450 +         in 
  10.451 +          if opr aconvc pow_tm andalso is_numeral r 
  10.452 +          then poly_pow (poly_of_term l) ((snd o HOLogic.dest_number o term_of) r)
  10.453 +          else if opr aconvc add_tm 
  10.454 +           then poly_add (poly_of_term l) (poly_of_term r)
  10.455 +          else if opr aconvc sub_tm 
  10.456 +           then poly_sub (poly_of_term l) (poly_of_term r)
  10.457 +          else if opr aconvc mul_tm 
  10.458 +           then poly_mul (poly_of_term l) (poly_of_term r)
  10.459 +          else if opr aconvc div_tm 
  10.460 +           then let 
  10.461 +                  val p = poly_of_term l 
  10.462 +                  val q = poly_of_term r 
  10.463 +                in if poly_isconst q then poly_cmul (Rat.inv (eval Ctermfunc.undefined q)) p
  10.464 +                   else error "poly_of_term: division by non-constant polynomial"
  10.465 +                end
  10.466 +          else poly_var tm
  10.467 + 
  10.468 +         end
  10.469 +         handle CTERM ("dest_comb",_) => poly_var tm)
  10.470 +   end
  10.471 +   handle CTERM ("dest_comb",_) => poly_var tm)
  10.472 +in
  10.473 +val poly_of_term = fn tm =>
  10.474 + if type_of (term_of tm) = @{typ real} then poly_of_term tm
  10.475 + else error "poly_of_term: term does not have real type"
  10.476 +end;
  10.477 +
  10.478 +(* String of vector (just a list of space-separated numbers).                *)
  10.479 +
  10.480 +fun sdpa_of_vector (v:vector) =
  10.481 + let 
  10.482 +  val n = dim v
  10.483 +  val strs = map (decimalize 20 o (fn i => Intfunc.tryapplyd (snd v) i rat_0)) (1 upto n) 
  10.484 + in fold1 (fn x => fn y => x ^ " " ^ y) strs ^ "\n"
  10.485 + end;
  10.486 +
  10.487 +fun increasing f ord (x,y) = ord (f x, f y);
  10.488 +fun triple_int_ord ((a,b,c),(a',b',c')) = 
  10.489 + prod_ord int_ord (prod_ord int_ord int_ord) 
  10.490 +    ((a,(b,c)),(a',(b',c')));
  10.491 +structure Inttriplefunc = FuncFun(type key = int*int*int val ord = triple_int_ord);
  10.492 +
  10.493 +(* String for block diagonal matrix numbered k.                              *)
  10.494 +
  10.495 +fun sdpa_of_blockdiagonal k m =
  10.496 + let 
  10.497 +  val pfx = string_of_int k ^" "
  10.498 +  val ents =
  10.499 +    Inttriplefunc.fold (fn ((b,i,j), c) => fn a => if i > j then a else ((b,i,j),c)::a) m []
  10.500 +  val entss = sort (increasing fst triple_int_ord ) ents
  10.501 + in  fold_rev (fn ((b,i,j),c) => fn a =>
  10.502 +     pfx ^ string_of_int b ^ " " ^ string_of_int i ^ " " ^ string_of_int j ^
  10.503 +     " " ^ decimalize 20 c ^ "\n" ^ a) entss ""
  10.504 + end;
  10.505 +
  10.506 +(* String for a matrix numbered k, in SDPA sparse format.                    *)
  10.507 +
  10.508 +fun sdpa_of_matrix k (m:matrix) =
  10.509 + let 
  10.510 +  val pfx = string_of_int k ^ " 1 "
  10.511 +  val ms = Intpairfunc.fold (fn ((i,j), c) => fn  a => if i > j then a else ((i,j),c)::a) (snd m) [] 
  10.512 +  val mss = sort (increasing fst (prod_ord int_ord int_ord)) ms 
  10.513 + in fold_rev (fn ((i,j),c) => fn a =>
  10.514 +     pfx ^ string_of_int i ^ " " ^ string_of_int j ^
  10.515 +     " " ^ decimalize 20 c ^ "\n" ^ a) mss ""
  10.516 + end;;
  10.517 +
  10.518 +(* ------------------------------------------------------------------------- *)
  10.519 +(* String in SDPA sparse format for standard SDP problem:                    *)
  10.520 +(*                                                                           *)
  10.521 +(*    X = v_1 * [M_1] + ... + v_m * [M_m] - [M_0] must be PSD                *)
  10.522 +(*    Minimize obj_1 * v_1 + ... obj_m * v_m                                 *)
  10.523 +(* ------------------------------------------------------------------------- *)
  10.524 +
  10.525 +fun sdpa_of_problem obj mats =
  10.526 + let 
  10.527 +  val m = length mats - 1
  10.528 +  val (n,_) = dimensions (hd mats) 
  10.529 + in
  10.530 +  string_of_int m ^ "\n" ^
  10.531 +  "1\n" ^
  10.532 +  string_of_int n ^ "\n" ^
  10.533 +  sdpa_of_vector obj ^
  10.534 +  fold_rev2 (fn k => fn m => fn a => sdpa_of_matrix (k - 1) m ^ a) (1 upto length mats) mats ""
  10.535 + end;
  10.536 +
  10.537 +fun index_char str chr pos =
  10.538 +  if pos >= String.size str then ~1
  10.539 +  else if String.sub(str,pos) = chr then pos
  10.540 +  else index_char str chr (pos + 1);
  10.541 +fun rat_of_quotient (a,b) = if b = 0 then rat_0 else Rat.rat_of_quotient (a,b);
  10.542 +fun rat_of_string s = 
  10.543 + let val n = index_char s #"/" 0 in
  10.544 +  if n = ~1 then s |> IntInf.fromString |> valOf |> Rat.rat_of_int
  10.545 +  else 
  10.546 +   let val SOME numer = IntInf.fromString(String.substring(s,0,n))
  10.547 +       val SOME den = IntInf.fromString (String.substring(s,n+1,String.size s - n - 1))
  10.548 +   in rat_of_quotient(numer, den)
  10.549 +   end
  10.550 + end;
  10.551 +
  10.552 +fun isspace x = x = " " ;
  10.553 +fun isnum x = x mem_string ["0","1","2","3","4","5","6","7","8","9"]
  10.554 +
  10.555 +(* More parser basics.                                                       *)
  10.556 +
  10.557 +local
  10.558 + open Scan
  10.559 +in 
  10.560 + val word = this_string
  10.561 + fun token s =
  10.562 +  repeat ($$ " ") |-- word s --| repeat ($$ " ")
  10.563 + val numeral = one isnum
  10.564 + val decimalint = bulk numeral >> (rat_of_string o implode)
  10.565 + val decimalfrac = bulk numeral
  10.566 +    >> (fn s => rat_of_string(implode s) // pow10 (length s))
  10.567 + val decimalsig =
  10.568 +    decimalint -- option (Scan.$$ "." |-- decimalfrac)
  10.569 +    >> (fn (h,NONE) => h | (h,SOME x) => h +/ x)
  10.570 + fun signed prs =
  10.571 +       $$ "-" |-- prs >> Rat.neg 
  10.572 +    || $$ "+" |-- prs
  10.573 +    || prs;
  10.574 +
  10.575 +fun emptyin def xs = if null xs then (def,xs) else Scan.fail xs
  10.576 +
  10.577 + val exponent = ($$ "e" || $$ "E") |-- signed decimalint;
  10.578 +
  10.579 + val decimal = signed decimalsig -- (emptyin rat_0|| exponent)
  10.580 +    >> (fn (h, x) => h */ pow10 (int_of_rat x));
  10.581 +end;
  10.582 +
  10.583 + fun mkparser p s =
  10.584 +  let val (x,rst) = p (explode s) 
  10.585 +  in if null rst then x 
  10.586 +     else error "mkparser: unparsed input"
  10.587 +  end;;
  10.588 +
  10.589 +(* Parse back csdp output.                                                      *)
  10.590 +
  10.591 + fun ignore inp = ((),[])
  10.592 + fun csdpoutput inp =  ((decimal -- Scan.bulk (Scan.$$ " " |-- Scan.option decimal) >> (fn (h,to) => map_filter I ((SOME h)::to))) --| ignore >> vector_of_list) inp
  10.593 + val parse_csdpoutput = mkparser csdpoutput
  10.594 +
  10.595 +(* Run prover on a problem in linear form.                       *)
  10.596 +
  10.597 +fun run_problem prover obj mats =
  10.598 +  parse_csdpoutput (prover (sdpa_of_problem obj mats))
  10.599 +
  10.600 +(* Try some apparently sensible scaling first. Note that this is purely to   *)
  10.601 +(* get a cleaner translation to floating-point, and doesn't affect any of    *)
  10.602 +(* the results, in principle. In practice it seems a lot better when there   *)
  10.603 +(* are extreme numbers in the original problem.                              *)
  10.604 +
  10.605 +  (* Version for (int*int) keys *)
  10.606 +local
  10.607 +  fun max_rat x y = if x </ y then y else x
  10.608 +  fun common_denominator fld amat acc =
  10.609 +      fld (fn (m,c) => fn a => lcm_rat (denominator_rat c) a) amat acc
  10.610 +  fun maximal_element fld amat acc =
  10.611 +    fld (fn (m,c) => fn maxa => max_rat maxa (abs_rat c)) amat acc 
  10.612 +fun float_of_rat x = let val (a,b) = Rat.quotient_of_rat x
  10.613 +                     in Real.fromLargeInt a / Real.fromLargeInt b end;
  10.614 +in
  10.615 +
  10.616 +fun pi_scale_then solver (obj:vector)  mats =
  10.617 + let 
  10.618 +  val cd1 = fold_rev (common_denominator Intpairfunc.fold) mats (rat_1)
  10.619 +  val cd2 = common_denominator Intfunc.fold (snd obj)  (rat_1) 
  10.620 +  val mats' = map (Intpairfunc.mapf (fn x => cd1 */ x)) mats
  10.621 +  val obj' = vector_cmul cd2 obj
  10.622 +  val max1 = fold_rev (maximal_element Intpairfunc.fold) mats' (rat_0)
  10.623 +  val max2 = maximal_element Intfunc.fold (snd obj') (rat_0) 
  10.624 +  val scal1 = pow2 (20 - trunc(Math.ln (float_of_rat max1) / Math.ln 2.0))
  10.625 +  val scal2 = pow2 (20 - trunc(Math.ln (float_of_rat max2) / Math.ln 2.0)) 
  10.626 +  val mats'' = map (Intpairfunc.mapf (fn x => x */ scal1)) mats'
  10.627 +  val obj'' = vector_cmul scal2 obj' 
  10.628 + in solver obj'' mats''
  10.629 +  end
  10.630 +end;
  10.631 +
  10.632 +(* Try some apparently sensible scaling first. Note that this is purely to   *)
  10.633 +(* get a cleaner translation to floating-point, and doesn't affect any of    *)
  10.634 +(* the results, in principle. In practice it seems a lot better when there   *)
  10.635 +(* are extreme numbers in the original problem.                              *)
  10.636 +
  10.637 +  (* Version for (int*int*int) keys *)
  10.638 +local
  10.639 +  fun max_rat x y = if x </ y then y else x
  10.640 +  fun common_denominator fld amat acc =
  10.641 +      fld (fn (m,c) => fn a => lcm_rat (denominator_rat c) a) amat acc
  10.642 +  fun maximal_element fld amat acc =
  10.643 +    fld (fn (m,c) => fn maxa => max_rat maxa (abs_rat c)) amat acc 
  10.644 +fun float_of_rat x = let val (a,b) = Rat.quotient_of_rat x
  10.645 +                     in Real.fromLargeInt a / Real.fromLargeInt b end;
  10.646 +fun int_of_float x = (trunc x handle Overflow => 0 | Domain => 0)
  10.647 +in
  10.648 +
  10.649 +fun tri_scale_then solver (obj:vector)  mats =
  10.650 + let 
  10.651 +  val cd1 = fold_rev (common_denominator Inttriplefunc.fold) mats (rat_1)
  10.652 +  val cd2 = common_denominator Intfunc.fold (snd obj)  (rat_1) 
  10.653 +  val mats' = map (Inttriplefunc.mapf (fn x => cd1 */ x)) mats
  10.654 +  val obj' = vector_cmul cd2 obj
  10.655 +  val max1 = fold_rev (maximal_element Inttriplefunc.fold) mats' (rat_0)
  10.656 +  val max2 = maximal_element Intfunc.fold (snd obj') (rat_0) 
  10.657 +  val scal1 = pow2 (20 - int_of_float(Math.ln (float_of_rat max1) / Math.ln 2.0))
  10.658 +  val scal2 = pow2 (20 - int_of_float(Math.ln (float_of_rat max2) / Math.ln 2.0)) 
  10.659 +  val mats'' = map (Inttriplefunc.mapf (fn x => x */ scal1)) mats'
  10.660 +  val obj'' = vector_cmul scal2 obj' 
  10.661 + in solver obj'' mats''
  10.662 +  end
  10.663 +end;
  10.664 +
  10.665 +(* Round a vector to "nice" rationals.                                       *)
  10.666 +
  10.667 +fun nice_rational n x = round_rat (n */ x) // n;;
  10.668 +fun nice_vector n ((d,v) : vector) = 
  10.669 + (d, Intfunc.fold (fn (i,c) => fn a => 
  10.670 +   let val y = nice_rational n c 
  10.671 +   in if c =/ rat_0 then a 
  10.672 +      else Intfunc.update (i,y) a end) v Intfunc.undefined):vector
  10.673 +
  10.674 +fun dest_ord f x = is_equal (f x);
  10.675 +
  10.676 +(* Stuff for "equations" ((int*int*int)->num functions).                         *)
  10.677 +
  10.678 +fun tri_equation_cmul c eq =
  10.679 +  if c =/ rat_0 then Inttriplefunc.undefined else Inttriplefunc.mapf (fn d => c */ d) eq;
  10.680 +
  10.681 +fun tri_equation_add eq1 eq2 = Inttriplefunc.combine (curry op +/) (fn x => x =/ rat_0) eq1 eq2;
  10.682 +
  10.683 +fun tri_equation_eval assig eq =
  10.684 + let fun value v = Inttriplefunc.apply assig v 
  10.685 + in Inttriplefunc.fold (fn (v, c) => fn a => a +/ value v */ c) eq rat_0
  10.686 + end;
  10.687 +
  10.688 +(* Eliminate among linear equations: return unconstrained variables and      *)
  10.689 +(* assignments for the others in terms of them. We give one pseudo-variable  *)
  10.690 +(* "one" that's used for a constant term.                                    *)
  10.691 +
  10.692 +local
  10.693 +  fun extract_first p l = case l of  (* FIXME : use find_first instead *)
  10.694 +   [] => error "extract_first"
  10.695 + | h::t => if p h then (h,t) else
  10.696 +          let val (k,s) = extract_first p t in (k,h::s) end
  10.697 +fun eliminate vars dun eqs = case vars of 
  10.698 +  [] => if forall Inttriplefunc.is_undefined eqs then dun
  10.699 +        else raise Unsolvable
  10.700 + | v::vs =>
  10.701 +  ((let 
  10.702 +    val (eq,oeqs) = extract_first (fn e => Inttriplefunc.defined e v) eqs 
  10.703 +    val a = Inttriplefunc.apply eq v
  10.704 +    val eq' = tri_equation_cmul ((Rat.neg rat_1) // a) (Inttriplefunc.undefine v eq)
  10.705 +    fun elim e =
  10.706 +     let val b = Inttriplefunc.tryapplyd e v rat_0 
  10.707 +     in if b =/ rat_0 then e else
  10.708 +        tri_equation_add e (tri_equation_cmul (Rat.neg b // a) eq)
  10.709 +     end
  10.710 +   in eliminate vs (Inttriplefunc.update (v,eq') (Inttriplefunc.mapf elim dun)) (map elim oeqs)
  10.711 +   end)
  10.712 +  handle Failure _ => eliminate vs dun eqs)
  10.713 +in
  10.714 +fun tri_eliminate_equations one vars eqs =
  10.715 + let 
  10.716 +  val assig = eliminate vars Inttriplefunc.undefined eqs
  10.717 +  val vs = Inttriplefunc.fold (fn (x, f) => fn a => remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig []
  10.718 +  in (distinct (dest_ord triple_int_ord) vs, assig)
  10.719 +  end
  10.720 +end;
  10.721 +
  10.722 +(* Eliminate all variables, in an essentially arbitrary order.               *)
  10.723 +
  10.724 +fun tri_eliminate_all_equations one =
  10.725 + let 
  10.726 +  fun choose_variable eq =
  10.727 +   let val (v,_) = Inttriplefunc.choose eq 
  10.728 +   in if is_equal (triple_int_ord(v,one)) then
  10.729 +      let val eq' = Inttriplefunc.undefine v eq 
  10.730 +      in if Inttriplefunc.is_undefined eq' then error "choose_variable" 
  10.731 +         else fst (Inttriplefunc.choose eq')
  10.732 +      end
  10.733 +    else v 
  10.734 +   end
  10.735 +  fun eliminate dun eqs = case eqs of 
  10.736 +    [] => dun
  10.737 +  | eq::oeqs =>
  10.738 +    if Inttriplefunc.is_undefined eq then eliminate dun oeqs else
  10.739 +    let val v = choose_variable eq
  10.740 +        val a = Inttriplefunc.apply eq v
  10.741 +        val eq' = tri_equation_cmul ((Rat.rat_of_int ~1) // a) 
  10.742 +                   (Inttriplefunc.undefine v eq)
  10.743 +        fun elim e =
  10.744 +         let val b = Inttriplefunc.tryapplyd e v rat_0 
  10.745 +         in if b =/ rat_0 then e 
  10.746 +            else tri_equation_add e (tri_equation_cmul (Rat.neg b // a) eq)
  10.747 +         end
  10.748 +    in eliminate (Inttriplefunc.update(v, eq') (Inttriplefunc.mapf elim dun)) 
  10.749 +                 (map elim oeqs) 
  10.750 +    end
  10.751 +in fn eqs =>
  10.752 + let 
  10.753 +  val assig = eliminate Inttriplefunc.undefined eqs
  10.754 +  val vs = Inttriplefunc.fold (fn (x, f) => fn a => remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig []
  10.755 + in (distinct (dest_ord triple_int_ord) vs,assig)
  10.756 + end
  10.757 +end;
  10.758 + 
  10.759 +(* Solve equations by assigning arbitrary numbers.                           *)
  10.760 +
  10.761 +fun tri_solve_equations one eqs =
  10.762 + let 
  10.763 +  val (vars,assigs) = tri_eliminate_all_equations one eqs
  10.764 +  val vfn = fold_rev (fn v => Inttriplefunc.update(v,rat_0)) vars 
  10.765 +            (Inttriplefunc.onefunc(one, Rat.rat_of_int ~1))
  10.766 +  val ass =
  10.767 +    Inttriplefunc.combine (curry op +/) (K false) 
  10.768 +    (Inttriplefunc.mapf (tri_equation_eval vfn) assigs) vfn 
  10.769 + in if forall (fn e => tri_equation_eval ass e =/ rat_0) eqs
  10.770 +    then Inttriplefunc.undefine one ass else raise Sanity
  10.771 + end;
  10.772 +
  10.773 +(* Multiply equation-parametrized poly by regular poly and add accumulator.  *)
  10.774 +
  10.775 +fun tri_epoly_pmul p q acc =
  10.776 + Monomialfunc.fold (fn (m1, c) => fn a =>
  10.777 +  Monomialfunc.fold (fn (m2,e) => fn b =>
  10.778 +   let val m =  monomial_mul m1 m2
  10.779 +       val es = Monomialfunc.tryapplyd b m Inttriplefunc.undefined 
  10.780 +   in Monomialfunc.update (m,tri_equation_add (tri_equation_cmul c e) es) b 
  10.781 +   end) q a) p acc ;
  10.782 +
  10.783 +(* Usual operations on equation-parametrized poly.                           *)
  10.784 +
  10.785 +fun tri_epoly_cmul c l =
  10.786 +  if c =/ rat_0 then Inttriplefunc.undefined else Inttriplefunc.mapf (tri_equation_cmul c) l;;
  10.787 +
  10.788 +val tri_epoly_neg = tri_epoly_cmul (Rat.rat_of_int ~1);
  10.789 +
  10.790 +val tri_epoly_add = Inttriplefunc.combine tri_equation_add Inttriplefunc.is_undefined;
  10.791 +
  10.792 +fun tri_epoly_sub p q = tri_epoly_add p (tri_epoly_neg q);;
  10.793 +
  10.794 +(* Stuff for "equations" ((int*int)->num functions).                         *)
  10.795 +
  10.796 +fun pi_equation_cmul c eq =
  10.797 +  if c =/ rat_0 then Inttriplefunc.undefined else Inttriplefunc.mapf (fn d => c */ d) eq;
  10.798 +
  10.799 +fun pi_equation_add eq1 eq2 = Inttriplefunc.combine (curry op +/) (fn x => x =/ rat_0) eq1 eq2;
  10.800 +
  10.801 +fun pi_equation_eval assig eq =
  10.802 + let fun value v = Inttriplefunc.apply assig v 
  10.803 + in Inttriplefunc.fold (fn (v, c) => fn a => a +/ value v */ c) eq rat_0
  10.804 + end;
  10.805 +
  10.806 +(* Eliminate among linear equations: return unconstrained variables and      *)
  10.807 +(* assignments for the others in terms of them. We give one pseudo-variable  *)
  10.808 +(* "one" that's used for a constant term.                                    *)
  10.809 +
  10.810 +local
  10.811 +fun extract_first p l = case l of 
  10.812 +   [] => error "extract_first"
  10.813 + | h::t => if p h then (h,t) else
  10.814 +          let val (k,s) = extract_first p t in (k,h::s) end
  10.815 +fun eliminate vars dun eqs = case vars of 
  10.816 +  [] => if forall Inttriplefunc.is_undefined eqs then dun
  10.817 +        else raise Unsolvable
  10.818 + | v::vs =>
  10.819 +   let 
  10.820 +    val (eq,oeqs) = extract_first (fn e => Inttriplefunc.defined e v) eqs 
  10.821 +    val a = Inttriplefunc.apply eq v
  10.822 +    val eq' = pi_equation_cmul ((Rat.neg rat_1) // a) (Inttriplefunc.undefine v eq)
  10.823 +    fun elim e =
  10.824 +     let val b = Inttriplefunc.tryapplyd e v rat_0 
  10.825 +     in if b =/ rat_0 then e else
  10.826 +        pi_equation_add e (pi_equation_cmul (Rat.neg b // a) eq)
  10.827 +     end
  10.828 +   in eliminate vs (Inttriplefunc.update (v,eq') (Inttriplefunc.mapf elim dun)) (map elim oeqs)
  10.829 +   end
  10.830 +  handle Failure _ => eliminate vs dun eqs
  10.831 +in
  10.832 +fun pi_eliminate_equations one vars eqs =
  10.833 + let 
  10.834 +  val assig = eliminate vars Inttriplefunc.undefined eqs
  10.835 +  val vs = Inttriplefunc.fold (fn (x, f) => fn a => remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig []
  10.836 +  in (distinct (dest_ord triple_int_ord) vs, assig)
  10.837 +  end
  10.838 +end;
  10.839 +
  10.840 +(* Eliminate all variables, in an essentially arbitrary order.               *)
  10.841 +
  10.842 +fun pi_eliminate_all_equations one =
  10.843 + let 
  10.844 +  fun choose_variable eq =
  10.845 +   let val (v,_) = Inttriplefunc.choose eq 
  10.846 +   in if is_equal (triple_int_ord(v,one)) then
  10.847 +      let val eq' = Inttriplefunc.undefine v eq 
  10.848 +      in if Inttriplefunc.is_undefined eq' then error "choose_variable" 
  10.849 +         else fst (Inttriplefunc.choose eq')
  10.850 +      end
  10.851 +    else v 
  10.852 +   end
  10.853 +  fun eliminate dun eqs = case eqs of 
  10.854 +    [] => dun
  10.855 +  | eq::oeqs =>
  10.856 +    if Inttriplefunc.is_undefined eq then eliminate dun oeqs else
  10.857 +    let val v = choose_variable eq
  10.858 +        val a = Inttriplefunc.apply eq v
  10.859 +        val eq' = pi_equation_cmul ((Rat.rat_of_int ~1) // a) 
  10.860 +                   (Inttriplefunc.undefine v eq)
  10.861 +        fun elim e =
  10.862 +         let val b = Inttriplefunc.tryapplyd e v rat_0 
  10.863 +         in if b =/ rat_0 then e 
  10.864 +            else pi_equation_add e (pi_equation_cmul (Rat.neg b // a) eq)
  10.865 +         end
  10.866 +    in eliminate (Inttriplefunc.update(v, eq') (Inttriplefunc.mapf elim dun)) 
  10.867 +                 (map elim oeqs) 
  10.868 +    end
  10.869 +in fn eqs =>
  10.870 + let 
  10.871 +  val assig = eliminate Inttriplefunc.undefined eqs
  10.872 +  val vs = Inttriplefunc.fold (fn (x, f) => fn a => remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig []
  10.873 + in (distinct (dest_ord triple_int_ord) vs,assig)
  10.874 + end
  10.875 +end;
  10.876 + 
  10.877 +(* Solve equations by assigning arbitrary numbers.                           *)
  10.878 +
  10.879 +fun pi_solve_equations one eqs =
  10.880 + let 
  10.881 +  val (vars,assigs) = pi_eliminate_all_equations one eqs
  10.882 +  val vfn = fold_rev (fn v => Inttriplefunc.update(v,rat_0)) vars 
  10.883 +            (Inttriplefunc.onefunc(one, Rat.rat_of_int ~1))
  10.884 +  val ass =
  10.885 +    Inttriplefunc.combine (curry op +/) (K false) 
  10.886 +    (Inttriplefunc.mapf (pi_equation_eval vfn) assigs) vfn 
  10.887 + in if forall (fn e => pi_equation_eval ass e =/ rat_0) eqs
  10.888 +    then Inttriplefunc.undefine one ass else raise Sanity
  10.889 + end;
  10.890 +
  10.891 +(* Multiply equation-parametrized poly by regular poly and add accumulator.  *)
  10.892 +
  10.893 +fun pi_epoly_pmul p q acc =
  10.894 + Monomialfunc.fold (fn (m1, c) => fn a =>
  10.895 +  Monomialfunc.fold (fn (m2,e) => fn b =>
  10.896 +   let val m =  monomial_mul m1 m2
  10.897 +       val es = Monomialfunc.tryapplyd b m Inttriplefunc.undefined 
  10.898 +   in Monomialfunc.update (m,pi_equation_add (pi_equation_cmul c e) es) b 
  10.899 +   end) q a) p acc ;
  10.900 +
  10.901 +(* Usual operations on equation-parametrized poly.                           *)
  10.902 +
  10.903 +fun pi_epoly_cmul c l =
  10.904 +  if c =/ rat_0 then Inttriplefunc.undefined else Inttriplefunc.mapf (pi_equation_cmul c) l;;
  10.905 +
  10.906 +val pi_epoly_neg = pi_epoly_cmul (Rat.rat_of_int ~1);
  10.907 +
  10.908 +val pi_epoly_add = Inttriplefunc.combine pi_equation_add Inttriplefunc.is_undefined;
  10.909 +
  10.910 +fun pi_epoly_sub p q = pi_epoly_add p (pi_epoly_neg q);;
  10.911 +
  10.912 +fun allpairs f l1 l2 =  fold_rev (fn x => (curry (op @)) (map (f x) l2)) l1 [];
  10.913 +
  10.914 +(* Hence produce the "relevant" monomials: those whose squares lie in the    *)
  10.915 +(* Newton polytope of the monomials in the input. (This is enough according  *)
  10.916 +(* to Reznik: "Extremal PSD forms with few terms", Duke Math. Journal,       *)
  10.917 +(* vol 45, pp. 363--374, 1978.                                               *)
  10.918 +(*                                                                           *)
  10.919 +(* These are ordered in sort of decreasing degree. In particular the         *)
  10.920 +(* constant monomial is last; this gives an order in diagonalization of the  *)
  10.921 +(* quadratic form that will tend to display constants.                       *)
  10.922 +
  10.923 +(* Diagonalize (Cholesky/LDU) the matrix corresponding to a quadratic form.  *)
  10.924 +
  10.925 +local
  10.926 +fun diagonalize n i m =
  10.927 + if Intpairfunc.is_undefined (snd m) then [] 
  10.928 + else
  10.929 +  let val a11 = Intpairfunc.tryapplyd (snd m) (i,i) rat_0 
  10.930 +  in if a11 </ rat_0 then raise Failure "diagonalize: not PSD"
  10.931 +    else if a11 =/ rat_0 then
  10.932 +          if Intfunc.is_undefined (snd (row i m)) then diagonalize n (i + 1) m
  10.933 +          else raise Failure "diagonalize: not PSD ___ "
  10.934 +    else
  10.935 +     let 
  10.936 +      val v = row i m
  10.937 +      val v' = (fst v, Intfunc.fold (fn (i, c) => fn a => 
  10.938 +       let val y = c // a11 
  10.939 +       in if y = rat_0 then a else Intfunc.update (i,y) a 
  10.940 +       end)  (snd v) Intfunc.undefined)
  10.941 +      fun upt0 x y a = if y = rat_0 then a else Intpairfunc.update (x,y) a
  10.942 +      val m' =
  10.943 +      ((n,n),
  10.944 +      iter (i+1,n) (fn j =>
  10.945 +          iter (i+1,n) (fn k =>
  10.946 +              (upt0 (j,k) (Intpairfunc.tryapplyd (snd m) (j,k) rat_0 -/ Intfunc.tryapplyd (snd v) j rat_0 */ Intfunc.tryapplyd (snd v') k rat_0))))
  10.947 +          Intpairfunc.undefined)
  10.948 +     in (a11,v')::diagonalize n (i + 1) m' 
  10.949 +     end
  10.950 +  end
  10.951 +in
  10.952 +fun diag m =
  10.953 + let 
  10.954 +   val nn = dimensions m 
  10.955 +   val n = fst nn 
  10.956 + in if snd nn <> n then error "diagonalize: non-square matrix" 
  10.957 +    else diagonalize n 1 m
  10.958 + end
  10.959 +end;
  10.960 +
  10.961 +fun gcd_rat a b = Rat.rat_of_int (Integer.gcd (int_of_rat a) (int_of_rat b));
  10.962 +
  10.963 +(* Adjust a diagonalization to collect rationals at the start.               *)
  10.964 +  (* FIXME : Potentially polymorphic keys, but here only: integers!! *)
  10.965 +local
  10.966 + fun upd0 x y a = if y =/ rat_0 then a else Intfunc.update(x,y) a;
  10.967 + fun mapa f (d,v) = 
  10.968 +  (d, Intfunc.fold (fn (i,c) => fn a => upd0 i (f c) a) v Intfunc.undefined)
  10.969 + fun adj (c,l) =
  10.970 + let val a = 
  10.971 +  Intfunc.fold (fn (i,c) => fn a => lcm_rat a (denominator_rat c)) 
  10.972 +    (snd l) rat_1 //
  10.973 +  Intfunc.fold (fn (i,c) => fn a => gcd_rat a (numerator_rat c)) 
  10.974 +    (snd l) rat_0
  10.975 +  in ((c // (a */ a)),mapa (fn x => a */ x) l)
  10.976 +  end
  10.977 +in
  10.978 +fun deration d = if null d then (rat_0,d) else
  10.979 + let val d' = map adj d
  10.980 +     val a = fold (lcm_rat o denominator_rat o fst) d' rat_1 //
  10.981 +          fold (gcd_rat o numerator_rat o fst) d' rat_0 
  10.982 + in ((rat_1 // a),map (fn (c,l) => (a */ c,l)) d')
  10.983 + end
  10.984 +end;
  10.985 + 
  10.986 +(* Enumeration of monomials with given multidegree bound.                    *)
  10.987 +
  10.988 +fun enumerate_monomials d vars = 
  10.989 + if d < 0 then []
  10.990 + else if d = 0 then [Ctermfunc.undefined]
  10.991 + else if null vars then [monomial_1] else
  10.992 + let val alts =
  10.993 +  map (fn k => let val oths = enumerate_monomials (d - k) (tl vars) 
  10.994 +               in map (fn ks => if k = 0 then ks else Ctermfunc.update (hd vars, k) ks) oths end) (0 upto d) 
  10.995 + in fold1 (curry op @) alts
  10.996 + end;
  10.997 +
  10.998 +(* Enumerate products of distinct input polys with degree <= d.              *)
  10.999 +(* We ignore any constant input polynomials.                                 *)
 10.1000 +(* Give the output polynomial and a record of how it was derived.            *)
 10.1001 +
 10.1002 +local
 10.1003 + open RealArith
 10.1004 +in
 10.1005 +fun enumerate_products d pols =
 10.1006 +if d = 0 then [(poly_const rat_1,Rational_lt rat_1)] 
 10.1007 +else if d < 0 then [] else
 10.1008 +case pols of 
 10.1009 +   [] => [(poly_const rat_1,Rational_lt rat_1)]
 10.1010 + | (p,b)::ps => 
 10.1011 +    let val e = multidegree p 
 10.1012 +    in if e = 0 then enumerate_products d ps else
 10.1013 +       enumerate_products d ps @
 10.1014 +       map (fn (q,c) => (poly_mul p q,Product(b,c)))
 10.1015 +         (enumerate_products (d - e) ps)
 10.1016 +    end
 10.1017 +end;
 10.1018 +
 10.1019 +(* Convert regular polynomial. Note that we treat (0,0,0) as -1.             *)
 10.1020 +
 10.1021 +fun epoly_of_poly p =
 10.1022 +  Monomialfunc.fold (fn (m,c) => fn a => Monomialfunc.update (m, Inttriplefunc.onefunc ((0,0,0), Rat.neg c)) a) p Monomialfunc.undefined;
 10.1023 +
 10.1024 +(* String for block diagonal matrix numbered k.                              *)
 10.1025 +
 10.1026 +fun sdpa_of_blockdiagonal k m =
 10.1027 + let 
 10.1028 +  val pfx = string_of_int k ^" "
 10.1029 +  val ents =
 10.1030 +    Inttriplefunc.fold 
 10.1031 +      (fn ((b,i,j),c) => fn a => if i > j then a else ((b,i,j),c)::a) 
 10.1032 +      m [] 
 10.1033 +  val entss = sort (increasing fst triple_int_ord) ents 
 10.1034 + in fold_rev (fn ((b,i,j),c) => fn a =>
 10.1035 +     pfx ^ string_of_int b ^ " " ^ string_of_int i ^ " " ^ string_of_int j ^
 10.1036 +     " " ^ decimalize 20 c ^ "\n" ^ a) entss ""
 10.1037 + end;
 10.1038 +
 10.1039 +(* SDPA for problem using block diagonal (i.e. multiple SDPs)                *)
 10.1040 +
 10.1041 +fun sdpa_of_blockproblem nblocks blocksizes obj mats =
 10.1042 + let val m = length mats - 1 
 10.1043 + in
 10.1044 +  string_of_int m ^ "\n" ^
 10.1045 +  string_of_int nblocks ^ "\n" ^
 10.1046 +  (fold1 (fn s => fn t => s^" "^t) (map string_of_int blocksizes)) ^
 10.1047 +  "\n" ^
 10.1048 +  sdpa_of_vector obj ^
 10.1049 +  fold_rev2 (fn k => fn m => fn a => sdpa_of_blockdiagonal (k - 1) m ^ a)
 10.1050 +    (1 upto length mats) mats ""
 10.1051 + end;
 10.1052 +
 10.1053 +(* Run prover on a problem in block diagonal form.                       *)
 10.1054 +
 10.1055 +fun run_blockproblem prover nblocks blocksizes obj mats=
 10.1056 +  parse_csdpoutput (prover (sdpa_of_blockproblem nblocks blocksizes obj mats))
 10.1057 +
 10.1058 +(* 3D versions of matrix operations to consider blocks separately.           *)
 10.1059 +
 10.1060 +val bmatrix_add = Inttriplefunc.combine (curry op +/) (fn x => x =/ rat_0);
 10.1061 +fun bmatrix_cmul c bm =
 10.1062 +  if c =/ rat_0 then Inttriplefunc.undefined
 10.1063 +  else Inttriplefunc.mapf (fn x => c */ x) bm;
 10.1064 +
 10.1065 +val bmatrix_neg = bmatrix_cmul (Rat.rat_of_int ~1);
 10.1066 +fun bmatrix_sub m1 m2 = bmatrix_add m1 (bmatrix_neg m2);;
 10.1067 +
 10.1068 +(* Smash a block matrix into components.                                     *)
 10.1069 +
 10.1070 +fun blocks blocksizes bm =
 10.1071 + map (fn (bs,b0) =>
 10.1072 +      let val m = Inttriplefunc.fold
 10.1073 +          (fn ((b,i,j),c) => fn a => if b = b0 then Intpairfunc.update ((i,j),c) a else a) bm Intpairfunc.undefined
 10.1074 +          val d = Intpairfunc.fold (fn ((i,j),c) => fn a => max a (max i j)) m 0 
 10.1075 +      in (((bs,bs),m):matrix) end)
 10.1076 + (blocksizes ~~ (1 upto length blocksizes));;
 10.1077 +
 10.1078 +(* FIXME : Get rid of this !!!*)
 10.1079 +local
 10.1080 +  fun tryfind_with msg f [] = raise Failure msg
 10.1081 +    | tryfind_with msg f (x::xs) = (f x handle Failure s => tryfind_with s f xs);
 10.1082 +in 
 10.1083 +  fun tryfind f = tryfind_with "tryfind" f
 10.1084 +end
 10.1085 +
 10.1086 +(*
 10.1087 +fun tryfind f [] = error "tryfind"
 10.1088 +  | tryfind f (x::xs) = (f x handle ERROR _ => tryfind f xs);
 10.1089 +*)
 10.1090 +
 10.1091 +(* Positiv- and Nullstellensatz. Flag "linf" forces a linear representation. *)
 10.1092 +
 10.1093 + 
 10.1094 +local
 10.1095 + open RealArith
 10.1096 +in
 10.1097 +fun real_positivnullstellensatz_general prover linf d eqs leqs pol =
 10.1098 +let 
 10.1099 + val vars = fold_rev (curry (gen_union (op aconvc)) o poly_variables) 
 10.1100 +              (pol::eqs @ map fst leqs) []
 10.1101 + val monoid = if linf then 
 10.1102 +      (poly_const rat_1,Rational_lt rat_1)::
 10.1103 +      (filter (fn (p,c) => multidegree p <= d) leqs)
 10.1104 +    else enumerate_products d leqs
 10.1105 + val nblocks = length monoid
 10.1106 + fun mk_idmultiplier k p =
 10.1107 +  let 
 10.1108 +   val e = d - multidegree p
 10.1109 +   val mons = enumerate_monomials e vars
 10.1110 +   val nons = mons ~~ (1 upto length mons) 
 10.1111 +  in (mons,
 10.1112 +      fold_rev (fn (m,n) => Monomialfunc.update(m,Inttriplefunc.onefunc((~k,~n,n),rat_1))) nons Monomialfunc.undefined)
 10.1113 +  end
 10.1114 +
 10.1115 + fun mk_sqmultiplier k (p,c) =
 10.1116 +  let 
 10.1117 +   val e = (d - multidegree p) div 2
 10.1118 +   val mons = enumerate_monomials e vars
 10.1119 +   val nons = mons ~~ (1 upto length mons) 
 10.1120 +  in (mons, 
 10.1121 +      fold_rev (fn (m1,n1) =>
 10.1122 +       fold_rev (fn (m2,n2) => fn  a =>
 10.1123 +        let val m = monomial_mul m1 m2 
 10.1124 +        in if n1 > n2 then a else
 10.1125 +          let val c = if n1 = n2 then rat_1 else rat_2
 10.1126 +              val e = Monomialfunc.tryapplyd a m Inttriplefunc.undefined 
 10.1127 +          in Monomialfunc.update(m, tri_equation_add (Inttriplefunc.onefunc((k,n1,n2), c)) e) a
 10.1128 +          end
 10.1129 +        end)  nons)
 10.1130 +       nons Monomialfunc.undefined)
 10.1131 +  end
 10.1132 +
 10.1133 +  val (sqmonlist,sqs) = split_list (map2 mk_sqmultiplier (1 upto length monoid) monoid)
 10.1134 +  val (idmonlist,ids) =  split_list(map2 mk_idmultiplier (1 upto length eqs) eqs)
 10.1135 +  val blocksizes = map length sqmonlist
 10.1136 +  val bigsum =
 10.1137 +    fold_rev2 (fn p => fn q => fn a => tri_epoly_pmul p q a) eqs ids
 10.1138 +            (fold_rev2 (fn (p,c) => fn s => fn a => tri_epoly_pmul p s a) monoid sqs
 10.1139 +                     (epoly_of_poly(poly_neg pol)))
 10.1140 +  val eqns = Monomialfunc.fold (fn (m,e) => fn a => e::a) bigsum []
 10.1141 +  val (pvs,assig) = tri_eliminate_all_equations (0,0,0) eqns
 10.1142 +  val qvars = (0,0,0)::pvs
 10.1143 +  val allassig = fold_rev (fn v => Inttriplefunc.update(v,(Inttriplefunc.onefunc(v,rat_1)))) pvs assig
 10.1144 +  fun mk_matrix v =
 10.1145 +    Inttriplefunc.fold (fn ((b,i,j), ass) => fn m => 
 10.1146 +        if b < 0 then m else
 10.1147 +         let val c = Inttriplefunc.tryapplyd ass v rat_0
 10.1148 +         in if c = rat_0 then m else
 10.1149 +            Inttriplefunc.update ((b,j,i), c) (Inttriplefunc.update ((b,i,j), c) m)
 10.1150 +         end)
 10.1151 +          allassig Inttriplefunc.undefined
 10.1152 +  val diagents = Inttriplefunc.fold
 10.1153 +    (fn ((b,i,j), e) => fn a => if b > 0 andalso i = j then tri_equation_add e a else a)
 10.1154 +    allassig Inttriplefunc.undefined
 10.1155 +
 10.1156 +  val mats = map mk_matrix qvars
 10.1157 +  val obj = (length pvs,
 10.1158 +            itern 1 pvs (fn v => fn i => Intfunc.updatep iszero (i,Inttriplefunc.tryapplyd diagents v rat_0))
 10.1159 +                        Intfunc.undefined)
 10.1160 +  val raw_vec = if null pvs then vector_0 0
 10.1161 +                else tri_scale_then (run_blockproblem prover nblocks blocksizes) obj mats
 10.1162 +  fun int_element (d,v) i = Intfunc.tryapplyd v i rat_0
 10.1163 +  fun cterm_element (d,v) i = Ctermfunc.tryapplyd v i rat_0
 10.1164 +
 10.1165 +  fun find_rounding d =
 10.1166 +   let 
 10.1167 +    val _ = if !debugging 
 10.1168 +          then writeln ("Trying rounding with limit "^Rat.string_of_rat d ^ "\n") 
 10.1169 +          else ()
 10.1170 +    val vec = nice_vector d raw_vec
 10.1171 +    val blockmat = iter (1,dim vec)
 10.1172 +     (fn i => fn a => bmatrix_add (bmatrix_cmul (int_element vec i) (nth mats i)) a)
 10.1173 +     (bmatrix_neg (nth mats 0))
 10.1174 +    val allmats = blocks blocksizes blockmat 
 10.1175 +   in (vec,map diag allmats)
 10.1176 +   end
 10.1177 +  val (vec,ratdias) =
 10.1178 +    if null pvs then find_rounding rat_1
 10.1179 +    else tryfind find_rounding (map Rat.rat_of_int (1 upto 31) @
 10.1180 +                                map pow2 (5 upto 66))
 10.1181 +  val newassigs =
 10.1182 +    fold_rev (fn k => Inttriplefunc.update (nth pvs (k - 1), int_element vec k))
 10.1183 +           (1 upto dim vec) (Inttriplefunc.onefunc ((0,0,0), Rat.rat_of_int ~1))
 10.1184 +  val finalassigs =
 10.1185 +    Inttriplefunc.fold (fn (v,e) => fn a => Inttriplefunc.update(v, tri_equation_eval newassigs e) a) allassig newassigs
 10.1186 +  fun poly_of_epoly p =
 10.1187 +    Monomialfunc.fold (fn (v,e) => fn a => Monomialfunc.updatep iszero (v,tri_equation_eval finalassigs e) a)
 10.1188 +          p Monomialfunc.undefined
 10.1189 +  fun  mk_sos mons =
 10.1190 +   let fun mk_sq (c,m) =
 10.1191 +    (c,fold_rev (fn k=> fn a => Monomialfunc.updatep iszero (nth mons (k - 1), int_element m k) a)
 10.1192 +                 (1 upto length mons) Monomialfunc.undefined)
 10.1193 +   in map mk_sq
 10.1194 +   end
 10.1195 +  val sqs = map2 mk_sos sqmonlist ratdias
 10.1196 +  val cfs = map poly_of_epoly ids
 10.1197 +  val msq = filter (fn (a,b) => not (null b)) (map2 pair monoid sqs)
 10.1198 +  fun eval_sq sqs = fold_rev (fn (c,q) => poly_add (poly_cmul c (poly_mul q q))) sqs poly_0
 10.1199 +  val sanity =
 10.1200 +    fold_rev (fn ((p,c),s) => poly_add (poly_mul p (eval_sq s))) msq
 10.1201 +           (fold_rev2 (fn p => fn q => poly_add (poly_mul p q)) cfs eqs
 10.1202 +                    (poly_neg pol))
 10.1203 +
 10.1204 +in if not(Monomialfunc.is_undefined sanity) then raise Sanity else
 10.1205 +  (cfs,map (fn (a,b) => (snd a,b)) msq)
 10.1206 + end
 10.1207 +
 10.1208 +
 10.1209 +end;
 10.1210 +
 10.1211 +(* Iterative deepening.                                                      *)
 10.1212 +
 10.1213 +fun deepen f n = 
 10.1214 +  (writeln ("Searching with depth limit " ^ string_of_int n) ; (f n handle Failure s => (writeln ("failed with message: " ^ s) ; deepen f (n+1))))
 10.1215 +
 10.1216 +(* The ordering so we can create canonical HOL polynomials.                  *)
 10.1217 +
 10.1218 +fun dest_monomial mon = sort (increasing fst cterm_ord) (Ctermfunc.graph mon);
 10.1219 +
 10.1220 +fun monomial_order (m1,m2) =
 10.1221 + if Ctermfunc.is_undefined m2 then LESS 
 10.1222 + else if Ctermfunc.is_undefined m1 then GREATER 
 10.1223 + else
 10.1224 +  let val mon1 = dest_monomial m1 
 10.1225 +      val mon2 = dest_monomial m2
 10.1226 +      val deg1 = fold (curry op + o snd) mon1 0
 10.1227 +      val deg2 = fold (curry op + o snd) mon2 0 
 10.1228 +  in if deg1 < deg2 then GREATER else if deg1 > deg2 then LESS
 10.1229 +     else list_ord (prod_ord cterm_ord int_ord) (mon1,mon2)
 10.1230 +  end;
 10.1231 +
 10.1232 +fun dest_poly p =
 10.1233 +  map (fn (m,c) => (c,dest_monomial m))
 10.1234 +      (sort (prod_ord monomial_order (K EQUAL)) (Monomialfunc.graph p));
 10.1235 +
 10.1236 +(* Map back polynomials and their composites to HOL.                         *)
 10.1237 +
 10.1238 +local
 10.1239 + open Thm Numeral RealArith
 10.1240 +in
 10.1241 +
 10.1242 +fun cterm_of_varpow x k = if k = 1 then x else capply (capply @{cterm "op ^ :: real => _"} x) 
 10.1243 +  (mk_cnumber @{ctyp nat} k)
 10.1244 +
 10.1245 +fun cterm_of_monomial m = 
 10.1246 + if Ctermfunc.is_undefined m then @{cterm "1::real"} 
 10.1247 + else 
 10.1248 +  let 
 10.1249 +   val m' = dest_monomial m
 10.1250 +   val vps = fold_rev (fn (x,k) => cons (cterm_of_varpow x k)) m' [] 
 10.1251 +  in fold1 (fn s => fn t => capply (capply @{cterm "op * :: real => _"} s) t) vps
 10.1252 +  end
 10.1253 +
 10.1254 +fun cterm_of_cmonomial (m,c) = if Ctermfunc.is_undefined m then cterm_of_rat c
 10.1255 +    else if c = Rat.one then cterm_of_monomial m
 10.1256 +    else capply (capply @{cterm "op *::real => _"} (cterm_of_rat c)) (cterm_of_monomial m);
 10.1257 +
 10.1258 +fun cterm_of_poly p = 
 10.1259 + if Monomialfunc.is_undefined p then @{cterm "0::real"} 
 10.1260 + else
 10.1261 +  let 
 10.1262 +   val cms = map cterm_of_cmonomial
 10.1263 +     (sort (prod_ord monomial_order (K EQUAL)) (Monomialfunc.graph p))
 10.1264 +  in fold1 (fn t1 => fn t2 => capply(capply @{cterm "op + :: real => _"} t1) t2) cms
 10.1265 +  end;
 10.1266 +
 10.1267 +fun cterm_of_sqterm (c,p) = Product(Rational_lt c,Square(cterm_of_poly p));
 10.1268 +
 10.1269 +fun cterm_of_sos (pr,sqs) = if null sqs then pr
 10.1270 +  else Product(pr,fold1 (fn a => fn b => Sum(a,b)) (map cterm_of_sqterm sqs));
 10.1271 +
 10.1272 +end
 10.1273 +
 10.1274 +(* Interface to HOL.                                                         *)
 10.1275 +local
 10.1276 +  open Thm Conv RealArith
 10.1277 +  val concl = dest_arg o cprop_of
 10.1278 +  fun simple_cterm_ord t u = TermOrd.fast_term_ord (term_of t, term_of u) = LESS
 10.1279 +in
 10.1280 +  (* FIXME: Replace tryfind by get_first !! *)
 10.1281 +fun real_nonlinear_prover prover ctxt =
 10.1282 + let 
 10.1283 +  val {add,mul,neg,pow,sub,main} =  Normalizer.semiring_normalizers_ord_wrapper ctxt
 10.1284 +      (valOf (NormalizerData.match ctxt @{cterm "(0::real) + 1"})) 
 10.1285 +     simple_cterm_ord
 10.1286 +  val (real_poly_add_conv,real_poly_mul_conv,real_poly_neg_conv,
 10.1287 +       real_poly_pow_conv,real_poly_sub_conv,real_poly_conv) = (add,mul,neg,pow,sub,main)
 10.1288 +  fun mainf  translator (eqs,les,lts) = 
 10.1289 +  let 
 10.1290 +   val eq0 = map (poly_of_term o dest_arg1 o concl) eqs
 10.1291 +   val le0 = map (poly_of_term o dest_arg o concl) les
 10.1292 +   val lt0 = map (poly_of_term o dest_arg o concl) lts
 10.1293 +   val eqp0 = map (fn (t,i) => (t,Axiom_eq i)) (eq0 ~~ (0 upto (length eq0 - 1)))
 10.1294 +   val lep0 = map (fn (t,i) => (t,Axiom_le i)) (le0 ~~ (0 upto (length le0 - 1)))
 10.1295 +   val ltp0 = map (fn (t,i) => (t,Axiom_lt i)) (lt0 ~~ (0 upto (length lt0 - 1)))
 10.1296 +   val (keq,eq) = List.partition (fn (p,_) => multidegree p = 0) eqp0
 10.1297 +   val (klep,lep) = List.partition (fn (p,_) => multidegree p = 0) lep0
 10.1298 +   val (kltp,ltp) = List.partition (fn (p,_) => multidegree p = 0) ltp0
 10.1299 +   fun trivial_axiom (p,ax) =
 10.1300 +    case ax of
 10.1301 +       Axiom_eq n => if eval Ctermfunc.undefined p <>/ Rat.zero then nth eqs n 
 10.1302 +                     else raise Failure "trivial_axiom: Not a trivial axiom"
 10.1303 +     | Axiom_le n => if eval Ctermfunc.undefined p </ Rat.zero then nth les n 
 10.1304 +                     else raise Failure "trivial_axiom: Not a trivial axiom"
 10.1305 +     | Axiom_lt n => if eval Ctermfunc.undefined p <=/ Rat.zero then nth lts n 
 10.1306 +                     else raise Failure "trivial_axiom: Not a trivial axiom"
 10.1307 +     | _ => error "trivial_axiom: Not a trivial axiom"
 10.1308 +   in 
 10.1309 +  ((let val th = tryfind trivial_axiom (keq @ klep @ kltp)
 10.1310 +   in fconv_rule (arg_conv (arg1_conv real_poly_conv) then_conv field_comp_conv) th end)
 10.1311 +   handle Failure _ => (
 10.1312 +    let 
 10.1313 +     val pol = fold_rev poly_mul (map fst ltp) (poly_const Rat.one)
 10.1314 +     val leq = lep @ ltp
 10.1315 +     fun tryall d =
 10.1316 +      let val e = multidegree pol
 10.1317 +          val k = if e = 0 then 0 else d div e
 10.1318 +          val eq' = map fst eq 
 10.1319 +      in tryfind (fn i => (d,i,real_positivnullstellensatz_general prover false d eq' leq
 10.1320 +                            (poly_neg(poly_pow pol i))))
 10.1321 +              (0 upto k)
 10.1322 +      end
 10.1323 +    val (d,i,(cert_ideal,cert_cone)) = deepen tryall 0
 10.1324 +    val proofs_ideal =
 10.1325 +      map2 (fn q => fn (p,ax) => Eqmul(cterm_of_poly q,ax)) cert_ideal eq
 10.1326 +    val proofs_cone = map cterm_of_sos cert_cone
 10.1327 +    val proof_ne = if null ltp then Rational_lt Rat.one else
 10.1328 +      let val p = fold1 (fn s => fn t => Product(s,t)) (map snd ltp) 
 10.1329 +      in  funpow i (fn q => Product(p,q)) (Rational_lt Rat.one)
 10.1330 +      end
 10.1331 +    val proof = fold1 (fn s => fn t => Sum(s,t))
 10.1332 +                           (proof_ne :: proofs_ideal @ proofs_cone) 
 10.1333 +    in writeln "Translating proof certificate to HOL";
 10.1334 +       translator (eqs,les,lts) proof
 10.1335 +    end))
 10.1336 +   end
 10.1337 + in mainf end
 10.1338 +end
 10.1339 +
 10.1340 +fun C f x y = f y x;
 10.1341 +  (* FIXME : This is very bad!!!*)
 10.1342 +fun subst_conv eqs t = 
 10.1343 + let 
 10.1344 +  val t' = fold (Thm.cabs o Thm.lhs_of) eqs t
 10.1345 + in Conv.fconv_rule (Thm.beta_conversion true) (fold (C combination) eqs (reflexive t'))
 10.1346 + end
 10.1347 +
 10.1348 +(* A wrapper that tries to substitute away variables first.                  *)
 10.1349 +
 10.1350 +local
 10.1351 + open Thm Conv RealArith
 10.1352 +  fun simple_cterm_ord t u = TermOrd.fast_term_ord (term_of t, term_of u) = LESS
 10.1353 + val concl = dest_arg o cprop_of
 10.1354 + val shuffle1 = 
 10.1355 +   fconv_rule (rewr_conv @{lemma "(a + x == y) == (x == y - (a::real))" by (atomize (full)) (simp add: ring_simps) })
 10.1356 + val shuffle2 =
 10.1357 +    fconv_rule (rewr_conv @{lemma "(x + a == y) ==  (x == y - (a::real))" by (atomize (full)) (simp add: ring_simps)})
 10.1358 + fun substitutable_monomial fvs tm = case term_of tm of
 10.1359 +    Free(_,@{typ real}) => if not (member (op aconvc) fvs tm) then (Rat.one,tm) 
 10.1360 +                           else raise Failure "substitutable_monomial"
 10.1361 +  | @{term "op * :: real => _"}$c$(t as Free _ ) => 
 10.1362 +     if is_ratconst (dest_arg1 tm) andalso not (member (op aconvc) fvs (dest_arg tm))
 10.1363 +         then (dest_ratconst (dest_arg1 tm),dest_arg tm) else raise Failure "substitutable_monomial"
 10.1364 +  | @{term "op + :: real => _"}$s$t => 
 10.1365 +       (substitutable_monomial (add_cterm_frees (dest_arg tm) fvs) (dest_arg1 tm)
 10.1366 +        handle Failure _ => substitutable_monomial (add_cterm_frees (dest_arg1 tm) fvs) (dest_arg tm))
 10.1367 +  | _ => raise Failure "substitutable_monomial"
 10.1368 +
 10.1369 +  fun isolate_variable v th = 
 10.1370 +   let val w = dest_arg1 (cprop_of th)
 10.1371 +   in if v aconvc w then th
 10.1372 +      else case term_of w of
 10.1373 +           @{term "op + :: real => _"}$s$t => 
 10.1374 +              if dest_arg1 w aconvc v then shuffle2 th 
 10.1375 +              else isolate_variable v (shuffle1 th)
 10.1376 +          | _ => error "isolate variable : This should not happen?"
 10.1377 +   end 
 10.1378 +in
 10.1379 +
 10.1380 +fun real_nonlinear_subst_prover prover ctxt =
 10.1381 + let 
 10.1382 +  val {add,mul,neg,pow,sub,main} =  Normalizer.semiring_normalizers_ord_wrapper ctxt
 10.1383 +      (valOf (NormalizerData.match ctxt @{cterm "(0::real) + 1"})) 
 10.1384 +     simple_cterm_ord
 10.1385 +
 10.1386 +  val (real_poly_add_conv,real_poly_mul_conv,real_poly_neg_conv,
 10.1387 +       real_poly_pow_conv,real_poly_sub_conv,real_poly_conv) = (add,mul,neg,pow,sub,main)
 10.1388 +
 10.1389 +  fun make_substitution th =
 10.1390 +   let 
 10.1391 +    val (c,v) = substitutable_monomial [] (dest_arg1(concl th))
 10.1392 +    val th1 = Drule.arg_cong_rule (capply @{cterm "op * :: real => _"} (cterm_of_rat (Rat.inv c))) (mk_meta_eq th)
 10.1393 +    val th2 = fconv_rule (binop_conv real_poly_mul_conv) th1
 10.1394 +   in fconv_rule (arg_conv real_poly_conv) (isolate_variable v th2)
 10.1395 +   end
 10.1396 +   fun oprconv cv ct = 
 10.1397 +    let val g = Thm.dest_fun2 ct
 10.1398 +    in if g aconvc @{cterm "op <= :: real => _"} 
 10.1399 +         orelse g aconvc @{cterm "op < :: real => _"} 
 10.1400 +       then arg_conv cv ct else arg1_conv cv ct
 10.1401 +    end
 10.1402 +  fun mainf translator =
 10.1403 +   let 
 10.1404 +    fun substfirst(eqs,les,lts) =
 10.1405 +      ((let 
 10.1406 +           val eth = tryfind make_substitution eqs
 10.1407 +           val modify = fconv_rule (arg_conv (oprconv(subst_conv [eth] then_conv real_poly_conv)))
 10.1408 +       in  substfirst
 10.1409 +             (filter_out (fn t => (Thm.dest_arg1 o Thm.dest_arg o cprop_of) t 
 10.1410 +                                   aconvc @{cterm "0::real"}) (map modify eqs),
 10.1411 +                                   map modify les,map modify lts)
 10.1412 +       end)
 10.1413 +       handle Failure  _ => real_nonlinear_prover prover ctxt translator (rev eqs, rev les, rev lts))
 10.1414 +    in substfirst
 10.1415 +   end
 10.1416 +
 10.1417 +
 10.1418 + in mainf
 10.1419 + end
 10.1420 +
 10.1421 +(* Overall function. *)
 10.1422 +
 10.1423 +fun real_sos prover ctxt t = gen_prover_real_arith ctxt (real_nonlinear_subst_prover prover ctxt) t;
 10.1424 +end;
 10.1425 +
 10.1426 +(* A tactic *)
 10.1427 +fun strip_all ct = 
 10.1428 + case term_of ct of 
 10.1429 +  Const("all",_) $ Abs (xn,xT,p) => 
 10.1430 +   let val (a,(v,t')) = (apsnd (Thm.dest_abs (SOME xn)) o Thm.dest_comb) ct
 10.1431 +   in apfst (cons v) (strip_all t')
 10.1432 +   end
 10.1433 +| _ => ([],ct)
 10.1434 +
 10.1435 +fun core_sos_conv prover ctxt t = Drule.arg_cong_rule @{cterm Trueprop} (real_sos prover ctxt (Thm.dest_arg t) RS @{thm Eq_TrueI})
 10.1436 +
 10.1437 +val known_sos_constants = 
 10.1438 +  [@{term "op ==>"}, @{term "Trueprop"}, 
 10.1439 +   @{term "op -->"}, @{term "op &"}, @{term "op |"}, 
 10.1440 +   @{term "Not"}, @{term "op = :: bool => _"}, 
 10.1441 +   @{term "All :: (real => _) => _"}, @{term "Ex :: (real => _) => _"}, 
 10.1442 +   @{term "op = :: real => _"}, @{term "op < :: real => _"}, 
 10.1443 +   @{term "op <= :: real => _"}, 
 10.1444 +   @{term "op + :: real => _"}, @{term "op - :: real => _"}, 
 10.1445 +   @{term "op * :: real => _"}, @{term "uminus :: real => _"}, 
 10.1446 +   @{term "op / :: real => _"}, @{term "inverse :: real => _"},
 10.1447 +   @{term "op ^ :: real => _"}, @{term "abs :: real => _"}, 
 10.1448 +   @{term "min :: real => _"}, @{term "max :: real => _"},
 10.1449 +   @{term "0::real"}, @{term "1::real"}, @{term "number_of :: int => real"},
 10.1450 +   @{term "number_of :: int => nat"},
 10.1451 +   @{term "Int.Bit0"}, @{term "Int.Bit1"}, 
 10.1452 +   @{term "Int.Pls"}, @{term "Int.Min"}];
 10.1453 +
 10.1454 +fun check_sos kcts ct = 
 10.1455 + let
 10.1456 +  val t = term_of ct
 10.1457 +  val _ = if not (null (Term.add_tfrees t []) 
 10.1458 +                  andalso null (Term.add_tvars t [])) 
 10.1459 +          then error "SOS: not sos. Additional type varables" else ()
 10.1460 +  val fs = Term.add_frees t []
 10.1461 +  val _ = if exists (fn ((_,T)) => not (T = @{typ "real"})) fs 
 10.1462 +          then error "SOS: not sos. Variables with type not real" else ()
 10.1463 +  val vs = Term.add_vars t []
 10.1464 +  val _ = if exists (fn ((_,T)) => not (T = @{typ "real"})) fs 
 10.1465 +          then error "SOS: not sos. Variables with type not real" else ()
 10.1466 +  val ukcs = subtract (fn (t,p) => Const p aconv t) kcts (Term.add_consts t [])
 10.1467 +  val _ = if  null ukcs then () 
 10.1468 +              else error ("SOSO: Unknown constants in Subgoal:" ^ commas (map fst ukcs))
 10.1469 +in () end
 10.1470 +
 10.1471 +fun core_sos_tac prover ctxt = CSUBGOAL (fn (ct, i) => 
 10.1472 +  let val _ = check_sos known_sos_constants ct
 10.1473 +      val (avs, p) = strip_all ct
 10.1474 +      val th = standard (fold_rev forall_intr avs (real_sos prover ctxt (Thm.dest_arg p)))
 10.1475 +  in rtac th i end);
 10.1476 +
 10.1477 +fun default_SOME f NONE v = SOME v
 10.1478 +  | default_SOME f (SOME v) _ = SOME v;
 10.1479 +
 10.1480 +fun lift_SOME f NONE a = f a
 10.1481 +  | lift_SOME f (SOME a) _ = SOME a;
 10.1482 +
 10.1483 +
 10.1484 +local
 10.1485 + val is_numeral = can (HOLogic.dest_number o term_of)
 10.1486 +in
 10.1487 +fun get_denom b ct = case term_of ct of
 10.1488 +  @{term "op / :: real => _"} $ _ $ _ => 
 10.1489 +     if is_numeral (Thm.dest_arg ct) then get_denom b (Thm.dest_arg1 ct)
 10.1490 +     else default_SOME (get_denom b) (get_denom b (Thm.dest_arg ct))   (Thm.dest_arg ct, b)
 10.1491 + | @{term "op < :: real => _"} $ _ $ _ => lift_SOME (get_denom true) (get_denom true (Thm.dest_arg ct)) (Thm.dest_arg1 ct)
 10.1492 + | @{term "op <= :: real => _"} $ _ $ _ => lift_SOME (get_denom true) (get_denom true (Thm.dest_arg ct)) (Thm.dest_arg1 ct)
 10.1493 + | _ $ _ => lift_SOME (get_denom b) (get_denom b (Thm.dest_fun ct)) (Thm.dest_arg ct)
 10.1494 + | _ => NONE
 10.1495 +end;
 10.1496 +
 10.1497 +fun elim_one_denom_tac ctxt = 
 10.1498 +CSUBGOAL (fn (P,i) => 
 10.1499 + case get_denom false P of 
 10.1500 +   NONE => no_tac
 10.1501 + | SOME (d,ord) => 
 10.1502 +     let 
 10.1503 +      val ss = simpset_of ctxt addsimps @{thms field_simps} 
 10.1504 +               addsimps [@{thm nonzero_power_divide}, @{thm power_divide}]
 10.1505 +      val th = instantiate' [] [SOME d, SOME (Thm.dest_arg P)] 
 10.1506 +         (if ord then @{lemma "(d=0 --> P) & (d>0 --> P) & (d<(0::real) --> P) ==> P" by auto}
 10.1507 +          else @{lemma "(d=0 --> P) & (d ~= (0::real) --> P) ==> P" by blast})
 10.1508 +     in (rtac th i THEN Simplifier.asm_full_simp_tac ss i) end);
 10.1509 +
 10.1510 +fun elim_denom_tac ctxt i = REPEAT (elim_one_denom_tac ctxt i);
 10.1511 +
 10.1512 +fun sos_tac prover ctxt = ObjectLogic.full_atomize_tac THEN' elim_denom_tac ctxt THEN' core_sos_tac prover ctxt
 10.1513 +
 10.1514 +
 10.1515 +end;
    11.1 --- a/src/HOL/Library/sos_wrapper.ML	Wed Aug 05 17:10:10 2009 +0200
    11.2 +++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
    11.3 @@ -1,158 +0,0 @@
    11.4 -(* Title:      sos_wrapper.ML
    11.5 -   Author:     Philipp Meyer, TU Muenchen
    11.6 -
    11.7 -Added functionality for sums of squares, e.g. calling a remote prover
    11.8 -*)
    11.9 -
   11.10 -signature SOS_WRAPPER =
   11.11 -sig
   11.12 -
   11.13 -  datatype prover_result = Success | PartialSuccess | Failure | Error
   11.14 -  type prover = string * (int -> string -> prover_result * string)
   11.15 -
   11.16 -  val setup: theory -> theory
   11.17 -end
   11.18 -
   11.19 -structure SosWrapper : SOS_WRAPPER=
   11.20 -struct
   11.21 -
   11.22 -datatype prover_result = Success | PartialSuccess | Failure | Error
   11.23 -type prover =
   11.24 -  string *                           (* command name *)
   11.25 -  (int -> string ->prover_result * string)   (* function to find failure from return value and output *)
   11.26 -
   11.27 -
   11.28 -(*** output control ***)
   11.29 -
   11.30 -fun debug s = Output.debug (fn () => s)
   11.31 -val answer = Output.priority
   11.32 -val write = Output.writeln
   11.33 -
   11.34 -(*** calling provers ***)
   11.35 -
   11.36 -val destdir = ref ""
   11.37 -
   11.38 -fun filename dir name =
   11.39 -  let
   11.40 -    val probfile = Path.basic (name ^ serial_string ())
   11.41 -  in
   11.42 -    if dir = "" then
   11.43 -      File.tmp_path probfile
   11.44 -    else
   11.45 -      if File.exists (Path.explode dir) then
   11.46 -        Path.append  (Path.explode dir) probfile
   11.47 -      else
   11.48 -        error ("No such directory: " ^ dir)
   11.49 -  end
   11.50 -
   11.51 -fun is_success Success = true
   11.52 -  | is_success PartialSuccess = true
   11.53 -  | is_success _ = false
   11.54 -fun str_of_status Success = "Success"
   11.55 -  | str_of_status PartialSuccess = "Partial Success"
   11.56 -  | str_of_status Failure= "Failure"
   11.57 -  | str_of_status Error= "Error"
   11.58 -
   11.59 -fun run_solver name (cmd, find_failure) input =
   11.60 -  let
   11.61 -    val _ = answer ("Calling solver: " ^ name)
   11.62 -
   11.63 -    (* create input file *)
   11.64 -    val dir = ! destdir
   11.65 -    val input_file = filename dir "sos_in" 
   11.66 -    val _ = File.write input_file input
   11.67 -
   11.68 -    val _ = debug "Solver input:"
   11.69 -    val _ = debug input
   11.70 -
   11.71 -    (* call solver *)
   11.72 -    val output_file = filename dir "sos_out"
   11.73 -    val (output, rv) = system_out (cmd ^ " " ^ (Path.implode input_file) ^
   11.74 -                         " " ^ (Path.implode output_file))
   11.75 - 
   11.76 -    (* read and analysize output *)
   11.77 -    val (res, res_msg) = find_failure rv output
   11.78 -    val result = if is_success res then File.read output_file else ""
   11.79 -
   11.80 -    (* remove temporary files *)
   11.81 -    val _ = if dir = "" then (File.rm input_file ; if File.exists output_file then File.rm output_file else ()) else ()
   11.82 -
   11.83 -    val _ = debug "Solver output:"
   11.84 -    val _ = debug output
   11.85 -    val _ = debug "Solver result:"
   11.86 -    val _ = debug result
   11.87 -
   11.88 -    val _ = answer (str_of_status res ^ ": " ^ res_msg)
   11.89 -
   11.90 -  in
   11.91 -    if is_success res then
   11.92 -      result
   11.93 -    else
   11.94 -      error ("Prover failed: " ^ res_msg)
   11.95 -  end
   11.96 -
   11.97 -(*** various provers ***)
   11.98 -
   11.99 -(* local csdp client *)
  11.100 -
  11.101 -fun find_csdp_run_failure rv _ =
  11.102 -  case rv of
  11.103 -    0 => (Success, "SDP solved")
  11.104 -  | 1 => (Failure, "SDP is primal infeasible")
  11.105 -  | 2 => (Failure, "SDP is dual infeasible")
  11.106 -  | 3 => (PartialSuccess, "SDP solved with reduced accuracy")
  11.107 -  | _ => (Failure, "return code is " ^ string_of_int rv)
  11.108 -
  11.109 -val csdp = ("csdp", find_csdp_run_failure)
  11.110 -
  11.111 -(* remote neos server *)
  11.112 -
  11.113 -fun find_neos_failure rv output =
  11.114 -  if rv = 2 then (Failure, "no solution") else
  11.115 -  if rv <> 0 then (Error, "return code is " ^ string_of_int rv) else
  11.116 -  let
  11.117 -    fun find_success str =
  11.118 -      if String.isPrefix "Success: " str then
  11.119 -        SOME (Success, unprefix "Success: " str)
  11.120 -      else if String.isPrefix "Partial Success: " str then
  11.121 -        SOME (PartialSuccess, unprefix "Partial Success: " str)
  11.122 -      else if String.isPrefix "Failure: " str then
  11.123 -        SOME (Failure, unprefix "Failure: " str)
  11.124 -      else
  11.125 -        NONE 
  11.126 -    val exit_line = get_first find_success (split_lines output)
  11.127 -  in
  11.128 -    case exit_line of
  11.129 -      SOME (status, msg) =>
  11.130 -        if String.isPrefix "SDP solved" msg then
  11.131 -          (status, msg)
  11.132 -        else (Failure, msg)
  11.133 -    | NONE => (Failure, "no success")
  11.134 -  end
  11.135 -
  11.136 -val neos_csdp = ("$ISABELLE_HOME/lib/scripts/neos/NeosCSDPClient.py", find_neos_failure)
  11.137 -
  11.138 -(* save provers in table *)
  11.139 -
  11.140 -val provers =
  11.141 -     Symtab.make [("remote_csdp", neos_csdp),("csdp", csdp)]
  11.142 -
  11.143 -fun get_prover name =
  11.144 -  case Symtab.lookup provers name of
  11.145 -    SOME prover => prover
  11.146 -  | NONE => error ("unknown prover: " ^ name)
  11.147 -
  11.148 -fun call_solver name =
  11.149 -    run_solver name (get_prover name)
  11.150 -
  11.151 -(* setup tactic *)
  11.152 -
  11.153 -val def_solver = "remote_csdp"
  11.154 -
  11.155 -fun sos_solver name = (SIMPLE_METHOD' o (Sos.sos_tac (call_solver name))) 
  11.156 -
  11.157 -val sos_method = Scan.optional (Scan.lift OuterParse.xname) def_solver >> sos_solver
  11.158 -
  11.159 -val setup = Method.setup @{binding sos} sos_method "Prove universal problems over the reals using sums of squares"
  11.160 -
  11.161 -end
    12.1 --- a/src/HOL/Library/sum_of_squares.ML	Wed Aug 05 17:10:10 2009 +0200
    12.2 +++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
    12.3 @@ -1,1754 +0,0 @@
    12.4 -(* Title:      sum_of_squares.ML
    12.5 -   Authors:    Amine Chaieb, University of Cambridge
    12.6 -               Philipp Meyer, TU Muenchen
    12.7 -
    12.8 -A tactic for proving nonlinear inequalities
    12.9 -*)
   12.10 -
   12.11 -signature SOS =
   12.12 -sig
   12.13 -
   12.14 -  val sos_tac : (string -> string) -> Proof.context -> int -> Tactical.tactic
   12.15 -
   12.16 -end
   12.17 -
   12.18 -structure Sos : SOS = 
   12.19 -struct
   12.20 -
   12.21 -
   12.22 -val rat_0 = Rat.zero;
   12.23 -val rat_1 = Rat.one;
   12.24 -val rat_2 = Rat.two;
   12.25 -val rat_10 = Rat.rat_of_int 10;
   12.26 -val rat_1_2 = rat_1 // rat_2;
   12.27 -val max = curry IntInf.max;
   12.28 -val min = curry IntInf.min;
   12.29 -
   12.30 -val denominator_rat = Rat.quotient_of_rat #> snd #> Rat.rat_of_int;
   12.31 -val numerator_rat = Rat.quotient_of_rat #> fst #> Rat.rat_of_int;
   12.32 -fun int_of_rat a = 
   12.33 -    case Rat.quotient_of_rat a of (i,1) => i | _ => error "int_of_rat: not an int";
   12.34 -fun lcm_rat x y = Rat.rat_of_int (Integer.lcm (int_of_rat x) (int_of_rat y));
   12.35 -
   12.36 -fun rat_pow r i = 
   12.37 - let fun pow r i = 
   12.38 -   if i = 0 then rat_1 else 
   12.39 -   let val d = pow r (i div 2)
   12.40 -   in d */ d */ (if i mod 2 = 0 then rat_1 else r)
   12.41 -   end
   12.42 - in if i < 0 then pow (Rat.inv r) (~ i) else pow r i end;
   12.43 -
   12.44 -fun round_rat r = 
   12.45 - let val (a,b) = Rat.quotient_of_rat (Rat.abs r)
   12.46 -     val d = a div b
   12.47 -     val s = if r </ rat_0 then (Rat.neg o Rat.rat_of_int) else Rat.rat_of_int
   12.48 -     val x2 = 2 * (a - (b * d))
   12.49 - in s (if x2 >= b then d + 1 else d) end
   12.50 -
   12.51 -val abs_rat = Rat.abs;
   12.52 -val pow2 = rat_pow rat_2;
   12.53 -val pow10 = rat_pow rat_10;
   12.54 -
   12.55 -val debugging = ref false;
   12.56 -
   12.57 -exception Sanity;
   12.58 -
   12.59 -exception Unsolvable;
   12.60 -
   12.61 -(* Turn a rational into a decimal string with d sig digits.                  *)
   12.62 -
   12.63 -local
   12.64 -fun normalize y =
   12.65 -  if abs_rat y </ (rat_1 // rat_10) then normalize (rat_10 */ y) - 1
   12.66 -  else if abs_rat y >=/ rat_1 then normalize (y // rat_10) + 1
   12.67 -  else 0 
   12.68 - in
   12.69 -fun decimalize d x =
   12.70 -  if x =/ rat_0 then "0.0" else
   12.71 -  let 
   12.72 -   val y = Rat.abs x
   12.73 -   val e = normalize y
   12.74 -   val z = pow10(~ e) */ y +/ rat_1
   12.75 -   val k = int_of_rat (round_rat(pow10 d */ z))
   12.76 -  in (if x </ rat_0 then "-0." else "0.") ^
   12.77 -     implode(tl(explode(string_of_int k))) ^
   12.78 -     (if e = 0 then "" else "e"^string_of_int e)
   12.79 -  end
   12.80 -end;
   12.81 -
   12.82 -(* Iterations over numbers, and lists indexed by numbers.                    *)
   12.83 -
   12.84 -fun itern k l f a =
   12.85 -  case l of
   12.86 -    [] => a
   12.87 -  | h::t => itern (k + 1) t f (f h k a);
   12.88 -
   12.89 -fun iter (m,n) f a =
   12.90 -  if n < m then a
   12.91 -  else iter (m+1,n) f (f m a);
   12.92 -
   12.93 -(* The main types.                                                           *)
   12.94 -
   12.95 -fun strict_ord ord (x,y) = case ord (x,y) of LESS => LESS | _ => GREATER
   12.96 -
   12.97 -structure Intpairfunc = FuncFun(type key = int*int val ord = prod_ord int_ord int_ord);
   12.98 -
   12.99 -type vector = int* Rat.rat Intfunc.T;
  12.100 -
  12.101 -type matrix = (int*int)*(Rat.rat Intpairfunc.T);
  12.102 -
  12.103 -type monomial = int Ctermfunc.T;
  12.104 -
  12.105 -val cterm_ord = (fn (s,t) => TermOrd.fast_term_ord(term_of s, term_of t))
  12.106 - fun monomial_ord (m1,m2) = list_ord (prod_ord cterm_ord int_ord) (Ctermfunc.graph m1, Ctermfunc.graph m2)
  12.107 -structure Monomialfunc = FuncFun(type key = monomial val ord = monomial_ord)
  12.108 -
  12.109 -type poly = Rat.rat Monomialfunc.T;
  12.110 -
  12.111 - fun iszero (k,r) = r =/ rat_0;
  12.112 -
  12.113 -fun fold_rev2 f l1 l2 b =
  12.114 -  case (l1,l2) of
  12.115 -    ([],[]) => b
  12.116 -  | (h1::t1,h2::t2) => f h1 h2 (fold_rev2 f t1 t2 b)
  12.117 -  | _ => error "fold_rev2";
  12.118 - 
  12.119 -(* Vectors. Conventionally indexed 1..n.                                     *)
  12.120 -
  12.121 -fun vector_0 n = (n,Intfunc.undefined):vector;
  12.122 -
  12.123 -fun dim (v:vector) = fst v;
  12.124 -
  12.125 -fun vector_const c n =
  12.126 -  if c =/ rat_0 then vector_0 n
  12.127 -  else (n,fold_rev (fn k => Intfunc.update (k,c)) (1 upto n) Intfunc.undefined) :vector;
  12.128 -
  12.129 -val vector_1 = vector_const rat_1;
  12.130 -
  12.131 -fun vector_cmul c (v:vector) =
  12.132 - let val n = dim v 
  12.133 - in if c =/ rat_0 then vector_0 n
  12.134 -    else (n,Intfunc.mapf (fn x => c */ x) (snd v))
  12.135 - end;
  12.136 -
  12.137 -fun vector_neg (v:vector) = (fst v,Intfunc.mapf Rat.neg (snd v)) :vector;
  12.138 -
  12.139 -fun vector_add (v1:vector) (v2:vector) =
  12.140 - let val m = dim v1  
  12.141 -     val n = dim v2 
  12.142 - in if m <> n then error "vector_add: incompatible dimensions"
  12.143 -    else (n,Intfunc.combine (curry op +/) (fn x => x =/ rat_0) (snd v1) (snd v2)) :vector 
  12.144 - end;
  12.145 -
  12.146 -fun vector_sub v1 v2 = vector_add v1 (vector_neg v2);
  12.147 -
  12.148 -fun vector_dot (v1:vector) (v2:vector) =
  12.149 - let val m = dim v1 
  12.150 -     val n = dim v2 
  12.151 - in if m <> n then error "vector_dot: incompatible dimensions" 
  12.152 -    else Intfunc.fold (fn (i,x) => fn a => x +/ a) 
  12.153 -        (Intfunc.combine (curry op */) (fn x => x =/ rat_0) (snd v1) (snd v2)) rat_0
  12.154 - end;
  12.155 -
  12.156 -fun vector_of_list l =
  12.157 - let val n = length l 
  12.158 - in (n,fold_rev2 (curry Intfunc.update) (1 upto n) l Intfunc.undefined) :vector
  12.159 - end;
  12.160 -
  12.161 -(* Matrices; again rows and columns indexed from 1.                          *)
  12.162 -
  12.163 -fun matrix_0 (m,n) = ((m,n),Intpairfunc.undefined):matrix;
  12.164 -
  12.165 -fun dimensions (m:matrix) = fst m;
  12.166 -
  12.167 -fun matrix_const c (mn as (m,n)) =
  12.168 -  if m <> n then error "matrix_const: needs to be square"
  12.169 -  else if c =/ rat_0 then matrix_0 mn
  12.170 -  else (mn,fold_rev (fn k => Intpairfunc.update ((k,k), c)) (1 upto n) Intpairfunc.undefined) :matrix;;
  12.171 -
  12.172 -val matrix_1 = matrix_const rat_1;
  12.173 -
  12.174 -fun matrix_cmul c (m:matrix) =
  12.175 - let val (i,j) = dimensions m 
  12.176 - in if c =/ rat_0 then matrix_0 (i,j)
  12.177 -    else ((i,j),Intpairfunc.mapf (fn x => c */ x) (snd m))
  12.178 - end;
  12.179 -
  12.180 -fun matrix_neg (m:matrix) = 
  12.181 -  (dimensions m, Intpairfunc.mapf Rat.neg (snd m)) :matrix;
  12.182 -
  12.183 -fun matrix_add (m1:matrix) (m2:matrix) =
  12.184 - let val d1 = dimensions m1 
  12.185 -     val d2 = dimensions m2 
  12.186 - in if d1 <> d2 
  12.187 -     then error "matrix_add: incompatible dimensions"
  12.188 -    else (d1,Intpairfunc.combine (curry op +/) (fn x => x =/ rat_0) (snd m1) (snd m2)) :matrix
  12.189 - end;;
  12.190 -
  12.191 -fun matrix_sub m1 m2 = matrix_add m1 (matrix_neg m2);
  12.192 -
  12.193 -fun row k (m:matrix) =
  12.194 - let val (i,j) = dimensions m 
  12.195 - in (j,
  12.196 -   Intpairfunc.fold (fn ((i,j), c) => fn a => if i = k then Intfunc.update (j,c) a else a) (snd m) Intfunc.undefined ) : vector
  12.197 - end;
  12.198 -
  12.199 -fun column k (m:matrix) =
  12.200 -  let val (i,j) = dimensions m 
  12.201 -  in (i,
  12.202 -   Intpairfunc.fold (fn ((i,j), c) => fn a => if j = k then Intfunc.update (i,c) a else a) (snd m)  Intfunc.undefined)
  12.203 -   : vector
  12.204 - end;
  12.205 -
  12.206 -fun transp (m:matrix) =
  12.207 -  let val (i,j) = dimensions m 
  12.208 -  in
  12.209 -  ((j,i),Intpairfunc.fold (fn ((i,j), c) => fn a => Intpairfunc.update ((j,i), c) a) (snd m) Intpairfunc.undefined) :matrix
  12.210 - end;
  12.211 -
  12.212 -fun diagonal (v:vector) =
  12.213 - let val n = dim v 
  12.214 - in ((n,n),Intfunc.fold (fn (i, c) => fn a => Intpairfunc.update ((i,i), c) a) (snd v) Intpairfunc.undefined) : matrix
  12.215 - end;
  12.216 -
  12.217 -fun matrix_of_list l =
  12.218 - let val m = length l 
  12.219 - in if m = 0 then matrix_0 (0,0) else
  12.220 -   let val n = length (hd l) 
  12.221 -   in ((m,n),itern 1 l (fn v => fn i => itern 1 v (fn c => fn j => Intpairfunc.update ((i,j), c))) Intpairfunc.undefined)
  12.222 -   end
  12.223 - end;
  12.224 -
  12.225 -(* Monomials.                                                                *)
  12.226 -
  12.227 -fun monomial_eval assig (m:monomial) =
  12.228 -  Ctermfunc.fold (fn (x, k) => fn a => a */ rat_pow (Ctermfunc.apply assig x) k)
  12.229 -        m rat_1;
  12.230 -val monomial_1 = (Ctermfunc.undefined:monomial);
  12.231 -
  12.232 -fun monomial_var x = Ctermfunc.onefunc (x, 1) :monomial;
  12.233 -
  12.234 -val (monomial_mul:monomial->monomial->monomial) =
  12.235 -  Ctermfunc.combine (curry op +) (K false);
  12.236 -
  12.237 -fun monomial_pow (m:monomial) k =
  12.238 -  if k = 0 then monomial_1
  12.239 -  else Ctermfunc.mapf (fn x => k * x) m;
  12.240 -
  12.241 -fun monomial_divides (m1:monomial) (m2:monomial) =
  12.242 -  Ctermfunc.fold (fn (x, k) => fn a => Ctermfunc.tryapplyd m2 x 0 >= k andalso a) m1 true;;
  12.243 -
  12.244 -fun monomial_div (m1:monomial) (m2:monomial) =
  12.245 - let val m = Ctermfunc.combine (curry op +) 
  12.246 -   (fn x => x = 0) m1 (Ctermfunc.mapf (fn x => ~ x) m2) 
  12.247 - in if Ctermfunc.fold (fn (x, k) => fn a => k >= 0 andalso a) m true then m
  12.248 -  else error "monomial_div: non-divisible"
  12.249 - end;
  12.250 -
  12.251 -fun monomial_degree x (m:monomial) = 
  12.252 -  Ctermfunc.tryapplyd m x 0;;
  12.253 -
  12.254 -fun monomial_lcm (m1:monomial) (m2:monomial) =
  12.255 -  fold_rev (fn x => Ctermfunc.update (x, max (monomial_degree x m1) (monomial_degree x m2)))
  12.256 -          (gen_union (is_equal o  cterm_ord) (Ctermfunc.dom m1, Ctermfunc.dom m2)) (Ctermfunc.undefined :monomial);
  12.257 -
  12.258 -fun monomial_multidegree (m:monomial) = 
  12.259 - Ctermfunc.fold (fn (x, k) => fn a => k + a) m 0;;
  12.260 -
  12.261 -fun monomial_variables m = Ctermfunc.dom m;;
  12.262 -
  12.263 -(* Polynomials.                                                              *)
  12.264 -
  12.265 -fun eval assig (p:poly) =
  12.266 -  Monomialfunc.fold (fn (m, c) => fn a => a +/ c */ monomial_eval assig m) p rat_0;
  12.267 -
  12.268 -val poly_0 = (Monomialfunc.undefined:poly);
  12.269 -
  12.270 -fun poly_isconst (p:poly) = 
  12.271 -  Monomialfunc.fold (fn (m, c) => fn a => Ctermfunc.is_undefined m andalso a) p true;
  12.272 -
  12.273 -fun poly_var x = Monomialfunc.onefunc (monomial_var x,rat_1) :poly;
  12.274 -
  12.275 -fun poly_const c =
  12.276 -  if c =/ rat_0 then poly_0 else Monomialfunc.onefunc(monomial_1, c);
  12.277 -
  12.278 -fun poly_cmul c (p:poly) =
  12.279 -  if c =/ rat_0 then poly_0
  12.280 -  else Monomialfunc.mapf (fn x => c */ x) p;
  12.281 -
  12.282 -fun poly_neg (p:poly) = (Monomialfunc.mapf Rat.neg p :poly);;
  12.283 -
  12.284 -fun poly_add (p1:poly) (p2:poly) =
  12.285 -  (Monomialfunc.combine (curry op +/) (fn x => x =/ rat_0) p1 p2 :poly);
  12.286 -
  12.287 -fun poly_sub p1 p2 = poly_add p1 (poly_neg p2);
  12.288 -
  12.289 -fun poly_cmmul (c,m) (p:poly) =
  12.290 - if c =/ rat_0 then poly_0
  12.291 - else if Ctermfunc.is_undefined m 
  12.292 -      then Monomialfunc.mapf (fn d => c */ d) p
  12.293 -      else Monomialfunc.fold (fn (m', d) => fn a => (Monomialfunc.update (monomial_mul m m', c */ d) a)) p poly_0;
  12.294 -
  12.295 -fun poly_mul (p1:poly) (p2:poly) =
  12.296 -  Monomialfunc.fold (fn (m, c) => fn a => poly_add (poly_cmmul (c,m) p2) a) p1 poly_0;
  12.297 -
  12.298 -fun poly_div (p1:poly) (p2:poly) =
  12.299 - if not(poly_isconst p2) 
  12.300 - then error "poly_div: non-constant" else
  12.301 - let val c = eval Ctermfunc.undefined p2 
  12.302 - in if c =/ rat_0 then error "poly_div: division by zero"
  12.303 -    else poly_cmul (Rat.inv c) p1
  12.304 - end;
  12.305 -
  12.306 -fun poly_square p = poly_mul p p;
  12.307 -
  12.308 -fun poly_pow p k =
  12.309 - if k = 0 then poly_const rat_1
  12.310 - else if k = 1 then p
  12.311 - else let val q = poly_square(poly_pow p (k div 2)) in
  12.312 -      if k mod 2 = 1 then poly_mul p q else q end;
  12.313 -
  12.314 -fun poly_exp p1 p2 =
  12.315 -  if not(poly_isconst p2) 
  12.316 -  then error "poly_exp: not a constant" 
  12.317 -  else poly_pow p1 (int_of_rat (eval Ctermfunc.undefined p2));
  12.318 -
  12.319 -fun degree x (p:poly) = 
  12.320 - Monomialfunc.fold (fn (m,c) => fn a => max (monomial_degree x m) a) p 0;
  12.321 -
  12.322 -fun multidegree (p:poly) =
  12.323 -  Monomialfunc.fold (fn (m, c) => fn a => max (monomial_multidegree m) a) p 0;
  12.324 -
  12.325 -fun poly_variables (p:poly) =
  12.326 -  sort cterm_ord (Monomialfunc.fold_rev (fn (m, c) => curry (gen_union (is_equal o  cterm_ord)) (monomial_variables m)) p []);;
  12.327 -
  12.328 -(* Order monomials for human presentation.                                   *)
  12.329 -
  12.330 -fun cterm_ord (t,t') = TermOrd.fast_term_ord (term_of t, term_of t');
  12.331 -
  12.332 -val humanorder_varpow = prod_ord cterm_ord (rev_order o int_ord);
  12.333 -
  12.334 -local
  12.335 - fun ord (l1,l2) = case (l1,l2) of
  12.336 -  (_,[]) => LESS 
  12.337 - | ([],_) => GREATER
  12.338 - | (h1::t1,h2::t2) => 
  12.339 -   (case humanorder_varpow (h1, h2) of 
  12.340 -     LESS => LESS
  12.341 -   | EQUAL => ord (t1,t2)
  12.342 -   | GREATER => GREATER)
  12.343 -in fun humanorder_monomial m1 m2 = 
  12.344 - ord (sort humanorder_varpow (Ctermfunc.graph m1),
  12.345 -  sort humanorder_varpow (Ctermfunc.graph m2))
  12.346 -end;
  12.347 -
  12.348 -fun fold1 f l =  case l of
  12.349 -   []     => error "fold1"
  12.350 - | [x]    => x
  12.351 - | (h::t) => f h (fold1 f t);
  12.352 -
  12.353 -(* Conversions to strings.                                                   *)
  12.354 -
  12.355 -fun string_of_vector min_size max_size (v:vector) =
  12.356 - let val n_raw = dim v 
  12.357 - in if n_raw = 0 then "[]" else
  12.358 -  let 
  12.359 -   val n = max min_size (min n_raw max_size) 
  12.360 -   val xs = map (Rat.string_of_rat o (fn i => Intfunc.tryapplyd (snd v) i rat_0)) (1 upto n) 
  12.361 -  in "[" ^ fold1 (fn s => fn t => s ^ ", " ^ t) xs ^
  12.362 -  (if n_raw > max_size then ", ...]" else "]")
  12.363 -  end
  12.364 - end;
  12.365 -
  12.366 -fun string_of_matrix max_size (m:matrix) =
  12.367 - let 
  12.368 -  val (i_raw,j_raw) = dimensions m
  12.369 -  val i = min max_size i_raw 
  12.370 -  val j = min max_size j_raw
  12.371 -  val rstr = map (fn k => string_of_vector j j (row k m)) (1 upto i) 
  12.372 - in "["^ fold1 (fn s => fn t => s^";\n "^t) rstr ^
  12.373 -  (if j > max_size then "\n ...]" else "]")
  12.374 - end;
  12.375 -
  12.376 -fun string_of_term t = 
  12.377 - case t of
  12.378 -   a$b => "("^(string_of_term a)^" "^(string_of_term b)^")"
  12.379 - | Abs x => 
  12.380 -    let val (xn, b) = Term.dest_abs x
  12.381 -    in "(\\"^xn^"."^(string_of_term b)^")"
  12.382 -    end
  12.383 - | Const(s,_) => s
  12.384 - | Free (s,_) => s
  12.385 - | Var((s,_),_) => s
  12.386 - | _ => error "string_of_term";
  12.387 -
  12.388 -val string_of_cterm = string_of_term o term_of;
  12.389 -
  12.390 -fun string_of_varpow x k =
  12.391 -  if k = 1 then string_of_cterm x 
  12.392 -  else string_of_cterm x^"^"^string_of_int k;
  12.393 -
  12.394 -fun string_of_monomial m =
  12.395 - if Ctermfunc.is_undefined m then "1" else
  12.396 - let val vps = fold_rev (fn (x,k) => fn a => string_of_varpow x k :: a)
  12.397 -  (sort humanorder_varpow (Ctermfunc.graph m)) [] 
  12.398 - in fold1 (fn s => fn t => s^"*"^t) vps
  12.399 - end;
  12.400 -
  12.401 -fun string_of_cmonomial (c,m) =
  12.402 - if Ctermfunc.is_undefined m then Rat.string_of_rat c
  12.403 - else if c =/ rat_1 then string_of_monomial m
  12.404 - else Rat.string_of_rat c ^ "*" ^ string_of_monomial m;;
  12.405 -
  12.406 -fun string_of_poly (p:poly) =
  12.407 - if Monomialfunc.is_undefined p then "<<0>>" else
  12.408 - let 
  12.409 -  val cms = sort (fn ((m1,_),(m2,_)) => humanorder_monomial m1  m2) (Monomialfunc.graph p)
  12.410 -  val s = fold (fn (m,c) => fn a =>
  12.411 -             if c </ rat_0 then a ^ " - " ^ string_of_cmonomial(Rat.neg c,m)
  12.412 -             else a ^ " + " ^ string_of_cmonomial(c,m))
  12.413 -          cms ""
  12.414 -  val s1 = String.substring (s, 0, 3)
  12.415 -  val s2 = String.substring (s, 3, String.size s - 3) 
  12.416 - in "<<" ^(if s1 = " + " then s2 else "-"^s2)^">>"
  12.417 - end;
  12.418 -
  12.419 -(* Conversion from HOL term.                                                 *)
  12.420 -
  12.421 -local
  12.422 - val neg_tm = @{cterm "uminus :: real => _"}
  12.423 - val add_tm = @{cterm "op + :: real => _"}
  12.424 - val sub_tm = @{cterm "op - :: real => _"}
  12.425 - val mul_tm = @{cterm "op * :: real => _"}
  12.426 - val inv_tm = @{cterm "inverse :: real => _"}
  12.427 - val div_tm = @{cterm "op / :: real => _"}
  12.428 - val pow_tm = @{cterm "op ^ :: real => _"}
  12.429 - val zero_tm = @{cterm "0:: real"}
  12.430 - val is_numeral = can (HOLogic.dest_number o term_of)
  12.431 - fun is_comb t = case t of _$_ => true | _ => false
  12.432 - fun poly_of_term tm =
  12.433 -  if tm aconvc zero_tm then poly_0
  12.434 -  else if RealArith.is_ratconst tm 
  12.435 -       then poly_const(RealArith.dest_ratconst tm)
  12.436 -  else 
  12.437 -  (let val (lop,r) = Thm.dest_comb tm
  12.438 -   in if lop aconvc neg_tm then poly_neg(poly_of_term r)
  12.439 -      else if lop aconvc inv_tm then
  12.440 -       let val p = poly_of_term r 
  12.441 -       in if poly_isconst p 
  12.442 -          then poly_const(Rat.inv (eval Ctermfunc.undefined p))
  12.443 -          else error "poly_of_term: inverse of non-constant polyomial"
  12.444 -       end
  12.445 -   else (let val (opr,l) = Thm.dest_comb lop
  12.446 -         in 
  12.447 -          if opr aconvc pow_tm andalso is_numeral r 
  12.448 -          then poly_pow (poly_of_term l) ((snd o HOLogic.dest_number o term_of) r)
  12.449 -          else if opr aconvc add_tm 
  12.450 -           then poly_add (poly_of_term l) (poly_of_term r)
  12.451 -          else if opr aconvc sub_tm 
  12.452 -           then poly_sub (poly_of_term l) (poly_of_term r)
  12.453 -          else if opr aconvc mul_tm 
  12.454 -           then poly_mul (poly_of_term l) (poly_of_term r)
  12.455 -          else if opr aconvc div_tm 
  12.456 -           then let 
  12.457 -                  val p = poly_of_term l 
  12.458 -                  val q = poly_of_term r 
  12.459 -                in if poly_isconst q then poly_cmul (Rat.inv (eval Ctermfunc.undefined q)) p
  12.460 -                   else error "poly_of_term: division by non-constant polynomial"
  12.461 -                end
  12.462 -          else poly_var tm
  12.463 - 
  12.464 -         end
  12.465 -         handle CTERM ("dest_comb",_) => poly_var tm)
  12.466 -   end
  12.467 -   handle CTERM ("dest_comb",_) => poly_var tm)
  12.468 -in
  12.469 -val poly_of_term = fn tm =>
  12.470 - if type_of (term_of tm) = @{typ real} then poly_of_term tm
  12.471 - else error "poly_of_term: term does not have real type"
  12.472 -end;
  12.473 -
  12.474 -(* String of vector (just a list of space-separated numbers).                *)
  12.475 -
  12.476 -fun sdpa_of_vector (v:vector) =
  12.477 - let 
  12.478 -  val n = dim v
  12.479 -  val strs = map (decimalize 20 o (fn i => Intfunc.tryapplyd (snd v) i rat_0)) (1 upto n) 
  12.480 - in fold1 (fn x => fn y => x ^ " " ^ y) strs ^ "\n"
  12.481 - end;
  12.482 -
  12.483 -fun increasing f ord (x,y) = ord (f x, f y);
  12.484 -fun triple_int_ord ((a,b,c),(a',b',c')) = 
  12.485 - prod_ord int_ord (prod_ord int_ord int_ord) 
  12.486 -    ((a,(b,c)),(a',(b',c')));
  12.487 -structure Inttriplefunc = FuncFun(type key = int*int*int val ord = triple_int_ord);
  12.488 -
  12.489 -(* String for block diagonal matrix numbered k.                              *)
  12.490 -
  12.491 -fun sdpa_of_blockdiagonal k m =
  12.492 - let 
  12.493 -  val pfx = string_of_int k ^" "
  12.494 -  val ents =
  12.495 -    Inttriplefunc.fold (fn ((b,i,j), c) => fn a => if i > j then a else ((b,i,j),c)::a) m []
  12.496 -  val entss = sort (increasing fst triple_int_ord ) ents
  12.497 - in  fold_rev (fn ((b,i,j),c) => fn a =>
  12.498 -     pfx ^ string_of_int b ^ " " ^ string_of_int i ^ " " ^ string_of_int j ^
  12.499 -     " " ^ decimalize 20 c ^ "\n" ^ a) entss ""
  12.500 - end;
  12.501 -
  12.502 -(* String for a matrix numbered k, in SDPA sparse format.                    *)
  12.503 -
  12.504 -fun sdpa_of_matrix k (m:matrix) =
  12.505 - let 
  12.506 -  val pfx = string_of_int k ^ " 1 "
  12.507 -  val ms = Intpairfunc.fold (fn ((i,j), c) => fn  a => if i > j then a else ((i,j),c)::a) (snd m) [] 
  12.508 -  val mss = sort (increasing fst (prod_ord int_ord int_ord)) ms 
  12.509 - in fold_rev (fn ((i,j),c) => fn a =>
  12.510 -     pfx ^ string_of_int i ^ " " ^ string_of_int j ^
  12.511 -     " " ^ decimalize 20 c ^ "\n" ^ a) mss ""
  12.512 - end;;
  12.513 -
  12.514 -(* ------------------------------------------------------------------------- *)
  12.515 -(* String in SDPA sparse format for standard SDP problem:                    *)
  12.516 -(*                                                                           *)
  12.517 -(*    X = v_1 * [M_1] + ... + v_m * [M_m] - [M_0] must be PSD                *)
  12.518 -(*    Minimize obj_1 * v_1 + ... obj_m * v_m                                 *)
  12.519 -(* ------------------------------------------------------------------------- *)
  12.520 -
  12.521 -fun sdpa_of_problem obj mats =
  12.522 - let 
  12.523 -  val m = length mats - 1
  12.524 -  val (n,_) = dimensions (hd mats) 
  12.525 - in
  12.526 -  string_of_int m ^ "\n" ^
  12.527 -  "1\n" ^
  12.528 -  string_of_int n ^ "\n" ^
  12.529 -  sdpa_of_vector obj ^
  12.530 -  fold_rev2 (fn k => fn m => fn a => sdpa_of_matrix (k - 1) m ^ a) (1 upto length mats) mats ""
  12.531 - end;
  12.532 -
  12.533 -fun index_char str chr pos =
  12.534 -  if pos >= String.size str then ~1
  12.535 -  else if String.sub(str,pos) = chr then pos
  12.536 -  else index_char str chr (pos + 1);
  12.537 -fun rat_of_quotient (a,b) = if b = 0 then rat_0 else Rat.rat_of_quotient (a,b);
  12.538 -fun rat_of_string s = 
  12.539 - let val n = index_char s #"/" 0 in
  12.540 -  if n = ~1 then s |> IntInf.fromString |> valOf |> Rat.rat_of_int
  12.541 -  else 
  12.542 -   let val SOME numer = IntInf.fromString(String.substring(s,0,n))
  12.543 -       val SOME den = IntInf.fromString (String.substring(s,n+1,String.size s - n - 1))
  12.544 -   in rat_of_quotient(numer, den)
  12.545 -   end
  12.546 - end;
  12.547 -
  12.548 -fun isspace x = x = " " ;
  12.549 -fun isnum x = x mem_string ["0","1","2","3","4","5","6","7","8","9"]
  12.550 -
  12.551 -(* More parser basics.                                                       *)
  12.552 -
  12.553 -local
  12.554 - open Scan
  12.555 -in 
  12.556 - val word = this_string
  12.557 - fun token s =
  12.558 -  repeat ($$ " ") |-- word s --| repeat ($$ " ")
  12.559 - val numeral = one isnum
  12.560 - val decimalint = bulk numeral >> (rat_of_string o implode)
  12.561 - val decimalfrac = bulk numeral
  12.562 -    >> (fn s => rat_of_string(implode s) // pow10 (length s))
  12.563 - val decimalsig =
  12.564 -    decimalint -- option (Scan.$$ "." |-- decimalfrac)
  12.565 -    >> (fn (h,NONE) => h | (h,SOME x) => h +/ x)
  12.566 - fun signed prs =
  12.567 -       $$ "-" |-- prs >> Rat.neg 
  12.568 -    || $$ "+" |-- prs
  12.569 -    || prs;
  12.570 -
  12.571 -fun emptyin def xs = if null xs then (def,xs) else Scan.fail xs
  12.572 -
  12.573 - val exponent = ($$ "e" || $$ "E") |-- signed decimalint;
  12.574 -
  12.575 - val decimal = signed decimalsig -- (emptyin rat_0|| exponent)
  12.576 -    >> (fn (h, x) => h */ pow10 (int_of_rat x));
  12.577 -end;
  12.578 -
  12.579 - fun mkparser p s =
  12.580 -  let val (x,rst) = p (explode s) 
  12.581 -  in if null rst then x 
  12.582 -     else error "mkparser: unparsed input"
  12.583 -  end;;
  12.584 -val parse_decimal = mkparser decimal;
  12.585 -
  12.586 -fun fix err prs = 
  12.587 -  prs || (fn x=> error err);
  12.588 -
  12.589 -fun listof prs sep err =
  12.590 -  prs -- Scan.bulk (sep |-- fix err prs) >> uncurry cons;
  12.591 -
  12.592 -(* Parse back a vector.                                                      *)
  12.593 -
  12.594 - val vector = 
  12.595 -    token "{" |-- listof decimal (token ",") "decimal" --| token "}"
  12.596 -               >> vector_of_list 
  12.597 - val parse_vector = mkparser vector
  12.598 - fun skipupto dscr prs inp =
  12.599 -   (dscr |-- prs 
  12.600 -    || Scan.one (K true) |-- skipupto dscr prs) inp 
  12.601 - fun ignore inp = ((),[])
  12.602 - fun sdpaoutput inp =  skipupto (word "xVec" -- token "=")
  12.603 -             (vector --| ignore) inp
  12.604 - fun csdpoutput inp =  ((decimal -- Scan.bulk (Scan.$$ " " |-- Scan.option decimal) >> (fn (h,to) => map_filter I ((SOME h)::to))) --| ignore >> vector_of_list) inp
  12.605 - val parse_sdpaoutput = mkparser sdpaoutput
  12.606 - val parse_csdpoutput = mkparser csdpoutput
  12.607 -
  12.608 -(* Run prover on a problem in linear form.                       *)
  12.609 -
  12.610 -fun run_problem prover obj mats =
  12.611 -  parse_csdpoutput (prover (sdpa_of_problem obj mats))
  12.612 -
  12.613 -(*
  12.614 -UNUSED
  12.615 -
  12.616 -(* Also parse the SDPA output to test success (CSDP yields a return code).   *)
  12.617 -
  12.618 -local
  12.619 - val prs = 
  12.620 -  skipupto (word "phase.value" -- token "=")
  12.621 -   (Scan.option (Scan.$$ "p") -- Scan.option (Scan.$$ "d") 
  12.622 -    -- (word "OPT" || word "FEAS")) 
  12.623 -in
  12.624 - fun sdpa_run_succeeded s = 
  12.625 -  (prs (explode s); true) handle _ => false
  12.626 -end;
  12.627 -
  12.628 -(* The default parameters. Unfortunately this goes to a fixed file.          *)
  12.629 -
  12.630 -val sdpa_default_parameters =
  12.631 -"100     unsigned int maxIteration; \n1.0E-7  double 0.0 < epsilonStar;\n1.0E2   double 0.0 < lambdaStar;\n2.0     double 1.0 < omegaStar;\n-1.0E5  double lowerBound;\n1.0E5   double upperBound;\n0.1     double 0.0 <= betaStar <  1.0;\n0.2     double 0.0 <= betaBar  <  1.0, betaStar <= betaBar;\n0.9     double 0.0 < gammaStar  <  1.0;\n1.0E-7  double 0.0 < epsilonDash;\n";;
  12.632 -
  12.633 -(* These were suggested by Makoto Yamashita for problems where we are        *)
  12.634 -(* right at the edge of the semidefinite cone, as sometimes happens.         *)
  12.635 -
  12.636 -val sdpa_alt_parameters =
  12.637 -"1000    unsigned int maxIteration;\n1.0E-7  double 0.0 < epsilonStar;\n1.0E4   double 0.0 < lambdaStar;\n2.0     double 1.0 < omegaStar;\n-1.0E5  double lowerBound;\n1.0E5   double upperBound;\n0.1     double 0.0 <= betaStar <  1.0;\n0.2     double 0.0 <= betaBar  <  1.0, betaStar <= betaBar;\n0.9     double 0.0 < gammaStar  <  1.0;\n1.0E-7  double 0.0 < epsilonDash;\n";;
  12.638 -
  12.639 -val sdpa_params = sdpa_alt_parameters;;
  12.640 -
  12.641 -(* CSDP parameters; so far I'm sticking with the defaults.                   *)
  12.642 -
  12.643 -val csdp_default_parameters =
  12.644 -"axtol=1.0e-8\natytol=1.0e-8\nobjtol=1.0e-8\npinftol=1.0e8\ndinftol=1.0e8\nmaxiter=100\nminstepfrac=0.9\nmaxstepfrac=0.97\nminstepp=1.0e-8\nminstepd=1.0e-8\nusexzgap=1\ntweakgap=0\naffine=0\nprintlevel=1\n";;
  12.645 -
  12.646 -val csdp_params = csdp_default_parameters;;
  12.647 -
  12.648 -fun tmp_file pre suf =
  12.649 - let val i = string_of_int (round (random()))
  12.650 -   val name = Path.append (Path.variable "ISABELLE_TMP") (Path.explode (pre ^ i ^ suf))
  12.651 - in 
  12.652 -   if File.exists name then tmp_file pre suf 
  12.653 -   else name 
  12.654 - end;
  12.655 -
  12.656 -(* Now call SDPA on a problem and parse back the output.                     *)
  12.657 -
  12.658 -fun run_sdpa dbg obj mats =
  12.659 - let 
  12.660 -  val input_file = tmp_file "sos" ".dat-s"
  12.661 -  val output_file = tmp_file "sos" ".out"
  12.662 -  val params_file = tmp_file "param" ".sdpa" 
  12.663 -  val current_dir = File.pwd()
  12.664 -  val _ = File.write input_file 
  12.665 -                         (sdpa_of_problem "" obj mats)
  12.666 -  val _ = File.write params_file sdpa_params
  12.667 -  val _ = File.cd (Path.variable "ISABELLE_TMP")
  12.668 -  val _ = File.system_command ("sdpa "^ (Path.implode input_file) ^ " " ^ 
  12.669 -                               (Path.implode output_file) ^
  12.670 -                               (if dbg then "" else "> /dev/null"))
  12.671 -  val opr = File.read output_file 
  12.672 - in if not(sdpa_run_succeeded opr) then error "sdpa: call failed" 
  12.673 -    else
  12.674 -      let val res = parse_sdpaoutput opr 
  12.675 -      in ((if dbg then ()
  12.676 -           else (File.rm input_file; File.rm output_file ; File.cd current_dir));
  12.677 -          res)
  12.678 -      end
  12.679 - end;
  12.680 -
  12.681 -fun sdpa obj mats = run_sdpa (!debugging) obj mats;
  12.682 -
  12.683 -(* The same thing with CSDP.                                                 *)
  12.684 -
  12.685 -fun run_csdp dbg obj mats =
  12.686 - let 
  12.687 -  val input_file = tmp_file "sos" ".dat-s"
  12.688 -  val output_file = tmp_file "sos" ".out"
  12.689 -  val params_file = tmp_file "param" ".csdp"
  12.690 -  val current_dir = File.pwd()
  12.691 -  val _ = File.write input_file (sdpa_of_problem "" obj mats)
  12.692 -  val _ = File.write params_file csdp_params
  12.693 -  val _ = File.cd (Path.variable "ISABELLE_TMP")
  12.694 -  val rv = system ("csdp "^(Path.implode input_file) ^ " " 
  12.695 -                   ^ (Path.implode output_file) ^
  12.696 -                   (if dbg then "" else "> /dev/null"))
  12.697 -  val  opr = File.read output_file 
  12.698 -  val res = parse_csdpoutput opr 
  12.699 - in
  12.700 -    ((if dbg then ()
  12.701 -      else (File.rm input_file; File.rm output_file ; File.cd current_dir));
  12.702 -     (rv,res))
  12.703 - end;
  12.704 -
  12.705 -fun csdp obj mats =
  12.706 - let 
  12.707 -  val (rv,res) = run_csdp (!debugging) obj mats 
  12.708 - in
  12.709 -   ((if rv = 1 orelse rv = 2 then error "csdp: Problem is infeasible"
  12.710 -    else if rv = 3 then writeln "csdp warning: Reduced accuracy"
  12.711 -    else if rv <> 0 then error ("csdp: error "^string_of_int rv)
  12.712 -    else ());
  12.713 -   res)
  12.714 - end;
  12.715 -
  12.716 -*)
  12.717 -
  12.718 -(* Try some apparently sensible scaling first. Note that this is purely to   *)
  12.719 -(* get a cleaner translation to floating-point, and doesn't affect any of    *)
  12.720 -(* the results, in principle. In practice it seems a lot better when there   *)
  12.721 -(* are extreme numbers in the original problem.                              *)
  12.722 -
  12.723 -  (* Version for (int*int) keys *)
  12.724 -local
  12.725 -  fun max_rat x y = if x </ y then y else x
  12.726 -  fun common_denominator fld amat acc =
  12.727 -      fld (fn (m,c) => fn a => lcm_rat (denominator_rat c) a) amat acc
  12.728 -  fun maximal_element fld amat acc =
  12.729 -    fld (fn (m,c) => fn maxa => max_rat maxa (abs_rat c)) amat acc 
  12.730 -fun float_of_rat x = let val (a,b) = Rat.quotient_of_rat x
  12.731 -                     in Real.fromLargeInt a / Real.fromLargeInt b end;
  12.732 -in
  12.733 -
  12.734 -fun pi_scale_then solver (obj:vector)  mats =
  12.735 - let 
  12.736 -  val cd1 = fold_rev (common_denominator Intpairfunc.fold) mats (rat_1)
  12.737 -  val cd2 = common_denominator Intfunc.fold (snd obj)  (rat_1) 
  12.738 -  val mats' = map (Intpairfunc.mapf (fn x => cd1 */ x)) mats
  12.739 -  val obj' = vector_cmul cd2 obj
  12.740 -  val max1 = fold_rev (maximal_element Intpairfunc.fold) mats' (rat_0)
  12.741 -  val max2 = maximal_element Intfunc.fold (snd obj') (rat_0) 
  12.742 -  val scal1 = pow2 (20 - trunc(Math.ln (float_of_rat max1) / Math.ln 2.0))
  12.743 -  val scal2 = pow2 (20 - trunc(Math.ln (float_of_rat max2) / Math.ln 2.0)) 
  12.744 -  val mats'' = map (Intpairfunc.mapf (fn x => x */ scal1)) mats'
  12.745 -  val obj'' = vector_cmul scal2 obj' 
  12.746 - in solver obj'' mats''
  12.747 -  end
  12.748 -end;
  12.749 -
  12.750 -(* Try some apparently sensible scaling first. Note that this is purely to   *)
  12.751 -(* get a cleaner translation to floating-point, and doesn't affect any of    *)
  12.752 -(* the results, in principle. In practice it seems a lot better when there   *)
  12.753 -(* are extreme numbers in the original problem.                              *)
  12.754 -
  12.755 -  (* Version for (int*int*int) keys *)
  12.756 -local
  12.757 -  fun max_rat x y = if x </ y then y else x
  12.758 -  fun common_denominator fld amat acc =
  12.759 -      fld (fn (m,c) => fn a => lcm_rat (denominator_rat c) a) amat acc
  12.760 -  fun maximal_element fld amat acc =
  12.761 -    fld (fn (m,c) => fn maxa => max_rat maxa (abs_rat c)) amat acc 
  12.762 -fun float_of_rat x = let val (a,b) = Rat.quotient_of_rat x
  12.763 -                     in Real.fromLargeInt a / Real.fromLargeInt b end;
  12.764 -fun int_of_float x = (trunc x handle Overflow => 0 | Domain => 0)
  12.765 -in
  12.766 -
  12.767 -fun tri_scale_then solver (obj:vector)  mats =
  12.768 - let 
  12.769 -  val cd1 = fold_rev (common_denominator Inttriplefunc.fold) mats (rat_1)
  12.770 -  val cd2 = common_denominator Intfunc.fold (snd obj)  (rat_1) 
  12.771 -  val mats' = map (Inttriplefunc.mapf (fn x => cd1 */ x)) mats
  12.772 -  val obj' = vector_cmul cd2 obj
  12.773 -  val max1 = fold_rev (maximal_element Inttriplefunc.fold) mats' (rat_0)
  12.774 -  val max2 = maximal_element Intfunc.fold (snd obj') (rat_0) 
  12.775 -  val scal1 = pow2 (20 - int_of_float(Math.ln (float_of_rat max1) / Math.ln 2.0))
  12.776 -  val scal2 = pow2 (20 - int_of_float(Math.ln (float_of_rat max2) / Math.ln 2.0)) 
  12.777 -  val mats'' = map (Inttriplefunc.mapf (fn x => x */ scal1)) mats'
  12.778 -  val obj'' = vector_cmul scal2 obj' 
  12.779 - in solver obj'' mats''
  12.780 -  end
  12.781 -end;
  12.782 -
  12.783 -(* Round a vector to "nice" rationals.                                       *)
  12.784 -
  12.785 -fun nice_rational n x = round_rat (n */ x) // n;;
  12.786 -fun nice_vector n ((d,v) : vector) = 
  12.787 - (d, Intfunc.fold (fn (i,c) => fn a => 
  12.788 -   let val y = nice_rational n c 
  12.789 -   in if c =/ rat_0 then a 
  12.790 -      else Intfunc.update (i,y) a end) v Intfunc.undefined):vector
  12.791 -
  12.792 -(*
  12.793 -UNUSED
  12.794 -
  12.795 -(* Reduce linear program to SDP (diagonal matrices) and test with CSDP. This *)
  12.796 -(* one tests A [-1;x1;..;xn] >= 0 (i.e. left column is negated constants).   *)
  12.797 -
  12.798 -fun linear_program_basic a =
  12.799 - let 
  12.800 -  val (m,n) = dimensions a
  12.801 -  val mats =  map (fn j => diagonal (column j a)) (1 upto n)
  12.802 -  val obj = vector_const rat_1 m 
  12.803 -  val (rv,res) = run_csdp false obj mats 
  12.804 - in if rv = 1 orelse rv = 2 then false
  12.805 -    else if rv = 0 then true
  12.806 -    else error "linear_program: An error occurred in the SDP solver"
  12.807 - end;
  12.808 -
  12.809 -(* Alternative interface testing A x >= b for matrix A, vector b.            *)
  12.810 -
  12.811 -fun linear_program a b =
  12.812 - let val (m,n) = dimensions a 
  12.813 - in if dim b <> m then error "linear_program: incompatible dimensions" 
  12.814 -    else
  12.815 -    let 
  12.816 -     val mats = diagonal b :: map (fn j => diagonal (column j a)) (1 upto n)
  12.817 -     val obj = vector_const rat_1 m 
  12.818 -     val (rv,res) = run_csdp false obj mats 
  12.819 -    in if rv = 1 orelse rv = 2 then false
  12.820 -       else if rv = 0 then true
  12.821 -       else error "linear_program: An error occurred in the SDP solver"
  12.822 -    end
  12.823 - end;
  12.824 -
  12.825 -(* Test whether a point is in the convex hull of others. Rather than use     *)
  12.826 -(* computational geometry, express as linear inequalities and call CSDP.     *)
  12.827 -(* This is a bit lazy of me, but it's easy and not such a bottleneck so far. *)
  12.828 -
  12.829 -fun in_convex_hull pts pt =
  12.830 - let 
  12.831 -  val pts1 = (1::pt) :: map (fn x => 1::x) pts 
  12.832 -  val pts2 = map (fn p => map (fn x => ~x) p @ p) pts1
  12.833 -  val n = length pts + 1
  12.834 -  val v = 2 * (length pt + 1)
  12.835 -  val m = v + n - 1 
  12.836 -  val mat = ((m,n),
  12.837 -  itern 1 pts2 (fn pts => fn j => itern 1 pts 
  12.838 -               (fn x => fn i => Intpairfunc.update ((i,j), Rat.rat_of_int x)))
  12.839 -  (iter (1,n) (fn i => Intpairfunc.update((v + i,i+1), rat_1)) 
  12.840 -      Intpairfunc.undefined))
  12.841 - in linear_program_basic mat
  12.842 - end;
  12.843 -
  12.844 -(* Filter down a set of points to a minimal set with the same convex hull.   *)
  12.845 -
  12.846 -local
  12.847 - fun augment1 (m::ms) = if in_convex_hull ms m then ms else ms@[m]
  12.848 - fun augment m ms = funpow 3 augment1 (m::ms)
  12.849 -in
  12.850 -fun minimal_convex_hull mons =
  12.851 - let val mons' = fold_rev augment (tl mons) [hd mons] 
  12.852 - in funpow (length mons') augment1 mons'
  12.853 - end
  12.854 -end;
  12.855 -
  12.856 -*)
  12.857 -
  12.858 -fun dest_ord f x = is_equal (f x);
  12.859 -
  12.860 -
  12.861 -
  12.862 -(* Stuff for "equations" ((int*int*int)->num functions).                         *)
  12.863 -
  12.864 -fun tri_equation_cmul c eq =
  12.865 -  if c =/ rat_0 then Inttriplefunc.undefined else Inttriplefunc.mapf (fn d => c */ d) eq;
  12.866 -
  12.867 -fun tri_equation_add eq1 eq2 = Inttriplefunc.combine (curry op +/) (fn x => x =/ rat_0) eq1 eq2;
  12.868 -
  12.869 -fun tri_equation_eval assig eq =
  12.870 - let fun value v = Inttriplefunc.apply assig v 
  12.871 - in Inttriplefunc.fold (fn (v, c) => fn a => a +/ value v */ c) eq rat_0
  12.872 - end;
  12.873 -
  12.874 -(* Eliminate among linear equations: return unconstrained variables and      *)
  12.875 -(* assignments for the others in terms of them. We give one pseudo-variable  *)
  12.876 -(* "one" that's used for a constant term.                                    *)
  12.877 -
  12.878 -local
  12.879 -  fun extract_first p l = case l of  (* FIXME : use find_first instead *)
  12.880 -   [] => error "extract_first"
  12.881 - | h::t => if p h then (h,t) else
  12.882 -          let val (k,s) = extract_first p t in (k,h::s) end
  12.883 -fun eliminate vars dun eqs = case vars of 
  12.884 -  [] => if forall Inttriplefunc.is_undefined eqs then dun
  12.885 -        else raise Unsolvable
  12.886 - | v::vs =>
  12.887 -  ((let 
  12.888 -    val (eq,oeqs) = extract_first (fn e => Inttriplefunc.defined e v) eqs 
  12.889 -    val a = Inttriplefunc.apply eq v
  12.890 -    val eq' = tri_equation_cmul ((Rat.neg rat_1) // a) (Inttriplefunc.undefine v eq)
  12.891 -    fun elim e =
  12.892 -     let val b = Inttriplefunc.tryapplyd e v rat_0 
  12.893 -     in if b =/ rat_0 then e else
  12.894 -        tri_equation_add e (tri_equation_cmul (Rat.neg b // a) eq)
  12.895 -     end
  12.896 -   in eliminate vs (Inttriplefunc.update (v,eq') (Inttriplefunc.mapf elim dun)) (map elim oeqs)
  12.897 -   end)
  12.898 -  handle ERROR _ => eliminate vs dun eqs)
  12.899 -in
  12.900 -fun tri_eliminate_equations one vars eqs =
  12.901 - let 
  12.902 -  val assig = eliminate vars Inttriplefunc.undefined eqs
  12.903 -  val vs = Inttriplefunc.fold (fn (x, f) => fn a => remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig []
  12.904 -  in (distinct (dest_ord triple_int_ord) vs, assig)
  12.905 -  end
  12.906 -end;
  12.907 -
  12.908 -(* Eliminate all variables, in an essentially arbitrary order.               *)
  12.909 -
  12.910 -fun tri_eliminate_all_equations one =
  12.911 - let 
  12.912 -  fun choose_variable eq =
  12.913 -   let val (v,_) = Inttriplefunc.choose eq 
  12.914 -   in if is_equal (triple_int_ord(v,one)) then
  12.915 -      let val eq' = Inttriplefunc.undefine v eq 
  12.916 -      in if Inttriplefunc.is_undefined eq' then error "choose_variable" 
  12.917 -         else fst (Inttriplefunc.choose eq')
  12.918 -      end
  12.919 -    else v 
  12.920 -   end
  12.921 -  fun eliminate dun eqs = case eqs of 
  12.922 -    [] => dun
  12.923 -  | eq::oeqs =>
  12.924 -    if Inttriplefunc.is_undefined eq then eliminate dun oeqs else
  12.925 -    let val v = choose_variable eq
  12.926 -        val a = Inttriplefunc.apply eq v
  12.927 -        val eq' = tri_equation_cmul ((Rat.rat_of_int ~1) // a) 
  12.928 -                   (Inttriplefunc.undefine v eq)
  12.929 -        fun elim e =
  12.930 -         let val b = Inttriplefunc.tryapplyd e v rat_0 
  12.931 -         in if b =/ rat_0 then e 
  12.932 -            else tri_equation_add e (tri_equation_cmul (Rat.neg b // a) eq)
  12.933 -         end
  12.934 -    in eliminate (Inttriplefunc.update(v, eq') (Inttriplefunc.mapf elim dun)) 
  12.935 -                 (map elim oeqs) 
  12.936 -    end
  12.937 -in fn eqs =>
  12.938 - let 
  12.939 -  val assig = eliminate Inttriplefunc.undefined eqs
  12.940 -  val vs = Inttriplefunc.fold (fn (x, f) => fn a => remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig []
  12.941 - in (distinct (dest_ord triple_int_ord) vs,assig)
  12.942 - end
  12.943 -end;
  12.944 - 
  12.945 -(* Solve equations by assigning arbitrary numbers.                           *)
  12.946 -
  12.947 -fun tri_solve_equations one eqs =
  12.948 - let 
  12.949 -  val (vars,assigs) = tri_eliminate_all_equations one eqs
  12.950 -  val vfn = fold_rev (fn v => Inttriplefunc.update(v,rat_0)) vars 
  12.951 -            (Inttriplefunc.onefunc(one, Rat.rat_of_int ~1))
  12.952 -  val ass =
  12.953 -    Inttriplefunc.combine (curry op +/) (K false) 
  12.954 -    (Inttriplefunc.mapf (tri_equation_eval vfn) assigs) vfn 
  12.955 - in if forall (fn e => tri_equation_eval ass e =/ rat_0) eqs
  12.956 -    then Inttriplefunc.undefine one ass else raise Sanity
  12.957 - end;
  12.958 -
  12.959 -(* Multiply equation-parametrized poly by regular poly and add accumulator.  *)
  12.960 -
  12.961 -fun tri_epoly_pmul p q acc =
  12.962 - Monomialfunc.fold (fn (m1, c) => fn a =>
  12.963 -  Monomialfunc.fold (fn (m2,e) => fn b =>
  12.964 -   let val m =  monomial_mul m1 m2
  12.965 -       val es = Monomialfunc.tryapplyd b m Inttriplefunc.undefined 
  12.966 -   in Monomialfunc.update (m,tri_equation_add (tri_equation_cmul c e) es) b 
  12.967 -   end) q a) p acc ;
  12.968 -
  12.969 -(* Usual operations on equation-parametrized poly.                           *)
  12.970 -
  12.971 -fun tri_epoly_cmul c l =
  12.972 -  if c =/ rat_0 then Inttriplefunc.undefined else Inttriplefunc.mapf (tri_equation_cmul c) l;;
  12.973 -
  12.974 -val tri_epoly_neg = tri_epoly_cmul (Rat.rat_of_int ~1);
  12.975 -
  12.976 -val tri_epoly_add = Inttriplefunc.combine tri_equation_add Inttriplefunc.is_undefined;
  12.977 -
  12.978 -fun tri_epoly_sub p q = tri_epoly_add p (tri_epoly_neg q);;
  12.979 -
  12.980 -(* Stuff for "equations" ((int*int)->num functions).                         *)
  12.981 -
  12.982 -fun pi_equation_cmul c eq =
  12.983 -  if c =/ rat_0 then Inttriplefunc.undefined else Inttriplefunc.mapf (fn d => c */ d) eq;
  12.984 -
  12.985 -fun pi_equation_add eq1 eq2 = Inttriplefunc.combine (curry op +/) (fn x => x =/ rat_0) eq1 eq2;
  12.986 -
  12.987 -fun pi_equation_eval assig eq =
  12.988 - let fun value v = Inttriplefunc.apply assig v 
  12.989 - in Inttriplefunc.fold (fn (v, c) => fn a => a +/ value v */ c) eq rat_0
  12.990 - end;
  12.991 -
  12.992 -(* Eliminate among linear equations: return unconstrained variables and      *)
  12.993 -(* assignments for the others in terms of them. We give one pseudo-variable  *)
  12.994 -(* "one" that's used for a constant term.                                    *)
  12.995 -
  12.996 -local
  12.997 -fun extract_first p l = case l of 
  12.998 -   [] => error "extract_first"
  12.999 - | h::t => if p h then (h,t) else
 12.1000 -          let val (k,s) = extract_first p t in (k,h::s) end
 12.1001 -fun eliminate vars dun eqs = case vars of 
 12.1002 -  [] => if forall Inttriplefunc.is_undefined eqs then dun
 12.1003 -        else raise Unsolvable
 12.1004 - | v::vs =>
 12.1005 -   let 
 12.1006 -    val (eq,oeqs) = extract_first (fn e => Inttriplefunc.defined e v) eqs 
 12.1007 -    val a = Inttriplefunc.apply eq v
 12.1008 -    val eq' = pi_equation_cmul ((Rat.neg rat_1) // a) (Inttriplefunc.undefine v eq)
 12.1009 -    fun elim e =
 12.1010 -     let val b = Inttriplefunc.tryapplyd e v rat_0 
 12.1011 -     in if b =/ rat_0 then e else
 12.1012 -        pi_equation_add e (pi_equation_cmul (Rat.neg b // a) eq)
 12.1013 -     end
 12.1014 -   in eliminate vs (Inttriplefunc.update (v,eq') (Inttriplefunc.mapf elim dun)) (map elim oeqs)
 12.1015 -   end
 12.1016 -  handle ERROR _ => eliminate vs dun eqs
 12.1017 -in
 12.1018 -fun pi_eliminate_equations one vars eqs =
 12.1019 - let 
 12.1020 -  val assig = eliminate vars Inttriplefunc.undefined eqs
 12.1021 -  val vs = Inttriplefunc.fold (fn (x, f) => fn a => remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig []
 12.1022 -  in (distinct (dest_ord triple_int_ord) vs, assig)
 12.1023 -  end
 12.1024 -end;
 12.1025 -
 12.1026 -(* Eliminate all variables, in an essentially arbitrary order.               *)
 12.1027 -
 12.1028 -fun pi_eliminate_all_equations one =
 12.1029 - let 
 12.1030 -  fun choose_variable eq =
 12.1031 -   let val (v,_) = Inttriplefunc.choose eq 
 12.1032 -   in if is_equal (triple_int_ord(v,one)) then
 12.1033 -      let val eq' = Inttriplefunc.undefine v eq 
 12.1034 -      in if Inttriplefunc.is_undefined eq' then error "choose_variable" 
 12.1035 -         else fst (Inttriplefunc.choose eq')
 12.1036 -      end
 12.1037 -    else v 
 12.1038 -   end
 12.1039 -  fun eliminate dun eqs = case eqs of 
 12.1040 -    [] => dun
 12.1041 -  | eq::oeqs =>
 12.1042 -    if Inttriplefunc.is_undefined eq then eliminate dun oeqs else
 12.1043 -    let val v = choose_variable eq
 12.1044 -        val a = Inttriplefunc.apply eq v
 12.1045 -        val eq' = pi_equation_cmul ((Rat.rat_of_int ~1) // a) 
 12.1046 -                   (Inttriplefunc.undefine v eq)
 12.1047 -        fun elim e =
 12.1048 -         let val b = Inttriplefunc.tryapplyd e v rat_0 
 12.1049 -         in if b =/ rat_0 then e 
 12.1050 -            else pi_equation_add e (pi_equation_cmul (Rat.neg b // a) eq)
 12.1051 -         end
 12.1052 -    in eliminate (Inttriplefunc.update(v, eq') (Inttriplefunc.mapf elim dun)) 
 12.1053 -                 (map elim oeqs) 
 12.1054 -    end
 12.1055 -in fn eqs =>
 12.1056 - let 
 12.1057 -  val assig = eliminate Inttriplefunc.undefined eqs
 12.1058 -  val vs = Inttriplefunc.fold (fn (x, f) => fn a => remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig []
 12.1059 - in (distinct (dest_ord triple_int_ord) vs,assig)
 12.1060 - end
 12.1061 -end;
 12.1062 - 
 12.1063 -(* Solve equations by assigning arbitrary numbers.                           *)
 12.1064 -
 12.1065 -fun pi_solve_equations one eqs =
 12.1066 - let 
 12.1067 -  val (vars,assigs) = pi_eliminate_all_equations one eqs
 12.1068 -  val vfn = fold_rev (fn v => Inttriplefunc.update(v,rat_0)) vars 
 12.1069 -            (Inttriplefunc.onefunc(one, Rat.rat_of_int ~1))
 12.1070 -  val ass =
 12.1071 -    Inttriplefunc.combine (curry op +/) (K false) 
 12.1072 -    (Inttriplefunc.mapf (pi_equation_eval vfn) assigs) vfn 
 12.1073 - in if forall (fn e => pi_equation_eval ass e =/ rat_0) eqs
 12.1074 -    then Inttriplefunc.undefine one ass else raise Sanity
 12.1075 - end;
 12.1076 -
 12.1077 -(* Multiply equation-parametrized poly by regular poly and add accumulator.  *)
 12.1078 -
 12.1079 -fun pi_epoly_pmul p q acc =
 12.1080 - Monomialfunc.fold (fn (m1, c) => fn a =>
 12.1081 -  Monomialfunc.fold (fn (m2,e) => fn b =>
 12.1082 -   let val m =  monomial_mul m1 m2
 12.1083 -       val es = Monomialfunc.tryapplyd b m Inttriplefunc.undefined 
 12.1084 -   in Monomialfunc.update (m,pi_equation_add (pi_equation_cmul c e) es) b 
 12.1085 -   end) q a) p acc ;
 12.1086 -
 12.1087 -(* Usual operations on equation-parametrized poly.                           *)
 12.1088 -
 12.1089 -fun pi_epoly_cmul c l =
 12.1090 -  if c =/ rat_0 then Inttriplefunc.undefined else Inttriplefunc.mapf (pi_equation_cmul c) l;;
 12.1091 -
 12.1092 -val pi_epoly_neg = pi_epoly_cmul (Rat.rat_of_int ~1);
 12.1093 -
 12.1094 -val pi_epoly_add = Inttriplefunc.combine pi_equation_add Inttriplefunc.is_undefined;
 12.1095 -
 12.1096 -fun pi_epoly_sub p q = pi_epoly_add p (pi_epoly_neg q);;
 12.1097 -
 12.1098 -fun allpairs f l1 l2 =  fold_rev (fn x => (curry (op @)) (map (f x) l2)) l1 [];
 12.1099 -
 12.1100 -(* Hence produce the "relevant" monomials: those whose squares lie in the    *)
 12.1101 -(* Newton polytope of the monomials in the input. (This is enough according  *)
 12.1102 -(* to Reznik: "Extremal PSD forms with few terms", Duke Math. Journal,       *)
 12.1103 -(* vol 45, pp. 363--374, 1978.                                               *)
 12.1104 -(*                                                                           *)
 12.1105 -(* These are ordered in sort of decreasing degree. In particular the         *)
 12.1106 -(* constant monomial is last; this gives an order in diagonalization of the  *)
 12.1107 -(* quadratic form that will tend to display constants.                       *)
 12.1108 -
 12.1109 -(*
 12.1110 -UNUSED
 12.1111 -
 12.1112 -fun newton_polytope pol =
 12.1113 - let 
 12.1114 -  val vars = poly_variables pol
 12.1115 -  val mons = map (fn m => map (fn x => monomial_degree x m) vars) 
 12.1116 -             (Monomialfunc.dom pol)
 12.1117 -  val ds = map (fn x => (degree x pol + 1) div 2) vars
 12.1118 -  val all = fold_rev (fn n => allpairs cons (0 upto n)) ds [[]]
 12.1119 -  val mons' = minimal_convex_hull mons
 12.1120 -  val all' =
 12.1121 -    filter (fn m => in_convex_hull mons' (map (fn x => 2 * x) m)) all 
 12.1122 - in map (fn m => fold_rev2 (fn v => fn i => fn a => if i = 0 then a else Ctermfunc.update (v,i) a)
 12.1123 -                        vars m monomial_1) (rev all')
 12.1124 - end;
 12.1125 -
 12.1126 -*)
 12.1127 -
 12.1128 -(* Diagonalize (Cholesky/LDU) the matrix corresponding to a quadratic form.  *)
 12.1129 -
 12.1130 -local
 12.1131 -fun diagonalize n i m =
 12.1132 - if Intpairfunc.is_undefined (snd m) then [] 
 12.1133 - else
 12.1134 -  let val a11 = Intpairfunc.tryapplyd (snd m) (i,i) rat_0 
 12.1135 -  in if a11 </ rat_0 then error "diagonalize: not PSD"
 12.1136 -    else if a11 =/ rat_0 then
 12.1137 -          if Intfunc.is_undefined (snd (row i m)) then diagonalize n (i + 1) m
 12.1138 -          else error "diagonalize: not PSD ___ "
 12.1139 -    else
 12.1140 -     let 
 12.1141 -      val v = row i m
 12.1142 -      val v' = (fst v, Intfunc.fold (fn (i, c) => fn a => 
 12.1143 -       let val y = c // a11 
 12.1144 -       in if y = rat_0 then a else Intfunc.update (i,y) a 
 12.1145 -       end)  (snd v) Intfunc.undefined)
 12.1146 -      fun upt0 x y a = if y = rat_0 then a else Intpairfunc.update (x,y) a
 12.1147 -      val m' =
 12.1148 -      ((n,n),
 12.1149 -      iter (i+1,n) (fn j =>
 12.1150 -          iter (i+1,n) (fn k =>
 12.1151 -              (upt0 (j,k) (Intpairfunc.tryapplyd (snd m) (j,k) rat_0 -/ Intfunc.tryapplyd (snd v) j rat_0 */ Intfunc.tryapplyd (snd v') k rat_0))))
 12.1152 -          Intpairfunc.undefined)
 12.1153 -     in (a11,v')::diagonalize n (i + 1) m' 
 12.1154 -     end
 12.1155 -  end
 12.1156 -in
 12.1157 -fun diag m =
 12.1158 - let 
 12.1159 -   val nn = dimensions m 
 12.1160 -   val n = fst nn 
 12.1161 - in if snd nn <> n then error "diagonalize: non-square matrix" 
 12.1162 -    else diagonalize n 1 m
 12.1163 - end
 12.1164 -end;
 12.1165 -
 12.1166 -fun gcd_rat a b = Rat.rat_of_int (Integer.gcd (int_of_rat a) (int_of_rat b));
 12.1167 -
 12.1168 -(* Adjust a diagonalization to collect rationals at the start.               *)
 12.1169 -  (* FIXME : Potentially polymorphic keys, but here only: integers!! *)
 12.1170 -local
 12.1171 - fun upd0 x y a = if y =/ rat_0 then a else Intfunc.update(x,y) a;
 12.1172 - fun mapa f (d,v) = 
 12.1173 -  (d, Intfunc.fold (fn (i,c) => fn a => upd0 i (f c) a) v Intfunc.undefined)
 12.1174 - fun adj (c,l) =
 12.1175 - let val a = 
 12.1176 -  Intfunc.fold (fn (i,c) => fn a => lcm_rat a (denominator_rat c)) 
 12.1177 -    (snd l) rat_1 //
 12.1178 -  Intfunc.fold (fn (i,c) => fn a => gcd_rat a (numerator_rat c)) 
 12.1179 -    (snd l) rat_0
 12.1180 -  in ((c // (a */ a)),mapa (fn x => a */ x) l)
 12.1181 -  end
 12.1182 -in
 12.1183 -fun deration d = if null d then (rat_0,d) else
 12.1184 - let val d' = map adj d
 12.1185 -     val a = fold (lcm_rat o denominator_rat o fst) d' rat_1 //
 12.1186 -          fold (gcd_rat o numerator_rat o fst) d' rat_0 
 12.1187 - in ((rat_1 // a),map (fn (c,l) => (a */ c,l)) d')
 12.1188 - end
 12.1189 -end;
 12.1190 - 
 12.1191 -(* Enumeration of monomials with given multidegree bound.                    *)
 12.1192 -
 12.1193 -fun enumerate_monomials d vars = 
 12.1194 - if d < 0 then []
 12.1195 - else if d = 0 then [Ctermfunc.undefined]
 12.1196 - else if null vars then [monomial_1] else
 12.1197 - let val alts =
 12.1198 -  map (fn k => let val oths = enumerate_monomials (d - k) (tl vars) 
 12.1199 -               in map (fn ks => if k = 0 then ks else Ctermfunc.update (hd vars, k) ks) oths end) (0 upto d) 
 12.1200 - in fold1 (curry op @) alts
 12.1201 - end;
 12.1202 -
 12.1203 -(* Enumerate products of distinct input polys with degree <= d.              *)
 12.1204 -(* We ignore any constant input polynomials.                                 *)
 12.1205 -(* Give the output polynomial and a record of how it was derived.            *)
 12.1206 -
 12.1207 -local
 12.1208 - open RealArith
 12.1209 -in
 12.1210 -fun enumerate_products d pols =
 12.1211 -if d = 0 then [(poly_const rat_1,Rational_lt rat_1)] 
 12.1212 -else if d < 0 then [] else
 12.1213 -case pols of 
 12.1214 -   [] => [(poly_const rat_1,Rational_lt rat_1)]
 12.1215 - | (p,b)::ps => 
 12.1216 -    let val e = multidegree p 
 12.1217 -    in if e = 0 then enumerate_products d ps else
 12.1218 -       enumerate_products d ps @
 12.1219 -       map (fn (q,c) => (poly_mul p q,Product(b,c)))
 12.1220 -         (enumerate_products (d - e) ps)
 12.1221 -    end
 12.1222 -end;
 12.1223 -
 12.1224 -(* Convert regular polynomial. Note that we treat (0,0,0) as -1.             *)
 12.1225 -
 12.1226 -fun epoly_of_poly p =
 12.1227 -  Monomialfunc.fold (fn (m,c) => fn a => Monomialfunc.update (m, Inttriplefunc.onefunc ((0,0,0), Rat.neg c)) a) p Monomialfunc.undefined;
 12.1228 -
 12.1229 -(* String for block diagonal matrix numbered k.                              *)
 12.1230 -
 12.1231 -fun sdpa_of_blockdiagonal k m =
 12.1232 - let 
 12.1233 -  val pfx = string_of_int k ^" "
 12.1234 -  val ents =
 12.1235 -    Inttriplefunc.fold 
 12.1236 -      (fn ((b,i,j),c) => fn a => if i > j then a else ((b,i,j),c)::a) 
 12.1237 -      m [] 
 12.1238 -  val entss = sort (increasing fst triple_int_ord) ents 
 12.1239 - in fold_rev (fn ((b,i,j),c) => fn a =>
 12.1240 -     pfx ^ string_of_int b ^ " " ^ string_of_int i ^ " " ^ string_of_int j ^
 12.1241 -     " " ^ decimalize 20 c ^ "\n" ^ a) entss ""
 12.1242 - end;
 12.1243 -
 12.1244 -(* SDPA for problem using block diagonal (i.e. multiple SDPs)                *)
 12.1245 -
 12.1246 -fun sdpa_of_blockproblem nblocks blocksizes obj mats =
 12.1247 - let val m = length mats - 1 
 12.1248 - in
 12.1249 -  string_of_int m ^ "\n" ^
 12.1250 -  string_of_int nblocks ^ "\n" ^
 12.1251 -  (fold1 (fn s => fn t => s^" "^t) (map string_of_int blocksizes)) ^
 12.1252 -  "\n" ^
 12.1253 -  sdpa_of_vector obj ^
 12.1254 -  fold_rev2 (fn k => fn m => fn a => sdpa_of_blockdiagonal (k - 1) m ^ a)
 12.1255 -    (1 upto length mats) mats ""
 12.1256 - end;
 12.1257 -
 12.1258 -(* Run prover on a problem in block diagonal form.                       *)
 12.1259 -
 12.1260 -fun run_blockproblem prover nblocks blocksizes obj mats=
 12.1261 -  parse_csdpoutput (prover (sdpa_of_blockproblem nblocks blocksizes obj mats))
 12.1262 -
 12.1263 -(*
 12.1264 -UNUSED
 12.1265 -
 12.1266 -(* Hence run CSDP on a problem in block diagonal form.                       *)
 12.1267 -
 12.1268 -fun run_csdp dbg nblocks blocksizes obj mats =
 12.1269 - let 
 12.1270 -  val input_file = tmp_file "sos" ".dat-s" 
 12.1271 -  val output_file = tmp_file "sos" ".out"
 12.1272 -  val params_file = tmp_file "param" ".csdp" 
 12.1273 -  val _ = File.write input_file
 12.1274 -   (sdpa_of_blockproblem "" nblocks blocksizes obj mats)
 12.1275 -  val _ = File.write params_file csdp_params
 12.1276 -  val current_dir = File.pwd()
 12.1277 -  val _ = File.cd (Path.variable "ISABELLE_TMP")
 12.1278 -  val rv = system ("csdp "^(Path.implode input_file) ^ " " 
 12.1279 -                   ^ (Path.implode output_file) ^
 12.1280 -                   (if dbg then "" else "> /dev/null"))
 12.1281 -  val  opr = File.read output_file 
 12.1282 -  val res = parse_csdpoutput opr 
 12.1283 - in
 12.1284 -   ((if dbg then ()
 12.1285 -     else (File.rm input_file ; File.rm output_file ; File.cd current_dir));
 12.1286 -    (rv,res))
 12.1287 - end;
 12.1288 -
 12.1289 -fun csdp nblocks blocksizes obj mats =
 12.1290 - let 
 12.1291 -  val (rv,res) = run_csdp (!debugging) nblocks blocksizes obj mats 
 12.1292 - in ((if rv = 1 orelse rv = 2 then error "csdp: Problem is infeasible"
 12.1293 -     else if rv = 3 then writeln "csdp warning: Reduced accuracy"
 12.1294 -     else if rv <> 0 then error ("csdp: error "^string_of_int rv)
 12.1295 -     else ());
 12.1296 -     res)
 12.1297 - end;
 12.1298 -*)
 12.1299 -
 12.1300 -(* 3D versions of matrix operations to consider blocks separately.           *)
 12.1301 -
 12.1302 -val bmatrix_add = Inttriplefunc.combine (curry op +/) (fn x => x =/ rat_0);
 12.1303 -fun bmatrix_cmul c bm =
 12.1304 -  if c =/ rat_0 then Inttriplefunc.undefined
 12.1305 -  else Inttriplefunc.mapf (fn x => c */ x) bm;
 12.1306 -
 12.1307 -val bmatrix_neg = bmatrix_cmul (Rat.rat_of_int ~1);
 12.1308 -fun bmatrix_sub m1 m2 = bmatrix_add m1 (bmatrix_neg m2);;
 12.1309 -
 12.1310 -(* Smash a block matrix into components.                                     *)
 12.1311 -
 12.1312 -fun blocks blocksizes bm =
 12.1313 - map (fn (bs,b0) =>
 12.1314 -      let val m = Inttriplefunc.fold
 12.1315 -          (fn ((b,i,j),c) => fn a => if b = b0 then Intpairfunc.update ((i,j),c) a else a) bm Intpairfunc.undefined
 12.1316 -          val d = Intpairfunc.fold (fn ((i,j),c) => fn a => max a (max i j)) m 0 
 12.1317 -      in (((bs,bs),m):matrix) end)
 12.1318 - (blocksizes ~~ (1 upto length blocksizes));;
 12.1319 -
 12.1320 -(* FIXME : Get rid of this !!!*)
 12.1321 -local
 12.1322 -  fun tryfind_with msg f [] = error msg
 12.1323 -    | tryfind_with msg f (x::xs) = (f x handle ERROR s => tryfind_with s f xs);
 12.1324 -in 
 12.1325 -  fun tryfind f = tryfind_with "tryfind" f
 12.1326 -end
 12.1327 -
 12.1328 -(*
 12.1329 -fun tryfind f [] = error "tryfind"
 12.1330 -  | tryfind f (x::xs) = (f x handle ERROR _ => tryfind f xs);
 12.1331 -*)
 12.1332 -
 12.1333 -(* Positiv- and Nullstellensatz. Flag "linf" forces a linear representation. *)
 12.1334 -
 12.1335 - 
 12.1336 -local
 12.1337 - open RealArith
 12.1338 -in
 12.1339 -fun real_positivnullstellensatz_general prover linf d eqs leqs pol =
 12.1340 -let 
 12.1341 - val vars = fold_rev (curry (gen_union (op aconvc)) o poly_variables) 
 12.1342 -              (pol::eqs @ map fst leqs) []
 12.1343 - val monoid = if linf then 
 12.1344 -      (poly_const rat_1,Rational_lt rat_1)::
 12.1345 -      (filter (fn (p,c) => multidegree p <= d) leqs)
 12.1346 -    else enumerate_products d leqs
 12.1347 - val nblocks = length monoid
 12.1348 - fun mk_idmultiplier k p =
 12.1349 -  let 
 12.1350 -   val e = d - multidegree p
 12.1351 -   val mons = enumerate_monomials e vars
 12.1352 -   val nons = mons ~~ (1 upto length mons) 
 12.1353 -  in (mons,
 12.1354 -      fold_rev (fn (m,n) => Monomialfunc.update(m,Inttriplefunc.onefunc((~k,~n,n),rat_1))) nons Monomialfunc.undefined)
 12.1355 -  end
 12.1356 -
 12.1357 - fun mk_sqmultiplier k (p,c) =
 12.1358 -  let 
 12.1359 -   val e = (d - multidegree p) div 2
 12.1360 -   val mons = enumerate_monomials e vars
 12.1361 -   val nons = mons ~~ (1 upto length mons) 
 12.1362 -  in (mons, 
 12.1363 -      fold_rev (fn (m1,n1) =>
 12.1364 -       fold_rev (fn (m2,n2) => fn  a =>
 12.1365 -        let val m = monomial_mul m1 m2 
 12.1366 -        in if n1 > n2 then a else
 12.1367 -          let val c = if n1 = n2 then rat_1 else rat_2
 12.1368 -              val e = Monomialfunc.tryapplyd a m Inttriplefunc.undefined 
 12.1369 -          in Monomialfunc.update(m, tri_equation_add (Inttriplefunc.onefunc((k,n1,n2), c)) e) a
 12.1370 -          end
 12.1371 -        end)  nons)
 12.1372 -       nons Monomialfunc.undefined)
 12.1373 -  end
 12.1374 -
 12.1375 -  val (sqmonlist,sqs) = split_list (map2 mk_sqmultiplier (1 upto length monoid) monoid)
 12.1376 -  val (idmonlist,ids) =  split_list(map2 mk_idmultiplier (1 upto length eqs) eqs)
 12.1377 -  val blocksizes = map length sqmonlist
 12.1378 -  val bigsum =
 12.1379 -    fold_rev2 (fn p => fn q => fn a => tri_epoly_pmul p q a) eqs ids
 12.1380 -            (fold_rev2 (fn (p,c) => fn s => fn a => tri_epoly_pmul p s a) monoid sqs
 12.1381 -                     (epoly_of_poly(poly_neg pol)))
 12.1382 -  val eqns = Monomialfunc.fold (fn (m,e) => fn a => e::a) bigsum []
 12.1383 -  val (pvs,assig) = tri_eliminate_all_equations (0,0,0) eqns
 12.1384 -  val qvars = (0,0,0)::pvs
 12.1385 -  val allassig = fold_rev (fn v => Inttriplefunc.update(v,(Inttriplefunc.onefunc(v,rat_1)))) pvs assig
 12.1386 -  fun mk_matrix v =
 12.1387 -    Inttriplefunc.fold (fn ((b,i,j), ass) => fn m => 
 12.1388 -        if b < 0 then m else
 12.1389 -         let val c = Inttriplefunc.tryapplyd ass v rat_0
 12.1390 -         in if c = rat_0 then m else
 12.1391 -            Inttriplefunc.update ((b,j,i), c) (Inttriplefunc.update ((b,i,j), c) m)
 12.1392 -         end)
 12.1393 -          allassig Inttriplefunc.undefined
 12.1394 -  val diagents = Inttriplefunc.fold
 12.1395 -    (fn ((b,i,j), e) => fn a => if b > 0 andalso i = j then tri_equation_add e a else a)
 12.1396 -    allassig Inttriplefunc.undefined
 12.1397 -
 12.1398 -  val mats = map mk_matrix qvars
 12.1399 -  val obj = (length pvs,
 12.1400 -            itern 1 pvs (fn v => fn i => Intfunc.updatep iszero (i,Inttriplefunc.tryapplyd diagents v rat_0))
 12.1401 -                        Intfunc.undefined)
 12.1402 -  val raw_vec = if null pvs then vector_0 0
 12.1403 -                else tri_scale_then (run_blockproblem prover nblocks blocksizes) obj mats
 12.1404 -  fun int_element (d,v) i = Intfunc.tryapplyd v i rat_0
 12.1405 -  fun cterm_element (d,v) i = Ctermfunc.tryapplyd v i rat_0
 12.1406 -
 12.1407 -  fun find_rounding d =
 12.1408 -   let 
 12.1409 -    val _ = if !debugging 
 12.1410 -          then writeln ("Trying rounding with limit "^Rat.string_of_rat d ^ "\n") 
 12.1411 -          else ()
 12.1412 -    val vec = nice_vector d raw_vec
 12.1413 -    val blockmat = iter (1,dim vec)
 12.1414 -     (fn i => fn a => bmatrix_add (bmatrix_cmul (int_element vec i) (nth mats i)) a)
 12.1415 -     (bmatrix_neg (nth mats 0))
 12.1416 -    val allmats = blocks blocksizes blockmat 
 12.1417 -   in (vec,map diag allmats)
 12.1418 -   end
 12.1419 -  val (vec,ratdias) =
 12.1420 -    if null pvs then find_rounding rat_1
 12.1421 -    else tryfind find_rounding (map Rat.rat_of_int (1 upto 31) @
 12.1422 -                                map pow2 (5 upto 66))
 12.1423 -  val newassigs =
 12.1424 -    fold_rev (fn k => Inttriplefunc.update (nth pvs (k - 1), int_element vec k))
 12.1425 -           (1 upto dim vec) (Inttriplefunc.onefunc ((0,0,0), Rat.rat_of_int ~1))
 12.1426 -  val finalassigs =
 12.1427 -    Inttriplefunc.fold (fn (v,e) => fn a => Inttriplefunc.update(v, tri_equation_eval newassigs e) a) allassig newassigs
 12.1428 -  fun poly_of_epoly p =
 12.1429 -    Monomialfunc.fold (fn (v,e) => fn a => Monomialfunc.updatep iszero (v,tri_equation_eval finalassigs e) a)
 12.1430 -          p Monomialfunc.undefined
 12.1431 -  fun  mk_sos mons =
 12.1432 -   let fun mk_sq (c,m) =
 12.1433 -    (c,fold_rev (fn k=> fn a => Monomialfunc.updatep iszero (nth mons (k - 1), int_element m k) a)
 12.1434 -                 (1 upto length mons) Monomialfunc.undefined)
 12.1435 -   in map mk_sq
 12.1436 -   end
 12.1437 -  val sqs = map2 mk_sos sqmonlist ratdias
 12.1438 -  val cfs = map poly_of_epoly ids
 12.1439 -  val msq = filter (fn (a,b) => not (null b)) (map2 pair monoid sqs)
 12.1440 -  fun eval_sq sqs = fold_rev (fn (c,q) => poly_add (poly_cmul c (poly_mul q q))) sqs poly_0
 12.1441 -  val sanity =
 12.1442 -    fold_rev (fn ((p,c),s) => poly_add (poly_mul p (eval_sq s))) msq
 12.1443 -           (fold_rev2 (fn p => fn q => poly_add (poly_mul p q)) cfs eqs
 12.1444 -                    (poly_neg pol))
 12.1445 -
 12.1446 -in if not(Monomialfunc.is_undefined sanity) then raise Sanity else
 12.1447 -  (cfs,map (fn (a,b) => (snd a,b)) msq)
 12.1448 - end
 12.1449 -
 12.1450 -
 12.1451 -end;
 12.1452 -
 12.1453 -(* Iterative deepening.                                                      *)
 12.1454 -
 12.1455 -fun deepen f n = 
 12.1456 -  (writeln ("Searching with depth limit " ^ string_of_int n) ; (f n handle ERROR s => (writeln ("failed with message: " ^ s) ; deepen f (n+1))))
 12.1457 -
 12.1458 -(* The ordering so we can create canonical HOL polynomials.                  *)
 12.1459 -
 12.1460 -fun dest_monomial mon = sort (increasing fst cterm_ord) (Ctermfunc.graph mon);
 12.1461 -
 12.1462 -fun monomial_order (m1,m2) =
 12.1463 - if Ctermfunc.is_undefined m2 then LESS 
 12.1464 - else if Ctermfunc.is_undefined m1 then GREATER 
 12.1465 - else
 12.1466 -  let val mon1 = dest_monomial m1 
 12.1467 -      val mon2 = dest_monomial m2
 12.1468 -      val deg1 = fold (curry op + o snd) mon1 0
 12.1469 -      val deg2 = fold (curry op + o snd) mon2 0 
 12.1470 -  in if deg1 < deg2 then GREATER else if deg1 > deg2 then LESS
 12.1471 -     else list_ord (prod_ord cterm_ord int_ord) (mon1,mon2)
 12.1472 -  end;
 12.1473 -
 12.1474 -fun dest_poly p =
 12.1475 -  map (fn (m,c) => (c,dest_monomial m))
 12.1476 -      (sort (prod_ord monomial_order (K EQUAL)) (Monomialfunc.graph p));
 12.1477 -
 12.1478 -(* Map back polynomials and their composites to HOL.                         *)
 12.1479 -
 12.1480 -local
 12.1481 - open Thm Numeral RealArith
 12.1482 -in
 12.1483 -
 12.1484 -fun cterm_of_varpow x k = if k = 1 then x else capply (capply @{cterm "op ^ :: real => _"} x) 
 12.1485 -  (mk_cnumber @{ctyp nat} k)
 12.1486 -
 12.1487 -fun cterm_of_monomial m = 
 12.1488 - if Ctermfunc.is_undefined m then @{cterm "1::real"} 
 12.1489 - else 
 12.1490 -  let 
 12.1491 -   val m' = dest_monomial m
 12.1492 -   val vps = fold_rev (fn (x,k) => cons (cterm_of_varpow x k)) m' [] 
 12.1493 -  in fold1 (fn s => fn t => capply (capply @{cterm "op * :: real => _"} s) t) vps
 12.1494 -  end
 12.1495 -
 12.1496 -fun cterm_of_cmonomial (m,c) = if Ctermfunc.is_undefined m then cterm_of_rat c
 12.1497 -    else if c = Rat.one then cterm_of_monomial m
 12.1498 -    else capply (capply @{cterm "op *::real => _"} (cterm_of_rat c)) (cterm_of_monomial m);
 12.1499 -
 12.1500 -fun cterm_of_poly p = 
 12.1501 - if Monomialfunc.is_undefined p then @{cterm "0::real"} 
 12.1502 - else
 12.1503 -  let 
 12.1504 -   val cms = map cterm_of_cmonomial
 12.1505 -     (sort (prod_ord monomial_order (K EQUAL)) (Monomialfunc.graph p))
 12.1506 -  in fold1 (fn t1 => fn t2 => capply(capply @{cterm "op + :: real => _"} t1) t2) cms
 12.1507 -  end;
 12.1508 -
 12.1509 -fun cterm_of_sqterm (c,p) = Product(Rational_lt c,Square(cterm_of_poly p));
 12.1510 -
 12.1511 -fun cterm_of_sos (pr,sqs) = if null sqs then pr
 12.1512 -  else Product(pr,fold1 (fn a => fn b => Sum(a,b)) (map cterm_of_sqterm sqs));
 12.1513 -
 12.1514 -end
 12.1515 -
 12.1516 -(* Interface to HOL.                                                         *)
 12.1517 -local
 12.1518 -  open Thm Conv RealArith
 12.1519 -  val concl = dest_arg o cprop_of
 12.1520 -  fun simple_cterm_ord t u = TermOrd.fast_term_ord (term_of t, term_of u) = LESS
 12.1521 -in
 12.1522 -  (* FIXME: Replace tryfind by get_first !! *)
 12.1523 -fun real_nonlinear_prover prover ctxt =
 12.1524 - let 
 12.1525 -  val {add,mul,neg,pow,sub,main} =  Normalizer.semiring_normalizers_ord_wrapper ctxt
 12.1526 -      (valOf (NormalizerData.match ctxt @{cterm "(0::real) + 1"})) 
 12.1527 -     simple_cterm_ord
 12.1528 -  val (real_poly_add_conv,real_poly_mul_conv,real_poly_neg_conv,
 12.1529 -       real_poly_pow_conv,real_poly_sub_conv,real_poly_conv) = (add,mul,neg,pow,sub,main)
 12.1530 -  fun mainf  translator (eqs,les,lts) = 
 12.1531 -  let 
 12.1532 -   val eq0 = map (poly_of_term o dest_arg1 o concl) eqs
 12.1533 -   val le0 = map (poly_of_term o dest_arg o concl) les
 12.1534 -   val lt0 = map (poly_of_term o dest_arg o concl) lts
 12.1535 -   val eqp0 = map (fn (t,i) => (t,Axiom_eq i)) (eq0 ~~ (0 upto (length eq0 - 1)))
 12.1536 -   val lep0 = map (fn (t,i) => (t,Axiom_le i)) (le0 ~~ (0 upto (length le0 - 1)))
 12.1537 -   val ltp0 = map (fn (t,i) => (t,Axiom_lt i)) (lt0 ~~ (0 upto (length lt0 - 1)))
 12.1538 -   val (keq,eq) = List.partition (fn (p,_) => multidegree p = 0) eqp0
 12.1539 -   val (klep,lep) = List.partition (fn (p,_) => multidegree p = 0) lep0
 12.1540 -   val (kltp,ltp) = List.partition (fn (p,_) => multidegree p = 0) ltp0
 12.1541 -   fun trivial_axiom (p,ax) =
 12.1542 -    case ax of
 12.1543 -       Axiom_eq n => if eval Ctermfunc.undefined p <>/ Rat.zero then nth eqs n 
 12.1544 -                     else error "trivial_axiom: Not a trivial axiom"
 12.1545 -     | Axiom_le n => if eval Ctermfunc.undefined p </ Rat.zero then nth les n 
 12.1546 -                     else error "trivial_axiom: Not a trivial axiom"
 12.1547 -     | Axiom_lt n => if eval Ctermfunc.undefined p <=/ Rat.zero then nth lts n 
 12.1548 -                     else error "trivial_axiom: Not a trivial axiom"
 12.1549 -     | _ => error "trivial_axiom: Not a trivial axiom"
 12.1550 -   in 
 12.1551 -  ((let val th = tryfind trivial_axiom (keq @ klep @ kltp)
 12.1552 -   in fconv_rule (arg_conv (arg1_conv real_poly_conv) then_conv field_comp_conv) th end)
 12.1553 -   handle ERROR _ => (
 12.1554 -    let 
 12.1555 -     val pol = fold_rev poly_mul (map fst ltp) (poly_const Rat.one)
 12.1556 -     val leq = lep @ ltp
 12.1557 -     fun tryall d =
 12.1558 -      let val e = multidegree pol
 12.1559 -          val k = if e = 0 then 0 else d div e
 12.1560 -          val eq' = map fst eq 
 12.1561 -      in tryfind (fn i => (d,i,real_positivnullstellensatz_general prover false d eq' leq
 12.1562 -                            (poly_neg(poly_pow pol i))))
 12.1563 -              (0 upto k)
 12.1564 -      end
 12.1565 -    val (d,i,(cert_ideal,cert_cone)) = deepen tryall 0
 12.1566 -    val proofs_ideal =
 12.1567 -      map2 (fn q => fn (p,ax) => Eqmul(cterm_of_poly q,ax)) cert_ideal eq
 12.1568 -    val proofs_cone = map cterm_of_sos cert_cone
 12.1569 -    val proof_ne = if null ltp then Rational_lt Rat.one else
 12.1570 -      let val p = fold1 (fn s => fn t => Product(s,t)) (map snd ltp) 
 12.1571 -      in  funpow i (fn q => Product(p,q)) (Rational_lt Rat.one)
 12.1572 -      end
 12.1573 -    val proof = fold1 (fn s => fn t => Sum(s,t))
 12.1574 -                           (proof_ne :: proofs_ideal @ proofs_cone) 
 12.1575 -    in writeln "Translating proof certificate to HOL";
 12.1576 -       translator (eqs,les,lts) proof
 12.1577 -    end))
 12.1578 -   end
 12.1579 - in mainf end
 12.1580 -end
 12.1581 -
 12.1582 -fun C f x y = f y x;
 12.1583 -  (* FIXME : This is very bad!!!*)
 12.1584 -fun subst_conv eqs t = 
 12.1585 - let 
 12.1586 -  val t' = fold (Thm.cabs o Thm.lhs_of) eqs t
 12.1587 - in Conv.fconv_rule (Thm.beta_conversion true) (fold (C combination) eqs (reflexive t'))
 12.1588 - end
 12.1589 -
 12.1590 -(* A wrapper that tries to substitute away variables first.                  *)
 12.1591 -
 12.1592 -local
 12.1593 - open Thm Conv RealArith
 12.1594 -  fun simple_cterm_ord t u = TermOrd.fast_term_ord (term_of t, term_of u) = LESS
 12.1595 - val concl = dest_arg o cprop_of
 12.1596 - val shuffle1 = 
 12.1597 -   fconv_rule (rewr_conv @{lemma "(a + x == y) == (x == y - (a::real))" by (atomize (full)) (simp add: ring_simps) })
 12.1598 - val shuffle2 =
 12.1599 -    fconv_rule (rewr_conv @{lemma "(x + a == y) ==  (x == y - (a::real))" by (atomize (full)) (simp add: ring_simps)})
 12.1600 - fun substitutable_monomial fvs tm = case term_of tm of
 12.1601 -    Free(_,@{typ real}) => if not (member (op aconvc) fvs tm) then (Rat.one,tm) 
 12.1602 -                           else error "substitutable_monomial"
 12.1603 -  | @{term "op * :: real => _"}$c$(t as Free _ ) => 
 12.1604 -     if is_ratconst (dest_arg1 tm) andalso not (member (op aconvc) fvs (dest_arg tm))
 12.1605 -         then (dest_ratconst (dest_arg1 tm),dest_arg tm) else error "substitutable_monomial"
 12.1606 -  | @{term "op + :: real => _"}$s$t => 
 12.1607 -       (substitutable_monomial (add_cterm_frees (dest_arg tm) fvs) (dest_arg1 tm)
 12.1608 -        handle ERROR _ => substitutable_monomial (add_cterm_frees (dest_arg1 tm) fvs) (dest_arg tm))
 12.1609 -  | _ => error "substitutable_monomial"
 12.1610 -
 12.1611 -  fun isolate_variable v th = 
 12.1612 -   let val w = dest_arg1 (cprop_of th)
 12.1613 -   in if v aconvc w then th
 12.1614 -      else case term_of w of
 12.1615 -           @{term "op + :: real => _"}$s$t => 
 12.1616 -              if dest_arg1 w aconvc v then shuffle2 th 
 12.1617 -              else isolate_variable v (shuffle1 th)
 12.1618 -          | _ => error "isolate variable : This should not happen?"
 12.1619 -   end 
 12.1620 -in
 12.1621 -
 12.1622 -fun real_nonlinear_subst_prover prover ctxt =
 12.1623 - let 
 12.1624 -  val {add,mul,neg,pow,sub,main} =  Normalizer.semiring_normalizers_ord_wrapper ctxt
 12.1625 -      (valOf (NormalizerData.match ctxt @{cterm "(0::real) + 1"})) 
 12.1626 -     simple_cterm_ord
 12.1627 -
 12.1628 -  val (real_poly_add_conv,real_poly_mul_conv,real_poly_neg_conv,
 12.1629 -       real_poly_pow_conv,real_poly_sub_conv,real_poly_conv) = (add,mul,neg,pow,sub,main)
 12.1630 -
 12.1631 -  fun make_substitution th =
 12.1632 -   let 
 12.1633 -    val (c,v) = substitutable_monomial [] (dest_arg1(concl th))
 12.1634 -    val th1 = Drule.arg_cong_rule (capply @{cterm "op * :: real => _"} (cterm_of_rat (Rat.inv c))) (mk_meta_eq th)
 12.1635 -    val th2 = fconv_rule (binop_conv real_poly_mul_conv) th1
 12.1636 -   in fconv_rule (arg_conv real_poly_conv) (isolate_variable v th2)
 12.1637 -   end
 12.1638 -   fun oprconv cv ct = 
 12.1639 -    let val g = Thm.dest_fun2 ct
 12.1640 -    in if g aconvc @{cterm "op <= :: real => _"} 
 12.1641 -         orelse g aconvc @{cterm "op < :: real => _"} 
 12.1642 -       then arg_conv cv ct else arg1_conv cv ct
 12.1643 -    end
 12.1644 -  fun mainf translator =
 12.1645 -   let 
 12.1646 -    fun substfirst(eqs,les,lts) =
 12.1647 -      ((let 
 12.1648 -           val eth = tryfind make_substitution eqs
 12.1649 -           val modify = fconv_rule (arg_conv (oprconv(subst_conv [eth] then_conv real_poly_conv)))
 12.1650 -       in  substfirst
 12.1651 -             (filter_out (fn t => (Thm.dest_arg1 o Thm.dest_arg o cprop_of) t 
 12.1652 -                                   aconvc @{cterm "0::real"}) (map modify eqs),
 12.1653 -                                   map modify les,map modify lts)
 12.1654 -       end)
 12.1655 -       handle ERROR  _ => real_nonlinear_prover prover ctxt translator (rev eqs, rev les, rev lts))
 12.1656 -    in substfirst
 12.1657 -   end
 12.1658 -
 12.1659 -
 12.1660 - in mainf
 12.1661 - end
 12.1662 -
 12.1663 -(* Overall function. *)
 12.1664 -
 12.1665 -fun real_sos prover ctxt t = gen_prover_real_arith ctxt (real_nonlinear_subst_prover prover ctxt) t;
 12.1666 -end;
 12.1667 -
 12.1668 -(* A tactic *)
 12.1669 -fun strip_all ct = 
 12.1670 - case term_of ct of 
 12.1671 -  Const("all",_) $ Abs (xn,xT,p) => 
 12.1672 -   let val (a,(v,t')) = (apsnd (Thm.dest_abs (SOME xn)) o Thm.dest_comb) ct
 12.1673 -   in apfst (cons v) (strip_all t')
 12.1674 -   end
 12.1675 -| _ => ([],ct)
 12.1676 -
 12.1677 -fun core_sos_conv prover ctxt t = Drule.arg_cong_rule @{cterm Trueprop} (real_sos prover ctxt (Thm.dest_arg t) RS @{thm Eq_TrueI})
 12.1678 -
 12.1679 -val known_sos_constants = 
 12.1680 -  [@{term "op ==>"}, @{term "Trueprop"}, 
 12.1681 -   @{term "op -->"}, @{term "op &"}, @{term "op |"}, 
 12.1682 -   @{term "Not"}, @{term "op = :: bool => _"}, 
 12.1683 -   @{term "All :: (real => _) => _"}, @{term "Ex :: (real => _) => _"}, 
 12.1684 -   @{term "op = :: real => _"}, @{term "op < :: real => _"}, 
 12.1685 -   @{term "op <= :: real => _"}, 
 12.1686 -   @{term "op + :: real => _"}, @{term "op - :: real => _"}, 
 12.1687 -   @{term "op * :: real => _"}, @{term "uminus :: real => _"}, 
 12.1688 -   @{term "op / :: real => _"}, @{term "inverse :: real => _"},
 12.1689 -   @{term "op ^ :: real => _"}, @{term "abs :: real => _"}, 
 12.1690 -   @{term "min :: real => _"}, @{term "max :: real => _"},
 12.1691 -   @{term "0::real"}, @{term "1::real"}, @{term "number_of :: int => real"},
 12.1692 -   @{term "number_of :: int => nat"},
 12.1693 -   @{term "Int.Bit0"}, @{term "Int.Bit1"}, 
 12.1694 -   @{term "Int.Pls"}, @{term "Int.Min"}];
 12.1695 -
 12.1696 -fun check_sos kcts ct = 
 12.1697 - let
 12.1698 -  val t = term_of ct
 12.1699 -  val _ = if not (null (Term.add_tfrees t []) 
 12.1700 -                  andalso null (Term.add_tvars t [])) 
 12.1701 -          then error "SOS: not sos. Additional type varables" else ()
 12.1702 -  val fs = Term.add_frees t []
 12.1703 -  val _ = if exists (fn ((_,T)) => not (T = @{typ "real"})) fs 
 12.1704 -          then error "SOS: not sos. Variables with type not real" else ()
 12.1705 -  val vs = Term.add_vars t []
 12.1706 -  val _ = if exists (fn ((_,T)) => not (T = @{typ "real"})) fs 
 12.1707 -          then error "SOS: not sos. Variables with type not real" else ()
 12.1708 -  val ukcs = subtract (fn (t,p) => Const p aconv t) kcts (Term.add_consts t [])
 12.1709 -  val _ = if  null ukcs then () 
 12.1710 -              else error ("SOSO: Unknown constants in Subgoal:" ^ commas (map fst ukcs))
 12.1711 -in () end
 12.1712 -
 12.1713 -fun core_sos_tac prover ctxt = CSUBGOAL (fn (ct, i) => 
 12.1714 -  let val _ = check_sos known_sos_constants ct
 12.1715 -      val (avs, p) = strip_all ct
 12.1716 -      val th = standard (fold_rev forall_intr avs (real_sos prover ctxt (Thm.dest_arg p)))
 12.1717 -  in rtac th i end);
 12.1718 -
 12.1719 -fun default_SOME f NONE v = SOME v
 12.1720 -  | default_SOME f (SOME v) _ = SOME v;
 12.1721 -
 12.1722 -fun lift_SOME f NONE a = f a
 12.1723 -  | lift_SOME f (SOME a) _ = SOME a;
 12.1724 -
 12.1725 -
 12.1726 -local
 12.1727 - val is_numeral = can (HOLogic.dest_number o term_of)
 12.1728 -in
 12.1729 -fun get_denom b ct = case term_of ct of
 12.1730 -  @{term "op / :: real => _"} $ _ $ _ => 
 12.1731 -     if is_numeral (Thm.dest_arg ct) then get_denom b (Thm.dest_arg1 ct)
 12.1732 -     else default_SOME (get_denom b) (get_denom b (Thm.dest_arg ct))   (Thm.dest_arg ct, b)
 12.1733 - | @{term "op < :: real => _"} $ _ $ _ => lift_SOME (get_denom true) (get_denom true (Thm.dest_arg ct)) (Thm.dest_arg1 ct)
 12.1734 - | @{term "op <= :: real => _"} $ _ $ _ => lift_SOME (get_denom true) (get_denom true (Thm.dest_arg ct)) (Thm.dest_arg1 ct)
 12.1735 - | _ $ _ => lift_SOME (get_denom b) (get_denom b (Thm.dest_fun ct)) (Thm.dest_arg ct)
 12.1736 - | _ => NONE
 12.1737 -end;
 12.1738 -
 12.1739 -fun elim_one_denom_tac ctxt = 
 12.1740 -CSUBGOAL (fn (P,i) => 
 12.1741 - case get_denom false P of 
 12.1742 -   NONE => no_tac
 12.1743 - | SOME (d,ord) => 
 12.1744 -     let 
 12.1745 -      val ss = simpset_of ctxt addsimps @{thms field_simps} 
 12.1746 -               addsimps [@{thm nonzero_power_divide}, @{thm power_divide}]
 12.1747 -      val th = instantiate' [] [SOME d, SOME (Thm.dest_arg P)] 
 12.1748 -         (if ord then @{lemma "(d=0 --> P) & (d>0 --> P) & (d<(0::real) --> P) ==> P" by auto}
 12.1749 -          else @{lemma "(d=0 --> P) & (d ~= (0::real) --> P) ==> P" by blast})
 12.1750 -     in (rtac th i THEN Simplifier.asm_full_simp_tac ss i) end);
 12.1751 -
 12.1752 -fun elim_denom_tac ctxt i = REPEAT (elim_one_denom_tac ctxt i);
 12.1753 -
 12.1754 -fun sos_tac prover ctxt = ObjectLogic.full_atomize_tac THEN' elim_denom_tac ctxt THEN' core_sos_tac prover ctxt
 12.1755 -
 12.1756 -
 12.1757 -end;