1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
|
(************************************************************************)
(* v * The Coq Proof Assistant / The Coq Development Team *)
(* <O___,, * CNRS-Ecole Polytechnique-INRIA Futurs-Universite Paris Sud *)
(* \VV/ **************************************************************)
(* // * This file is distributed under the terms of the *)
(* * GNU Lesser General Public License Version 2.1 *)
(************************************************************************)
(* *)
(* Micromega: A reflexive tactic using the Positivstellensatz *)
(* *)
(* Frédéric Besson (Irisa/Inria) 2006-2008 *)
(* *)
(************************************************************************)
open Big_int
open Num
open Sos
module Mc = Micromega
module Ml2C = Mutils.CamlToCoq
module C2Ml = Mutils.CoqToCaml
let debug = false
module M =
struct
open Mc
let rec expr_to_term = function
| PEc z -> Const (C2Ml.q_to_num z)
| PEX v -> Var ("x"^(string_of_int (C2Ml.index v)))
| PEmul(p1,p2) ->
let p1 = expr_to_term p1 in
let p2 = expr_to_term p2 in
let res = Mul(p1,p2) in res
| PEadd(p1,p2) -> Add(expr_to_term p1, expr_to_term p2)
| PEsub(p1,p2) -> Sub(expr_to_term p1, expr_to_term p2)
| PEpow(p,n) -> Pow(expr_to_term p , C2Ml.n n)
| PEopp p -> Opp (expr_to_term p)
(* let term_to_expr e =
let e' = term_to_expr e in
if debug
then Printf.printf "term_to_expr : %s - %s\n"
(string_of_poly (poly_of_term e))
(string_of_poly (poly_of_term (expr_to_term e')));
e' *)
end
open M
open List
open Mutils
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 ;"
(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"
let partition_expr l =
let rec f i = function
| [] -> ([],[],[])
| (e,k)::l ->
let (eq,ge,neq) = f (i+1) l in
match k with
| Mc.Equal -> ((e,i)::eq,ge,neq)
| Mc.NonStrict -> (eq,(e,Axiom_le i)::ge,neq)
| Mc.Strict -> (* e > 0 == e >= 0 /\ e <> 0 *)
(eq, (e,Axiom_lt i)::ge,(e,Axiom_lt i)::neq)
| Mc.NonEqual -> (eq,ge,(e,Axiom_eq i)::neq)
(* Not quite sure -- Coq interface has changed *)
in f 0 l
let rec sets_of_list l =
match l with
| [] -> [[]]
| e::l -> let s = sets_of_list l in
s@(List.map (fun s0 -> e::s0) s)
(* The exploration is probably not complete - for simple cases, it works... *)
let real_nonlinear_prover d l =
try
let (eq,ge,neq) = partition_expr l in
let rec elim_const = function
[] -> []
| (x,y)::l -> let p = poly_of_term (expr_to_term x) in
if poly_isconst p
then elim_const l
else (p,y)::(elim_const l) in
let eq = elim_const eq in
let peq = List.map fst eq in
let pge = List.map
(fun (e,psatz) -> poly_of_term (expr_to_term e),psatz) ge in
let monoids = List.map (fun m -> (List.fold_right (fun (p,kd) y ->
let p = poly_of_term (expr_to_term p) in
match kd with
| Axiom_lt i -> poly_mul p y
| Axiom_eq i -> poly_mul (poly_pow p 2) y
| _ -> failwith "monoids") m (poly_const (Int 1)) , map snd m))
(sets_of_list neq) in
let (cert_ideal, cert_cone,monoid) = deepen_until d (fun d ->
list_try_find (fun m -> let (ci,cc) =
real_positivnullstellensatz_general false d peq pge (poly_neg (fst m) ) in
(ci,cc,snd m)) monoids) 0 in
let proofs_ideal = map2 (fun q i -> Eqmul(term_of_poly q,Axiom_eq i))
cert_ideal (List.map snd eq) in
let proofs_cone = map term_of_sos cert_cone in
let proof_ne =
let (neq , lt) = List.partition
(function Axiom_eq _ -> true | _ -> false ) monoid in
let sq = match
(List.map (function Axiom_eq i -> i | _ -> failwith "error") neq)
with
| [] -> Rational_lt (Int 1)
| l -> Monoid l in
List.fold_right (fun x y -> Product(x,y)) lt sq in
let proof = list_fold_right_elements
(fun s t -> Sum(s,t)) (proof_ne :: proofs_ideal @ proofs_cone) in
Some proof
with
| Sos.TooDeep -> None
(* This is somewhat buggy, over Z, strict inequality vanish... *)
let pure_sos l =
(* If there is no strict inequality,
I should nonetheless be able to try something - over Z > is equivalent to -1 >= *)
try
let l = List.combine l (interval 0 (length l -1)) in
let (lt,i) = try (List.find (fun (x,_) -> snd x = Mc.Strict) l)
with Not_found -> List.hd l in
let plt = poly_neg (poly_of_term (expr_to_term (fst lt))) in
let (n,polys) = sumofsquares plt in (* n * (ci * pi^2) *)
let pos = Product (Rational_lt n,
List.fold_right (fun (c,p) rst -> Sum (Product (Rational_lt c, Square
(term_of_poly p)), rst))
polys (Rational_lt (Int 0))) in
let proof = Sum(Axiom_lt i, pos) in
(* let s,proof' = scale_certificate proof in
let cert = snd (cert_of_pos proof') in *)
Some proof
with
| Not_found -> (* This is no strict inequality *) None
| x -> None
type micromega_polys = (Micromega.q Mc.pExpr * Mc.op1) list
type csdp_certificate = Sos.positivstellensatz option
type provername = string * int option
let main () =
if Array.length Sys.argv <> 3 then
(Printf.printf "Usage: csdpcert inputfile outputfile\n"; exit 1);
let input_file = Sys.argv.(1) in
let output_file = Sys.argv.(2) in
let inch = open_in input_file in
let (prover,poly) = (input_value inch : provername * micromega_polys) in
close_in inch;
let cert =
match prover with
| "real_nonlinear_prover", Some d -> real_nonlinear_prover d poly
| "pure_sos", None -> pure_sos poly
| prover, _ -> (Printf.printf "unknown prover: %s\n" prover; exit 1) in
let outch = open_out output_file in
output_value outch (cert:csdp_certificate);
close_out outch;
exit 0;;
let _ = main () in ()
|
