diff options
Diffstat (limited to 'plugins/micromega/sos.ml')
| -rw-r--r-- | plugins/micromega/sos.ml | 1859 |
1 files changed, 1859 insertions, 0 deletions
diff --git a/plugins/micromega/sos.ml b/plugins/micromega/sos.ml new file mode 100644 index 00000000..3029496b --- /dev/null +++ b/plugins/micromega/sos.ml @@ -0,0 +1,1859 @@ +(* ========================================================================= *) +(* - 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 *) +(* - Frédéric Besson (fbesson@irisa.fr) is using it to feed micromega *) +(* ========================================================================= *) + +(* ========================================================================= *) +(* Nonlinear universal reals procedure using SOS decomposition. *) +(* ========================================================================= *) +open Num;; +open List;; +open Sos_types;; +open Sos_lib;; + +(* +prioritize_real();; +*) + +let debugging = ref false;; + +exception Sanity;; + +exception Unsolvable;; + +(* ------------------------------------------------------------------------- *) +(* Turn a rational into a decimal string with d sig digits. *) +(* ------------------------------------------------------------------------- *) + +let decimalize = + let rec normalize y = + if abs_num y </ Int 1 // Int 10 then normalize (Int 10 */ y) - 1 + else if abs_num y >=/ Int 1 then normalize (y // Int 10) + 1 + else 0 in + fun d x -> + if x =/ Int 0 then "0.0" else + let y = abs_num x in + let e = normalize y in + let z = pow10(-e) */ y +/ Int 1 in + let k = round_num(pow10 d */ z) in + (if x </ Int 0 then "-0." else "0.") ^ + 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. *) +(* ------------------------------------------------------------------------- *) + +let rec itern k l f a = + match l with + [] -> a + | h::t -> itern (k + 1) t f (f h k a);; + +let rec iter (m,n) f a = + if n < m then a + else iter (m+1,n) f (f m a);; + +(* ------------------------------------------------------------------------- *) +(* The main types. *) +(* ------------------------------------------------------------------------- *) + +type vector = int*(int,num)func;; + +type matrix = (int*int)*(int*int,num)func;; + +type monomial = (vname,int)func;; + +type poly = (monomial,num)func;; + +(* ------------------------------------------------------------------------- *) +(* Assignment avoiding zeros. *) +(* ------------------------------------------------------------------------- *) + +let (|-->) x y a = if y =/ Int 0 then a else (x |-> y) a;; + +(* ------------------------------------------------------------------------- *) +(* This can be generic. *) +(* ------------------------------------------------------------------------- *) + +let element (d,v) i = tryapplyd v i (Int 0);; + +let mapa f (d,v) = + d,foldl (fun a i c -> (i |--> f(c)) a) undefined v;; + +let is_zero (d,v) = + match v with + Empty -> true + | _ -> false;; + +(* ------------------------------------------------------------------------- *) +(* Vectors. Conventionally indexed 1..n. *) +(* ------------------------------------------------------------------------- *) + +let vector_0 n = (n,undefined:vector);; + +let dim (v:vector) = fst v;; + +let vector_const c n = + if c =/ Int 0 then vector_0 n + else (n,itlist (fun k -> k |-> c) (1--n) undefined :vector);; + +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) + +let vector_neg (v:vector) = (fst v,mapf minus_num (snd v) :vector);; + +let vector_add (v1:vector) (v2:vector) = + let m = dim v1 and n = dim v2 in + if m <> n then failwith "vector_add: incompatible dimensions" else + (n,combine (+/) (fun x -> x =/ Int 0) (snd v1) (snd 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);; + +(* ------------------------------------------------------------------------- *) +(* Matrices; again rows and columns indexed from 1. *) +(* ------------------------------------------------------------------------- *) + +let matrix_0 (m,n) = ((m,n),undefined:matrix);; + +let dimensions (m:matrix) = fst m;; + +let matrix_const c (m,n as mn) = + if m <> n then failwith "matrix_const: needs to be square" + else if c =/ Int 0 then matrix_0 mn + else (mn,itlist (fun k -> (k,k) |-> c) (1--n) undefined :matrix);; + +let matrix_1 = matrix_const (Int 1);; + +let matrix_cmul c (m:matrix) = + let (i,j) = dimensions m in + if c =/ Int 0 then matrix_0 (i,j) + else (i,j),mapf (fun x -> c */ x) (snd m);; + +let matrix_neg (m:matrix) = (dimensions m,mapf minus_num (snd m) :matrix);; + +let matrix_add (m1:matrix) (m2:matrix) = + let d1 = dimensions m1 and d2 = dimensions m2 in + if d1 <> d2 then failwith "matrix_add: incompatible dimensions" + else (d1,combine (+/) (fun x -> x =/ Int 0) (snd m1) (snd m2) :matrix);; + +let matrix_sub m1 m2 = matrix_add m1 (matrix_neg m2);; + +let row k (m:matrix) = + let i,j = dimensions m in + (j, + foldl (fun a (i,j) c -> if i = k then (j |-> c) a else a) undefined (snd m) + : vector);; + +let column k (m:matrix) = + let i,j = dimensions m in + (i, + foldl (fun a (i,j) c -> if j = k then (i |-> c) a else a) undefined (snd m) + : vector);; + +let transp (m:matrix) = + let i,j = dimensions m in + ((j,i),foldl (fun a (i,j) c -> ((j,i) |-> c) a) undefined (snd m) :matrix);; + +let diagonal (v:vector) = + let n = dim v in + ((n,n),foldl (fun a i c -> ((i,i) |-> c) a) undefined (snd v) : matrix);; + +let matrix_of_list l = + let m = length l in + if m = 0 then matrix_0 (0,0) else + let n = length (hd l) in + (m,n),itern 1 l (fun v i -> itern 1 v (fun c j -> (i,j) |-> c)) undefined;; + +(* ------------------------------------------------------------------------- *) +(* Monomials. *) +(* ------------------------------------------------------------------------- *) + +let monomial_eval assig (m:monomial) = + foldl (fun a x k -> a */ power_num (apply assig x) (Int k)) + (Int 1) m;; + +let monomial_1 = (undefined:monomial);; + +let monomial_var x = (x |=> 1 :monomial);; + +let (monomial_mul:monomial->monomial->monomial) = + combine (+) (fun x -> false);; + +let monomial_pow (m:monomial) k = + if k = 0 then monomial_1 + else mapf (fun x -> k * x) m;; + +let monomial_divides (m1:monomial) (m2:monomial) = + foldl (fun a x k -> tryapplyd m2 x 0 >= k & a) true m1;; + +let monomial_div (m1:monomial) (m2:monomial) = + let m = combine (+) (fun x -> x = 0) m1 (mapf (fun x -> -x) m2) in + if foldl (fun a x k -> k >= 0 & a) true m then m + else failwith "monomial_div: non-divisible";; + +let monomial_degree x (m:monomial) = tryapplyd m x 0;; + +let monomial_lcm (m1:monomial) (m2:monomial) = + (itlist (fun x -> x |-> max (monomial_degree x m1) (monomial_degree x m2)) + (union (dom m1) (dom m2)) undefined :monomial);; + +let monomial_multidegree (m:monomial) = foldl (fun a x k -> k + a) 0 m;; + +let monomial_variables m = dom m;; + +(* ------------------------------------------------------------------------- *) +(* Polynomials. *) +(* ------------------------------------------------------------------------- *) + +let eval assig (p:poly) = + foldl (fun a m c -> a +/ c */ monomial_eval assig m) (Int 0) p;; + +let poly_0 = (undefined:poly);; + +let poly_isconst (p:poly) = foldl (fun a m c -> m = monomial_1 & a) true p;; + +let poly_var x = ((monomial_var x) |=> Int 1 :poly);; + +let poly_const c = + if c =/ Int 0 then poly_0 else (monomial_1 |=> c);; + +let poly_cmul c (p:poly) = + if c =/ Int 0 then poly_0 + else mapf (fun x -> c */ x) p;; + +let poly_neg (p:poly) = (mapf minus_num p :poly);; + +let poly_add (p1:poly) (p2:poly) = + (combine (+/) (fun x -> x =/ Int 0) p1 p2 :poly);; + +let poly_sub p1 p2 = poly_add p1 (poly_neg p2);; + +let poly_cmmul (c,m) (p:poly) = + if c =/ Int 0 then poly_0 + else if m = monomial_1 then mapf (fun d -> c */ d) p + else foldl (fun a m' d -> (monomial_mul m m' |-> c */ d) a) poly_0 p;; + +let poly_mul (p1:poly) (p2:poly) = + foldl (fun a m c -> poly_add (poly_cmmul (c,m) p2) a) poly_0 p1;; + +let poly_div (p1:poly) (p2:poly) = + if not(poly_isconst p2) then failwith "poly_div: non-constant" else + let c = eval undefined p2 in + if c =/ Int 0 then failwith "poly_div: division by zero" + else poly_cmul (Int 1 // c) p1;; + +let poly_square p = poly_mul p p;; + +let rec poly_pow p k = + if k = 0 then poly_const (Int 1) + else if k = 1 then p + else let q = poly_square(poly_pow p (k / 2)) in + if k mod 2 = 1 then poly_mul p q else q;; + +let poly_exp p1 p2 = + if not(poly_isconst p2) then failwith "poly_exp: not a constant" else + poly_pow p1 (Num.int_of_num (eval undefined p2));; + +let degree x (p:poly) = foldl (fun a m c -> max (monomial_degree x m) a) 0 p;; + +let multidegree (p:poly) = + foldl (fun a m c -> max (monomial_multidegree m) a) 0 p;; + +let poly_variables (p:poly) = + foldr (fun m c -> union (monomial_variables m)) p [];; + +(* ------------------------------------------------------------------------- *) +(* Order monomials for human presentation. *) +(* ------------------------------------------------------------------------- *) + +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 + fun m1 m2 -> m1 = m2 or + ord (sort humanorder_varpow (graph m1)) + (sort humanorder_varpow (graph m2));; + +(* ------------------------------------------------------------------------- *) +(* Conversions to strings. *) +(* ------------------------------------------------------------------------- *) + +let string_of_vector min_size max_size (v:vector) = + let n_raw = dim v in + if n_raw = 0 then "[]" else + let n = max min_size (min n_raw max_size) in + let xs = map ((o) string_of_num (element v)) (1--n) in + "[" ^ end_itlist (fun s t -> s ^ ", " ^ t) xs ^ + (if n_raw > max_size then ", ...]" else "]");; + +let string_of_matrix max_size (m:matrix) = + let i_raw,j_raw = dimensions m in + let i = min max_size i_raw and j = min max_size j_raw in + let rstr = map (fun k -> string_of_vector j j (row k m)) (1--i) in + "["^end_itlist(fun s t -> s^";\n "^t) rstr ^ + (if j > max_size then "\n ...]" else "]");; + +let string_of_vname (v:vname): string = (v: string);; + +let rec string_of_term t = + match t with + Opp t1 -> "(- " ^ string_of_term t1 ^ ")" +| Add (t1, t2) -> + "(" ^ (string_of_term t1) ^ " + " ^ (string_of_term t2) ^ ")" +| Sub (t1, t2) -> + "(" ^ (string_of_term t1) ^ " - " ^ (string_of_term t2) ^ ")" +| Mul (t1, t2) -> + "(" ^ (string_of_term t1) ^ " * " ^ (string_of_term t2) ^ ")" +| Inv t1 -> "(/ " ^ string_of_term t1 ^ ")" +| Div (t1, t2) -> + "(" ^ (string_of_term t1) ^ " / " ^ (string_of_term t2) ^ ")" +| Pow (t1, n1) -> + "(" ^ (string_of_term t1) ^ " ^ " ^ (string_of_int n1) ^ ")" +| Zero -> "0" +| Var v -> "x" ^ (string_of_vname v) +| Const x -> string_of_num x;; + + +let string_of_varpow x k = + if k = 1 then string_of_vname x else string_of_vname x^"^"^string_of_int k;; + +let string_of_monomial m = + if m = monomial_1 then "1" else + let vps = List.fold_right (fun (x,k) a -> string_of_varpow x k :: a) + (sort humanorder_varpow (graph m)) [] in + end_itlist (fun s t -> s^"*"^t) vps;; + +let string_of_cmonomial (c,m) = + if m = monomial_1 then string_of_num c + else if c =/ Int 1 then string_of_monomial m + else string_of_num c ^ "*" ^ string_of_monomial m;; + +let string_of_poly (p:poly) = + if p = poly_0 then "<<0>>" else + let cms = sort (fun (m1,_) (m2,_) -> humanorder_monomial m1 m2) (graph p) in + let s = + List.fold_left (fun a (m,c) -> + if c </ Int 0 then a ^ " - " ^ string_of_cmonomial(minus_num c,m) + else a ^ " + " ^ string_of_cmonomial(c,m)) + "" cms in + let s1 = String.sub s 0 3 + and s2 = String.sub s 3 (String.length s - 3) in + "<<" ^(if s1 = " + " then s2 else "-"^s2)^">>";; + +(* ------------------------------------------------------------------------- *) +(* Printers. *) +(* ------------------------------------------------------------------------- *) + +let print_vector v = Format.print_string(string_of_vector 0 20 v);; + +let print_matrix m = Format.print_string(string_of_matrix 20 m);; + +let print_monomial m = Format.print_string(string_of_monomial m);; + +let print_poly m = Format.print_string(string_of_poly m);; + +(* +#install_printer print_vector;; +#install_printer print_matrix;; +#install_printer print_monomial;; +#install_printer print_poly;; +*) + +(* ------------------------------------------------------------------------- *) +(* Conversion from term. *) +(* ------------------------------------------------------------------------- *) + +let rec poly_of_term t = match t with + Zero -> poly_0 +| Const n -> poly_const n +| Var x -> poly_var x +| Opp t1 -> poly_neg (poly_of_term t1) +| Inv t1 -> + let p = poly_of_term t1 in + if poly_isconst p then poly_const(Int 1 // eval undefined p) + else failwith "poly_of_term: inverse of non-constant polyomial" +| Add (l, r) -> poly_add (poly_of_term l) (poly_of_term r) +| Sub (l, r) -> poly_sub (poly_of_term l) (poly_of_term r) +| Mul (l, r) -> poly_mul (poly_of_term l) (poly_of_term r) +| Div (l, r) -> + let p = poly_of_term l and q = poly_of_term r in + if poly_isconst q then poly_cmul (Int 1 // eval undefined q) p + else failwith "poly_of_term: division by non-constant polynomial" +| Pow (t, n) -> + poly_pow (poly_of_term t) n;; + +(* ------------------------------------------------------------------------- *) +(* String of vector (just a list of space-separated numbers). *) +(* ------------------------------------------------------------------------- *) + +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";; + +(* ------------------------------------------------------------------------- *) +(* String for block diagonal matrix numbered k. *) +(* ------------------------------------------------------------------------- *) + +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 "";; + +(* ------------------------------------------------------------------------- *) +(* String for a matrix numbered k, in SDPA sparse format. *) +(* ------------------------------------------------------------------------- *) + +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 "";; + +(* ------------------------------------------------------------------------- *) +(* String in SDPA sparse format for standard SDP problem: *) +(* *) +(* X = v_1 * [M_1] + ... + v_m * [M_m] - [M_0] must be PSD *) +(* Minimize obj_1 * v_1 + ... obj_m * v_m *) +(* ------------------------------------------------------------------------- *) + +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 "";; + +(* ------------------------------------------------------------------------- *) +(* More parser basics. *) +(* ------------------------------------------------------------------------- *) + +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 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,[x]) -> h +/ x | (h,_) -> h) 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,[x]) -> h */ power_num (Int 10) x | (h,_) -> h);; + +let mkparser p s = + let x,rst = p(explode s) in + if rst = [] then x else failwith "mkparser: unparsed input";; + +let parse_decimal = mkparser decimal;; + +(* ------------------------------------------------------------------------- *) +(* Parse back a vector. *) +(* ------------------------------------------------------------------------- *) + +let parse_sdpaoutput,parse_csdpoutput = + let vector = + token "{" ++ listof decimal (token ",") "decimal" ++ token "}" + >> (fun ((_,v),_) -> vector_of_list v) 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 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 ignore (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. *) +(* ------------------------------------------------------------------------- *) + +let csdp_default_parameters = +"axtol=1.0e-8 +atytol=1.0e-8 +objtol=1.0e-8 +pinftol=1.0e8 +dinftol=1.0e8 +maxiter=100 +minstepfrac=0.9 +maxstepfrac=0.97 +minstepp=1.0e-8 +minstepd=1.0e-8 +usexzgap=1 +tweakgap=0 +affine=0 +printlevel=1 +";; + +let csdp_params = csdp_default_parameters;; + +(* ------------------------------------------------------------------------- *) +(* Now call CSDP on a problem and parse back the output. *) +(* ------------------------------------------------------------------------- *) + +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;; + +(* ------------------------------------------------------------------------- *) +(* Try some apparently sensible scaling first. Note that this is purely to *) +(* get a cleaner translation to floating-point, and doesn't affect any of *) +(* the results, in principle. In practice it seems a lot better when there *) +(* are extreme numbers in the original problem. *) +(* ------------------------------------------------------------------------- *) + +let scale_then = + let common_denominator amat acc = + foldl (fun a m c -> lcm_num (denominator c) a) acc amat + and maximal_element amat acc = + foldl (fun maxa m c -> max_num maxa (abs_num c)) acc amat in + fun solver obj mats -> + let cd1 = itlist common_denominator mats (Int 1) + and cd2 = common_denominator (snd obj) (Int 1) in + let mats' = map (mapf (fun x -> cd1 */ x)) mats + and obj' = vector_cmul cd2 obj in + let max1 = itlist maximal_element mats' (Int 0) + and max2 = maximal_element (snd obj') (Int 0) in + let scal1 = pow2 (20-int_of_float(log(float_of_num max1) /. log 2.0)) + and scal2 = pow2 (20-int_of_float(log(float_of_num max2) /. log 2.0)) in + let mats'' = map (mapf (fun x -> x */ scal1)) mats' + and obj'' = vector_cmul scal2 obj' in + solver obj'' mats'';; + +(* ------------------------------------------------------------------------- *) +(* Round a vector to "nice" rationals. *) +(* ------------------------------------------------------------------------- *) + +let nice_rational n x = round_num (n */ x) // n;; + +let nice_vector n = mapa (nice_rational n);; + +(* ------------------------------------------------------------------------- *) +(* Reduce linear program to SDP (diagonal matrices) and test with CSDP. This *) +(* one tests A [-1;x1;..;xn] >= 0 (i.e. left column is negated constants). *) +(* ------------------------------------------------------------------------- *) + +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 + 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 *) +(* computational geometry, express as linear inequalities and call CSDP. *) +(* This is a bit lazy of me, but it's easy and not such a bottleneck so far. *) +(* ------------------------------------------------------------------------- *) + +let in_convex_hull pts pt = + let pts1 = (1::pt) :: map (fun x -> 1::x) pts in + let pts2 = map (fun p -> map (fun x -> -x) p @ p) pts1 in + let n = length pts + 1 + and v = 2 * (length pt + 1) in + let m = v + n - 1 in + let mat = + (m,n), + itern 1 pts2 (fun pts j -> itern 1 pts (fun x i -> (i,j) |-> Int x)) + (iter (1,n) (fun i -> (v + i,i+1) |-> Int 1) undefined) in + linear_program_basic mat;; + +(* ------------------------------------------------------------------------- *) +(* Filter down a set of points to a minimal set with the same convex hull. *) +(* ------------------------------------------------------------------------- *) + +let minimal_convex_hull = + let augment1 = function + | [] -> assert false + | (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 + funpow (length mons') augment1 mons';; + +(* ------------------------------------------------------------------------- *) +(* Stuff for "equations" (generic A->num functions). *) +(* ------------------------------------------------------------------------- *) + +let equation_cmul c eq = + if c =/ Int 0 then Empty else mapf (fun d -> c */ d) eq;; + +let equation_add eq1 eq2 = combine (+/) (fun x -> x =/ Int 0) eq1 eq2;; + +let equation_eval assig eq = + let value v = apply assig v in + foldl (fun a v c -> a +/ value(v) */ c) (Int 0) eq;; + +(* ------------------------------------------------------------------------- *) +(* Eliminate among linear equations: return unconstrained variables and *) +(* assignments for the others in terms of them. We give one pseudo-variable *) +(* "one" that's used for a constant term. *) +(* ------------------------------------------------------------------------- *) + +let failstore = ref [];; + +let eliminate_equations = + let rec extract_first p l = + match l with + [] -> failwith "extract_first" + | h::t -> if p(h) then h,t else + let k,s = extract_first p t in + k,h::s in + let rec eliminate vars dun eqs = + match vars with + [] -> if forall is_undefined eqs then dun + 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 + let eq' = equation_cmul (Int(-1) // a) (undefine v eq) in + let elim e = + let b = tryapplyd e v (Int 0) in + if b =/ Int 0 then e else + equation_add e (equation_cmul (minus_num b // a) eq) in + eliminate vs ((v |-> eq') (mapf elim dun)) (map elim oeqs) + with Failure _ -> eliminate vs dun eqs in + fun one vars eqs -> + let assig = eliminate vars undefined eqs in + let vs = foldl (fun a x f -> subtract (dom f) [one] @ a) [] assig in + setify vs,assig;; + +(* ------------------------------------------------------------------------- *) +(* Eliminate all variables, in an essentially arbitrary order. *) +(* ------------------------------------------------------------------------- *) + +let eliminate_all_equations one = + let choose_variable eq = + let (v,_) = choose eq in + if v = one then + let eq' = undefine v eq in + if is_undefined eq' then failwith "choose_variable" else + let (w,_) = choose eq' in w + else v in + let rec eliminate dun eqs = + match eqs with + [] -> dun + | eq::oeqs -> + if is_undefined eq then eliminate dun oeqs else + let v = choose_variable eq in + let a = apply eq v in + let eq' = equation_cmul (Int(-1) // a) (undefine v eq) in + let elim e = + let b = tryapplyd e v (Int 0) in + if b =/ Int 0 then e else + equation_add e (equation_cmul (minus_num b // a) eq) in + eliminate ((v |-> eq') (mapf elim dun)) (map elim oeqs) in + fun eqs -> + let assig = eliminate undefined eqs in + let vs = foldl (fun a x f -> subtract (dom f) [one] @ a) [] assig in + setify vs,assig;; + +(* ------------------------------------------------------------------------- *) +(* Solve equations by assigning arbitrary numbers. *) +(* ------------------------------------------------------------------------- *) + +let solve_equations one eqs = + let vars,assigs = eliminate_all_equations one eqs in + let vfn = itlist (fun v -> (v |-> Int 0)) vars (one |=> Int(-1)) in + let ass = + combine (+/) (fun c -> false) (mapf (equation_eval vfn) assigs) vfn in + if forall (fun e -> equation_eval ass e =/ Int 0) eqs + then undefine one ass else raise Sanity;; + +(* ------------------------------------------------------------------------- *) +(* Hence produce the "relevant" monomials: those whose squares lie in the *) +(* Newton polytope of the monomials in the input. (This is enough according *) +(* to Reznik: "Extremal PSD forms with few terms", Duke Math. Journal, *) +(* vol 45, pp. 363--374, 1978. *) +(* *) +(* These are ordered in sort of decreasing degree. In particular the *) +(* constant monomial is last; this gives an order in diagonalization of the *) +(* quadratic form that will tend to display constants. *) +(* ------------------------------------------------------------------------- *) + +let newton_polytope pol = + let vars = poly_variables pol in + let mons = map (fun m -> map (fun x -> monomial_degree x m) vars) (dom pol) + and ds = map (fun x -> (degree x pol + 1) / 2) vars in + let all = itlist (fun n -> allpairs (fun h t -> h::t) (0--n)) ds [[]] + and mons' = minimal_convex_hull mons in + let all' = + filter (fun m -> in_convex_hull mons' (map (fun x -> 2 * x) m)) all in + map (fun m -> itlist2 (fun v i a -> if i = 0 then a else (v |-> i) a) + vars m monomial_1) (rev all');; + +(* ------------------------------------------------------------------------- *) +(* Diagonalize (Cholesky/LDU) the matrix corresponding to a quadratic form. *) +(* ------------------------------------------------------------------------- *) + +let diag m = + let nn = dimensions m in + let n = fst nn in + if snd nn <> n then failwith "diagonalize: non-square matrix" else + let rec diagonalize i m = + if is_zero m then [] else + let a11 = element m (i,i) in + if a11 </ Int 0 then failwith "diagonalize: not PSD" + else if a11 =/ Int 0 then + if is_zero(row i m) then diagonalize (i + 1) m + else failwith "diagonalize: not PSD" + else + let v = row i m in + let v' = mapa (fun a1k -> a1k // a11) v in + let m' = + (n,n), + iter (i+1,n) (fun j -> + iter (i+1,n) (fun k -> + ((j,k) |--> (element m (j,k) -/ element v j */ element v' k)))) + undefined in + (a11,v')::diagonalize (i + 1) m' in + diagonalize 1 m;; + +(* ------------------------------------------------------------------------- *) +(* Adjust a diagonalization to collect rationals at the start. *) +(* ------------------------------------------------------------------------- *) + +let deration d = + if d = [] then Int 0,d else + let adj(c,l) = + let a = foldl (fun a i c -> lcm_num a (denominator c)) (Int 1) (snd l) // + 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 + (Int 1 // a),map (fun (c,l) -> (a */ c,l)) d';; + +(* ------------------------------------------------------------------------- *) +(* Enumeration of monomials with given multidegree bound. *) +(* ------------------------------------------------------------------------- *) + +let rec enumerate_monomials d vars = + if d < 0 then [] + else if d = 0 then [undefined] + else if vars = [] then [monomial_1] else + let alts = + map (fun k -> let oths = enumerate_monomials (d - k) (tl vars) in + map (fun ks -> if k = 0 then ks else (hd vars |-> k) ks) oths) + (0--d) in + end_itlist (@) alts;; + +(* ------------------------------------------------------------------------- *) +(* Enumerate products of distinct input polys with degree <= d. *) +(* We ignore any constant input polynomials. *) +(* Give the output polynomial and a record of how it was derived. *) +(* ------------------------------------------------------------------------- *) + +let rec enumerate_products d pols = + if d = 0 then [poly_const num_1,Rational_lt num_1] else if d < 0 then [] else + match pols with + [] -> [poly_const num_1,Rational_lt num_1] + | (p,b)::ps -> let e = multidegree p in + if e = 0 then enumerate_products d ps else + enumerate_products d ps @ + map (fun (q,c) -> poly_mul p q,Product(b,c)) + (enumerate_products (d - e) ps);; + +(* ------------------------------------------------------------------------- *) +(* Multiply equation-parametrized poly by regular poly and add accumulator. *) +(* ------------------------------------------------------------------------- *) + +let epoly_pmul p q acc = + foldl (fun a m1 c -> + foldl (fun b m2 e -> + let m = monomial_mul m1 m2 in + let es = tryapplyd b m undefined in + (m |-> equation_add (equation_cmul c e) es) b) + a q) acc p;; + +(* ------------------------------------------------------------------------- *) +(* Usual operations on equation-parametrized poly. *) +(* ------------------------------------------------------------------------- *) + +let epoly_cmul c l = + if c =/ Int 0 then undefined else mapf (equation_cmul c) l;; + +let epoly_neg = epoly_cmul (Int(-1));; + +let epoly_add = combine equation_add is_undefined;; + +let epoly_sub p q = epoly_add p (epoly_neg q);; + +(* ------------------------------------------------------------------------- *) +(* Convert regular polynomial. Note that we treat (0,0,0) as -1. *) +(* ------------------------------------------------------------------------- *) + +let epoly_of_poly p = + foldl (fun a m c -> (m |-> ((0,0,0) |=> minus_num c)) a) undefined p;; + +(* ------------------------------------------------------------------------- *) +(* String for block diagonal matrix numbered k. *) +(* ------------------------------------------------------------------------- *) + +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 "";; + +(* ------------------------------------------------------------------------- *) +(* SDPA for problem using block diagonal (i.e. multiple SDPs) *) +(* ------------------------------------------------------------------------- *) + +let sdpa_of_blockproblem comment nblocks blocksizes obj mats = + let m = length mats - 1 in + "\"" ^ comment ^ "\"\n" ^ + string_of_int m ^ "\n" ^ + string_of_int nblocks ^ "\n" ^ + (end_itlist (fun s t -> s^" "^t) (map string_of_int blocksizes)) ^ + "\n" ^ + sdpa_of_vector obj ^ + itlist2 (fun k m a -> sdpa_of_blockdiagonal (k - 1) m ^ a) + (1--length mats) mats "";; + +(* ------------------------------------------------------------------------- *) +(* Hence run CSDP on a problem in block diagonal form. *) +(* ------------------------------------------------------------------------- *) + +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. *) +(* ------------------------------------------------------------------------- *) + +let bmatrix_add = combine (+/) (fun x -> x =/ Int 0);; + +let bmatrix_cmul c bm = + if c =/ Int 0 then undefined + else mapf (fun x -> c */ x) bm;; + +let bmatrix_neg = bmatrix_cmul (Int(-1));; + +let bmatrix_sub m1 m2 = bmatrix_add m1 (bmatrix_neg m2);; + +(* ------------------------------------------------------------------------- *) +(* Smash a block matrix into components. *) +(* ------------------------------------------------------------------------- *) + +let blocks blocksizes bm = + map (fun (bs,b0) -> + let m = foldl + (fun a (b,i,j) c -> if b = b0 then ((i,j) |-> c) a else a) + undefined bm in + (((bs,bs),m):matrix)) + (zip blocksizes (1--length blocksizes));; + +(* ------------------------------------------------------------------------- *) +(* Positiv- and Nullstellensatz. Flag "linf" forces a linear representation. *) +(* ------------------------------------------------------------------------- *) + +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 + (poly_const num_1,Rational_lt num_1):: + (filter (fun (p,c) -> multidegree p <= d) leqs) + else enumerate_products d leqs in + let nblocks = length monoid in + let mk_idmultiplier k p = + let e = d - multidegree p in + let mons = enumerate_monomials e vars in + let nons = zip mons (1--length mons) in + mons, + itlist (fun (m,n) -> (m |-> ((-k,-n,n) |=> Int 1))) nons undefined in + let mk_sqmultiplier k (p,c) = + let e = (d - multidegree p) / 2 in + let mons = enumerate_monomials e vars in + let nons = zip mons (1--length mons) in + mons, + itlist (fun (m1,n1) -> + itlist (fun (m2,n2) a -> + let m = monomial_mul m1 m2 in + if n1 > n2 then a else + let c = if n1 = n2 then Int 1 else Int 2 in + let e = tryapplyd a m undefined in + (m |-> equation_add ((k,n1,n2) |=> c) e) a) + nons) + nons undefined in + let sqmonlist,sqs = unzip(map2 mk_sqmultiplier (1--length monoid) monoid) + and idmonlist,ids = unzip(map2 mk_idmultiplier (1--length eqs) eqs) in + let blocksizes = map length sqmonlist in + let bigsum = + itlist2 (fun p q a -> epoly_pmul p q a) eqs ids + (itlist2 (fun (p,c) s a -> epoly_pmul p s a) monoid sqs + (epoly_of_poly(poly_neg pol))) in + let eqns = foldl (fun a m e -> e::a) [] bigsum in + let pvs,assig = eliminate_all_equations (0,0,0) eqns in + let qvars = (0,0,0)::pvs in + let allassig = itlist (fun v -> (v |-> (v |=> Int 1))) pvs assig in + let mk_matrix v = + foldl (fun m (b,i,j) ass -> if b < 0 then m else + let c = tryapplyd ass v (Int 0) in + if c =/ Int 0 then m else + ((b,j,i) |-> c) (((b,i,j) |-> c) m)) + undefined allassig in + let diagents = foldl + (fun a (b,i,j) e -> if b > 0 & i = j then equation_add e a else a) + undefined allassig in + let mats = map mk_matrix qvars + and obj = length pvs, + 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 nblocks blocksizes) obj mats in + let find_rounding d = + (if !debugging then + (Format.print_string("Trying rounding with limit "^string_of_num d); + Format.print_newline()) + else ()); + let vec = nice_vector d raw_vec in + let blockmat = iter (1,dim vec) + (fun i a -> bmatrix_add (bmatrix_cmul (element vec i) (el i mats)) a) + (bmatrix_neg (el 0 mats)) in + let allmats = blocks blocksizes blockmat in + vec,map diag allmats in + let vec,ratdias = + if pvs = [] then find_rounding num_1 + else tryfind find_rounding (map Num.num_of_int (1--31) @ + map pow2 (5--66)) in + let newassigs = + itlist (fun k -> el (k - 1) pvs |-> element vec k) + (1--dim vec) ((0,0,0) |=> Int(-1)) in + let finalassigs = + foldl (fun a v e -> (v |-> equation_eval newassigs e) a) newassigs + allassig in + let poly_of_epoly p = + foldl (fun a v e -> (v |--> equation_eval finalassigs e) a) + undefined p in + let mk_sos mons = + let mk_sq (c,m) = + c,itlist (fun k a -> (el (k - 1) mons |--> element m k) a) + (1--length mons) undefined in + map mk_sq in + let sqs = map2 mk_sos sqmonlist ratdias + and cfs = map poly_of_epoly ids in + let msq = filter (fun (a,b) -> b <> []) (map2 (fun a b -> a,b) monoid sqs) in + let eval_sq sqs = itlist + (fun (c,q) -> poly_add (poly_cmul c (poly_mul q q))) sqs poly_0 in + let sanity = + itlist (fun ((p,c),s) -> poly_add (poly_mul p (eval_sq s))) msq + (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;; + +(* ------------------------------------------------------------------------- *) +(* 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);; + +(* ------------------------------------------------------------------------- *) +(* The ordering so we can create canonical HOL polynomials. *) +(* ------------------------------------------------------------------------- *) + +let dest_monomial mon = sort (increasing fst) (graph mon);; + +let monomial_order = + let rec lexorder l1 l2 = + match (l1,l2) with + [],[] -> true + | vps,[] -> false + | [],vps -> true + | ((x1,n1)::vs1),((x2,n2)::vs2) -> + if x1 < x2 then true + else if x2 < x1 then false + else if n1 < n2 then false + else if n2 < n1 then true + else lexorder vs1 vs2 in + fun m1 m2 -> + if m2 = monomial_1 then true else if m1 = monomial_1 then false else + let mon1 = dest_monomial m1 and mon2 = dest_monomial m2 in + let deg1 = itlist ((o) (+) snd) mon1 0 + and deg2 = itlist ((o) (+) snd) mon2 0 in + if deg1 < deg2 then false else if deg1 > deg2 then true + else lexorder mon1 mon2;; + +let dest_poly p = + map (fun (m,c) -> c,dest_monomial m) + (sort (fun (m1,_) (m2,_) -> monomial_order m1 m2) (graph p));; + +(* ------------------------------------------------------------------------- *) +(* Map back polynomials and their composites to HOL. *) +(* ------------------------------------------------------------------------- *) + +let term_of_varpow = + fun x k -> + if k = 1 then Var x else Pow (Var x, k);; + +let term_of_monomial = + fun m -> if m = monomial_1 then Const num_1 else + let m' = dest_monomial m in + let vps = itlist (fun (x,k) a -> term_of_varpow x k :: a) m' [] in + end_itlist (fun s t -> Mul (s,t)) vps;; + +let term_of_cmonomial = + fun (m,c) -> + if m = monomial_1 then Const c + else if c =/ num_1 then term_of_monomial m + else Mul (Const c,term_of_monomial m);; + +let term_of_poly = + fun p -> + if p = poly_0 then Zero else + let cms = map term_of_cmonomial + (sort (fun (m1,_) (m2,_) -> monomial_order m1 m2) (graph p)) in + end_itlist (fun t1 t2 -> Add (t1,t2)) cms;; + +let term_of_sqterm (c,p) = + Product(Rational_lt c,Square(term_of_poly p));; + +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));; + +(* ------------------------------------------------------------------------- *) +(* 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 </ num_0 -> 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. *) +(* ------------------------------------------------------------------------- *) + +let rec allpermutations l = + if l = [] then [[]] else + itlist (fun h acc -> map (fun t -> h::t) + (allpermutations (subtract l [h])) @ acc) l [];; + +let allvarorders l = + map (fun vlis x -> index x vlis) (allpermutations l);; + +let changevariables_monomial zoln (m:monomial) = + foldl (fun a x k -> (assoc x zoln |-> k) a) monomial_1 m;; + +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_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. *) +(* ------------------------------------------------------------------------- *) + +let sumofsquares_general_symmetry tool pol = + let vars = poly_variables pol + and lpps = newton_polytope pol in + let n = length lpps in + let sym_eqs = + let invariants = filter + (fun vars' -> + is_undefined(poly_sub pol (changevariables (zip vars vars') pol))) + (allpermutations vars) in + let lpns = zip lpps (1--length lpps) in + let lppcs = + filter (fun (m,(n1,n2)) -> n1 <= n2) + (allpairs + (fun (m1,n1) (m2,n2) -> (m1,m2),(n1,n2)) lpns lpns) in + let clppcs = end_itlist (@) + (map (fun ((m1,m2),(n1,n2)) -> + map (fun vars' -> + (changevariables_monomial (zip vars vars') m1, + changevariables_monomial (zip vars vars') m2),(n1,n2)) + invariants) + lppcs) in + let clppcs_dom = setify(map fst clppcs) in + let clppcs_cls = map (fun d -> filter (fun (e,_) -> e = d) clppcs) + clppcs_dom in + let eqvcls = map (o setify (map snd)) clppcs_cls in + let mk_eq cls acc = + match cls with + [] -> raise Sanity + | [h] -> acc + | h::t -> map (fun k -> (k |-> Int(-1)) (h |=> Int 1)) t @ acc in + itlist mk_eq eqvcls [] in + let eqs = foldl (fun a x y -> y::a) [] + (itern 1 lpps (fun m1 n1 -> + itern 1 lpps (fun m2 n2 f -> + let m = monomial_mul m1 m2 in + if n1 > n2 then f else + let c = if n1 = n2 then Int 1 else Int 2 in + (m |-> ((n1,n2) |-> c) (tryapplyd f m undefined)) f)) + (foldl (fun a m c -> (m |-> ((0,0)|=>c)) a) + undefined pol)) @ + sym_eqs in + let pvs,assig = eliminate_all_equations (0,0) eqs in + let allassig = itlist (fun v -> (v |-> (v |=> Int 1))) pvs assig in + let qvars = (0,0)::pvs in + let diagents = + end_itlist equation_add (map (fun i -> apply allassig (i,i)) (1--n)) in + let mk_matrix v = + ((n,n), + foldl (fun m (i,j) ass -> let c = tryapplyd ass v (Int 0) in + if c =/ Int 0 then m else + ((j,i) |-> c) (((i,j) |-> c) m)) + undefined allassig :matrix) in + let mats = map mk_matrix qvars + and obj = length pvs, + itern 1 pvs (fun v i -> (i |--> tryapplyd diagents v (Int 0))) + undefined in + let raw_vec = if pvs = [] then vector_0 0 else tool obj mats in + let find_rounding d = + (if !debugging then + (Format.print_string("Trying rounding with limit "^string_of_num d); + Format.print_newline()) + else ()); + let vec = nice_vector d raw_vec in + let mat = iter (1,dim vec) + (fun i a -> matrix_add (matrix_cmul (element vec i) (el i mats)) a) + (matrix_neg (el 0 mats)) in + deration(diag mat) in + let rat,dia = + if pvs = [] then + let mat = matrix_neg (el 0 mats) in + deration(diag mat) + else + tryfind find_rounding (map Num.num_of_int (1--31) @ + map pow2 (5--66)) in + let poly_of_lin(d,v) = + d,foldl(fun a i c -> (el (i - 1) lpps |-> c) a) undefined (snd v) in + let lins = map poly_of_lin dia in + let sqs = map (fun (d,l) -> poly_mul (poly_const d) (poly_pow l 2)) lins in + 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 = 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`;; + +*****) |