summaryrefslogtreecommitdiff
path: root/contrib/micromega/certificate.ml
diff options
context:
space:
mode:
Diffstat (limited to 'contrib/micromega/certificate.ml')
-rw-r--r--contrib/micromega/certificate.ml618
1 files changed, 618 insertions, 0 deletions
diff --git a/contrib/micromega/certificate.ml b/contrib/micromega/certificate.ml
new file mode 100644
index 00000000..88e882e6
--- /dev/null
+++ b/contrib/micromega/certificate.ml
@@ -0,0 +1,618 @@
+(************************************************************************)
+(* 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 *)
+(* *)
+(************************************************************************)
+
+(* We take as input a list of polynomials [p1...pn] and return an unfeasibility
+ certificate polynomial. *)
+
+(*open Micromega.Polynomial*)
+open Big_int
+open Num
+
+module Mc = Micromega
+module Ml2C = Mutils.CamlToCoq
+module C2Ml = Mutils.CoqToCaml
+
+let (<+>) = add_num
+let (<->) = minus_num
+let (<*>) = mult_num
+
+type var = Mc.positive
+
+module Monomial :
+sig
+ type t
+ val const : t
+ val var : var -> t
+ val find : var -> t -> int
+ val mult : var -> t -> t
+ val prod : t -> t -> t
+ val compare : t -> t -> int
+ val pp : out_channel -> t -> unit
+ val fold : (var -> int -> 'a -> 'a) -> t -> 'a -> 'a
+end
+ =
+struct
+ (* A monomial is represented by a multiset of variables *)
+ module Map = Map.Make(struct type t = var let compare = Pervasives.compare end)
+ open Map
+
+ type t = int Map.t
+
+ (* The monomial that corresponds to a constant *)
+ let const = Map.empty
+
+ (* The monomial 'x' *)
+ let var x = Map.add x 1 Map.empty
+
+ (* Get the degre of a variable in a monomial *)
+ let find x m = try find x m with Not_found -> 0
+
+ (* Multiply a monomial by a variable *)
+ let mult x m = add x ( (find x m) + 1) m
+
+ (* Product of monomials *)
+ let prod m1 m2 = Map.fold (fun k d m -> add k ((find k m) + d) m) m1 m2
+
+ (* Total ordering of monomials *)
+ let compare m1 m2 = Map.compare Pervasives.compare m1 m2
+
+ let pp o m = Map.iter (fun k v ->
+ if v = 1 then Printf.fprintf o "x%i." (C2Ml.index k)
+ else Printf.fprintf o "x%i^%i." (C2Ml.index k) v) m
+
+ let fold = fold
+
+end
+
+
+module Poly :
+ (* A polynomial is a map of monomials *)
+ (*
+ This is probably a naive implementation
+ (expected to be fast enough - Coq is probably the bottleneck)
+ *The new ring contribution is using a sparse Horner representation.
+ *)
+sig
+ type t
+ val get : Monomial.t -> t -> num
+ val variable : var -> t
+ val add : Monomial.t -> num -> t -> t
+ val constant : num -> t
+ val mult : Monomial.t -> num -> t -> t
+ val product : t -> t -> t
+ val addition : t -> t -> t
+ val uminus : t -> t
+ val fold : (Monomial.t -> num -> 'a -> 'a) -> t -> 'a -> 'a
+ val pp : out_channel -> t -> unit
+ val compare : t -> t -> int
+end =
+struct
+ (*normalisation bug : 0*x ... *)
+ module P = Map.Make(Monomial)
+ open P
+
+ type t = num P.t
+
+ let pp o p = P.iter (fun k v ->
+ if compare_num v (Int 0) <> 0
+ then
+ if Monomial.compare Monomial.const k = 0
+ then Printf.fprintf o "%s " (string_of_num v)
+ else Printf.fprintf o "%s*%a " (string_of_num v) Monomial.pp k) p
+
+ (* Get the coefficient of monomial mn *)
+ let get : Monomial.t -> t -> num =
+ fun mn p -> try find mn p with Not_found -> (Int 0)
+
+
+ (* The polynomial 1.x *)
+ let variable : var -> t =
+ fun x -> add (Monomial.var x) (Int 1) empty
+
+ (*The constant polynomial *)
+ let constant : num -> t =
+ fun c -> add (Monomial.const) c empty
+
+ (* The addition of a monomial *)
+
+ let add : Monomial.t -> num -> t -> t =
+ fun mn v p ->
+ let vl = (get mn p) <+> v in
+ add mn vl p
+
+
+ (** Design choice: empty is not a polynomial
+ I do not remember why ....
+ **)
+
+ (* The product by a monomial *)
+ let mult : Monomial.t -> num -> t -> t =
+ fun mn v p ->
+ fold (fun mn' v' res -> P.add (Monomial.prod mn mn') (v<*>v') res) p empty
+
+
+ let addition : t -> t -> t =
+ fun p1 p2 -> fold (fun mn v p -> add mn v p) p1 p2
+
+
+ let product : t -> t -> t =
+ fun p1 p2 ->
+ fold (fun mn v res -> addition (mult mn v p2) res ) p1 empty
+
+
+ let uminus : t -> t =
+ fun p -> map (fun v -> minus_num v) p
+
+ let fold = P.fold
+
+ let compare = compare compare_num
+end
+
+open Mutils
+type 'a number_spec = {
+ bigint_to_number : big_int -> 'a;
+ number_to_num : 'a -> num;
+ zero : 'a;
+ unit : 'a;
+ mult : 'a -> 'a -> 'a;
+ eqb : 'a -> 'a -> Mc.bool
+}
+
+let z_spec = {
+ bigint_to_number = Ml2C.bigint ;
+ number_to_num = (fun x -> Big_int (C2Ml.z_big_int x));
+ zero = Mc.Z0;
+ unit = Mc.Zpos Mc.XH;
+ mult = Mc.zmult;
+ eqb = Mc.zeq_bool
+}
+
+
+let q_spec = {
+ bigint_to_number = (fun x -> {Mc.qnum = Ml2C.bigint x; Mc.qden = Mc.XH});
+ number_to_num = C2Ml.q_to_num;
+ zero = {Mc.qnum = Mc.Z0;Mc.qden = Mc.XH};
+ unit = {Mc.qnum = (Mc.Zpos Mc.XH) ; Mc.qden = Mc.XH};
+ mult = Mc.qmult;
+ eqb = Mc.qeq_bool
+}
+
+let r_spec = z_spec
+
+
+
+
+let dev_form n_spec p =
+ let rec dev_form p =
+ match p with
+ | Mc.PEc z -> Poly.constant (n_spec.number_to_num z)
+ | Mc.PEX v -> Poly.variable v
+ | Mc.PEmul(p1,p2) ->
+ let p1 = dev_form p1 in
+ let p2 = dev_form p2 in
+ Poly.product p1 p2
+ | Mc.PEadd(p1,p2) -> Poly.addition (dev_form p1) (dev_form p2)
+ | Mc.PEopp p -> Poly.uminus (dev_form p)
+ | Mc.PEsub(p1,p2) -> Poly.addition (dev_form p1) (Poly.uminus (dev_form p2))
+ | Mc.PEpow(p,n) ->
+ let p = dev_form p in
+ let n = C2Ml.n n in
+ let rec pow n =
+ if n = 0
+ then Poly.constant (n_spec.number_to_num n_spec.unit)
+ else Poly.product p (pow (n-1)) in
+ pow n in
+ dev_form p
+
+
+let monomial_to_polynomial mn =
+ Monomial.fold
+ (fun v i acc ->
+ let mn = if i = 1 then Mc.PEX v else Mc.PEpow (Mc.PEX v ,Ml2C.n i) in
+ if acc = Mc.PEc (Mc.Zpos Mc.XH)
+ then mn
+ else Mc.PEmul(mn,acc))
+ mn
+ (Mc.PEc (Mc.Zpos Mc.XH))
+
+let list_to_polynomial vars l =
+ assert (List.for_all (fun x -> ceiling_num x =/ x) l);
+ let var x = monomial_to_polynomial (List.nth vars x) in
+ let rec xtopoly p i = function
+ | [] -> p
+ | c::l -> if c =/ (Int 0) then xtopoly p (i+1) l
+ else let c = Mc.PEc (Ml2C.bigint (numerator c)) in
+ let mn =
+ if c = Mc.PEc (Mc.Zpos Mc.XH)
+ then var i
+ else Mc.PEmul (c,var i) in
+ let p' = if p = Mc.PEc Mc.Z0 then mn else
+ Mc.PEadd (mn, p) in
+ xtopoly p' (i+1) l in
+
+ xtopoly (Mc.PEc Mc.Z0) 0 l
+
+let rec fixpoint f x =
+ let y' = f x in
+ if y' = x then y'
+ else fixpoint f y'
+
+
+
+
+
+
+
+
+let rec_simpl_cone n_spec e =
+ let simpl_cone =
+ Mc.simpl_cone n_spec.zero n_spec.unit n_spec.mult n_spec.eqb in
+
+ let rec rec_simpl_cone = function
+ | Mc.S_Mult(t1, t2) ->
+ simpl_cone (Mc.S_Mult (rec_simpl_cone t1, rec_simpl_cone t2))
+ | Mc.S_Add(t1,t2) ->
+ simpl_cone (Mc.S_Add (rec_simpl_cone t1, rec_simpl_cone t2))
+ | x -> simpl_cone x in
+ rec_simpl_cone e
+
+
+let simplify_cone n_spec c = fixpoint (rec_simpl_cone n_spec) c
+
+type cone_prod =
+ Const of cone
+ | Ideal of cone *cone
+ | Mult of cone * cone
+ | Other of cone
+and cone = Mc.zWitness
+
+
+
+let factorise_linear_cone c =
+
+ let rec cone_list c l =
+ match c with
+ | Mc.S_Add (x,r) -> cone_list r (x::l)
+ | _ -> c :: l in
+
+ let factorise c1 c2 =
+ match c1 , c2 with
+ | Mc.S_Ideal(x,y) , Mc.S_Ideal(x',y') ->
+ if x = x' then Some (Mc.S_Ideal(x, Mc.S_Add(y,y'))) else None
+ | Mc.S_Mult(x,y) , Mc.S_Mult(x',y') ->
+ if x = x' then Some (Mc.S_Mult(x, Mc.S_Add(y,y'))) else None
+ | _ -> None in
+
+ let rec rebuild_cone l pending =
+ match l with
+ | [] -> (match pending with
+ | None -> Mc.S_Z
+ | Some p -> p
+ )
+ | e::l ->
+ (match pending with
+ | None -> rebuild_cone l (Some e)
+ | Some p -> (match factorise p e with
+ | None -> Mc.S_Add(p, rebuild_cone l (Some e))
+ | Some f -> rebuild_cone l (Some f) )
+ ) in
+
+ (rebuild_cone (List.sort Pervasives.compare (cone_list c [])) None)
+
+
+
+(* The binding with Fourier might be a bit obsolete
+ -- how does it handle equalities ? *)
+
+(* Certificates are elements of the cone such that P = 0 *)
+
+(* To begin with, we search for certificates of the form:
+ a1.p1 + ... an.pn + b1.q1 +... + bn.qn + c = 0
+ where pi >= 0 qi > 0
+ ai >= 0
+ bi >= 0
+ Sum bi + c >= 1
+ This is a linear problem: each monomial is considered as a variable.
+ Hence, we can use fourier.
+
+ The variable c is at index 0
+*)
+
+open Mfourier
+ (*module Fourier = Fourier(Vector.VList)(SysSet(Vector.VList))*)
+ (*module Fourier = Fourier(Vector.VSparse)(SysSetAlt(Vector.VSparse))*)
+module Fourier = Mfourier.Fourier(Vector.VSparse)(*(SysSetAlt(Vector.VMap))*)
+
+module Vect = Fourier.Vect
+open Fourier.Cstr
+
+(* fold_left followed by a rev ! *)
+
+let constrain_monomial mn l =
+ let coeffs = List.fold_left (fun acc p -> (Poly.get mn p)::acc) [] l in
+ if mn = Monomial.const
+ then
+ { coeffs = Vect.from_list ((Big_int unit_big_int):: (List.rev coeffs)) ;
+ op = Eq ;
+ cst = Big_int zero_big_int }
+ else
+ { coeffs = Vect.from_list ((Big_int zero_big_int):: (List.rev coeffs)) ;
+ op = Eq ;
+ cst = Big_int zero_big_int }
+
+
+let positivity l =
+ let rec xpositivity i l =
+ match l with
+ | [] -> []
+ | (_,Mc.Equal)::l -> xpositivity (i+1) l
+ | (_,_)::l ->
+ {coeffs = Vect.update (i+1) (fun _ -> Int 1) Vect.null ;
+ op = Ge ;
+ cst = Int 0 } :: (xpositivity (i+1) l)
+ in
+ xpositivity 0 l
+
+
+let string_of_op = function
+ | Mc.Strict -> "> 0"
+ | Mc.NonStrict -> ">= 0"
+ | Mc.Equal -> "= 0"
+ | Mc.NonEqual -> "<> 0"
+
+
+
+(* If the certificate includes at least one strict inequality,
+ the obtained polynomial can also be 0 *)
+let build_linear_system l =
+
+ (* Gather the monomials: HINT add up of the polynomials *)
+ let l' = List.map fst l in
+ let monomials =
+ List.fold_left (fun acc p -> Poly.addition p acc) (Poly.constant (Int 0)) l'
+ in (* For each monomial, compute a constraint *)
+ let s0 =
+ Poly.fold (fun mn _ res -> (constrain_monomial mn l')::res) monomials [] in
+ (* I need at least something strictly positive *)
+ let strict = {
+ coeffs = Vect.from_list ((Big_int unit_big_int)::
+ (List.map (fun (x,y) ->
+ match y with Mc.Strict ->
+ Big_int unit_big_int
+ | _ -> Big_int zero_big_int) l));
+ op = Ge ; cst = Big_int unit_big_int } in
+ (* Add the positivity constraint *)
+ {coeffs = Vect.from_list ([Big_int unit_big_int]) ;
+ op = Ge ;
+ cst = Big_int zero_big_int}::(strict::(positivity l)@s0)
+
+
+let big_int_to_z = Ml2C.bigint
+
+(* For Q, this is a pity that the certificate has been scaled
+ -- at a lower layer, certificates are using nums... *)
+let make_certificate n_spec cert li =
+ let bint_to_cst = n_spec.bigint_to_number in
+ match cert with
+ | [] -> None
+ | e::cert' ->
+ let cst = match compare_big_int e zero_big_int with
+ | 0 -> Mc.S_Z
+ | 1 -> Mc.S_Pos (bint_to_cst e)
+ | _ -> failwith "positivity error"
+ in
+ let rec scalar_product cert l =
+ match cert with
+ | [] -> Mc.S_Z
+ | c::cert -> match l with
+ | [] -> failwith "make_certificate(1)"
+ | i::l ->
+ let r = scalar_product cert l in
+ match compare_big_int c zero_big_int with
+ | -1 -> Mc.S_Add (
+ Mc.S_Ideal (Mc.PEc ( bint_to_cst c), Mc.S_In (Ml2C.nat i)),
+ r)
+ | 0 -> r
+ | _ -> Mc.S_Add (
+ Mc.S_Mult (Mc.S_Pos (bint_to_cst c), Mc.S_In (Ml2C.nat i)),
+ r) in
+
+ Some ((factorise_linear_cone
+ (simplify_cone n_spec (Mc.S_Add (cst, scalar_product cert' li)))))
+
+
+exception Found of Monomial.t
+
+let raw_certificate l =
+ let sys = build_linear_system l in
+ try
+ match Fourier.find_point sys with
+ | None -> None
+ | Some cert -> Some (rats_to_ints (Vect.to_list cert))
+ (* should not use rats_to_ints *)
+ with x ->
+ if debug
+ then (Printf.printf "raw certificate %s" (Printexc.to_string x);
+ flush stdout) ;
+ None
+
+
+let simple_linear_prover to_constant l =
+ let (lc,li) = List.split l in
+ match raw_certificate lc with
+ | None -> None (* No certificate *)
+ | Some cert -> make_certificate to_constant cert li
+
+
+
+let linear_prover n_spec l =
+ let li = List.combine l (interval 0 (List.length l -1)) in
+ let (l1,l') = List.partition
+ (fun (x,_) -> if snd' x = Mc.NonEqual then true else false) li in
+ let l' = List.map
+ (fun (c,i) -> let (Mc.Pair(x,y)) = c in
+ match y with
+ Mc.NonEqual -> failwith "cannot happen"
+ | y -> ((dev_form n_spec x, y),i)) l' in
+
+ simple_linear_prover n_spec l'
+
+
+let linear_prover n_spec l =
+ try linear_prover n_spec l with
+ x -> (print_string (Printexc.to_string x); None)
+
+(* zprover.... *)
+
+(* I need to gather the set of variables --->
+ Then go for fold
+ Once I have an interval, I need a certificate : 2 other fourier elims.
+ (I could probably get the certificate directly
+ as it is done in the fourier contrib.)
+*)
+
+let make_linear_system l =
+ let l' = List.map fst l in
+ let monomials = List.fold_left (fun acc p -> Poly.addition p acc)
+ (Poly.constant (Int 0)) l' in
+ let monomials = Poly.fold
+ (fun mn _ l -> if mn = Monomial.const then l else mn::l) monomials [] in
+ (List.map (fun (c,op) ->
+ {coeffs = Vect.from_list (List.map (fun mn -> (Poly.get mn c)) monomials) ;
+ op = op ;
+ cst = minus_num ( (Poly.get Monomial.const c))}) l
+ ,monomials)
+
+
+open Interval
+let pplus x y = Mc.PEadd(x,y)
+let pmult x y = Mc.PEmul(x,y)
+let pconst x = Mc.PEc x
+let popp x = Mc.PEopp x
+
+let debug = false
+
+(* keep track of enumerated vectors *)
+let rec mem p x l =
+ match l with [] -> false | e::l -> if p x e then true else mem p x l
+
+let rec remove_assoc p x l =
+ match l with [] -> [] | e::l -> if p x (fst e) then
+ remove_assoc p x l else e::(remove_assoc p x l)
+
+let eq x y = Vect.compare x y = 0
+
+(* Beurk... this code is a shame *)
+
+let rec zlinear_prover sys = xzlinear_prover [] sys
+
+and xzlinear_prover enum l : (Mc.proofTerm option) =
+ match linear_prover z_spec l with
+ | Some prf -> Some (Mc.RatProof prf)
+ | None ->
+ let ll = List.fold_right (fun (Mc.Pair(e,k)) r -> match k with
+ Mc.NonEqual -> r
+ | k -> (dev_form z_spec e ,
+ match k with
+ | Mc.Strict | Mc.NonStrict -> Ge
+ (* Loss of precision -- weakness of fourier*)
+ | Mc.Equal -> Eq
+ | Mc.NonEqual -> failwith "Cannot happen") :: r) l [] in
+
+ let (sys,var) = make_linear_system ll in
+ let res =
+ match Fourier.find_Q_interval sys with
+ | Some(i,x,j) -> if i =/ j
+ then Some(i,Vect.set x (Int 1) Vect.null,i) else None
+ | None -> None in
+ let res = match res with
+ | None ->
+ begin
+ let candidates = List.fold_right
+ (fun cstr acc ->
+ let gcd = Big_int (Vect.gcd cstr.coeffs) in
+ let vect = Vect.mul (Int 1 // gcd) cstr.coeffs in
+ if mem eq vect enum then acc
+ else ((vect,Fourier.optimise vect sys)::acc)) sys [] in
+ let candidates = List.fold_left (fun l (x,i) ->
+ match i with
+ None -> (x,Empty)::l
+ | Some i -> (x,i)::l) [] (candidates) in
+ match List.fold_left (fun (x1,i1) (x2,i2) ->
+ if smaller_itv i1 i2
+ then (x1,i1) else (x2,i2)) (Vect.null,Itv(None,None)) candidates
+ with
+ | (i,Empty) -> None
+ | (x,Itv(Some i, Some j)) -> Some(i,x,j)
+ | (x,Point n) -> Some(n,x,n)
+ | x -> match Fourier.find_Q_interval sys with
+ | None -> None
+ | Some(i,x,j) ->
+ if i =/ j
+ then Some(i,Vect.set x (Int 1) Vect.null,i)
+ else None
+ end
+ | _ -> res in
+
+ match res with
+ | Some (lb,e,ub) ->
+ let (lbn,lbd) =
+ (Ml2C.bigint (sub_big_int (numerator lb) unit_big_int),
+ Ml2C.bigint (denominator lb)) in
+ let (ubn,ubd) =
+ (Ml2C.bigint (add_big_int unit_big_int (numerator ub)) ,
+ Ml2C.bigint (denominator ub)) in
+ let expr = list_to_polynomial var (Vect.to_list e) in
+ (match
+ (*x <= ub -> x > ub *)
+ linear_prover z_spec
+ (Mc.Pair(pplus (pmult (pconst ubd) expr) (popp (pconst ubn)),
+ Mc.NonStrict) :: l),
+ (* lb <= x -> lb > x *)
+ linear_prover z_spec
+ (Mc.Pair( pplus (popp (pmult (pconst lbd) expr)) (pconst lbn) ,
+ Mc.NonStrict)::l)
+ with
+ | Some cub , Some clb ->
+ (match zlinear_enum (e::enum) expr
+ (ceiling_num lb) (floor_num ub) l
+ with
+ | None -> None
+ | Some prf ->
+ Some (Mc.EnumProof(Ml2C.q lb,expr,Ml2C.q ub,clb,cub,prf)))
+ | _ -> None
+ )
+ | _ -> None
+and xzlinear_enum enum expr clb cub l =
+ if clb >/ cub
+ then Some Mc.Nil
+ else
+ let pexpr = pplus (popp (pconst (Ml2C.bigint (numerator clb)))) expr in
+ let sys' = (Mc.Pair(pexpr, Mc.Equal))::l in
+ match xzlinear_prover enum sys' with
+ | None -> if debug then print_string "zlp?"; None
+ | Some prf -> if debug then print_string "zlp!";
+ match zlinear_enum enum expr (clb +/ (Int 1)) cub l with
+ | None -> None
+ | Some prfl -> Some (Mc.Cons(prf,prfl))
+
+and zlinear_enum enum expr clb cub l =
+ let res = xzlinear_enum enum expr clb cub l in
+ if debug then Printf.printf "zlinear_enum %s %s -> %s\n"
+ (string_of_num clb)
+ (string_of_num cub)
+ (match res with
+ | None -> "None"
+ | Some r -> "Some") ; res
+
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-18 19:50:26 +0000