1
(************************************************************************)
2
(* v * The Coq Proof Assistant / The Coq Development Team *)
3
(* <O___,, * CNRS-Ecole Polytechnique-INRIA Futurs-Universite Paris Sud *)
4
(* \VV/ **************************************************************)
5
(* // * This file is distributed under the terms of the *)
6
(* * GNU Lesser General Public License Version 2.1 *)
7
(************************************************************************)
9
(* Micromega: A reflexive tactic using the Positivstellensatz *)
11
(* Frédéric Besson (Irisa/Inria) 2006-2008 *)
13
(************************************************************************)
15
(* We take as input a list of polynomials [p1...pn] and return an unfeasibility
16
certificate polynomial. *)
18
(*open Micromega.Polynomial*)
23
module Ml2C = Mutils.CamlToCoq
24
module C2Ml = Mutils.CoqToCaml
30
type var = Mc.positive
37
val find : var -> t -> int
38
val mult : var -> t -> t
39
val prod : t -> t -> t
40
val compare : t -> t -> int
41
val pp : out_channel -> t -> unit
42
val fold : (var -> int -> 'a -> 'a) -> t -> 'a -> 'a
46
(* A monomial is represented by a multiset of variables *)
47
module Map = Map.Make(struct type t = var let compare = Pervasives.compare end)
52
(* The monomial that corresponds to a constant *)
55
(* The monomial 'x' *)
56
let var x = Map.add x 1 Map.empty
58
(* Get the degre of a variable in a monomial *)
59
let find x m = try find x m with Not_found -> 0
61
(* Multiply a monomial by a variable *)
62
let mult x m = add x ( (find x m) + 1) m
64
(* Product of monomials *)
65
let prod m1 m2 = Map.fold (fun k d m -> add k ((find k m) + d) m) m1 m2
67
(* Total ordering of monomials *)
68
let compare m1 m2 = Map.compare Pervasives.compare m1 m2
70
let pp o m = Map.iter (fun k v ->
71
if v = 1 then Printf.fprintf o "x%i." (C2Ml.index k)
72
else Printf.fprintf o "x%i^%i." (C2Ml.index k) v) m
80
(* A polynomial is a map of monomials *)
82
This is probably a naive implementation
83
(expected to be fast enough - Coq is probably the bottleneck)
84
*The new ring contribution is using a sparse Horner representation.
88
val get : Monomial.t -> t -> num
89
val variable : var -> t
90
val add : Monomial.t -> num -> t -> t
91
val constant : num -> t
92
val mult : Monomial.t -> num -> t -> t
93
val product : t -> t -> t
94
val addition : t -> t -> t
96
val fold : (Monomial.t -> num -> 'a -> 'a) -> t -> 'a -> 'a
97
val pp : out_channel -> t -> unit
98
val compare : t -> t -> int
101
(*normalisation bug : 0*x ... *)
102
module P = Map.Make(Monomial)
107
let pp o p = P.iter (fun k v ->
108
if compare_num v (Int 0) <> 0
110
if Monomial.compare Monomial.const k = 0
111
then Printf.fprintf o "%s " (string_of_num v)
112
else Printf.fprintf o "%s*%a " (string_of_num v) Monomial.pp k) p
114
(* Get the coefficient of monomial mn *)
115
let get : Monomial.t -> t -> num =
116
fun mn p -> try find mn p with Not_found -> (Int 0)
119
(* The polynomial 1.x *)
120
let variable : var -> t =
121
fun x -> add (Monomial.var x) (Int 1) empty
123
(*The constant polynomial *)
124
let constant : num -> t =
125
fun c -> add (Monomial.const) c empty
127
(* The addition of a monomial *)
129
let add : Monomial.t -> num -> t -> t =
131
let vl = (get mn p) <+> v in
135
(** Design choice: empty is not a polynomial
136
I do not remember why ....
139
(* The product by a monomial *)
140
let mult : Monomial.t -> num -> t -> t =
142
fold (fun mn' v' res -> P.add (Monomial.prod mn mn') (v<*>v') res) p empty
145
let addition : t -> t -> t =
146
fun p1 p2 -> fold (fun mn v p -> add mn v p) p1 p2
149
let product : t -> t -> t =
151
fold (fun mn v res -> addition (mult mn v p2) res ) p1 empty
154
let uminus : t -> t =
155
fun p -> map (fun v -> minus_num v) p
159
let compare = compare compare_num
163
type 'a number_spec = {
164
bigint_to_number : big_int -> 'a;
165
number_to_num : 'a -> num;
168
mult : 'a -> 'a -> 'a;
169
eqb : 'a -> 'a -> Mc.bool
173
bigint_to_number = Ml2C.bigint ;
174
number_to_num = (fun x -> Big_int (C2Ml.z_big_int x));
176
unit = Mc.Zpos Mc.XH;
183
bigint_to_number = (fun x -> {Mc.qnum = Ml2C.bigint x; Mc.qden = Mc.XH});
184
number_to_num = C2Ml.q_to_num;
185
zero = {Mc.qnum = Mc.Z0;Mc.qden = Mc.XH};
186
unit = {Mc.qnum = (Mc.Zpos Mc.XH) ; Mc.qden = Mc.XH};
196
let dev_form n_spec p =
199
| Mc.PEc z -> Poly.constant (n_spec.number_to_num z)
200
| Mc.PEX v -> Poly.variable v
202
let p1 = dev_form p1 in
203
let p2 = dev_form p2 in
205
| Mc.PEadd(p1,p2) -> Poly.addition (dev_form p1) (dev_form p2)
206
| Mc.PEopp p -> Poly.uminus (dev_form p)
207
| Mc.PEsub(p1,p2) -> Poly.addition (dev_form p1) (Poly.uminus (dev_form p2))
209
let p = dev_form p in
213
then Poly.constant (n_spec.number_to_num n_spec.unit)
214
else Poly.product p (pow (n-1)) in
219
let monomial_to_polynomial mn =
222
let mn = if i = 1 then Mc.PEX v else Mc.PEpow (Mc.PEX v ,Ml2C.n i) in
223
if acc = Mc.PEc (Mc.Zpos Mc.XH)
225
else Mc.PEmul(mn,acc))
227
(Mc.PEc (Mc.Zpos Mc.XH))
229
let list_to_polynomial vars l =
230
assert (List.for_all (fun x -> ceiling_num x =/ x) l);
231
let var x = monomial_to_polynomial (List.nth vars x) in
232
let rec xtopoly p i = function
234
| c::l -> if c =/ (Int 0) then xtopoly p (i+1) l
235
else let c = Mc.PEc (Ml2C.bigint (numerator c)) in
237
if c = Mc.PEc (Mc.Zpos Mc.XH)
239
else Mc.PEmul (c,var i) in
240
let p' = if p = Mc.PEc Mc.Z0 then mn else
242
xtopoly p' (i+1) l in
244
xtopoly (Mc.PEc Mc.Z0) 0 l
246
let rec fixpoint f x =
258
let rec_simpl_cone n_spec e =
260
Mc.simpl_cone n_spec.zero n_spec.unit n_spec.mult n_spec.eqb in
262
let rec rec_simpl_cone = function
263
| Mc.S_Mult(t1, t2) ->
264
simpl_cone (Mc.S_Mult (rec_simpl_cone t1, rec_simpl_cone t2))
266
simpl_cone (Mc.S_Add (rec_simpl_cone t1, rec_simpl_cone t2))
267
| x -> simpl_cone x in
271
let simplify_cone n_spec c = fixpoint (rec_simpl_cone n_spec) c
275
| Ideal of cone *cone
276
| Mult of cone * cone
278
and cone = Mc.zWitness
282
let factorise_linear_cone c =
284
let rec cone_list c l =
286
| Mc.S_Add (x,r) -> cone_list r (x::l)
289
let factorise c1 c2 =
291
| Mc.S_Ideal(x,y) , Mc.S_Ideal(x',y') ->
292
if x = x' then Some (Mc.S_Ideal(x, Mc.S_Add(y,y'))) else None
293
| Mc.S_Mult(x,y) , Mc.S_Mult(x',y') ->
294
if x = x' then Some (Mc.S_Mult(x, Mc.S_Add(y,y'))) else None
297
let rec rebuild_cone l pending =
299
| [] -> (match pending with
305
| None -> rebuild_cone l (Some e)
306
| Some p -> (match factorise p e with
307
| None -> Mc.S_Add(p, rebuild_cone l (Some e))
308
| Some f -> rebuild_cone l (Some f) )
311
(rebuild_cone (List.sort Pervasives.compare (cone_list c [])) None)
315
(* The binding with Fourier might be a bit obsolete
316
-- how does it handle equalities ? *)
318
(* Certificates are elements of the cone such that P = 0 *)
320
(* To begin with, we search for certificates of the form:
321
a1.p1 + ... an.pn + b1.q1 +... + bn.qn + c = 0
326
This is a linear problem: each monomial is considered as a variable.
327
Hence, we can use fourier.
329
The variable c is at index 0
333
(*module Fourier = Fourier(Vector.VList)(SysSet(Vector.VList))*)
334
(*module Fourier = Fourier(Vector.VSparse)(SysSetAlt(Vector.VSparse))*)
335
module Fourier = Mfourier.Fourier(Vector.VSparse)(*(SysSetAlt(Vector.VMap))*)
337
module Vect = Fourier.Vect
340
(* fold_left followed by a rev ! *)
342
let constrain_monomial mn l =
343
let coeffs = List.fold_left (fun acc p -> (Poly.get mn p)::acc) [] l in
344
if mn = Monomial.const
346
{ coeffs = Vect.from_list ((Big_int unit_big_int):: (List.rev coeffs)) ;
348
cst = Big_int zero_big_int }
350
{ coeffs = Vect.from_list ((Big_int zero_big_int):: (List.rev coeffs)) ;
352
cst = Big_int zero_big_int }
356
let rec xpositivity i l =
359
| (_,Mc.Equal)::l -> xpositivity (i+1) l
361
{coeffs = Vect.update (i+1) (fun _ -> Int 1) Vect.null ;
363
cst = Int 0 } :: (xpositivity (i+1) l)
368
let string_of_op = function
370
| Mc.NonStrict -> ">= 0"
372
| Mc.NonEqual -> "<> 0"
376
(* If the certificate includes at least one strict inequality,
377
the obtained polynomial can also be 0 *)
378
let build_linear_system l =
380
(* Gather the monomials: HINT add up of the polynomials *)
381
let l' = List.map fst l in
383
List.fold_left (fun acc p -> Poly.addition p acc) (Poly.constant (Int 0)) l'
384
in (* For each monomial, compute a constraint *)
386
Poly.fold (fun mn _ res -> (constrain_monomial mn l')::res) monomials [] in
387
(* I need at least something strictly positive *)
389
coeffs = Vect.from_list ((Big_int unit_big_int)::
390
(List.map (fun (x,y) ->
391
match y with Mc.Strict ->
393
| _ -> Big_int zero_big_int) l));
394
op = Ge ; cst = Big_int unit_big_int } in
395
(* Add the positivity constraint *)
396
{coeffs = Vect.from_list ([Big_int unit_big_int]) ;
398
cst = Big_int zero_big_int}::(strict::(positivity l)@s0)
401
let big_int_to_z = Ml2C.bigint
403
(* For Q, this is a pity that the certificate has been scaled
404
-- at a lower layer, certificates are using nums... *)
405
let make_certificate n_spec cert li =
406
let bint_to_cst = n_spec.bigint_to_number in
410
let cst = match compare_big_int e zero_big_int with
412
| 1 -> Mc.S_Pos (bint_to_cst e)
413
| _ -> failwith "positivity error"
415
let rec scalar_product cert l =
418
| c::cert -> match l with
419
| [] -> failwith "make_certificate(1)"
421
let r = scalar_product cert l in
422
match compare_big_int c zero_big_int with
424
Mc.S_Ideal (Mc.PEc ( bint_to_cst c), Mc.S_In (Ml2C.nat i)),
428
Mc.S_Mult (Mc.S_Pos (bint_to_cst c), Mc.S_In (Ml2C.nat i)),
431
Some ((factorise_linear_cone
432
(simplify_cone n_spec (Mc.S_Add (cst, scalar_product cert' li)))))
435
exception Found of Monomial.t
437
let raw_certificate l =
438
let sys = build_linear_system l in
440
match Fourier.find_point sys with
442
| Some cert -> Some (rats_to_ints (Vect.to_list cert))
443
(* should not use rats_to_ints *)
446
then (Printf.printf "raw certificate %s" (Printexc.to_string x);
451
let simple_linear_prover to_constant l =
452
let (lc,li) = List.split l in
453
match raw_certificate lc with
454
| None -> None (* No certificate *)
455
| Some cert -> make_certificate to_constant cert li
459
let linear_prover n_spec l =
460
let li = List.combine l (interval 0 (List.length l -1)) in
461
let (l1,l') = List.partition
462
(fun (x,_) -> if snd' x = Mc.NonEqual then true else false) li in
464
(fun (c,i) -> let (Mc.Pair(x,y)) = c in
466
Mc.NonEqual -> failwith "cannot happen"
467
| y -> ((dev_form n_spec x, y),i)) l' in
469
simple_linear_prover n_spec l'
472
let linear_prover n_spec l =
473
try linear_prover n_spec l with
474
x -> (print_string (Printexc.to_string x); None)
478
(* I need to gather the set of variables --->
480
Once I have an interval, I need a certificate : 2 other fourier elims.
481
(I could probably get the certificate directly
482
as it is done in the fourier contrib.)
485
let make_linear_system l =
486
let l' = List.map fst l in
487
let monomials = List.fold_left (fun acc p -> Poly.addition p acc)
488
(Poly.constant (Int 0)) l' in
489
let monomials = Poly.fold
490
(fun mn _ l -> if mn = Monomial.const then l else mn::l) monomials [] in
491
(List.map (fun (c,op) ->
492
{coeffs = Vect.from_list (List.map (fun mn -> (Poly.get mn c)) monomials) ;
494
cst = minus_num ( (Poly.get Monomial.const c))}) l
499
let pplus x y = Mc.PEadd(x,y)
500
let pmult x y = Mc.PEmul(x,y)
501
let pconst x = Mc.PEc x
502
let popp x = Mc.PEopp x
506
(* keep track of enumerated vectors *)
508
match l with [] -> false | e::l -> if p x e then true else mem p x l
510
let rec remove_assoc p x l =
511
match l with [] -> [] | e::l -> if p x (fst e) then
512
remove_assoc p x l else e::(remove_assoc p x l)
514
let eq x y = Vect.compare x y = 0
516
let remove e l = List.fold_left (fun l x -> if eq x e then l else x::l) [] l
519
(* The prover is (probably) incomplete --
520
only searching for naive cutting planes *)
523
let ll = List.fold_right (
524
fun (Mc.Pair(e,k)) r ->
526
| Mc.NonStrict -> (dev_form z_spec e , Ge)::r
527
| Mc.Equal -> (dev_form z_spec e , Eq)::r
528
(* we already know the bound -- don't compute it again *)
529
| _ -> failwith "Cannot happen candidates") sys [] in
531
let (sys,var_mn) = make_linear_system ll in
532
let vars = mapi (fun _ i -> Vect.set i (Int 1) Vect.null) var_mn in
533
(List.fold_left (fun l cstr ->
534
let gcd = Big_int (Vect.gcd cstr.coeffs) in
535
if gcd =/ (Int 1) && cstr.op = Eq
537
else (Vect.mul (Int 1 // gcd) cstr.coeffs)::l) [] sys) @ vars
540
let rec xzlinear_prover planes sys =
541
match linear_prover z_spec sys with
542
| Some prf -> Some (Mc.RatProof prf)
543
| None -> (* find the candidate with the smallest range *)
544
(* Grrr - linear_prover is also calling 'make_linear_system' *)
545
let ll = List.fold_right (fun (Mc.Pair(e,k)) r -> match k with
547
| k -> (dev_form z_spec e ,
551
| Mc.Strict | Mc.NonEqual -> failwith "Cannot happen") :: r) sys [] in
552
let (ll,var) = make_linear_system ll in
553
let candidates = List.fold_left (fun acc vect ->
554
match Fourier.optimise vect ll with
557
(* Printf.printf "%s in %s\n" (Vect.string vect) (string_of_intrvl i) ; *)
559
(vect,i) ::acc) [] planes in
561
let smallest_interval =
562
match List.fold_left (fun (x1,i1) (x2,i2) ->
564
then (x1,i1) else (x2,i2)) (Vect.null,Itv(None,None)) candidates
566
| (x,Itv(Some i, Some j)) -> Some(i,x,j)
567
| (x,Point n) -> Some(n,x,n)
568
| x -> None (* This might be a cutting plane *)
570
match smallest_interval with
573
(Ml2C.bigint (sub_big_int (numerator lb) unit_big_int),
574
Ml2C.bigint (denominator lb)) in
576
(Ml2C.bigint (add_big_int unit_big_int (numerator ub)) ,
577
Ml2C.bigint (denominator ub)) in
578
let expr = list_to_polynomial var (Vect.to_list e) in
580
(*x <= ub -> x > ub *)
582
(Mc.Pair(pplus (pmult (pconst ubd) expr) (popp (pconst ubn)),
583
Mc.NonStrict) :: sys),
584
(* lb <= x -> lb > x *)
586
(Mc.Pair( pplus (popp (pmult (pconst lbd) expr)) (pconst lbn) ,
589
| Some cub , Some clb ->
590
(match zlinear_enum (remove e planes) expr
591
(ceiling_num lb) (floor_num ub) sys
595
Some (Mc.EnumProof(Ml2C.q lb,expr,Ml2C.q ub,clb,cub,prf)))
599
and zlinear_enum planes expr clb cub l =
603
let pexpr = pplus (popp (pconst (Ml2C.bigint (numerator clb)))) expr in
604
let sys' = (Mc.Pair(pexpr, Mc.Equal))::l in
606
match xzlinear_prover planes sys' with
607
| None -> if debug then print_string "zlp?"; None
608
| Some prf -> if debug then print_string "zlp!";
609
match zlinear_enum planes expr (clb +/ (Int 1)) cub l with
611
| Some prfl -> Some (Mc.Cons(prf,prfl))
613
let zlinear_prover sys =
614
let candidates = candidates sys in
615
(* Printf.printf "candidates %d" (List.length candidates) ; *)
616
xzlinear_prover candidates sys
620
let rec scale_term t =
622
| Zero -> unit_big_int , Zero
623
| Const n -> (denominator n) , Const (Big_int (numerator n))
624
| Var n -> unit_big_int , Var n
625
| Inv _ -> failwith "scale_term : not implemented"
626
| Opp t -> let s, t = scale_term t in s, Opp t
627
| Add(t1,t2) -> let s1,y1 = scale_term t1 and s2,y2 = scale_term t2 in
628
let g = gcd_big_int s1 s2 in
629
let s1' = div_big_int s1 g in
630
let s2' = div_big_int s2 g in
631
let e = mult_big_int g (mult_big_int s1' s2') in
632
if (compare_big_int e unit_big_int) = 0
633
then (unit_big_int, Add (y1,y2))
634
else e, Add (Mul(Const (Big_int s2'), y1),
635
Mul (Const (Big_int s1'), y2))
636
| Sub _ -> failwith "scale term: not implemented"
637
| Mul(y,z) -> let s1,y1 = scale_term y and s2,y2 = scale_term z in
638
mult_big_int s1 s2 , Mul (y1, y2)
639
| Pow(t,n) -> let s,t = scale_term t in
640
power_big_int_positive_int s n , Pow(t,n)
641
| _ -> failwith "scale_term : not implemented"
644
let (s,t') = scale_term t in
648
let get_index_of_ith_match f i l =
649
let rec get j res l =
651
| [] -> failwith "bad index"
654
(if j = i then res else get (j+1) (res+1) l )
655
else get j (res+1) l in
659
let rec scale_certificate pos = match pos with
660
| Axiom_eq i -> unit_big_int , Axiom_eq i
661
| Axiom_le i -> unit_big_int , Axiom_le i
662
| Axiom_lt i -> unit_big_int , Axiom_lt i
663
| Monoid l -> unit_big_int , Monoid l
664
| Rational_eq n -> (denominator n) , Rational_eq (Big_int (numerator n))
665
| Rational_le n -> (denominator n) , Rational_le (Big_int (numerator n))
666
| Rational_lt n -> (denominator n) , Rational_lt (Big_int (numerator n))
667
| Square t -> let s,t' = scale_term t in
668
mult_big_int s s , Square t'
669
| Eqmul (t, y) -> let s1,y1 = scale_term t and s2,y2 = scale_certificate y in
670
mult_big_int s1 s2 , Eqmul (y1,y2)
671
| Sum (y, z) -> let s1,y1 = scale_certificate y
672
and s2,y2 = scale_certificate z in
673
let g = gcd_big_int s1 s2 in
674
let s1' = div_big_int s1 g in
675
let s2' = div_big_int s2 g in
676
mult_big_int g (mult_big_int s1' s2'),
677
Sum (Product(Rational_le (Big_int s2'), y1),
678
Product (Rational_le (Big_int s1'), y2))
680
let s1,y1 = scale_certificate y and s2,y2 = scale_certificate z in
681
mult_big_int s1 s2 , Product (y1,y2)
685
let rec term_to_q_expr = function
686
| Const n -> PEc (Ml2C.q n)
687
| Zero -> PEc ( Ml2C.q (Int 0))
688
| Var s -> PEX (Ml2C.index
689
(int_of_string (String.sub s 1 (String.length s - 1))))
690
| Mul(p1,p2) -> PEmul(term_to_q_expr p1, term_to_q_expr p2)
691
| Add(p1,p2) -> PEadd(term_to_q_expr p1, term_to_q_expr p2)
692
| Opp p -> PEopp (term_to_q_expr p)
693
| Pow(t,n) -> PEpow (term_to_q_expr t,Ml2C.n n)
694
| Sub(t1,t2) -> PEsub (term_to_q_expr t1, term_to_q_expr t2)
695
| _ -> failwith "term_to_q_expr: not implemented"
697
let q_cert_of_pos pos =
698
let rec _cert_of_pos = function
699
Axiom_eq i -> Mc.S_In (Ml2C.nat i)
700
| Axiom_le i -> Mc.S_In (Ml2C.nat i)
701
| Axiom_lt i -> Mc.S_In (Ml2C.nat i)
702
| Monoid l -> Mc.S_Monoid (Ml2C.list Ml2C.nat l)
703
| Rational_eq n | Rational_le n | Rational_lt n ->
704
if compare_num n (Int 0) = 0 then Mc.S_Z else
706
| Square t -> Mc.S_Square (term_to_q_expr t)
707
| Eqmul (t, y) -> Mc.S_Ideal(term_to_q_expr t, _cert_of_pos y)
708
| Sum (y, z) -> Mc.S_Add (_cert_of_pos y, _cert_of_pos z)
709
| Product (y, z) -> Mc.S_Mult (_cert_of_pos y, _cert_of_pos z) in
710
simplify_cone q_spec (_cert_of_pos pos)
713
let rec term_to_z_expr = function
714
| Const n -> PEc (Ml2C.bigint (big_int_of_num n))
716
| Var s -> PEX (Ml2C.index
717
(int_of_string (String.sub s 1 (String.length s - 1))))
718
| Mul(p1,p2) -> PEmul(term_to_z_expr p1, term_to_z_expr p2)
719
| Add(p1,p2) -> PEadd(term_to_z_expr p1, term_to_z_expr p2)
720
| Opp p -> PEopp (term_to_z_expr p)
721
| Pow(t,n) -> PEpow (term_to_z_expr t,Ml2C.n n)
722
| Sub(t1,t2) -> PEsub (term_to_z_expr t1, term_to_z_expr t2)
723
| _ -> failwith "term_to_z_expr: not implemented"
725
let z_cert_of_pos pos =
726
let s,pos = (scale_certificate pos) in
727
let rec _cert_of_pos = function
728
Axiom_eq i -> Mc.S_In (Ml2C.nat i)
729
| Axiom_le i -> Mc.S_In (Ml2C.nat i)
730
| Axiom_lt i -> Mc.S_In (Ml2C.nat i)
731
| Monoid l -> Mc.S_Monoid (Ml2C.list Ml2C.nat l)
732
| Rational_eq n | Rational_le n | Rational_lt n ->
733
if compare_num n (Int 0) = 0 then Mc.S_Z else
734
Mc.S_Pos (Ml2C.bigint (big_int_of_num n))
735
| Square t -> Mc.S_Square (term_to_z_expr t)
736
| Eqmul (t, y) -> Mc.S_Ideal(term_to_z_expr t, _cert_of_pos y)
737
| Sum (y, z) -> Mc.S_Add (_cert_of_pos y, _cert_of_pos z)
738
| Product (y, z) -> Mc.S_Mult (_cert_of_pos y, _cert_of_pos z) in
739
simplify_cone z_spec (_cert_of_pos pos)
added
removed
contrib/micromega/certificate.ml
