From 3e55afd7a92e8a58f278d94fe459fda273d2e78d Mon Sep 17 00:00:00 2001 From: fbesson Date: Tue, 25 Aug 2009 20:02:48 +0000 Subject: git-svn-id: svn+ssh://scm.gforge.inria.fr/svn/coq/trunk@12294 85f007b7-540e-0410-9357-904b9bb8a0f7 --- Makefile.common | 3 +- doc/refman/Micromega.tex | 8 +- plugins/micromega/certificate.ml | 4 +- plugins/micromega/coq_micromega.ml | 6 +- plugins/micromega/csdpcert.ml | 45 +- plugins/micromega/micromega_plugin.mllib | 1 + plugins/micromega/mutils.ml | 14 +- plugins/micromega/sos.ml | 1526 +++++++++++++++--------------- plugins/micromega/sos.mli | 32 +- plugins/micromega/sos_lib.ml | 621 ++++++++++++ plugins/micromega/sos_types.ml | 68 ++ test-suite/csdp.cache | Bin 748 -> 44878 bytes test-suite/micromega/csdp.cache | Bin 48254 -> 36967 bytes 13 files changed, 1528 insertions(+), 800 deletions(-) create mode 100644 plugins/micromega/sos_lib.ml create mode 100644 plugins/micromega/sos_types.ml diff --git a/Makefile.common b/Makefile.common index ccf79331d..1a6424c95 100644 --- a/Makefile.common +++ b/Makefile.common @@ -255,7 +255,8 @@ PARSERREQUIRESCMX:=$(LINKCMX) $(PARSERCMA:.cma=.cmxa) CSDPCERTCMO:=$(addprefix plugins/micromega/, \ mutils.cmo micromega.cmo mfourier.cmo certificate.cmo \ - sos.cmo csdpcert.cmo ) + sos_types.cmo sos_lib.cmo sos.cmo csdpcert.cmo ) + CSDPCERTCMX:= $(CSDPCERTCMO:.cmo=.cmx) DEBUGPRINTERS:=dev/top_printers.cmo dev/vm_printers.cmo dev/printers.cma diff --git a/doc/refman/Micromega.tex b/doc/refman/Micromega.tex index 00e7c0413..a43cd15b0 100644 --- a/doc/refman/Micromega.tex +++ b/doc/refman/Micromega.tex @@ -8,7 +8,8 @@ Section~\ref{sec:psatz-back} presents some background explaining the proof princ % Section~\ref{sec:lia} explains how to get a complete procedure for linear integer arithmetic. -\section{The {\tt psatz} tactic in a hurry} +\asection{The {\tt psatz} tactic in a hurry} +\tacindex{psatz} \label{sec:psatz-hurry} Load the {\tt Psatz} module ({\tt Require Psatz}.). This module defines the tactics: {\tt lia}, {\tt psatzl D}, %{\tt sos D} @@ -54,7 +55,7 @@ The following table details for each domain $D \in \{\mathbb{Z},\mathbb{Q},\math \end{array} \] -\section{\emph{Positivstellensatz} refutations} +\asection{\emph{Positivstellensatz} refutations} \label{sec:psatz-back} The name {\tt psatz} is an abbreviation for \emph{positivstellensatz} -- literally positivity theorem -- which @@ -121,7 +122,8 @@ refutation. %% This amounts to searching for $p$ in the cone without generators \emph{i.e.}, $Cone(\{\})$. % -\section{ {\tt lia} : the linear integer arithmetic tactic } +\asection{ {\tt lia} : the linear integer arithmetic tactic } +\tacindex{lia} \label{sec:lia} The tactic {\tt lia} offers an alternative to the {\tt omega} and {\tt romega} tactic (see diff --git a/plugins/micromega/certificate.ml b/plugins/micromega/certificate.ml index 229b1d0e1..2a1c2fe22 100644 --- a/plugins/micromega/certificate.ml +++ b/plugins/micromega/certificate.ml @@ -18,6 +18,7 @@ (*open Micromega.Polynomial*) open Big_int open Num +open Sos_lib module Mc = Micromega module Ml2C = Mutils.CamlToCoq @@ -671,7 +672,8 @@ let zlinear_prover sys = let res = xzlinear_prover candidates sys in (*Printf.printf "Time prover : %f" (Sys.time () -. t0) ;*) res -open Sos +open Sos_types +open Mutils let rec scale_term t = match t with diff --git a/plugins/micromega/coq_micromega.ml b/plugins/micromega/coq_micromega.ml index 85045c661..5e13db1b6 100644 --- a/plugins/micromega/coq_micromega.ml +++ b/plugins/micromega/coq_micromega.ml @@ -1323,7 +1323,7 @@ let lift_ratproof prover l = | Some c -> Some (Mc.RatProof( c,Mc.DoneProof)) type micromega_polys = (Micromega.q Mc.pol * Mc.op1) list -type csdp_certificate = S of Sos.positivstellensatz option | F of string +type csdp_certificate = S of Sos_types.positivstellensatz option | F of string type provername = string * int option open Persistent_cache @@ -1336,14 +1336,14 @@ end) let csdp_cache = "csdp.cache" -let really_call_csdpcert : provername -> micromega_polys -> Sos.positivstellensatz option = +let really_call_csdpcert : provername -> micromega_polys -> Sos_types.positivstellensatz option = fun provername poly -> let cmdname = List.fold_left Filename.concat (Envars.coqlib ()) ["plugins"; "micromega"; "csdpcert" ^ Coq_config.exec_extension] in - match command cmdname [| cmdname |] (provername,poly) with + match ((command cmdname [| cmdname |] (provername,poly)) : csdp_certificate) with | F str -> failwith str | S res -> res diff --git a/plugins/micromega/csdpcert.ml b/plugins/micromega/csdpcert.ml index 96b369acf..78087c070 100644 --- a/plugins/micromega/csdpcert.ml +++ b/plugins/micromega/csdpcert.ml @@ -15,17 +15,24 @@ open Big_int open Num open Sos +open Sos_types +open Sos_lib + module Mc = Micromega module Ml2C = Mutils.CamlToCoq module C2Ml = Mutils.CoqToCaml type micromega_polys = (Micromega.q Mc.pol * Mc.op1) list -type csdp_certificate = S of Sos.positivstellensatz option | F of string +type csdp_certificate = S of Sos_types.positivstellensatz option | F of string type provername = string * int option -let debug = false +let debug = true +let flags = [Open_append;Open_binary;Open_creat] + +let chan = open_out_gen flags 0o666 "trace" + module M = struct @@ -58,16 +65,16 @@ let rec canonical_sum_to_string = function s -> failwith "not implemented" let print_canonical_sum m = Format.print_string (canonical_sum_to_string m) -let print_list_term l = - print_string "print_list_term\n"; - List.iter (fun (e,k) -> Printf.printf "q: %s %s ;" +let print_list_term o l = + output_string o "print_list_term\n"; + List.iter (fun (e,k) -> Printf.fprintf o "q: %s %s ;" (string_of_poly (poly_of_term (expr_to_term e))) (match k with Mc.Equal -> "= " | Mc.Strict -> "> " | Mc.NonStrict -> ">= " - | _ -> failwith "not_implemented")) l ; - print_string "\n" + | _ -> failwith "not_implemented")) (List.map (fun (e, o) -> Mc.denorm e , o) l) ; + output_string o "\n" let partition_expr l = @@ -142,7 +149,7 @@ let real_nonlinear_prover d l = (fun s t -> Sum(s,t)) (proof_ne :: proofs_ideal @ proofs_cone) in S (Some proof) with - | Sos.TooDeep -> S None + | Sos_lib.TooDeep -> S None | x -> F (Printexc.to_string x) (* This is somewhat buggy, over Z, strict inequality vanish... *) @@ -166,8 +173,8 @@ let pure_sos l = let cert = snd (cert_of_pos proof') in *) S (Some proof) with - | Not_found -> (* This is no strict inequality *) F "pure_sos cannot solve this goal" - | x -> F (Printf.sprintf "pure_sos : %s" (Printexc.to_string x)) +(* | Sos.CsdpNotFound -> F "Sos.CsdpNotFound" *) + | x -> (* May be that could be refined *) S None @@ -178,11 +185,27 @@ let run_prover prover pb = | prover, _ -> (Printf.printf "unknown prover: %s\n" prover; exit 1) +let output_csdp_certificate o = function + | S None -> output_string o "S None" + | S (Some p) -> Printf.fprintf o "S (Some %a)" output_psatz p + | F s -> Printf.fprintf o "F %s" s + + let main () = + try let (prover,poly) = (input_value stdin : provername * micromega_polys) in let cert = run_prover prover poly in +(* Printf.fprintf chan "%a -> %a" print_list_term poly output_csdp_certificate cert ; + close_out chan ; *) + output_value stdout (cert:csdp_certificate); - exit 0 ;; + flush stdout ; + Marshal.to_channel chan (cert:csdp_certificate) [] ; + flush chan ; + exit 0 + with x -> (Printf.fprintf chan "error %s" (Printexc.to_string x) ; exit 1) + +;; let _ = main () in () diff --git a/plugins/micromega/micromega_plugin.mllib b/plugins/micromega/micromega_plugin.mllib index 2bd79ce04..debc296e3 100644 --- a/plugins/micromega/micromega_plugin.mllib +++ b/plugins/micromega/micromega_plugin.mllib @@ -1,3 +1,4 @@ +Sos_types Mutils Micromega Mfourier diff --git a/plugins/micromega/mutils.ml b/plugins/micromega/mutils.ml index 9cee84946..a0158b156 100644 --- a/plugins/micromega/mutils.ml +++ b/plugins/micromega/mutils.ml @@ -360,21 +360,18 @@ module TagSet = Set.Make(Tag) let command exe_path args vl = - (* creating pipes for stdin, stdout, stderr *) let (stdin_read,stdin_write) = Unix.pipe () and (stdout_read,stdout_write) = Unix.pipe () and (stderr_read,stderr_write) = Unix.pipe () in - (* Creating channels from file descr *) - let outch = Unix.out_channel_of_descr stdin_write in - let inch = Unix.in_channel_of_descr stdout_read in (* Create the process *) let pid = Unix.create_process exe_path args stdin_read stdout_write stderr_write in (* Write the data on the stdin of the created process *) - Marshal.to_channel outch vl [Marshal.No_sharing] ; + let outch = Unix.out_channel_of_descr stdin_write in + output_value outch vl ; flush outch ; (* Wait for its completion *) @@ -384,14 +381,15 @@ let command exe_path args vl = (fun () -> (* Recover the result *) match status with - | Unix.WEXITED 0 -> Marshal.from_channel inch + | Unix.WEXITED 0 -> + let inch = Unix.in_channel_of_descr stdout_read in + begin try Marshal.from_channel inch with x -> failwith (Printf.sprintf "command \"%s\" exited %s" exe_path (Printexc.to_string x)) end | Unix.WEXITED i -> failwith (Printf.sprintf "command \"%s\" exited %i" exe_path i) | Unix.WSIGNALED i -> failwith (Printf.sprintf "command \"%s\" killed %i" exe_path i) | Unix.WSTOPPED i -> failwith (Printf.sprintf "command \"%s\" stopped %i" exe_path i)) (fun () -> (* Cleanup *) - close_out outch ; close_in inch ; - List.iter (fun x -> try Unix.close x with _ -> ()) [stdin_read; stdin_write; stderr_write; stdout_write] + List.iter (fun x -> try Unix.close x with _ -> ()) [stdin_read; stdin_write; stdout_read ; stdout_write ; stderr_read; stderr_write] ) diff --git a/plugins/micromega/sos.ml b/plugins/micromega/sos.ml index c75648af2..87e55c9e1 100644 --- a/plugins/micromega/sos.ml +++ b/plugins/micromega/sos.ml @@ -1,5 +1,5 @@ (* ========================================================================= *) -(* - This code originates from John Harrison's HOL LIGHT 2.20 *) +(* - This code originates from John Harrison's HOL LIGHT 2.30 *) (* (see file LICENSE.sos for license, copyright and disclaimer) *) (* - Laurent Théry (thery@sophia.inria.fr) has isolated the HOL *) (* independent bits *) @@ -9,9 +9,14 @@ (* ========================================================================= *) (* Nonlinear universal reals procedure using SOS decomposition. *) (* ========================================================================= *) - open Num;; open List;; +open Sos_types;; +open Sos_lib;; + +(* +prioritize_real();; +*) let debugging = ref false;; @@ -19,522 +24,6 @@ exception Sanity;; exception Unsolvable;; -(* ------------------------------------------------------------------------- *) -(* Comparisons that are reflexive on NaN and also short-circuiting. *) -(* ------------------------------------------------------------------------- *) - -let (=?) = fun x y -> Pervasives.compare x y = 0;; -let ( Pervasives.compare x y < 0;; -let (<=?) = fun x y -> Pervasives.compare x y <= 0;; -let (>?) = fun x y -> Pervasives.compare x y > 0;; -let (>=?) = fun x y -> Pervasives.compare x y >= 0;; - -(* ------------------------------------------------------------------------- *) -(* Combinators. *) -(* ------------------------------------------------------------------------- *) - -let (o) = fun f g x -> f(g x);; - -(* ------------------------------------------------------------------------- *) -(* Some useful functions on "num" type. *) -(* ------------------------------------------------------------------------- *) - - -let num_0 = Int 0 -and num_1 = Int 1 -and num_2 = Int 2 -and num_10 = Int 10;; - -let pow2 n = power_num num_2 (Int n);; -let pow10 n = power_num num_10 (Int n);; - -let numdom r = - let r' = Ratio.normalize_ratio (ratio_of_num r) in - num_of_big_int(Ratio.numerator_ratio r'), - num_of_big_int(Ratio.denominator_ratio r');; - -let numerator = (o) fst numdom -and denominator = (o) snd numdom;; - -let gcd_num n1 n2 = - num_of_big_int(Big_int.gcd_big_int (big_int_of_num n1) (big_int_of_num n2));; - -let lcm_num x y = - if x =/ num_0 & y =/ num_0 then num_0 - else abs_num((x */ y) // gcd_num x y);; - - -(* ------------------------------------------------------------------------- *) -(* List basics. *) -(* ------------------------------------------------------------------------- *) - -let rec el n l = - if n = 0 then hd l else el (n - 1) (tl l);; - - -(* ------------------------------------------------------------------------- *) -(* Various versions of list iteration. *) -(* ------------------------------------------------------------------------- *) - -let rec itlist f l b = - match l with - [] -> b - | (h::t) -> f h (itlist f t b);; - -let rec end_itlist f l = - match l with - [] -> failwith "end_itlist" - | [x] -> x - | (h::t) -> f h (end_itlist f t);; - -let rec itlist2 f l1 l2 b = - match (l1,l2) with - ([],[]) -> b - | (h1::t1,h2::t2) -> f h1 h2 (itlist2 f t1 t2 b) - | _ -> failwith "itlist2";; - -(* ------------------------------------------------------------------------- *) -(* All pairs arising from applying a function over two lists. *) -(* ------------------------------------------------------------------------- *) - -let rec allpairs f l1 l2 = - match l1 with - h1::t1 -> itlist (fun x a -> f h1 x :: a) l2 (allpairs f t1 l2) - | [] -> [];; - -(* ------------------------------------------------------------------------- *) -(* String operations (surely there is a better way...) *) -(* ------------------------------------------------------------------------- *) - -let implode l = itlist (^) l "";; - -let explode s = - let rec exap n l = - if n < 0 then l else - exap (n - 1) ((String.sub s n 1)::l) in - exap (String.length s - 1) [];; - - -(* ------------------------------------------------------------------------- *) -(* Attempting function or predicate applications. *) -(* ------------------------------------------------------------------------- *) - -let can f x = try (f x; true) with Failure _ -> false;; - - -(* ------------------------------------------------------------------------- *) -(* Repetition of a function. *) -(* ------------------------------------------------------------------------- *) - -let rec funpow n f x = - if n < 1 then x else funpow (n-1) f (f x);; - - -(* ------------------------------------------------------------------------- *) -(* term?? *) -(* ------------------------------------------------------------------------- *) - -type vname = string;; - -type term = -| Zero -| Const of Num.num -| Var of vname -| Inv of term -| Opp of term -| Add of (term * term) -| Sub of (term * term) -| Mul of (term * term) -| Div of (term * term) -| Pow of (term * int);; - - -(* ------------------------------------------------------------------------- *) -(* Data structure for Positivstellensatz refutations. *) -(* ------------------------------------------------------------------------- *) - -type positivstellensatz = - Axiom_eq of int - | Axiom_le of int - | Axiom_lt of int - | Rational_eq of num - | Rational_le of num - | Rational_lt of num - | Square of term - | Monoid of int list - | Eqmul of term * positivstellensatz - | Sum of positivstellensatz * positivstellensatz - | Product of positivstellensatz * positivstellensatz;; - - - -(* ------------------------------------------------------------------------- *) -(* Replication and sequences. *) -(* ------------------------------------------------------------------------- *) - -let rec replicate x n = - if n < 1 then [] - else x::(replicate x (n - 1));; - -let rec (--) = fun m n -> if m > n then [] else m::((m + 1) -- n);; - -(* ------------------------------------------------------------------------- *) -(* Various useful list operations. *) -(* ------------------------------------------------------------------------- *) - -let rec forall p l = - match l with - [] -> true - | h::t -> p(h) & forall p t;; - -let rec tryfind f l = - match l with - [] -> failwith "tryfind" - | (h::t) -> try f h with Failure _ -> tryfind f t;; - -let index x = - let rec ind n l = - match l with - [] -> failwith "index" - | (h::t) -> if x =? h then n else ind (n + 1) t in - ind 0;; - -(* ------------------------------------------------------------------------- *) -(* "Set" operations on lists. *) -(* ------------------------------------------------------------------------- *) - -let rec mem x lis = - match lis with - [] -> false - | (h::t) -> x =? h or mem x t;; - -let insert x l = - if mem x l then l else x::l;; - -let union l1 l2 = itlist insert l1 l2;; - -let subtract l1 l2 = filter (fun x -> not (mem x l2)) l1;; - -(* ------------------------------------------------------------------------- *) -(* Merging and bottom-up mergesort. *) -(* ------------------------------------------------------------------------- *) - -let rec merge ord l1 l2 = - match l1 with - [] -> l2 - | h1::t1 -> match l2 with - [] -> l1 - | h2::t2 -> if ord h1 h2 then h1::(merge ord t1 l2) - else h2::(merge ord l1 t2);; - - -(* ------------------------------------------------------------------------- *) -(* Common measure predicates to use with "sort". *) -(* ------------------------------------------------------------------------- *) - -let increasing f x y = f x ? f y;; - -(* ------------------------------------------------------------------------- *) -(* Zipping, unzipping etc. *) -(* ------------------------------------------------------------------------- *) - -let rec zip l1 l2 = - match (l1,l2) with - ([],[]) -> [] - | (h1::t1,h2::t2) -> (h1,h2)::(zip t1 t2) - | _ -> failwith "zip";; - -let rec unzip = - function [] -> [],[] - | ((a,b)::rest) -> let alist,blist = unzip rest in - (a::alist,b::blist);; - -(* ------------------------------------------------------------------------- *) -(* Iterating functions over lists. *) -(* ------------------------------------------------------------------------- *) - -let rec do_list f l = - match l with - [] -> () - | (h::t) -> (f h; do_list f t);; - -(* ------------------------------------------------------------------------- *) -(* Sorting. *) -(* ------------------------------------------------------------------------- *) - -let rec sort cmp lis = - match lis with - [] -> [] - | piv::rest -> - let r,l = partition (cmp piv) rest in - (sort cmp l) @ (piv::(sort cmp r));; - -(* ------------------------------------------------------------------------- *) -(* Removing adjacent (NB!) equal elements from list. *) -(* ------------------------------------------------------------------------- *) - -let rec uniq l = - match l with - x::(y::_ as t) -> let t' = uniq t in - if x =? y then t' else - if t'==t then l else x::t' - | _ -> l;; - -(* ------------------------------------------------------------------------- *) -(* Convert list into set by eliminating duplicates. *) -(* ------------------------------------------------------------------------- *) - -let setify s = uniq (sort (<=?) s);; - -(* ------------------------------------------------------------------------- *) -(* Polymorphic finite partial functions via Patricia trees. *) -(* *) -(* The point of this strange representation is that it is canonical (equal *) -(* functions have the same encoding) yet reasonably efficient on average. *) -(* *) -(* Idea due to Diego Olivier Fernandez Pons (OCaml list, 2003/11/10). *) -(* ------------------------------------------------------------------------- *) - -type ('a,'b)func = - Empty - | Leaf of int * ('a*'b)list - | Branch of int * int * ('a,'b)func * ('a,'b)func;; - -(* ------------------------------------------------------------------------- *) -(* Undefined function. *) -(* ------------------------------------------------------------------------- *) - -let undefined = Empty;; - -(* ------------------------------------------------------------------------- *) -(* In case of equality comparison worries, better use this. *) -(* ------------------------------------------------------------------------- *) - -let is_undefined f = - match f with - Empty -> true - | _ -> false;; - -(* ------------------------------------------------------------------------- *) -(* Operation analagous to "map" for lists. *) -(* ------------------------------------------------------------------------- *) - -let mapf = - let rec map_list f l = - match l with - [] -> [] - | (x,y)::t -> (x,f(y))::(map_list f t) in - let rec mapf f t = - match t with - Empty -> Empty - | Leaf(h,l) -> Leaf(h,map_list f l) - | Branch(p,b,l,r) -> Branch(p,b,mapf f l,mapf f r) in - mapf;; - -(* ------------------------------------------------------------------------- *) -(* Operations analogous to "fold" for lists. *) -(* ------------------------------------------------------------------------- *) - -let foldl = - let rec foldl_list f a l = - match l with - [] -> a - | (x,y)::t -> foldl_list f (f a x y) t in - let rec foldl f a t = - match t with - Empty -> a - | Leaf(h,l) -> foldl_list f a l - | Branch(p,b,l,r) -> foldl f (foldl f a l) r in - foldl;; - -let foldr = - let rec foldr_list f l a = - match l with - [] -> a - | (x,y)::t -> f x y (foldr_list f t a) in - let rec foldr f t a = - match t with - Empty -> a - | Leaf(h,l) -> foldr_list f l a - | Branch(p,b,l,r) -> foldr f l (foldr f r a) in - foldr;; - -(* ------------------------------------------------------------------------- *) -(* Redefinition and combination. *) -(* ------------------------------------------------------------------------- *) - -let (|->),combine = - let ldb x y = let z = x lxor y in z land (-z) in - let newbranch p1 t1 p2 t2 = - let b = ldb p1 p2 in - let p = p1 land (b - 1) in - if p1 land b = 0 then Branch(p,b,t1,t2) - else Branch(p,b,t2,t1) in - let rec define_list (x,y as xy) l = - match l with - (a,b as ab)::t -> - if x =? a then xy::t - else if x [xy] - and combine_list op z l1 l2 = - match (l1,l2) with - [],_ -> l2 - | _,[] -> l1 - | ((x1,y1 as xy1)::t1,(x2,y2 as xy2)::t2) -> - if x1 ) x y = - let k = Hashtbl.hash x in - let rec upd t = - match t with - Empty -> Leaf (k,[x,y]) - | Leaf(h,l) -> - if h = k then Leaf(h,define_list (x,y) l) - else newbranch h t k (Leaf(k,[x,y])) - | Branch(p,b,l,r) -> - if k land (b - 1) <> p then newbranch p t k (Leaf(k,[x,y])) - else if k land b = 0 then Branch(p,b,upd l,r) - else Branch(p,b,l,upd r) in - upd in - let rec combine op z t1 t2 = - match (t1,t2) with - Empty,_ -> t2 - | _,Empty -> t1 - | Leaf(h1,l1),Leaf(h2,l2) -> - if h1 = h2 then - let l = combine_list op z l1 l2 in - if l = [] then Empty else Leaf(h1,l) - else newbranch h1 t1 h2 t2 - | (Leaf(k,lis) as lf),(Branch(p,b,l,r) as br) | - (Branch(p,b,l,r) as br),(Leaf(k,lis) as lf) -> - if k land (b - 1) = p then - if k land b = 0 then - let l' = combine op z lf l in - if is_undefined l' then r else Branch(p,b,l',r) - else - let r' = combine op z lf r in - if is_undefined r' then l else Branch(p,b,l,r') - else - newbranch k lf p br - | Branch(p1,b1,l1,r1),Branch(p2,b2,l2,r2) -> - if b1 < b2 then - if p2 land (b1 - 1) <> p1 then newbranch p1 t1 p2 t2 - else if p2 land b1 = 0 then - let l = combine op z l1 t2 in - if is_undefined l then r1 else Branch(p1,b1,l,r1) - else - let r = combine op z r1 t2 in - if is_undefined r then l1 else Branch(p1,b1,l1,r) - else if b2 < b1 then - if p1 land (b2 - 1) <> p2 then newbranch p1 t1 p2 t2 - else if p1 land b2 = 0 then - let l = combine op z t1 l2 in - if is_undefined l then r2 else Branch(p2,b2,l,r2) - else - let r = combine op z t1 r2 in - if is_undefined r then l2 else Branch(p2,b2,l2,r) - else if p1 = p2 then - let l = combine op z l1 l2 and r = combine op z r1 r2 in - if is_undefined l then r - else if is_undefined r then l else Branch(p1,b1,l,r) - else - newbranch p1 t1 p2 t2 in - (|->),combine;; - -(* ------------------------------------------------------------------------- *) -(* Special case of point function. *) -(* ------------------------------------------------------------------------- *) - -let (|=>) = fun x y -> (x |-> y) undefined;; - - -(* ------------------------------------------------------------------------- *) -(* Grab an arbitrary element. *) -(* ------------------------------------------------------------------------- *) - -let rec choose t = - match t with - Empty -> failwith "choose: completely undefined function" - | Leaf(h,l) -> hd l - | Branch(b,p,t1,t2) -> choose t1;; - -(* ------------------------------------------------------------------------- *) -(* Application. *) -(* ------------------------------------------------------------------------- *) - -let applyd = - let rec apply_listd l d x = - match l with - (a,b)::t -> if x =? a then b - else if x >? a then apply_listd t d x else d x - | [] -> d x in - fun f d x -> - let k = Hashtbl.hash x in - let rec look t = - match t with - Leaf(h,l) when h = k -> apply_listd l d x - | Branch(p,b,l,r) -> look (if k land b = 0 then l else r) - | _ -> d x in - look f;; - -let apply f = applyd f (fun x -> failwith "apply");; - -let tryapplyd f a d = applyd f (fun x -> d) a;; - -let defined f x = try apply f x; true with Failure _ -> false;; - -(* ------------------------------------------------------------------------- *) -(* Undefinition. *) -(* ------------------------------------------------------------------------- *) - -let undefine = - let rec undefine_list x l = - match l with - (a,b as ab)::t -> - if x =? a then t - else if x [] in - fun x -> - let k = Hashtbl.hash x in - let rec und t = - match t with - Leaf(h,l) when h = k -> - let l' = undefine_list x l in - if l' == l then t - else if l' = [] then Empty - else Leaf(h,l') - | Branch(p,b,l,r) when k land (b - 1) = p -> - if k land b = 0 then - let l' = und l in - if l' == l then t - else if is_undefined l' then r - else Branch(p,b,l',r) - else - let r' = und r in - if r' == r then t - else if is_undefined r' then l - else Branch(p,b,l,r') - | _ -> t in - und;; - - -(* ------------------------------------------------------------------------- *) -(* Mapping to sorted-list representation of the graph, domain and range. *) -(* ------------------------------------------------------------------------- *) - -let graph f = setify (foldl (fun a x y -> (x,y)::a) [] f);; - -let dom f = setify(foldl (fun a x y -> x::a) [] f);; - -let ran f = setify(foldl (fun a x y -> y::a) [] f);; - (* ------------------------------------------------------------------------- *) (* Turn a rational into a decimal string with d sig digits. *) (* ------------------------------------------------------------------------- *) @@ -554,7 +43,6 @@ let decimalize = implode(tl(explode(string_of_num k))) ^ (if e = 0 then "" else "e"^string_of_int e);; - (* ------------------------------------------------------------------------- *) (* Iterations over numbers, and lists indexed by numbers. *) (* ------------------------------------------------------------------------- *) @@ -617,7 +105,7 @@ let vector_1 = vector_const (Int 1);; let vector_cmul c (v:vector) = let n = dim v in if c =/ Int 0 then vector_0 n - else n,mapf (fun x -> c */ x) (snd v);; + else n,mapf (fun x -> c */ x) (snd v) let vector_neg (v:vector) = (fst v,mapf minus_num (snd v) :vector);; @@ -628,6 +116,12 @@ let vector_add (v1:vector) (v2:vector) = let vector_sub v1 v2 = vector_add v1 (vector_neg v2);; +let vector_dot (v1:vector) (v2:vector) = + let m = dim v1 and n = dim v2 in + if m <> n then failwith "vector_add: incompatible dimensions" else + foldl (fun a i x -> x +/ a) (Int 0) + (combine ( */ ) (fun x -> x =/ Int 0) (snd v1) (snd v2));; + let vector_of_list l = let n = length l in (n,itlist2 (|->) (1--n) l undefined :vector);; @@ -789,13 +283,13 @@ let poly_variables (p:poly) = (* Order monomials for human presentation. *) (* ------------------------------------------------------------------------- *) -let humanorder_varpow (x1,k1) (x2,k2) = x1 < x2 or (x1 = x2 & k1 > k2);; +let humanorder_varpow (x1,k1) (x2,k2) = x1 < x2 or x1 = x2 & k1 > k2;; let humanorder_monomial = let rec ord l1 l2 = match (l1,l2) with _,[] -> true | [],_ -> false - | h1::t1,h2::t2 -> humanorder_varpow h1 h2 or (h1 = h2 & ord t1 t2) in + | h1::t1,h2::t2 -> humanorder_varpow h1 h2 or h1 = h2 & ord t1 t2 in fun m1 m2 -> m1 = m2 or ord (sort humanorder_varpow (graph m1)) (sort humanorder_varpow (graph m2));; @@ -886,7 +380,7 @@ let print_poly m = Format.print_string(string_of_poly m);; *) (* ------------------------------------------------------------------------- *) -(* Conversion from term. *) +(* Conversion from term. *) (* ------------------------------------------------------------------------- *) let rec poly_of_term t = match t with @@ -914,7 +408,7 @@ let rec poly_of_term t = match t with let sdpa_of_vector (v:vector) = let n = dim v in - let strs = map (o (decimalize 20) (element v)) (1--n) in + let strs = map (o (decimalize 20) (element v)) (1--n) in end_itlist (fun x y -> x ^ " " ^ y) strs ^ "\n";; (* ------------------------------------------------------------------------- *) @@ -965,103 +459,28 @@ let sdpa_of_problem comment obj mats = (* More parser basics. *) (* ------------------------------------------------------------------------- *) -exception Noparse;; - - -let isspace,issep,isbra,issymb,isalpha,isnum,isalnum = - let charcode s = Char.code(String.get s 0) in - let spaces = " \t\n\r" - and separators = ",;" - and brackets = "()[]{}" - and symbs = "\\!@#$%^&*-+|\\<=>/?~.:" - and alphas = "'abcdefghijklmnopqrstuvwxyz_ABCDEFGHIJKLMNOPQRSTUVWXYZ" - and nums = "0123456789" in - let allchars = spaces^separators^brackets^symbs^alphas^nums in - let csetsize = itlist ((o) max charcode) (explode allchars) 256 in - let ctable = Array.make csetsize 0 in - do_list (fun c -> Array.set ctable (charcode c) 1) (explode spaces); - do_list (fun c -> Array.set ctable (charcode c) 2) (explode separators); - do_list (fun c -> Array.set ctable (charcode c) 4) (explode brackets); - do_list (fun c -> Array.set ctable (charcode c) 8) (explode symbs); - do_list (fun c -> Array.set ctable (charcode c) 16) (explode alphas); - do_list (fun c -> Array.set ctable (charcode c) 32) (explode nums); - let isspace c = Array.get ctable (charcode c) = 1 - and issep c = Array.get ctable (charcode c) = 2 - and isbra c = Array.get ctable (charcode c) = 4 - and issymb c = Array.get ctable (charcode c) = 8 - and isalpha c = Array.get ctable (charcode c) = 16 - and isnum c = Array.get ctable (charcode c) = 32 - and isalnum c = Array.get ctable (charcode c) >= 16 in - isspace,issep,isbra,issymb,isalpha,isnum,isalnum;; - -let (||) parser1 parser2 input = - try parser1 input - with Noparse -> parser2 input;; - -let (++) parser1 parser2 input = - let result1,rest1 = parser1 input in - let result2,rest2 = parser2 rest1 in - (result1,result2),rest2;; - -let rec many prs input = - try let result,next = prs input in - let results,rest = many prs next in - (result::results),rest - with Noparse -> [],input;; - -let (>>) prs treatment input = - let result,rest = prs input in - treatment(result),rest;; - -let fix err prs input = - try prs input - with Noparse -> failwith (err ^ " expected");; - -let rec listof prs sep err = - prs ++ many (sep ++ fix err prs >> snd) >> (fun (h,t) -> h::t);; - -let possibly prs input = - try let x,rest = prs input in [x],rest - with Noparse -> [],input;; - -let some p = - function - [] -> raise Noparse - | (h::t) -> if p h then (h,t) else raise Noparse;; - -let a tok = some (fun item -> item = tok);; - -let rec atleast n prs i = - (if n <= 0 then many prs - else prs ++ atleast (n - 1) prs >> (fun (h,t) -> h::t)) i;; - -let finished input = - if input = [] then 0,input else failwith "Unparsed input";; - let word s = end_itlist (fun p1 p2 -> (p1 ++ p2) >> (fun (s,t) -> s^t)) (map a (explode s));; - let token s = many (some isspace) ++ word s ++ many (some isspace) >> (fun ((_,t),_) -> t);; let decimal = let numeral = some isnum in - let decimalint = atleast 1 numeral >> ((o) Num.num_of_string implode) in + let decimalint = atleast 1 numeral >> ((o) Num.num_of_string implode) in let decimalfrac = atleast 1 numeral >> (fun s -> Num.num_of_string(implode s) // pow10 (length s)) in let decimalsig = decimalint ++ possibly (a "." ++ decimalfrac >> snd) - >> (function (h,[]) -> h | (h,[x]) -> h +/ x | _ -> failwith "decimalsig") in + >> (function (h,[]) -> h | (h,[x]) -> h +/ x) in let signed prs = a "-" ++ prs >> ((o) minus_num snd) || a "+" ++ prs >> snd || prs in let exponent = (a "e" || a "E") ++ signed decimalint >> snd in signed decimalsig ++ possibly exponent - >> (function (h,[]) -> h | (h,[x]) -> h */ power_num (Int 10) x | _ -> - failwith "exponent");; + >> (function (h,[]) -> h | (h,[x]) -> h */ power_num (Int 10) x);; let mkparser p s = let x,rst = p(explode s) in @@ -1073,15 +492,72 @@ let parse_decimal = mkparser decimal;; (* Parse back a vector. *) (* ------------------------------------------------------------------------- *) -let parse_csdpoutput = +let parse_sdpaoutput,parse_csdpoutput = + let vector = + token "{" ++ listof decimal (token ",") "decimal" ++ token "}" + >> (fun ((_,v),_) -> vector_of_list v) in + let parse_vector = mkparser vector in let rec skipupto dscr prs inp = (dscr ++ prs >> snd || some (fun c -> true) ++ skipupto dscr prs >> snd) inp in let ignore inp = (),[] in + let sdpaoutput = + skipupto (word "xVec" ++ token "=") + (vector ++ ignore >> fst) in let csdpoutput = (decimal ++ many(a " " ++ decimal >> snd) >> (fun (h,t) -> h::t)) ++ (a " " ++ a "\n" ++ ignore) >> ((o) vector_of_list fst) in - mkparser csdpoutput;; + mkparser sdpaoutput,mkparser csdpoutput;; + +(* ------------------------------------------------------------------------- *) +(* Also parse the SDPA output to test success (CSDP yields a return code). *) +(* ------------------------------------------------------------------------- *) + +let sdpa_run_succeeded = + let rec skipupto dscr prs inp = + (dscr ++ prs >> snd + || some (fun c -> true) ++ skipupto dscr prs >> snd) inp in + let prs = skipupto (word "phase.value" ++ token "=") + (possibly (a "p") ++ possibly (a "d") ++ + (word "OPT" || word "FEAS")) in + fun s -> try prs (explode s); true with Noparse -> false;; + +(* ------------------------------------------------------------------------- *) +(* The default parameters. Unfortunately this goes to a fixed file. *) +(* ------------------------------------------------------------------------- *) + +let sdpa_default_parameters = +"100 unsigned int maxIteration; +1.0E-7 double 0.0 < epsilonStar; +1.0E2 double 0.0 < lambdaStar; +2.0 double 1.0 < omegaStar; +-1.0E5 double lowerBound; +1.0E5 double upperBound; +0.1 double 0.0 <= betaStar < 1.0; +0.2 double 0.0 <= betaBar < 1.0, betaStar <= betaBar; +0.9 double 0.0 < gammaStar < 1.0; +1.0E-7 double 0.0 < epsilonDash; +";; + +(* ------------------------------------------------------------------------- *) +(* These were suggested by Makoto Yamashita for problems where we are *) +(* right at the edge of the semidefinite cone, as sometimes happens. *) +(* ------------------------------------------------------------------------- *) + +let sdpa_alt_parameters = +"1000 unsigned int maxIteration; +1.0E-7 double 0.0 < epsilonStar; +1.0E4 double 0.0 < lambdaStar; +2.0 double 1.0 < omegaStar; +-1.0E5 double lowerBound; +1.0E5 double upperBound; +0.1 double 0.0 <= betaStar < 1.0; +0.2 double 0.0 <= betaBar < 1.0, betaStar <= betaBar; +0.9 double 0.0 < gammaStar < 1.0; +1.0E-7 double 0.0 < epsilonDash; +";; + +let sdpa_params = sdpa_alt_parameters;; (* ------------------------------------------------------------------------- *) (* CSDP parameters; so far I'm sticking with the defaults. *) @@ -1107,47 +583,57 @@ printlevel=1 let csdp_params = csdp_default_parameters;; (* ------------------------------------------------------------------------- *) -(* The same thing with CSDP. *) +(* Now call SDPA on a problem and parse back the output. *) (* ------------------------------------------------------------------------- *) -let buffer_add_line buff line = - Buffer.add_string buff line; Buffer.add_char buff '\n' - -let string_of_file filename = - let fd = open_in filename in - let buff = Buffer.create 16 in - try while true do buffer_add_line buff (input_line fd) done; failwith "" - with End_of_file -> (close_in fd; Buffer.contents buff) - -let file_of_string filename s = - let fd = Pervasives.open_out filename in - output_string fd s; close_out fd - +let run_sdpa dbg obj mats = + let input_file = Filename.temp_file "sos" ".dat-s" in + let output_file = + String.sub input_file 0 (String.length input_file - 6) ^ ".out" + and params_file = Filename.concat (!temp_path) "param.sdpa" in + file_of_string input_file (sdpa_of_problem "" obj mats); + file_of_string params_file sdpa_params; + Sys.command("cd "^ !temp_path ^ + "; sdpa "^input_file ^ " " ^ output_file ^ + (if dbg then "" else "> /dev/null")); + let op = string_of_file output_file in + if not(sdpa_run_succeeded op) then failwith "sdpa: call failed" else + let res = parse_sdpaoutput op in + ((if dbg then () + else (Sys.remove input_file; Sys.remove output_file)); + res);; + +let sdpa obj mats = run_sdpa (!debugging) obj mats;; -exception CsdpInfeasible -exception CsdpNotFound +(* ------------------------------------------------------------------------- *) +(* The same thing with CSDP. *) +(* ------------------------------------------------------------------------- *) -let run_csdp dbg string_problem = +let run_csdp dbg obj mats = let input_file = Filename.temp_file "sos" ".dat-s" in - let output_file = Filename.temp_file "sos" ".dat-s" in - let temp_path = Filename.dirname input_file in - let params_file = Filename.concat temp_path "param.csdp" in - file_of_string input_file string_problem; + let output_file = + String.sub input_file 0 (String.length input_file - 6) ^ ".out" + and params_file = Filename.concat (!temp_path) "param.csdp" in + file_of_string input_file (sdpa_of_problem "" obj mats); file_of_string params_file csdp_params; - let rv = Sys.command("cd "^temp_path^"; csdp "^input_file^" "^output_file^ + let rv = Sys.command("cd "^(!temp_path)^"; csdp "^input_file ^ + " " ^ output_file ^ (if dbg then "" else "> /dev/null")) in - if rv = 1 or rv = 2 then raise CsdpInfeasible; - if rv = 127 then raise CsdpNotFound ; - if rv <> 0 & rv <> 3 (* reduced accuracy *) then - failwith("csdp: error "^string_of_int rv); - let string_result = string_of_file output_file in - if not dbg then - (Sys.remove input_file; Sys.remove output_file; Sys.remove params_file); - string_result + let op = string_of_file output_file in + let res = parse_csdpoutput op in + ((if dbg then () + else (Sys.remove input_file; Sys.remove output_file)); + rv,res);; let csdp obj mats = - try parse_csdpoutput (run_csdp !debugging (sdpa_of_problem "" obj mats)) - with CsdpInfeasible -> failwith "csdp: Problem is infeasible" + let rv,res = run_csdp (!debugging) obj mats in + (if rv = 1 or rv = 2 then failwith "csdp: Problem is infeasible" + else if rv = 3 then () + (* Format.print_string "csdp warning: Reduced accuracy"; + Format.print_newline() *) + else if rv <> 0 then failwith("csdp: error "^string_of_int rv) + else ()); + res;; (* ------------------------------------------------------------------------- *) (* Try some apparently sensible scaling first. Note that this is purely to *) @@ -1191,8 +677,24 @@ let linear_program_basic a = let m,n = dimensions a in let mats = map (fun j -> diagonal (column j a)) (1--n) and obj = vector_const (Int 1) m in - try ignore (run_csdp false (sdpa_of_problem "" obj mats)); true - with CsdpInfeasible -> false + let rv,res = run_csdp false obj mats in + if rv = 1 or rv = 2 then false + else if rv = 0 then true + else failwith "linear_program: An error occurred in the SDP solver";; + +(* ------------------------------------------------------------------------- *) +(* Alternative interface testing A x >= b for matrix A, vector b. *) +(* ------------------------------------------------------------------------- *) + +let linear_program a b = + let m,n = dimensions a in + if dim b <> m then failwith "linear_program: incompatible dimensions" else + let mats = diagonal b :: map (fun j -> diagonal (column j a)) (1--n) + and obj = vector_const (Int 1) m in + let rv,res = run_csdp false obj mats in + if rv = 1 or rv = 2 then false + else if rv = 0 then true + else failwith "linear_program: An error occurred in the SDP solver";; (* ------------------------------------------------------------------------- *) (* Test whether a point is in the convex hull of others. Rather than use *) @@ -1217,9 +719,7 @@ let in_convex_hull pts pt = (* ------------------------------------------------------------------------- *) let minimal_convex_hull = - let augment1 = function (m::ms) -> if in_convex_hull ms m then ms else ms@[m] - | _ -> failwith "augment1" - in + let augment1 (m::ms) = if in_convex_hull ms m then ms else ms@[m] in let augment m ms = funpow 3 augment1 (m::ms) in fun mons -> let mons' = itlist augment (tl mons) [hd mons] in @@ -1244,6 +744,7 @@ let equation_eval assig eq = (* "one" that's used for a constant term. *) (* ------------------------------------------------------------------------- *) +let failstore = ref [];; let eliminate_equations = let rec extract_first p l = @@ -1255,7 +756,7 @@ let eliminate_equations = let rec eliminate vars dun eqs = match vars with [] -> if forall is_undefined eqs then dun - else (raise Unsolvable) + else (failstore := [vars,dun,eqs]; raise Unsolvable) | v::vs -> try let eq,oeqs = extract_first (fun e -> defined e v) eqs in let a = apply eq v in @@ -1373,8 +874,8 @@ let deration d = foldl (fun a i c -> gcd_num a (numerator c)) (Int 0) (snd l) in (c // (a */ a)),mapa (fun x -> a */ x) l in let d' = map adj d in - let a = itlist ((o) lcm_num ((o) denominator fst)) d' (Int 1) // - itlist ((o) gcd_num ((o) numerator fst)) d' (Int 0) in + let a = itlist ((o) lcm_num ( (o) denominator fst)) d' (Int 1) // + itlist ((o) gcd_num ( (o) numerator fst)) d' (Int 0) in (Int 1 // a),map (fun (c,l) -> (a */ c,l)) d';; (* ------------------------------------------------------------------------- *) @@ -1426,9 +927,9 @@ let epoly_pmul p q acc = let epoly_cmul c l = if c =/ Int 0 then undefined else mapf (equation_cmul c) l;; -let epoly_neg x = epoly_cmul (Int(-1)) x;; +let epoly_neg = epoly_cmul (Int(-1));; -let epoly_add x = combine equation_add is_undefined x;; +let epoly_add = combine equation_add is_undefined;; let epoly_sub p q = epoly_add p (epoly_neg q);; @@ -1471,12 +972,33 @@ let sdpa_of_blockproblem comment nblocks blocksizes obj mats = (* Hence run CSDP on a problem in block diagonal form. *) (* ------------------------------------------------------------------------- *) -let csdp_blocks nblocks blocksizes obj mats = - let string_problem = sdpa_of_blockproblem "" nblocks blocksizes obj mats in - try parse_csdpoutput (run_csdp !debugging string_problem) - with - | CsdpInfeasible -> failwith "csdp: Problem is infeasible" - +let run_csdp dbg nblocks blocksizes obj mats = + let input_file = Filename.temp_file "sos" ".dat-s" in + let output_file = + String.sub input_file 0 (String.length input_file - 6) ^ ".out" + and params_file = Filename.concat (!temp_path) "param.csdp" in + file_of_string input_file + (sdpa_of_blockproblem "" nblocks blocksizes obj mats); + file_of_string params_file csdp_params; + let rv = Sys.command("cd "^(!temp_path)^"; csdp "^input_file ^ + " " ^ output_file ^ + (if dbg then "" else "> /dev/null")) in + let op = string_of_file output_file in + let res = parse_csdpoutput op in + ((if dbg then () + else (Sys.remove input_file; Sys.remove output_file)); + rv,res);; + +let csdp nblocks blocksizes obj mats = + let rv,res = run_csdp (!debugging) nblocks blocksizes obj mats in + (if rv = 1 or rv = 2 then failwith "csdp: Problem is infeasible" + else if rv = 3 then () + (*Format.print_string "csdp warning: Reduced accuracy"; + Format.print_newline() *) + else if rv <> 0 then failwith("csdp: error "^string_of_int rv) + else ()); + res;; + (* ------------------------------------------------------------------------- *) (* 3D versions of matrix operations to consider blocks separately. *) (* ------------------------------------------------------------------------- *) @@ -1500,7 +1022,7 @@ let blocks blocksizes bm = let m = foldl (fun a (b,i,j) c -> if b = b0 then ((i,j) |-> c) a else a) undefined bm in - (*let d = foldl (fun a (i,j) c -> max a (max i j)) 0 m in*) + let d = foldl (fun a (i,j) c -> max a (max i j)) 0 m in (((bs,bs),m):matrix)) (zip blocksizes (1--length blocksizes));; @@ -1508,9 +1030,7 @@ let blocks blocksizes bm = (* Positiv- and Nullstellensatz. Flag "linf" forces a linear representation. *) (* ------------------------------------------------------------------------- *) -let real_positivnullstellensatz_general linf d eqs leqs pol - : poly list * (positivstellensatz * (num * poly) list) list = - +let real_positivnullstellensatz_general linf d eqs leqs pol = let vars = itlist ((o) union poly_variables) (pol::eqs @ map fst leqs) [] in let monoid = if linf then @@ -1563,7 +1083,7 @@ let real_positivnullstellensatz_general linf d eqs leqs pol itern 1 pvs (fun v i -> (i |--> tryapplyd diagents v (Int 0))) undefined in let raw_vec = if pvs = [] then vector_0 0 - else scale_then (csdp_blocks nblocks blocksizes) obj mats in + else scale_then (csdp nblocks blocksizes) obj mats in let find_rounding d = (if !debugging then (Format.print_string("Trying rounding with limit "^string_of_num d); @@ -1603,24 +1123,20 @@ let real_positivnullstellensatz_general linf d eqs leqs pol (itlist2 (fun p q -> poly_add (poly_mul p q)) cfs eqs (poly_neg pol)) in if not(is_undefined sanity) then raise Sanity else - cfs,map (fun (a,b) -> snd a,b) msq;; + cfs,map (fun (a,b) -> snd a,b) msq;; +(* ------------------------------------------------------------------------- *) +(* Iterative deepening. *) +(* ------------------------------------------------------------------------- *) -let term_of_monoid l1 m = itlist (fun i m -> Mul (nth l1 i,m)) m (Const num_1) - -let rec term_of_pos l1 x = match x with - Axiom_eq i -> failwith "term_of_pos" - | Axiom_le i -> nth l1 i - | Axiom_lt i -> nth l1 i - | Monoid m -> term_of_monoid l1 m - | Rational_eq n -> Const n - | Rational_le n -> Const n - | Rational_lt n -> Const n - | Square t -> Pow (t, 2) - | Eqmul (t, y) -> Mul (t, term_of_pos l1 y) - | Sum (y, z) -> Add (term_of_pos l1 y, term_of_pos l1 z) - | Product (y, z) -> Mul (term_of_pos l1 y, term_of_pos l1 z);; +let rec deepen f n = + try print_string "Searching with depth limit "; + print_int n; print_newline(); f n + with Failure _ -> deepen f (n + 1);; +(* ------------------------------------------------------------------------- *) +(* The ordering so we can create canonical HOL polynomials. *) +(* ------------------------------------------------------------------------- *) let dest_monomial mon = sort (increasing fst) (graph mon);; @@ -1649,7 +1165,7 @@ let dest_poly p = (sort (fun (m1,_) (m2,_) -> monomial_order m1 m2) (graph p));; (* ------------------------------------------------------------------------- *) -(* Map back polynomials and their composites to term. *) +(* Map back polynomials and their composites to HOL. *) (* ------------------------------------------------------------------------- *) let term_of_varpow = @@ -1682,74 +1198,196 @@ let term_of_sos (pr,sqs) = if sqs = [] then pr else Product(pr,end_itlist (fun a b -> Sum(a,b)) (map term_of_sqterm sqs));; -let rec deepen f n = - try (*print_string "Searching with depth limit "; - print_int n; print_newline();*) f n - with Failure _ -> deepen f (n + 1);; - - - - - -exception TooDeep - -let deepen_until limit f n = - match compare limit 0 with - | 0 -> raise TooDeep - | -1 -> deepen f n - | _ -> - let rec d_until f n = - try if !debugging - then (print_string "Searching with depth limit "; - print_int n; print_newline()) ; f n - with Failure x -> - if !debugging then (Printf.printf "solver error : %s\n" x) ; - if n = limit then raise TooDeep else d_until f (n + 1) in - d_until f n - - -(* patch to remove zero polynomials from equalities. - In this case, hol light loops *) - -let real_nonlinear_prover depthmax eqs les lts = - let eq = map poly_of_term eqs - and le = map poly_of_term les - and lt = map poly_of_term lts in - let pol = itlist poly_mul lt (poly_const num_1) - and lep = map (fun (t,i) -> t,Axiom_le i) (zip le (0--(length le - 1))) - and ltp = map (fun (t,i) -> t,Axiom_lt i) (zip lt (0--(length lt - 1))) - and eqp = itlist2 (fun t i res -> - if t = undefined then res else (t,Axiom_eq i)::res) eq (0--(length eq - 1)) [] - in - - let proof = - let leq = lep @ ltp in - let eq = List.map fst eqp in - let tryall d = - let e = multidegree pol (*and pol' = poly_neg pol*) in - let k = if e = 0 then 1 else d / e in - tryfind (fun i -> d,i, - real_positivnullstellensatz_general false d eq leq - (poly_neg(poly_pow pol i))) - (0--k) in - let d,i,(cert_ideal,cert_cone) = deepen_until depthmax tryall 0 in - let proofs_ideal = - map2 (fun q i -> Eqmul(term_of_poly q,i)) - cert_ideal (List.map snd eqp) - and proofs_cone = map term_of_sos cert_cone - and proof_ne = - if lt = [] then Rational_lt num_1 else - let p = end_itlist (fun s t -> Product(s,t)) (map snd ltp) in - funpow i (fun q -> Product(p,q)) (Rational_lt num_1) in - end_itlist (fun s t -> Sum(s,t)) (proof_ne :: proofs_ideal @ proofs_cone) in - if !debugging then (print_string("Translating proof certificate to Coq"); print_newline()); - proof;; - - +(* ------------------------------------------------------------------------- *) +(* Interface to HOL. *) +(* ------------------------------------------------------------------------- *) +(* +let REAL_NONLINEAR_PROVER translator (eqs,les,lts) = + let eq0 = map (poly_of_term o lhand o concl) eqs + and le0 = map (poly_of_term o lhand o concl) les + and lt0 = map (poly_of_term o lhand o concl) lts in + let eqp0 = map (fun (t,i) -> t,Axiom_eq i) (zip eq0 (0--(length eq0 - 1))) + and lep0 = map (fun (t,i) -> t,Axiom_le i) (zip le0 (0--(length le0 - 1))) + and ltp0 = map (fun (t,i) -> t,Axiom_lt i) (zip lt0 (0--(length lt0 - 1))) in + let keq,eq = partition (fun (p,_) -> multidegree p = 0) eqp0 + and klep,lep = partition (fun (p,_) -> multidegree p = 0) lep0 + and kltp,ltp = partition (fun (p,_) -> multidegree p = 0) ltp0 in + let trivial_axiom (p,ax) = + match ax with + Axiom_eq n when eval undefined p <>/ num_0 -> el n eqs + | Axiom_le n when eval undefined p el n les + | Axiom_lt n when eval undefined p <=/ num_0 -> el n lts + | _ -> failwith "not a trivial axiom" in + try let th = tryfind trivial_axiom (keq @ klep @ kltp) in + CONV_RULE (LAND_CONV REAL_POLY_CONV THENC REAL_RAT_RED_CONV) th + with Failure _ -> + let pol = itlist poly_mul (map fst ltp) (poly_const num_1) in + let leq = lep @ ltp in + let tryall d = + let e = multidegree pol in + let k = if e = 0 then 0 else d / e in + let eq' = map fst eq in + tryfind (fun i -> d,i,real_positivnullstellensatz_general false d eq' leq + (poly_neg(poly_pow pol i))) + (0--k) in + let d,i,(cert_ideal,cert_cone) = deepen tryall 0 in + let proofs_ideal = + map2 (fun q (p,ax) -> Eqmul(term_of_poly q,ax)) cert_ideal eq + and proofs_cone = map term_of_sos cert_cone + and proof_ne = + if ltp = [] then Rational_lt num_1 else + let p = end_itlist (fun s t -> Product(s,t)) (map snd ltp) in + funpow i (fun q -> Product(p,q)) (Rational_lt num_1) in + let proof = end_itlist (fun s t -> Sum(s,t)) + (proof_ne :: proofs_ideal @ proofs_cone) in + print_string("Translating proof certificate to HOL"); + print_newline(); + translator (eqs,les,lts) proof;; +*) +(* ------------------------------------------------------------------------- *) +(* A wrapper that tries to substitute away variables first. *) +(* ------------------------------------------------------------------------- *) +(* +let REAL_NONLINEAR_SUBST_PROVER = + let zero = `&0:real` + and mul_tm = `( * ):real->real->real` + and shuffle1 = + CONV_RULE(REWR_CONV(REAL_ARITH `a + x = (y:real) <=> x = y - a`)) + and shuffle2 = + CONV_RULE(REWR_CONV(REAL_ARITH `x + a = (y:real) <=> x = y - a`)) in + let rec substitutable_monomial fvs tm = + match tm with + Var(_,Tyapp("real",[])) when not (mem tm fvs) -> Int 1,tm + | Comb(Comb(Const("real_mul",_),c),(Var(_,_) as t)) + when is_ratconst c & not (mem t fvs) + -> rat_of_term c,t + | Comb(Comb(Const("real_add",_),s),t) -> + (try substitutable_monomial (union (frees t) fvs) s + with Failure _ -> substitutable_monomial (union (frees s) fvs) t) + | _ -> failwith "substitutable_monomial" + and isolate_variable v th = + match lhs(concl th) with + x when x = v -> th + | Comb(Comb(Const("real_add",_),(Var(_,Tyapp("real",[])) as x)),t) + when x = v -> shuffle2 th + | Comb(Comb(Const("real_add",_),s),t) -> + isolate_variable v(shuffle1 th) in + let make_substitution th = + let (c,v) = substitutable_monomial [] (lhs(concl th)) in + let th1 = AP_TERM (mk_comb(mul_tm,term_of_rat(Int 1 // c))) th in + let th2 = CONV_RULE(BINOP_CONV REAL_POLY_MUL_CONV) th1 in + CONV_RULE (RAND_CONV REAL_POLY_CONV) (isolate_variable v th2) in + fun translator -> + let rec substfirst(eqs,les,lts) = + try let eth = tryfind make_substitution eqs in + let modify = + CONV_RULE(LAND_CONV(SUBS_CONV[eth] THENC REAL_POLY_CONV)) in + substfirst(filter (fun t -> lhand(concl t) <> zero) (map modify eqs), + map modify les,map modify lts) + with Failure _ -> REAL_NONLINEAR_PROVER translator (eqs,les,lts) in + substfirst;; +*) +(* ------------------------------------------------------------------------- *) +(* Overall function. *) +(* ------------------------------------------------------------------------- *) +(* +let REAL_SOS = + let init = GEN_REWRITE_CONV ONCE_DEPTH_CONV [DECIMAL] + and pure = GEN_REAL_ARITH REAL_NONLINEAR_SUBST_PROVER in + fun tm -> let th = init tm in EQ_MP (SYM th) (pure(rand(concl th)));; +*) +(* ------------------------------------------------------------------------- *) +(* Add hacks for division. *) +(* ------------------------------------------------------------------------- *) +(* +let REAL_SOSFIELD = + let inv_tm = `inv:real->real` in + let prenex_conv = + TOP_DEPTH_CONV BETA_CONV THENC + PURE_REWRITE_CONV[FORALL_SIMP; EXISTS_SIMP; real_div; + REAL_INV_INV; REAL_INV_MUL; GSYM REAL_POW_INV] THENC + NNFC_CONV THENC DEPTH_BINOP_CONV `(/\)` CONDS_CELIM_CONV THENC + PRENEX_CONV + and setup_conv = NNF_CONV THENC WEAK_CNF_CONV THENC CONJ_CANON_CONV + and core_rule t = + try REAL_ARITH t + with Failure _ -> try REAL_RING t + with Failure _ -> REAL_SOS t + and is_inv = + let is_div = is_binop `(/):real->real->real` in + fun tm -> (is_div tm or (is_comb tm & rator tm = inv_tm)) & + not(is_ratconst(rand tm)) in + let BASIC_REAL_FIELD tm = + let is_freeinv t = is_inv t & free_in t tm in + let itms = setify(map rand (find_terms is_freeinv tm)) in + let hyps = map (fun t -> SPEC t REAL_MUL_RINV) itms in + let tm' = itlist (fun th t -> mk_imp(concl th,t)) hyps tm in + let itms' = map (curry mk_comb inv_tm) itms in + let gvs = map (genvar o type_of) itms' in + let tm'' = subst (zip gvs itms') tm' in + let th1 = setup_conv tm'' in + let cjs = conjuncts(rand(concl th1)) in + let ths = map core_rule cjs in + let th2 = EQ_MP (SYM th1) (end_itlist CONJ ths) in + rev_itlist (C MP) hyps (INST (zip itms' gvs) th2) in + fun tm -> + let th0 = prenex_conv tm in + let tm0 = rand(concl th0) in + let avs,bod = strip_forall tm0 in + let th1 = setup_conv bod in + let ths = map BASIC_REAL_FIELD (conjuncts(rand(concl th1))) in + EQ_MP (SYM th0) (GENL avs (EQ_MP (SYM th1) (end_itlist CONJ ths)));; +*) +(* ------------------------------------------------------------------------- *) +(* Integer version. *) +(* ------------------------------------------------------------------------- *) +(* +let INT_SOS = + let atom_CONV = + let pth = prove + (`(~(x <= y) <=> y + &1 <= x:int) /\ + (~(x < y) <=> y <= x) /\ + (~(x = y) <=> x + &1 <= y \/ y + &1 <= x) /\ + (x < y <=> x + &1 <= y)`, + REWRITE_TAC[INT_NOT_LE; INT_NOT_LT; INT_NOT_EQ; INT_LT_DISCRETE]) in + GEN_REWRITE_CONV I [pth] + and bub_CONV = GEN_REWRITE_CONV TOP_SWEEP_CONV + [int_eq; int_le; int_lt; int_ge; int_gt; + int_of_num_th; int_neg_th; int_add_th; int_mul_th; + int_sub_th; int_pow_th; int_abs_th; int_max_th; int_min_th] in + let base_CONV = TRY_CONV atom_CONV THENC bub_CONV in + let NNF_NORM_CONV = GEN_NNF_CONV false + (base_CONV,fun t -> base_CONV t,base_CONV(mk_neg t)) in + let init_CONV = + GEN_REWRITE_CONV DEPTH_CONV [FORALL_SIMP; EXISTS_SIMP] THENC + GEN_REWRITE_CONV DEPTH_CONV [INT_GT; INT_GE] THENC + CONDS_ELIM_CONV THENC NNF_NORM_CONV in + let p_tm = `p:bool` + and not_tm = `(~)` in + let pth = TAUT(mk_eq(mk_neg(mk_neg p_tm),p_tm)) in + fun tm -> + let th0 = INST [tm,p_tm] pth + and th1 = NNF_NORM_CONV(mk_neg tm) in + let th2 = REAL_SOS(mk_neg(rand(concl th1))) in + EQ_MP th0 (EQ_MP (AP_TERM not_tm (SYM th1)) th2);; +*) +(* ------------------------------------------------------------------------- *) +(* Natural number version. *) +(* ------------------------------------------------------------------------- *) +(* +let SOS_RULE tm = + let avs = frees tm in + let tm' = list_mk_forall(avs,tm) in + let th1 = NUM_TO_INT_CONV tm' in + let th2 = INT_SOS (rand(concl th1)) in + SPECL avs (EQ_MP (SYM th1) th2);; +*) (* ------------------------------------------------------------------------- *) (* Now pure SOS stuff. *) (* ------------------------------------------------------------------------- *) +(*prioritize_real();;*) + (* ------------------------------------------------------------------------- *) (* Some combinatorial helper functions. *) (* ------------------------------------------------------------------------- *) @@ -1769,6 +1407,88 @@ let changevariables zoln pol = foldl (fun a m c -> (changevariables_monomial zoln m |-> c) a) poly_0 pol;; +(* ------------------------------------------------------------------------- *) +(* Return to original non-block matrices. *) +(* ------------------------------------------------------------------------- *) + +let sdpa_of_vector (v:vector) = + let n = dim v in + let strs = map (o (decimalize 20) (element v)) (1--n) in + end_itlist (fun x y -> x ^ " " ^ y) strs ^ "\n";; + +let sdpa_of_blockdiagonal k m = + let pfx = string_of_int k ^" " in + let ents = + foldl (fun a (b,i,j) c -> if i > j then a else ((b,i,j),c)::a) [] m in + let entss = sort (increasing fst) ents in + itlist (fun ((b,i,j),c) a -> + pfx ^ string_of_int b ^ " " ^ string_of_int i ^ " " ^ string_of_int j ^ + " " ^ decimalize 20 c ^ "\n" ^ a) entss "";; + +let sdpa_of_matrix k (m:matrix) = + let pfx = string_of_int k ^ " 1 " in + let ms = foldr (fun (i,j) c a -> if i > j then a else ((i,j),c)::a) + (snd m) [] in + let mss = sort (increasing fst) ms in + itlist (fun ((i,j),c) a -> + pfx ^ string_of_int i ^ " " ^ string_of_int j ^ + " " ^ decimalize 20 c ^ "\n" ^ a) mss "";; + +let sdpa_of_problem comment obj mats = + let m = length mats - 1 + and n,_ = dimensions (hd mats) in + "\"" ^ comment ^ "\"\n" ^ + string_of_int m ^ "\n" ^ + "1\n" ^ + string_of_int n ^ "\n" ^ + sdpa_of_vector obj ^ + itlist2 (fun k m a -> sdpa_of_matrix (k - 1) m ^ a) + (1--length mats) mats "";; + +let run_sdpa dbg obj mats = + let input_file = Filename.temp_file "sos" ".dat-s" in + let output_file = + String.sub input_file 0 (String.length input_file - 6) ^ ".out" + and params_file = Filename.concat (!temp_path) "param.sdpa" in + file_of_string input_file (sdpa_of_problem "" obj mats); + file_of_string params_file sdpa_params; + Sys.command("cd "^(!temp_path)^"; sdpa "^input_file ^ " " ^ output_file ^ + (if dbg then "" else "> /dev/null")); + let op = string_of_file output_file in + if not(sdpa_run_succeeded op) then failwith "sdpa: call failed" else + let res = parse_sdpaoutput op in + ((if dbg then () + else (Sys.remove input_file; Sys.remove output_file)); + res);; + +let sdpa obj mats = run_sdpa (!debugging) obj mats;; + +let run_csdp dbg obj mats = + let input_file = Filename.temp_file "sos" ".dat-s" in + let output_file = + String.sub input_file 0 (String.length input_file - 6) ^ ".out" + and params_file = Filename.concat (!temp_path) "param.csdp" in + file_of_string input_file (sdpa_of_problem "" obj mats); + file_of_string params_file csdp_params; + let rv = Sys.command("cd "^(!temp_path)^"; csdp "^input_file ^ + " " ^ output_file ^ + (if dbg then "" else "> /dev/null")) in + let op = string_of_file output_file in + let res = parse_csdpoutput op in + ((if dbg then () + else (Sys.remove input_file; Sys.remove output_file)); + rv,res);; + +let csdp obj mats = + let rv,res = run_csdp (!debugging) obj mats in + (if rv = 1 or rv = 2 then failwith "csdp: Problem is infeasible" + else if rv = 3 then () +(* (Format.print_string "csdp warning: Reduced accuracy"; + Format.print_newline()) *) + else if rv <> 0 then failwith("csdp: error "^string_of_int rv) + else ()); + res;; + (* ------------------------------------------------------------------------- *) (* Sum-of-squares function with some lowbrow symmetry reductions. *) (* ------------------------------------------------------------------------- *) @@ -1782,11 +1502,11 @@ let sumofsquares_general_symmetry tool pol = (fun vars' -> is_undefined(poly_sub pol (changevariables (zip vars vars') pol))) (allpermutations vars) in -(* let lpps2 = allpairs monomial_mul lpps lpps in*) -(* let lpp2_classes = + let lpps2 = allpairs monomial_mul lpps lpps in + let lpp2_classes = setify(map (fun m -> setify(map (fun vars' -> changevariables_monomial (zip vars vars') m) - invariants)) lpps2) in *) + invariants)) lpps2) in let lpns = zip lpps (1--length lpps) in let lppcs = filter (fun (m,(n1,n2)) -> n1 <= n2) @@ -1859,5 +1579,327 @@ let sumofsquares_general_symmetry tool pol = let sos = poly_cmul rat (end_itlist poly_add sqs) in if is_undefined(poly_sub sos pol) then rat,lins else raise Sanity;; -let (sumofsquares: poly -> Num.num * (( Num.num * poly) list)) = -sumofsquares_general_symmetry csdp;; +let sumofsquares = sumofsquares_general_symmetry csdp;; + +(* ------------------------------------------------------------------------- *) +(* Pure HOL SOS conversion. *) +(* ------------------------------------------------------------------------- *) +(* +let SOS_CONV = + let mk_square = + let pow_tm = `(pow)` and two_tm = `2` in + fun tm -> mk_comb(mk_comb(pow_tm,tm),two_tm) + and mk_prod = mk_binop `( * )` + and mk_sum = mk_binop `(+)` in + fun tm -> + let k,sos = sumofsquares(poly_of_term tm) in + let mk_sqtm(c,p) = + mk_prod (term_of_rat(k */ c)) (mk_square(term_of_poly p)) in + let tm' = end_itlist mk_sum (map mk_sqtm sos) in + let th = REAL_POLY_CONV tm and th' = REAL_POLY_CONV tm' in + TRANS th (SYM th');; +*) +(* ------------------------------------------------------------------------- *) +(* Attempt to prove &0 <= x by direct SOS decomposition. *) +(* ------------------------------------------------------------------------- *) +(* +let PURE_SOS_TAC = + let tac = + MATCH_ACCEPT_TAC(REWRITE_RULE[GSYM REAL_POW_2] REAL_LE_SQUARE) ORELSE + MATCH_ACCEPT_TAC REAL_LE_SQUARE ORELSE + (MATCH_MP_TAC REAL_LE_ADD THEN CONJ_TAC) ORELSE + (MATCH_MP_TAC REAL_LE_MUL THEN CONJ_TAC) ORELSE + CONV_TAC(RAND_CONV REAL_RAT_REDUCE_CONV THENC REAL_RAT_LE_CONV) in + REPEAT GEN_TAC THEN REWRITE_TAC[real_ge] THEN + GEN_REWRITE_TAC I [GSYM REAL_SUB_LE] THEN + CONV_TAC(RAND_CONV SOS_CONV) THEN + REPEAT tac THEN NO_TAC;; + +let PURE_SOS tm = prove(tm,PURE_SOS_TAC);; +*) +(* ------------------------------------------------------------------------- *) +(* Examples. *) +(* ------------------------------------------------------------------------- *) + +(***** + +time REAL_SOS + `a1 >= &0 /\ a2 >= &0 /\ + (a1 * a1 + a2 * a2 = b1 * b1 + b2 * b2 + &2) /\ + (a1 * b1 + a2 * b2 = &0) + ==> a1 * a2 - b1 * b2 >= &0`;; + +time REAL_SOS `&3 * x + &7 * a < &4 /\ &3 < &2 * x ==> a < &0`;; + +time REAL_SOS + `b pow 2 < &4 * a * c ==> ~(a * x pow 2 + b * x + c = &0)`;; + +time REAL_SOS + `(a * x pow 2 + b * x + c = &0) ==> b pow 2 >= &4 * a * c`;; + +time REAL_SOS + `&0 <= x /\ x <= &1 /\ &0 <= y /\ y <= &1 + ==> x pow 2 + y pow 2 < &1 \/ + (x - &1) pow 2 + y pow 2 < &1 \/ + x pow 2 + (y - &1) pow 2 < &1 \/ + (x - &1) pow 2 + (y - &1) pow 2 < &1`;; + +time REAL_SOS + `&0 <= b /\ &0 <= c /\ &0 <= x /\ &0 <= y /\ + (x pow 2 = c) /\ (y pow 2 = a pow 2 * c + b) + ==> a * c <= y * x`;; + +time REAL_SOS + `&0 <= x /\ &0 <= y /\ &0 <= z /\ x + y + z <= &3 + ==> x * y + x * z + y * z >= &3 * x * y * z`;; + +time REAL_SOS + `(x pow 2 + y pow 2 + z pow 2 = &1) ==> (x + y + z) pow 2 <= &3`;; + +time REAL_SOS + `(w pow 2 + x pow 2 + y pow 2 + z pow 2 = &1) + ==> (w + x + y + z) pow 2 <= &4`;; + +time REAL_SOS + `x >= &1 /\ y >= &1 ==> x * y >= x + y - &1`;; + +time REAL_SOS + `x > &1 /\ y > &1 ==> x * y > x + y - &1`;; + +time REAL_SOS + `abs(x) <= &1 + ==> abs(&64 * x pow 7 - &112 * x pow 5 + &56 * x pow 3 - &7 * x) <= &1`;; + +time REAL_SOS + `abs(x - z) <= e /\ abs(y - z) <= e /\ &0 <= u /\ &0 <= v /\ (u + v = &1) + ==> abs((u * x + v * y) - z) <= e`;; + +(* ------------------------------------------------------------------------- *) +(* One component of denominator in dodecahedral example. *) +(* ------------------------------------------------------------------------- *) + +time REAL_SOS + `&2 <= x /\ x <= &125841 / &50000 /\ + &2 <= y /\ y <= &125841 / &50000 /\ + &2 <= z /\ z <= &125841 / &50000 + ==> &2 * (x * z + x * y + y * z) - (x * x + y * y + z * z) >= &0`;; + +(* ------------------------------------------------------------------------- *) +(* Over a larger but simpler interval. *) +(* ------------------------------------------------------------------------- *) + +time REAL_SOS + `&2 <= x /\ x <= &4 /\ &2 <= y /\ y <= &4 /\ &2 <= z /\ z <= &4 + ==> &0 <= &2 * (x * z + x * y + y * z) - (x * x + y * y + z * z)`;; + +(* ------------------------------------------------------------------------- *) +(* We can do 12. I think 12 is a sharp bound; see PP's certificate. *) +(* ------------------------------------------------------------------------- *) + +time REAL_SOS + `&2 <= x /\ x <= &4 /\ &2 <= y /\ y <= &4 /\ &2 <= z /\ z <= &4 + ==> &12 <= &2 * (x * z + x * y + y * z) - (x * x + y * y + z * z)`;; + +(* ------------------------------------------------------------------------- *) +(* Gloptipoly example. *) +(* ------------------------------------------------------------------------- *) + +(*** This works but normalization takes minutes + +time REAL_SOS + `(x - y - &2 * x pow 4 = &0) /\ &0 <= x /\ x <= &2 /\ &0 <= y /\ y <= &3 + ==> y pow 2 - &7 * y - &12 * x + &17 >= &0`;; + + ***) + +(* ------------------------------------------------------------------------- *) +(* Inequality from sci.math (see "Leon-Sotelo, por favor"). *) +(* ------------------------------------------------------------------------- *) + +time REAL_SOS + `&0 <= x /\ &0 <= y /\ (x * y = &1) + ==> x + y <= x pow 2 + y pow 2`;; + +time REAL_SOS + `&0 <= x /\ &0 <= y /\ (x * y = &1) + ==> x * y * (x + y) <= x pow 2 + y pow 2`;; + +time REAL_SOS + `&0 <= x /\ &0 <= y ==> x * y * (x + y) pow 2 <= (x pow 2 + y pow 2) pow 2`;; + +(* ------------------------------------------------------------------------- *) +(* Some examples over integers and natural numbers. *) +(* ------------------------------------------------------------------------- *) + +time SOS_RULE `!m n. 2 * m + n = (n + m) + m`;; +time SOS_RULE `!n. ~(n = 0) ==> (0 MOD n = 0)`;; +time SOS_RULE `!m n. m < n ==> (m DIV n = 0)`;; +time SOS_RULE `!n:num. n <= n * n`;; +time SOS_RULE `!m n. n * (m DIV n) <= m`;; +time SOS_RULE `!n. ~(n = 0) ==> (0 DIV n = 0)`;; +time SOS_RULE `!m n p. ~(p = 0) /\ m <= n ==> m DIV p <= n DIV p`;; +time SOS_RULE `!a b n. ~(a = 0) ==> (n <= b DIV a <=> a * n <= b)`;; + +(* ------------------------------------------------------------------------- *) +(* This is particularly gratifying --- cf hideous manual proof in arith.ml *) +(* ------------------------------------------------------------------------- *) + +(*** This doesn't now seem to work as well as it did; what changed? + +time SOS_RULE + `!a b c d. ~(b = 0) /\ b * c < (a + 1) * d ==> c DIV d <= a DIV b`;; + + ***) + +(* ------------------------------------------------------------------------- *) +(* Key lemma for injectivity of Cantor-type pairing functions. *) +(* ------------------------------------------------------------------------- *) + +time SOS_RULE + `!x1 y1 x2 y2. ((x1 + y1) EXP 2 + x1 + 1 = (x2 + y2) EXP 2 + x2 + 1) + ==> (x1 + y1 = x2 + y2)`;; + +time SOS_RULE + `!x1 y1 x2 y2. ((x1 + y1) EXP 2 + x1 + 1 = (x2 + y2) EXP 2 + x2 + 1) /\ + (x1 + y1 = x2 + y2) + ==> (x1 = x2) /\ (y1 = y2)`;; + +time SOS_RULE + `!x1 y1 x2 y2. + (((x1 + y1) EXP 2 + 3 * x1 + y1) DIV 2 = + ((x2 + y2) EXP 2 + 3 * x2 + y2) DIV 2) + ==> (x1 + y1 = x2 + y2)`;; + +time SOS_RULE + `!x1 y1 x2 y2. + (((x1 + y1) EXP 2 + 3 * x1 + y1) DIV 2 = + ((x2 + y2) EXP 2 + 3 * x2 + y2) DIV 2) /\ + (x1 + y1 = x2 + y2) + ==> (x1 = x2) /\ (y1 = y2)`;; + +(* ------------------------------------------------------------------------- *) +(* Reciprocal multiplication (actually just ARITH_RULE does these). *) +(* ------------------------------------------------------------------------- *) + +time SOS_RULE `x <= 127 ==> ((86 * x) DIV 256 = x DIV 3)`;; + +time SOS_RULE `x < 2 EXP 16 ==> ((104858 * x) DIV (2 EXP 20) = x DIV 10)`;; + +(* ------------------------------------------------------------------------- *) +(* This is more impressive since it's really nonlinear. See REMAINDER_DECODE *) +(* ------------------------------------------------------------------------- *) + +time SOS_RULE `0 < m /\ m < n ==> ((m * ((n * x) DIV m + 1)) DIV n = x)`;; + +(* ------------------------------------------------------------------------- *) +(* Some conversion examples. *) +(* ------------------------------------------------------------------------- *) + +time SOS_CONV + `&2 * x pow 4 + &2 * x pow 3 * y - x pow 2 * y pow 2 + &5 * y pow 4`;; + +time SOS_CONV + `x pow 4 - (&2 * y * z + &1) * x pow 2 + + (y pow 2 * z pow 2 + &2 * y * z + &2)`;; + +time SOS_CONV `&4 * x pow 4 + + &4 * x pow 3 * y - &7 * x pow 2 * y pow 2 - &2 * x * y pow 3 + + &10 * y pow 4`;; + +time SOS_CONV `&4 * x pow 4 * y pow 6 + x pow 2 - x * y pow 2 + y pow 2`;; + +time SOS_CONV + `&4096 * (x pow 4 + x pow 2 + z pow 6 - &3 * x pow 2 * z pow 2) + &729`;; + +time SOS_CONV + `&120 * x pow 2 - &63 * x pow 4 + &10 * x pow 6 + + &30 * x * y - &120 * y pow 2 + &120 * y pow 4 + &31`;; + +time SOS_CONV + `&9 * x pow 2 * y pow 4 + &9 * x pow 2 * z pow 4 + &36 * x pow 2 * y pow 3 + + &36 * x pow 2 * y pow 2 - &48 * x * y * z pow 2 + &4 * y pow 4 + + &4 * z pow 4 - &16 * y pow 3 + &16 * y pow 2`;; + +time SOS_CONV + `(x pow 2 + y pow 2 + z pow 2) * + (x pow 4 * y pow 2 + x pow 2 * y pow 4 + + z pow 6 - &3 * x pow 2 * y pow 2 * z pow 2)`;; + +time SOS_CONV + `x pow 4 + y pow 4 + z pow 4 - &4 * x * y * z + x + y + z + &3`;; + +(*** I think this will work, but normalization is slow + +time SOS_CONV + `&100 * (x pow 4 + y pow 4 + z pow 4 - &4 * x * y * z + x + y + z) + &212`;; + + ***) + +time SOS_CONV + `&100 * ((&2 * x - &2) pow 2 + (x pow 3 - &8 * x - &2) pow 2) - &588`;; + +time SOS_CONV + `x pow 2 * (&120 - &63 * x pow 2 + &10 * x pow 4) + &30 * x * y + + &30 * y pow 2 * (&4 * y pow 2 - &4) + &31`;; + +(* ------------------------------------------------------------------------- *) +(* Example of basic rule. *) +(* ------------------------------------------------------------------------- *) + +time PURE_SOS + `!x. x pow 4 + y pow 4 + z pow 4 - &4 * x * y * z + x + y + z + &3 + >= &1 / &7`;; + +time PURE_SOS + `&0 <= &98 * x pow 12 + + -- &980 * x pow 10 + + &3038 * x pow 8 + + -- &2968 * x pow 6 + + &1022 * x pow 4 + + -- &84 * x pow 2 + + &2`;; + +time PURE_SOS + `!x. &0 <= &2 * x pow 14 + + -- &84 * x pow 12 + + &1022 * x pow 10 + + -- &2968 * x pow 8 + + &3038 * x pow 6 + + -- &980 * x pow 4 + + &98 * x pow 2`;; + +(* ------------------------------------------------------------------------- *) +(* From Zeng et al, JSC vol 37 (2004), p83-99. *) +(* All of them work nicely with pure SOS_CONV, except (maybe) the one noted. *) +(* ------------------------------------------------------------------------- *) + +PURE_SOS + `x pow 6 + y pow 6 + z pow 6 - &3 * x pow 2 * y pow 2 * z pow 2 >= &0`;; + +PURE_SOS `x pow 4 + y pow 4 + z pow 4 + &1 - &4*x*y*z >= &0`;; + +PURE_SOS `x pow 4 + &2*x pow 2*z + x pow 2 - &2*x*y*z + &2*y pow 2*z pow 2 + +&2*y*z pow 2 + &2*z pow 2 - &2*x + &2* y*z + &1 >= &0`;; + +(**** This is harder. Interestingly, this fails the pure SOS test, it seems. + Yet only on rounding(!?) Poor Newton polytope optimization or something? + But REAL_SOS does finally converge on the second run at level 12! + +REAL_SOS +`x pow 4*y pow 4 - &2*x pow 5*y pow 3*z pow 2 + x pow 6*y pow 2*z pow 4 + &2*x +pow 2*y pow 3*z - &4* x pow 3*y pow 2*z pow 3 + &2*x pow 4*y*z pow 5 + z pow +2*y pow 2 - &2*z pow 4*y*x + z pow 6*x pow 2 >= &0`;; + + ****) + +PURE_SOS +`x pow 4 + &4*x pow 2*y pow 2 + &2*x*y*z pow 2 + &2*x*y*w pow 2 + y pow 4 + z +pow 4 + w pow 4 + &2*z pow 2*w pow 2 + &2*x pow 2*w + &2*y pow 2*w + &2*x*y + +&3*w pow 2 + &2*z pow 2 + &1 >= &0`;; + +PURE_SOS +`w pow 6 + &2*z pow 2*w pow 3 + x pow 4 + y pow 4 + z pow 4 + &2*x pow 2*w + +&2*x pow 2*z + &3*x pow 2 + w pow 2 + &2*z*w + z pow 2 + &2*z + &2*w + &1 >= +&0`;; + +*****) diff --git a/plugins/micromega/sos.mli b/plugins/micromega/sos.mli index 31c9518c8..42e22ffec 100644 --- a/plugins/micromega/sos.mli +++ b/plugins/micromega/sos.mli @@ -6,33 +6,7 @@ (* * GNU Lesser General Public License Version 2.1 *) (************************************************************************) - -type vname = string;; - -type term = -| Zero -| Const of Num.num -| Var of vname -| Inv of term -| Opp of term -| Add of (term * term) -| Sub of (term * term) -| Mul of (term * term) -| Div of (term * term) -| Pow of (term * int) - -type positivstellensatz = - Axiom_eq of int - | Axiom_le of int - | Axiom_lt of int - | Rational_eq of Num.num - | Rational_le of Num.num - | Rational_lt of Num.num - | Square of term - | Monoid of int list - | Eqmul of term * positivstellensatz - | Sum of positivstellensatz * positivstellensatz - | Product of positivstellensatz * positivstellensatz +open Sos_types type poly @@ -55,10 +29,6 @@ val term_of_sos : positivstellensatz * (Num.num * poly) list -> val string_of_poly : poly -> string -exception TooDeep - -val deepen_until : int -> (int -> 'a) -> int -> 'a - val real_positivnullstellensatz_general : bool -> int -> poly list -> (poly * positivstellensatz) list -> poly -> poly list * (positivstellensatz * (Num.num * poly) list) list diff --git a/plugins/micromega/sos_lib.ml b/plugins/micromega/sos_lib.ml new file mode 100644 index 000000000..a9228365e --- /dev/null +++ b/plugins/micromega/sos_lib.ml @@ -0,0 +1,621 @@ +(* ========================================================================= *) +(* - This code originates from John Harrison's HOL LIGHT 2.30 *) +(* (see file LICENSE.sos for license, copyright and disclaimer) *) +(* This code is the HOL LIGHT library code used by sos.ml *) +(* - Laurent Théry (thery@sophia.inria.fr) has isolated the HOL *) +(* independent bits *) +(* - Frédéric Besson (fbesson@irisa.fr) is using it to feed micromega *) +(* ========================================================================= *) +open Sos_types +open Num +open List + +let debugging = ref false;; + +(* ------------------------------------------------------------------------- *) +(* Comparisons that are reflexive on NaN and also short-circuiting. *) +(* ------------------------------------------------------------------------- *) + +let (=?) = fun x y -> Pervasives.compare x y = 0;; +let ( Pervasives.compare x y < 0;; +let (<=?) = fun x y -> Pervasives.compare x y <= 0;; +let (>?) = fun x y -> Pervasives.compare x y > 0;; +let (>=?) = fun x y -> Pervasives.compare x y >= 0;; + +(* ------------------------------------------------------------------------- *) +(* Combinators. *) +(* ------------------------------------------------------------------------- *) + +let (o) = fun f g x -> f(g x);; + +(* ------------------------------------------------------------------------- *) +(* Some useful functions on "num" type. *) +(* ------------------------------------------------------------------------- *) + + +let num_0 = Int 0 +and num_1 = Int 1 +and num_2 = Int 2 +and num_10 = Int 10;; + +let pow2 n = power_num num_2 (Int n);; +let pow10 n = power_num num_10 (Int n);; + +let numdom r = + let r' = Ratio.normalize_ratio (ratio_of_num r) in + num_of_big_int(Ratio.numerator_ratio r'), + num_of_big_int(Ratio.denominator_ratio r');; + +let numerator = (o) fst numdom +and denominator = (o) snd numdom;; + +let gcd_num n1 n2 = + num_of_big_int(Big_int.gcd_big_int (big_int_of_num n1) (big_int_of_num n2));; + +let lcm_num x y = + if x =/ num_0 & y =/ num_0 then num_0 + else abs_num((x */ y) // gcd_num x y);; + + +(* ------------------------------------------------------------------------- *) +(* List basics. *) +(* ------------------------------------------------------------------------- *) + +let rec el n l = + if n = 0 then hd l else el (n - 1) (tl l);; + + +(* ------------------------------------------------------------------------- *) +(* Various versions of list iteration. *) +(* ------------------------------------------------------------------------- *) + +let rec itlist f l b = + match l with + [] -> b + | (h::t) -> f h (itlist f t b);; + +let rec end_itlist f l = + match l with + [] -> failwith "end_itlist" + | [x] -> x + | (h::t) -> f h (end_itlist f t);; + +let rec itlist2 f l1 l2 b = + match (l1,l2) with + ([],[]) -> b + | (h1::t1,h2::t2) -> f h1 h2 (itlist2 f t1 t2 b) + | _ -> failwith "itlist2";; + +(* ------------------------------------------------------------------------- *) +(* All pairs arising from applying a function over two lists. *) +(* ------------------------------------------------------------------------- *) + +let rec allpairs f l1 l2 = + match l1 with + h1::t1 -> itlist (fun x a -> f h1 x :: a) l2 (allpairs f t1 l2) + | [] -> [];; + +(* ------------------------------------------------------------------------- *) +(* String operations (surely there is a better way...) *) +(* ------------------------------------------------------------------------- *) + +let implode l = itlist (^) l "";; + +let explode s = + let rec exap n l = + if n < 0 then l else + exap (n - 1) ((String.sub s n 1)::l) in + exap (String.length s - 1) [];; + + +(* ------------------------------------------------------------------------- *) +(* Attempting function or predicate applications. *) +(* ------------------------------------------------------------------------- *) + +let can f x = try (f x; true) with Failure _ -> false;; + + +(* ------------------------------------------------------------------------- *) +(* Repetition of a function. *) +(* ------------------------------------------------------------------------- *) + +let rec funpow n f x = + if n < 1 then x else funpow (n-1) f (f x);; + + + +(* ------------------------------------------------------------------------- *) +(* Replication and sequences. *) +(* ------------------------------------------------------------------------- *) + +let rec replicate x n = + if n < 1 then [] + else x::(replicate x (n - 1));; + +let rec (--) = fun m n -> if m > n then [] else m::((m + 1) -- n);; + +(* ------------------------------------------------------------------------- *) +(* Various useful list operations. *) +(* ------------------------------------------------------------------------- *) + +let rec forall p l = + match l with + [] -> true + | h::t -> p(h) & forall p t;; + +let rec tryfind f l = + match l with + [] -> failwith "tryfind" + | (h::t) -> try f h with Failure _ -> tryfind f t;; + +let index x = + let rec ind n l = + match l with + [] -> failwith "index" + | (h::t) -> if x =? h then n else ind (n + 1) t in + ind 0;; + +(* ------------------------------------------------------------------------- *) +(* "Set" operations on lists. *) +(* ------------------------------------------------------------------------- *) + +let rec mem x lis = + match lis with + [] -> false + | (h::t) -> x =? h or mem x t;; + +let insert x l = + if mem x l then l else x::l;; + +let union l1 l2 = itlist insert l1 l2;; + +let subtract l1 l2 = filter (fun x -> not (mem x l2)) l1;; + +(* ------------------------------------------------------------------------- *) +(* Merging and bottom-up mergesort. *) +(* ------------------------------------------------------------------------- *) + +let rec merge ord l1 l2 = + match l1 with + [] -> l2 + | h1::t1 -> match l2 with + [] -> l1 + | h2::t2 -> if ord h1 h2 then h1::(merge ord t1 l2) + else h2::(merge ord l1 t2);; + + +(* ------------------------------------------------------------------------- *) +(* Common measure predicates to use with "sort". *) +(* ------------------------------------------------------------------------- *) + +let increasing f x y = f x ? f y;; + +(* ------------------------------------------------------------------------- *) +(* Zipping, unzipping etc. *) +(* ------------------------------------------------------------------------- *) + +let rec zip l1 l2 = + match (l1,l2) with + ([],[]) -> [] + | (h1::t1,h2::t2) -> (h1,h2)::(zip t1 t2) + | _ -> failwith "zip";; + +let rec unzip = + function [] -> [],[] + | ((a,b)::rest) -> let alist,blist = unzip rest in + (a::alist,b::blist);; + +(* ------------------------------------------------------------------------- *) +(* Iterating functions over lists. *) +(* ------------------------------------------------------------------------- *) + +let rec do_list f l = + match l with + [] -> () + | (h::t) -> (f h; do_list f t);; + +(* ------------------------------------------------------------------------- *) +(* Sorting. *) +(* ------------------------------------------------------------------------- *) + +let rec sort cmp lis = + match lis with + [] -> [] + | piv::rest -> + let r,l = partition (cmp piv) rest in + (sort cmp l) @ (piv::(sort cmp r));; + +(* ------------------------------------------------------------------------- *) +(* Removing adjacent (NB!) equal elements from list. *) +(* ------------------------------------------------------------------------- *) + +let rec uniq l = + match l with + x::(y::_ as t) -> let t' = uniq t in + if x =? y then t' else + if t'==t then l else x::t' + | _ -> l;; + +(* ------------------------------------------------------------------------- *) +(* Convert list into set by eliminating duplicates. *) +(* ------------------------------------------------------------------------- *) + +let setify s = uniq (sort (<=?) s);; + +(* ------------------------------------------------------------------------- *) +(* Polymorphic finite partial functions via Patricia trees. *) +(* *) +(* The point of this strange representation is that it is canonical (equal *) +(* functions have the same encoding) yet reasonably efficient on average. *) +(* *) +(* Idea due to Diego Olivier Fernandez Pons (OCaml list, 2003/11/10). *) +(* ------------------------------------------------------------------------- *) + +type ('a,'b)func = + Empty + | Leaf of int * ('a*'b)list + | Branch of int * int * ('a,'b)func * ('a,'b)func;; + +(* ------------------------------------------------------------------------- *) +(* Undefined function. *) +(* ------------------------------------------------------------------------- *) + +let undefined = Empty;; + +(* ------------------------------------------------------------------------- *) +(* In case of equality comparison worries, better use this. *) +(* ------------------------------------------------------------------------- *) + +let is_undefined f = + match f with + Empty -> true + | _ -> false;; + +(* ------------------------------------------------------------------------- *) +(* Operation analagous to "map" for lists. *) +(* ------------------------------------------------------------------------- *) + +let mapf = + let rec map_list f l = + match l with + [] -> [] + | (x,y)::t -> (x,f(y))::(map_list f t) in + let rec mapf f t = + match t with + Empty -> Empty + | Leaf(h,l) -> Leaf(h,map_list f l) + | Branch(p,b,l,r) -> Branch(p,b,mapf f l,mapf f r) in + mapf;; + +(* ------------------------------------------------------------------------- *) +(* Operations analogous to "fold" for lists. *) +(* ------------------------------------------------------------------------- *) + +let foldl = + let rec foldl_list f a l = + match l with + [] -> a + | (x,y)::t -> foldl_list f (f a x y) t in + let rec foldl f a t = + match t with + Empty -> a + | Leaf(h,l) -> foldl_list f a l + | Branch(p,b,l,r) -> foldl f (foldl f a l) r in + foldl;; + +let foldr = + let rec foldr_list f l a = + match l with + [] -> a + | (x,y)::t -> f x y (foldr_list f t a) in + let rec foldr f t a = + match t with + Empty -> a + | Leaf(h,l) -> foldr_list f l a + | Branch(p,b,l,r) -> foldr f l (foldr f r a) in + foldr;; + +(* ------------------------------------------------------------------------- *) +(* Redefinition and combination. *) +(* ------------------------------------------------------------------------- *) + +let (|->),combine = + let ldb x y = let z = x lxor y in z land (-z) in + let newbranch p1 t1 p2 t2 = + let b = ldb p1 p2 in + let p = p1 land (b - 1) in + if p1 land b = 0 then Branch(p,b,t1,t2) + else Branch(p,b,t2,t1) in + let rec define_list (x,y as xy) l = + match l with + (a,b as ab)::t -> + if x =? a then xy::t + else if x [xy] + and combine_list op z l1 l2 = + match (l1,l2) with + [],_ -> l2 + | _,[] -> l1 + | ((x1,y1 as xy1)::t1,(x2,y2 as xy2)::t2) -> + if x1 ) x y = + let k = Hashtbl.hash x in + let rec upd t = + match t with + Empty -> Leaf (k,[x,y]) + | Leaf(h,l) -> + if h = k then Leaf(h,define_list (x,y) l) + else newbranch h t k (Leaf(k,[x,y])) + | Branch(p,b,l,r) -> + if k land (b - 1) <> p then newbranch p t k (Leaf(k,[x,y])) + else if k land b = 0 then Branch(p,b,upd l,r) + else Branch(p,b,l,upd r) in + upd in + let rec combine op z t1 t2 = + match (t1,t2) with + Empty,_ -> t2 + | _,Empty -> t1 + | Leaf(h1,l1),Leaf(h2,l2) -> + if h1 = h2 then + let l = combine_list op z l1 l2 in + if l = [] then Empty else Leaf(h1,l) + else newbranch h1 t1 h2 t2 + | (Leaf(k,lis) as lf),(Branch(p,b,l,r) as br) | + (Branch(p,b,l,r) as br),(Leaf(k,lis) as lf) -> + if k land (b - 1) = p then + if k land b = 0 then + let l' = combine op z lf l in + if is_undefined l' then r else Branch(p,b,l',r) + else + let r' = combine op z lf r in + if is_undefined r' then l else Branch(p,b,l,r') + else + newbranch k lf p br + | Branch(p1,b1,l1,r1),Branch(p2,b2,l2,r2) -> + if b1 < b2 then + if p2 land (b1 - 1) <> p1 then newbranch p1 t1 p2 t2 + else if p2 land b1 = 0 then + let l = combine op z l1 t2 in + if is_undefined l then r1 else Branch(p1,b1,l,r1) + else + let r = combine op z r1 t2 in + if is_undefined r then l1 else Branch(p1,b1,l1,r) + else if b2 < b1 then + if p1 land (b2 - 1) <> p2 then newbranch p1 t1 p2 t2 + else if p1 land b2 = 0 then + let l = combine op z t1 l2 in + if is_undefined l then r2 else Branch(p2,b2,l,r2) + else + let r = combine op z t1 r2 in + if is_undefined r then l2 else Branch(p2,b2,l2,r) + else if p1 = p2 then + let l = combine op z l1 l2 and r = combine op z r1 r2 in + if is_undefined l then r + else if is_undefined r then l else Branch(p1,b1,l,r) + else + newbranch p1 t1 p2 t2 in + (|->),combine;; + +(* ------------------------------------------------------------------------- *) +(* Special case of point function. *) +(* ------------------------------------------------------------------------- *) + +let (|=>) = fun x y -> (x |-> y) undefined;; + + +(* ------------------------------------------------------------------------- *) +(* Grab an arbitrary element. *) +(* ------------------------------------------------------------------------- *) + +let rec choose t = + match t with + Empty -> failwith "choose: completely undefined function" + | Leaf(h,l) -> hd l + | Branch(b,p,t1,t2) -> choose t1;; + +(* ------------------------------------------------------------------------- *) +(* Application. *) +(* ------------------------------------------------------------------------- *) + +let applyd = + let rec apply_listd l d x = + match l with + (a,b)::t -> if x =? a then b + else if x >? a then apply_listd t d x else d x + | [] -> d x in + fun f d x -> + let k = Hashtbl.hash x in + let rec look t = + match t with + Leaf(h,l) when h = k -> apply_listd l d x + | Branch(p,b,l,r) -> look (if k land b = 0 then l else r) + | _ -> d x in + look f;; + +let apply f = applyd f (fun x -> failwith "apply");; + +let tryapplyd f a d = applyd f (fun x -> d) a;; + +let defined f x = try apply f x; true with Failure _ -> false;; + +(* ------------------------------------------------------------------------- *) +(* Undefinition. *) +(* ------------------------------------------------------------------------- *) + +let undefine = + let rec undefine_list x l = + match l with + (a,b as ab)::t -> + if x =? a then t + else if x [] in + fun x -> + let k = Hashtbl.hash x in + let rec und t = + match t with + Leaf(h,l) when h = k -> + let l' = undefine_list x l in + if l' == l then t + else if l' = [] then Empty + else Leaf(h,l') + | Branch(p,b,l,r) when k land (b - 1) = p -> + if k land b = 0 then + let l' = und l in + if l' == l then t + else if is_undefined l' then r + else Branch(p,b,l',r) + else + let r' = und r in + if r' == r then t + else if is_undefined r' then l + else Branch(p,b,l,r') + | _ -> t in + und;; + + +(* ------------------------------------------------------------------------- *) +(* Mapping to sorted-list representation of the graph, domain and range. *) +(* ------------------------------------------------------------------------- *) + +let graph f = setify (foldl (fun a x y -> (x,y)::a) [] f);; + +let dom f = setify(foldl (fun a x y -> x::a) [] f);; + +let ran f = setify(foldl (fun a x y -> y::a) [] f);; + +(* ------------------------------------------------------------------------- *) +(* More parser basics. *) +(* ------------------------------------------------------------------------- *) + +exception Noparse;; + + +let isspace,issep,isbra,issymb,isalpha,isnum,isalnum = + let charcode s = Char.code(String.get s 0) in + let spaces = " \t\n\r" + and separators = ",;" + and brackets = "()[]{}" + and symbs = "\\!@#$%^&*-+|\\<=>/?~.:" + and alphas = "'abcdefghijklmnopqrstuvwxyz_ABCDEFGHIJKLMNOPQRSTUVWXYZ" + and nums = "0123456789" in + let allchars = spaces^separators^brackets^symbs^alphas^nums in + let csetsize = itlist ((o) max charcode) (explode allchars) 256 in + let ctable = Array.make csetsize 0 in + do_list (fun c -> Array.set ctable (charcode c) 1) (explode spaces); + do_list (fun c -> Array.set ctable (charcode c) 2) (explode separators); + do_list (fun c -> Array.set ctable (charcode c) 4) (explode brackets); + do_list (fun c -> Array.set ctable (charcode c) 8) (explode symbs); + do_list (fun c -> Array.set ctable (charcode c) 16) (explode alphas); + do_list (fun c -> Array.set ctable (charcode c) 32) (explode nums); + let isspace c = Array.get ctable (charcode c) = 1 + and issep c = Array.get ctable (charcode c) = 2 + and isbra c = Array.get ctable (charcode c) = 4 + and issymb c = Array.get ctable (charcode c) = 8 + and isalpha c = Array.get ctable (charcode c) = 16 + and isnum c = Array.get ctable (charcode c) = 32 + and isalnum c = Array.get ctable (charcode c) >= 16 in + isspace,issep,isbra,issymb,isalpha,isnum,isalnum;; + +let (||) parser1 parser2 input = + try parser1 input + with Noparse -> parser2 input;; + +let (++) parser1 parser2 input = + let result1,rest1 = parser1 input in + let result2,rest2 = parser2 rest1 in + (result1,result2),rest2;; + +let rec many prs input = + try let result,next = prs input in + let results,rest = many prs next in + (result::results),rest + with Noparse -> [],input;; + +let (>>) prs treatment input = + let result,rest = prs input in + treatment(result),rest;; + +let fix err prs input = + try prs input + with Noparse -> failwith (err ^ " expected");; + +let rec listof prs sep err = + prs ++ many (sep ++ fix err prs >> snd) >> (fun (h,t) -> h::t);; + +let possibly prs input = + try let x,rest = prs input in [x],rest + with Noparse -> [],input;; + +let some p = + function + [] -> raise Noparse + | (h::t) -> if p h then (h,t) else raise Noparse;; + +let a tok = some (fun item -> item = tok);; + +let rec atleast n prs i = + (if n <= 0 then many prs + else prs ++ atleast (n - 1) prs >> (fun (h,t) -> h::t)) i;; + +let finished input = + if input = [] then 0,input else failwith "Unparsed input";; + +(* ------------------------------------------------------------------------- *) + +let temp_path = ref Filename.temp_dir_name;; + +(* ------------------------------------------------------------------------- *) +(* Convenient conversion between files and (lists of) strings. *) +(* ------------------------------------------------------------------------- *) + +let strings_of_file filename = + let fd = try Pervasives.open_in filename + with Sys_error _ -> + failwith("strings_of_file: can't open "^filename) in + let rec suck_lines acc = + try let l = Pervasives.input_line fd in + suck_lines (l::acc) + with End_of_file -> rev acc in + let data = suck_lines [] in + (Pervasives.close_in fd; data);; + +let string_of_file filename = + end_itlist (fun s t -> s^"\n"^t) (strings_of_file filename);; + +let file_of_string filename s = + let fd = Pervasives.open_out filename in + output_string fd s; close_out fd;; + + +(* ------------------------------------------------------------------------- *) +(* Iterative deepening. *) +(* ------------------------------------------------------------------------- *) + +let rec deepen f n = + try (*print_string "Searching with depth limit "; + print_int n; print_newline();*) f n + with Failure _ -> deepen f (n + 1);; + +exception TooDeep + +let deepen_until limit f n = + match compare limit 0 with + | 0 -> raise TooDeep + | -1 -> deepen f n + | _ -> + let rec d_until f n = + try(* if !debugging + then (print_string "Searching with depth limit "; + print_int n; print_newline()) ;*) f n + with Failure x -> + (*if !debugging then (Printf.printf "solver error : %s\n" x) ; *) + if n = limit then raise TooDeep else d_until f (n + 1) in + d_until f n diff --git a/plugins/micromega/sos_types.ml b/plugins/micromega/sos_types.ml new file mode 100644 index 000000000..fe481ecc3 --- /dev/null +++ b/plugins/micromega/sos_types.ml @@ -0,0 +1,68 @@ +(************************************************************************) +(* v * The Coq Proof Assistant / The Coq Development Team *) +(* output_string o "0" + | Const n -> output_string o (string_of_num n) + | Var n -> Printf.fprintf o "v%s" n + | Inv t -> Printf.fprintf o "1/(%a)" output_term t + | Opp t -> Printf.fprintf o "- (%a)" output_term t + | Add(t1,t2) -> Printf.fprintf o "(%a)+(%a)" output_term t1 output_term t2 + | Sub(t1,t2) -> Printf.fprintf o "(%a)-(%a)" output_term t1 output_term t2 + | Mul(t1,t2) -> Printf.fprintf o "(%a)*(%a)" output_term t1 output_term t2 + | Div(t1,t2) -> Printf.fprintf o "(%a)/(%a)" output_term t1 output_term t2 + | Pow(t1,i) -> Printf.fprintf o "(%a)^(%i)" output_term t1 i +(* ------------------------------------------------------------------------- *) +(* Data structure for Positivstellensatz refutations. *) +(* ------------------------------------------------------------------------- *) + +type positivstellensatz = + Axiom_eq of int + | Axiom_le of int + | Axiom_lt of int + | Rational_eq of num + | Rational_le of num + | Rational_lt of num + | Square of term + | Monoid of int list + | Eqmul of term * positivstellensatz + | Sum of positivstellensatz * positivstellensatz + | Product of positivstellensatz * positivstellensatz;; + + +let rec output_psatz o = function + | Axiom_eq i -> Printf.fprintf o "Aeq(%i)" i + | Axiom_le i -> Printf.fprintf o "Ale(%i)" i + | Axiom_lt i -> Printf.fprintf o "Alt(%i)" i + | Rational_eq n -> Printf.fprintf o "eq(%s)" (string_of_num n) + | Rational_le n -> Printf.fprintf o "le(%s)" (string_of_num n) + | Rational_lt n -> Printf.fprintf o "lt(%s)" (string_of_num n) + | Square t -> Printf.fprintf o "(%a)^2" output_term t + | Monoid l -> Printf.fprintf o "monoid" + | Eqmul (t,ps) -> Printf.fprintf o "%a * %a" output_term t output_psatz ps + | Sum (t1,t2) -> Printf.fprintf o "%a + %a" output_psatz t1 output_psatz t2 + | Product (t1,t2) -> Printf.fprintf o "%a * %a" output_psatz t1 output_psatz t2 diff --git a/test-suite/csdp.cache b/test-suite/csdp.cache index f6f068c04..d3a700c04 100644 Binary files a/test-suite/csdp.cache and b/test-suite/csdp.cache differ diff --git a/test-suite/micromega/csdp.cache b/test-suite/micromega/csdp.cache index 9e620f734..961d159bf 100644 Binary files a/test-suite/micromega/csdp.cache and b/test-suite/micromega/csdp.cache differ -- cgit v1.2.3