diff options
| author | Emilio Jesus Gallego Arias | 2020-03-04 14:29:10 -0500 |
|---|---|---|
| committer | Emilio Jesus Gallego Arias | 2020-03-04 21:17:46 -0500 |
| commit | 8af9dbdcc27934deda35f6c8fbdecdfe869b09c5 (patch) | |
| tree | e0e4f6ece4b2bbfc01b7662d43519ff49a27143a /plugins/micromega | |
| parent | 33ab70ac3a8d08afb67d9602d7c23da7133ad0f4 (diff) | |
[micromega] Add numerical compatibility layer.
Only significant change is in gcd / lcm which now are typed in `Z.t`
Diffstat (limited to 'plugins/micromega')
25 files changed, 853 insertions, 667 deletions
diff --git a/plugins/micromega/certificate.ml b/plugins/micromega/certificate.ml index e946ffd8bc..c788c7f147 100644 --- a/plugins/micromega/certificate.ml +++ b/plugins/micromega/certificate.ml @@ -19,12 +19,13 @@ let debug = false -open Big_int -open Num open Polynomial module Mc = Micromega module Ml2C = Mutils.CamlToCoq module C2Ml = Mutils.CoqToCaml +open NumCompat +open Q.Notations +open Mutils let use_simplex = ref true @@ -32,11 +33,9 @@ type ('prf, 'model) res = Prf of 'prf | Model of 'model | Unknown type zres = (Mc.zArithProof, int * Mc.z list) res type qres = (Mc.q Mc.psatz, int * Mc.q list) res -open Mutils - type 'a number_spec = - { bigint_to_number : big_int -> 'a - ; number_to_num : 'a -> num + { bigint_to_number : Z.t -> 'a + ; number_to_num : 'a -> Q.t ; zero : 'a ; unit : 'a ; mult : 'a -> 'a -> 'a @@ -44,7 +43,7 @@ type 'a number_spec = let z_spec = { bigint_to_number = Ml2C.bigint - ; number_to_num = (fun x -> Big_int (C2Ml.z_big_int x)) + ; number_to_num = (fun x -> Q.of_bigint (C2Ml.z_big_int x)) ; zero = Mc.Z0 ; unit = Mc.Zpos Mc.XH ; mult = Mc.Z.mul @@ -124,17 +123,16 @@ let constrain_variable v l = let coeffs = List.fold_left (fun acc p -> Vect.get v p.coeffs :: acc) [] l in { coeffs = Vect.from_list - (Big_int zero_big_int :: Big_int zero_big_int :: List.rev coeffs) + (Q.of_bigint Z.zero :: Q.of_bigint Z.zero :: List.rev coeffs) ; op = Eq - ; cst = Big_int zero_big_int } + ; cst = Q.of_bigint Z.zero } let constrain_constant l = - let coeffs = List.fold_left (fun acc p -> minus_num p.cst :: acc) [] l in + let coeffs = List.fold_left (fun acc p -> Q.neg p.cst :: acc) [] l in { coeffs = - Vect.from_list - (Big_int zero_big_int :: Big_int unit_big_int :: List.rev coeffs) + Vect.from_list (Q.of_bigint Z.zero :: Q.of_bigint Z.one :: List.rev coeffs) ; op = Eq - ; cst = Big_int zero_big_int } + ; cst = Q.of_bigint Z.zero } let positivity l = let rec xpositivity i l = @@ -144,16 +142,16 @@ let positivity l = match c.op with | Eq -> xpositivity (i + 1) l | _ -> - { coeffs = Vect.update (i + 1) (fun _ -> Int 1) Vect.null + { coeffs = Vect.update (i + 1) (fun _ -> Q.one) Vect.null ; op = Ge - ; cst = Int 0 } + ; cst = Q.zero } :: xpositivity (i + 1) l ) in xpositivity 1 l let cstr_of_poly (p, o) = let c, l = Vect.decomp_cst p in - {coeffs = l; op = o; cst = minus_num c} + {coeffs = l; op = o; cst = Q.neg c} let variables_of_cstr c = Vect.variables c.coeffs @@ -175,25 +173,23 @@ let build_dual_linear_system l = let strict = { coeffs = Vect.from_list - ( Big_int zero_big_int :: Big_int unit_big_int + ( Q.of_bigint Z.zero :: Q.of_bigint Z.one :: List.map (fun c -> - if is_strict c then Big_int unit_big_int - else Big_int zero_big_int) + if is_strict c then Q.of_bigint Z.one else Q.of_bigint Z.zero) l ) ; op = Ge - ; cst = Big_int unit_big_int } + ; cst = Q.of_bigint Z.one } in (* Add the positivity constraint *) - { coeffs = Vect.from_list [Big_int zero_big_int; Big_int unit_big_int] + { coeffs = Vect.from_list [Q.of_bigint Z.zero; Q.of_bigint Z.one] ; op = Ge - ; cst = Big_int zero_big_int } + ; cst = Q.of_bigint Z.zero } :: ((strict :: positivity l) @ (c :: s0)) -open Util - (** [direct_linear_prover l] does not handle strict inegalities *) let fourier_linear_prover l = + let open Util in match Mfourier.Fourier.find_point l with | Inr prf -> if debug then Printf.printf "AProof : %a\n" Mfourier.pp_proof prf; @@ -211,6 +207,7 @@ let direct_linear_prover l = else fourier_linear_prover l let find_point l = + let open Util in if !use_simplex then Simplex.find_point l else match Mfourier.Fourier.find_point l with @@ -237,8 +234,8 @@ let dual_raw_certificate l = match Vect.choose cert with | None -> failwith "dual_raw_certificate: empty_certificate" | Some _ -> - (*Some (rats_to_ints (Vect.to_list (Vect.decr_var 2 (Vect.set 1 (Int 0) cert))))*) - Some (Vect.normalise (Vect.decr_var 2 (Vect.set 1 (Int 0) cert))) ) + (*Some (rats_to_ints (Vect.to_list (Vect.decr_var 2 (Vect.set 1 Q.zero cert))))*) + Some (Vect.normalise (Vect.decr_var 2 (Vect.set 1 Q.zero cert))) ) (* should not use rats_to_ints *) with x when CErrors.noncritical x -> if debug then ( @@ -306,14 +303,14 @@ exception FoundProof of ProofFormat.prf_rule let check_int_sat (cstr, prf) = let {coeffs; op; cst} = cstr in match Vect.choose coeffs with - | None -> if eval_op op (Int 0) cst then Tauto else Unsat prf + | None -> if eval_op op Q.zero cst then Tauto else Unsat prf | _ -> ( let gcdi = Vect.gcd coeffs in - let gcd = Big_int gcdi in - if eq_num gcd (Int 1) then Normalise (cstr, prf) - else if Int.equal (sign_num (mod_num cst gcd)) 0 then begin + let gcd = Q.of_bigint gcdi in + if gcd =/ Q.one then Normalise (cstr, prf) + else if Int.equal (Q.sign (Q.mod_ cst gcd)) 0 then begin (* We can really normalise *) - assert (sign_num gcd >= 1); + assert (Q.sign gcd >= 1); let cstr = {coeffs = Vect.div gcd coeffs; op; cst = cst // gcd} in Normalise (cstr, ProofFormat.Gcd (gcdi, prf)) (* Normalise(cstr,CutPrf prf)*) @@ -323,7 +320,7 @@ let check_int_sat (cstr, prf) = | Eq -> Unsat (ProofFormat.CutPrf prf) | Ge -> let cstr = - {coeffs = Vect.div gcd coeffs; op; cst = ceiling_num (cst // gcd)} + {coeffs = Vect.div gcd coeffs; op; cst = Q.ceiling (cst // gcd)} in Cut (cstr, ProofFormat.CutPrf prf) | Gt -> failwith "check_sat : Unexpected operator" ) @@ -351,7 +348,7 @@ let is_linear_for v pc = *) let is_linear_substitution sys ((p, o), prf) = - let pred v = v =/ Int 1 || v =/ Int (-1) in + let pred v = v =/ Q.one || v =/ Q.neg_one in match o with | Eq -> ( match @@ -521,28 +518,31 @@ open Sos_types let rec scale_term t = match t with - | Zero -> (unit_big_int, Zero) - | Const n -> (denominator n, Const (Big_int (numerator n))) - | Var n -> (unit_big_int, Var n) + | Zero -> (Z.one, Zero) + | Const n -> (Q.den n, Const (Q.of_bigint (Q.num n))) + | Var n -> (Z.one, Var n) | Opp t -> let s, t = scale_term t in (s, Opp t) | Add (t1, t2) -> let s1, y1 = scale_term t1 and s2, y2 = scale_term t2 in - let g = gcd_big_int s1 s2 in - let s1' = div_big_int s1 g in - let s2' = div_big_int s2 g in - let e = mult_big_int g (mult_big_int s1' s2') in - if Int.equal (compare_big_int e unit_big_int) 0 then - (unit_big_int, Add (y1, y2)) - else (e, Add (Mul (Const (Big_int s2'), y1), Mul (Const (Big_int s1'), y2))) + let g = Z.gcd s1 s2 in + let s1' = Z.div s1 g in + let s2' = Z.div s2 g in + let e = Z.mul g (Z.mul s1' s2') in + if Int.equal (Z.compare e Z.one) 0 then (Z.one, Add (y1, y2)) + else + ( e + , Add + (Mul (Const (Q.of_bigint s2'), y1), Mul (Const (Q.of_bigint s1'), y2)) + ) | Sub _ -> failwith "scale term: not implemented" | Mul (y, z) -> let s1, y1 = scale_term y and s2, y2 = scale_term z in - (mult_big_int s1 s2, Mul (y1, y2)) + (Z.mul s1 s2, Mul (y1, y2)) | Pow (t, n) -> let s, t = scale_term t in - (power_big_int_positive_int s n, Pow (t, n)) + (Z.power_int s n, Pow (t, n)) let scale_term t = let s, t' = scale_term t in @@ -550,37 +550,38 @@ let scale_term t = let rec scale_certificate pos = match pos with - | Axiom_eq i -> (unit_big_int, Axiom_eq i) - | Axiom_le i -> (unit_big_int, Axiom_le i) - | Axiom_lt i -> (unit_big_int, Axiom_lt i) - | Monoid l -> (unit_big_int, Monoid l) - | Rational_eq n -> (denominator n, Rational_eq (Big_int (numerator n))) - | Rational_le n -> (denominator n, Rational_le (Big_int (numerator n))) - | Rational_lt n -> (denominator n, Rational_lt (Big_int (numerator n))) + | Axiom_eq i -> (Z.one, Axiom_eq i) + | Axiom_le i -> (Z.one, Axiom_le i) + | Axiom_lt i -> (Z.one, Axiom_lt i) + | Monoid l -> (Z.one, Monoid l) + | Rational_eq n -> (Q.den n, Rational_eq (Q.of_bigint (Q.num n))) + | Rational_le n -> (Q.den n, Rational_le (Q.of_bigint (Q.num n))) + | Rational_lt n -> (Q.den n, Rational_lt (Q.of_bigint (Q.num n))) | Square t -> let s, t' = scale_term t in - (mult_big_int s s, Square t') + (Z.mul s s, Square t') | Eqmul (t, y) -> let s1, y1 = scale_term t and s2, y2 = scale_certificate y in - (mult_big_int s1 s2, Eqmul (y1, y2)) + (Z.mul s1 s2, Eqmul (y1, y2)) | Sum (y, z) -> let s1, y1 = scale_certificate y and s2, y2 = scale_certificate z in - let g = gcd_big_int s1 s2 in - let s1' = div_big_int s1 g in - let s2' = div_big_int s2 g in - ( mult_big_int g (mult_big_int s1' s2') + let g = Z.gcd s1 s2 in + let s1' = Z.div s1 g in + let s2' = Z.div s2 g in + ( Z.mul g (Z.mul s1' s2') , Sum - ( Product (Rational_le (Big_int s2'), y1) - , Product (Rational_le (Big_int s1'), y2) ) ) + ( Product (Rational_le (Q.of_bigint s2'), y1) + , Product (Rational_le (Q.of_bigint s1'), y2) ) ) | Product (y, z) -> let s1, y1 = scale_certificate y and s2, y2 = scale_certificate z in - (mult_big_int s1 s2, Product (y1, y2)) + (Z.mul s1 s2, Product (y1, y2)) +module Z_ = Z open Micromega let rec term_to_q_expr = function | Const n -> PEc (Ml2C.q n) - | Zero -> PEc (Ml2C.q (Int 0)) + | Zero -> PEc (Ml2C.q Q.zero) | Var s -> PEX (Ml2C.index (int_of_string (String.sub s 1 (String.length s - 1)))) | Mul (p1, p2) -> PEmul (term_to_q_expr p1, term_to_q_expr p2) @@ -590,8 +591,8 @@ let rec term_to_q_expr = function | Sub (t1, t2) -> PEsub (term_to_q_expr t1, term_to_q_expr t2) let term_to_q_pol e = - Mc.norm_aux (Ml2C.q (Int 0)) (Ml2C.q (Int 1)) Mc.qplus Mc.qmult Mc.qminus - Mc.qopp Mc.qeq_bool (term_to_q_expr e) + Mc.norm_aux (Ml2C.q Q.zero) (Ml2C.q Q.one) Mc.qplus Mc.qmult Mc.qminus Mc.qopp + Mc.qeq_bool (term_to_q_expr e) let rec product l = match l with @@ -606,7 +607,7 @@ let q_cert_of_pos pos = | Axiom_lt i -> Mc.PsatzIn (Ml2C.nat i) | Monoid l -> product l | Rational_eq n | Rational_le n | Rational_lt n -> - if Int.equal (compare_num n (Int 0)) 0 then Mc.PsatzZ + if Int.equal (Q.compare n Q.zero) 0 then Mc.PsatzZ else Mc.PsatzC (Ml2C.q n) | Square t -> Mc.PsatzSquare (term_to_q_pol t) | Eqmul (t, y) -> Mc.PsatzMulC (term_to_q_pol t, _cert_of_pos y) @@ -616,7 +617,7 @@ let q_cert_of_pos pos = simplify_cone q_spec (_cert_of_pos pos) let rec term_to_z_expr = function - | Const n -> PEc (Ml2C.bigint (big_int_of_num n)) + | Const n -> PEc (Ml2C.bigint (Q.to_bigint n)) | Zero -> PEc Z0 | Var s -> PEX (Ml2C.index (int_of_string (String.sub s 1 (String.length s - 1)))) @@ -638,11 +639,11 @@ let z_cert_of_pos pos = | Axiom_lt i -> Mc.PsatzIn (Ml2C.nat i) | Monoid l -> product l | Rational_eq n | Rational_le n | Rational_lt n -> - if Int.equal (compare_num n (Int 0)) 0 then Mc.PsatzZ - else Mc.PsatzC (Ml2C.bigint (big_int_of_num n)) + if Int.equal (Q.compare n Q.zero) 0 then Mc.PsatzZ + else Mc.PsatzC (Ml2C.bigint (Q.to_bigint n)) | Square t -> Mc.PsatzSquare (term_to_z_pol t) | Eqmul (t, y) -> - let is_unit = match t with Const n -> n =/ Int 1 | _ -> false in + let is_unit = match t with Const n -> n =/ Q.one | _ -> false in if is_unit then _cert_of_pos y else Mc.PsatzMulC (term_to_z_pol t, _cert_of_pos y) | Sum (y, z) -> Mc.PsatzAdd (_cert_of_pos y, _cert_of_pos z) @@ -655,8 +656,6 @@ open Mutils Given a constraint, all the coefficients are always integers. *) -open Num -open Big_int open Polynomial type prf_sys = (cstr * ProofFormat.prf_rule) list @@ -674,19 +673,18 @@ let pivot v (c1, p1) (c2, p2) = (ProofFormat.mul_cst_proof cv1 p1) (ProofFormat.mul_cst_proof cv2 p2) ) in - match (Vect.get v v1, Vect.get v v2) with - | Int 0, _ | _, Int 0 -> None - | a, b -> - if Int.equal (sign_num a * sign_num b) (-1) then - let cv1 = abs_num b and cv2 = abs_num a in - Some (xpivot cv1 cv2) - else if op1 == Eq then - let cv1 = minus_num (b */ Int (sign_num a)) and cv2 = abs_num a in - Some (xpivot cv1 cv2) - else if op2 == Eq then - let cv1 = abs_num b and cv2 = minus_num (a */ Int (sign_num b)) in - Some (xpivot cv1 cv2) - else None + let a, b = (Vect.get v v1, Vect.get v v2) in + if a =/ Q.zero || b =/ Q.zero then None + else if Int.equal (Q.sign a * Q.sign b) (-1) then + let cv1 = Q.abs b and cv2 = Q.abs a in + Some (xpivot cv1 cv2) + else if op1 == Eq then + let cv1 = Q.neg (b */ Q.of_int (Q.sign a)) and cv2 = Q.abs a in + Some (xpivot cv1 cv2) + else if op2 == Eq then + let cv1 = Q.abs b and cv2 = Q.neg (a */ Q.of_int (Q.sign b)) in + Some (xpivot cv1 cv2) + else None (* op2 could be Eq ... this might happen *) @@ -705,21 +703,17 @@ let simpl_sys sys = Source: http://en.wikipedia.org/wiki/Extended_Euclidean_algorithm *) let rec ext_gcd a b = - if Int.equal (sign_big_int b) 0 then (unit_big_int, zero_big_int) + if Int.equal (Z_.sign b) 0 then (Z_.one, Z_.zero) else - let q, r = quomod_big_int a b in + let q, r = Z_.quomod a b in let s, t = ext_gcd b r in - (t, sub_big_int s (mult_big_int q t)) + (t, Z_.sub s (Z_.mul q t)) let extract_coprime (c1, p1) (c2, p2) = if c1.op == Eq && c2.op == Eq then Vect.exists2 (fun n1 n2 -> - Int.equal - (compare_big_int - (gcd_big_int (numerator n1) (numerator n2)) - unit_big_int) - 0) + Int.equal (Z_.compare (Z_.gcd (Q.num n1) (Q.num n2)) Z_.one) 0) c1.coeffs c2.coeffs else None @@ -742,8 +736,8 @@ let reduce_coprime psys = match oeq with | None -> None (* Nothing to do *) | Some ((v, n1, n2), (c1, p1), (c2, p2)) -> - let l1, l2 = ext_gcd (numerator n1) (numerator n2) in - let l1' = Big_int l1 and l2' = Big_int l2 in + let l1, l2 = ext_gcd (Q.num n1) (Q.num n2) in + let l1' = Q.of_bigint l1 and l2' = Q.of_bigint l2 in let cstr = { coeffs = Vect.add (Vect.mul l1' c1.coeffs) (Vect.mul l2' c2.coeffs) ; op = Eq @@ -761,7 +755,7 @@ let reduce_unary psys = let is_unary_equation (cstr, prf) = if cstr.op == Eq then Vect.find - (fun v n -> if n =/ Int 1 || n =/ Int (-1) then Some v else None) + (fun v n -> if n =/ Q.one || n =/ Q.neg_one then Some v else None) cstr.coeffs else None in @@ -775,13 +769,12 @@ let reduce_var_change psys = match Vect.choose vect with | None -> None | Some (x, v, vect) -> ( - let v = numerator v in + let v = Q.num v in match Vect.find (fun x' v' -> - let v' = numerator v' in - if eq_big_int (gcd_big_int v v') unit_big_int then Some (x', v') - else None) + let v' = Q.num v' in + if Z_.equal (Z_.gcd v v') Z_.one then Some (x', v') else None) vect with | Some (x', v') -> Some ((x, v), (x', v')) @@ -795,12 +788,12 @@ let reduce_var_change psys = | None -> None | Some (((x, v), (x', v')), (c, p)) -> let l1, l2 = ext_gcd v v' in - let l1, l2 = (Big_int l1, Big_int l2) in + let l1, l2 = (Q.of_bigint l1, Q.of_bigint l2) in let pivot_eq (c', p') = let {coeffs; op; cst} = c' in let vx = Vect.get x coeffs in let vx' = Vect.get x' coeffs in - let m = minus_num ((vx */ l1) +/ (vx' */ l2)) in + let m = Q.neg ((vx */ l1) +/ (vx' */ l2)) in Some ( { coeffs = Vect.add (Vect.mul m c.coeffs) coeffs ; op @@ -818,7 +811,7 @@ let reduction_equations psys = (** [get_bound sys] returns upon success an interval (lb,e,ub) with proofs *) let get_bound sys = let is_small (v, i) = - match Itv.range i with None -> false | Some i -> i <=/ Int 1 + match Itv.range i with None -> false | Some i -> i <=/ Q.one in let select_best (x1, i1) (x2, i2) = if Itv.smaller_itv i1 i2 then (x1, i1) else (x2, i2) @@ -858,18 +851,20 @@ let get_bound sys = in match smallest_interval with | Some (lb, e, ub) -> ( - let lbn, lbd = (sub_big_int (numerator lb) unit_big_int, denominator lb) in - let ubn, ubd = (add_big_int unit_big_int (numerator ub), denominator ub) in + let lbn, lbd = (Z_.sub (Q.num lb) Z_.one, Q.den lb) in + let ubn, ubd = (Z_.add Z_.one (Q.num ub), Q.den ub) in (* x <= ub -> x > ub *) match ( direct_linear_prover - ( {coeffs = Vect.mul (Big_int ubd) e; op = Ge; cst = Big_int ubn} + ( { coeffs = Vect.mul (Q.of_bigint ubd) e + ; op = Ge + ; cst = Q.of_bigint ubn } :: sys ) , (* lb <= x -> lb > x *) direct_linear_prover - ( { coeffs = Vect.mul (minus_num (Big_int lbd)) e + ( { coeffs = Vect.mul (Q.neg (Q.of_bigint lbd)) e ; op = Ge - ; cst = minus_num (Big_int lbn) } + ; cst = Q.neg (Q.of_bigint lbn) } :: sys ) ) with | Some cub, Some clb -> @@ -879,7 +874,7 @@ let get_bound sys = let check_sys sys = List.for_all - (fun (c, p) -> Vect.for_all (fun _ n -> sign_num n <> 0) c.coeffs) + (fun (c, p) -> Vect.for_all (fun _ n -> Q.sign n <> 0) c.coeffs) sys open ProofFormat @@ -896,8 +891,8 @@ let xlia (can_enum : bool) reduction_equations sys = | Some (prf1, (lb, e, ub), prf2) -> ( if debug then Printf.printf "Found interval: %a in [%s;%s] -> " Vect.pp e - (string_of_num lb) (string_of_num ub); - match start_enum id e (ceiling_num lb) (floor_num ub) sys with + (Q.to_string lb) (Q.to_string ub); + match start_enum id e (Q.ceiling lb) (Q.floor ub) sys with | Prf prfl -> Prf (ProofFormat.Enum @@ -916,7 +911,7 @@ let xlia (can_enum : bool) reduction_equations sys = match aux_lia (id + 1) ((eq, ProofFormat.Def id) :: sys) with | Unknown | Model _ -> Unknown | Prf prf -> ( - match start_enum id e (clb +/ Int 1) cub sys with + match start_enum id e (clb +/ Q.one) cub sys with | Prf l -> Prf (prf :: l) | _ -> Unknown ) and aux_lia (id : int) (sys : prf_sys) = @@ -964,7 +959,7 @@ let xlia (can_enum : bool) reduction_equations sys = if Mc.zChecker sys' prf then Some prf else raise Certificate.BadCertificate with Failure s -> (Printf.printf "%s" s ; Some prf) - *) + *) Prf prf let xlia_simplex env red sys = diff --git a/plugins/micromega/coq_micromega.ml b/plugins/micromega/coq_micromega.ml index 4b656f8e61..c3f59b4208 100644 --- a/plugins/micromega/coq_micromega.ml +++ b/plugins/micromega/coq_micromega.ml @@ -518,7 +518,7 @@ module M = struct | Mc.Zneg p -> EConstr.mkApp (Lazy.force coq_NEG, [|dump_positive p|]) let pp_z o x = - Printf.fprintf o "%s" (Big_int.string_of_big_int (CoqToCaml.z_big_int x)) + Printf.fprintf o "%s" (NumCompat.Z.to_string (CoqToCaml.z_big_int x)) let dump_q q = EConstr.mkApp @@ -636,14 +636,14 @@ module M = struct in pp_pol o e - (* let pp_clause pp_c o (f: 'cst clause) = - List.iter (fun ((p,_),(t,_)) -> Printf.fprintf o "(%a @%a)" (pp_pol pp_c) p Tag.pp t) f *) + (* let pp_clause pp_c o (f: 'cst clause) = + List.iter (fun ((p,_),(t,_)) -> Printf.fprintf o "(%a @%a)" (pp_pol pp_c) p Tag.pp t) f *) let pp_clause_tag o (f : 'cst clause) = List.iter (fun ((p, _), (t, _)) -> Printf.fprintf o "(_ @%a)" Tag.pp t) f - (* let pp_cnf pp_c o (f:'cst cnf) = - List.iter (fun l -> Printf.fprintf o "[%a]" (pp_clause pp_c) l) f *) + (* let pp_cnf pp_c o (f:'cst cnf) = + List.iter (fun l -> Printf.fprintf o "[%a]" (pp_clause pp_c) l) f *) let pp_cnf_tag o (f : 'cst cnf) = List.iter (fun l -> Printf.fprintf o "[%a]" pp_clause_tag l) f @@ -819,16 +819,16 @@ module M = struct let elements env = env.vars - (* let string_of_env gl env = - let rec string_of_env i env acc = - match env with - | [] -> acc - | e::env -> string_of_env (i+1) env - (IMap.add i - (Pp.string_of_ppcmds - (Printer.pr_econstr_env gl.env gl.sigma e)) acc) in - string_of_env 1 env IMap.empty - *) + (* let string_of_env gl env = + let rec string_of_env i env acc = + match env with + | [] -> acc + | e::env -> string_of_env (i+1) env + (IMap.add i + (Pp.string_of_ppcmds + (Printer.pr_econstr_env gl.env gl.sigma e)) acc) in + string_of_env 1 env IMap.empty + *) let pp gl env = let ppl = List.mapi @@ -951,7 +951,7 @@ module M = struct (* NB: R is a different story. Because it is axiomatised, reducing would not be effective. Therefore, there is a specific parser for constant over R - *) + *) let rconst_assoc = [ (coq_Rplus, fun x y -> Mc.CPlus (x, y)) @@ -1613,14 +1613,14 @@ let compact_proofs (cnf_ff : 'cst cnf) res (cnf_ff' : 'cst cnf) = in List.assoc formula new_cl in - (* if debug then - begin - Printf.printf "\ncompact_proof : %a %a %a" - (pp_ml_list prover.pp_f) (List.map fst old_cl) - prover.pp_prf prf - (pp_ml_list prover.pp_f) (List.map fst new_cl) ; - flush stdout - end ; *) + (* if debug then + begin + Printf.printf "\ncompact_proof : %a %a %a" + (pp_ml_list prover.pp_f) (List.map fst old_cl) + prover.pp_prf prf + (pp_ml_list prover.pp_f) (List.map fst new_cl) ; + flush stdout + end ; *) let res = try prover.compact prf remap with x when CErrors.noncritical x -> ( @@ -1790,14 +1790,14 @@ let micromega_tauto pre_process cnf spec prover env flush stdout end; (* Even if it does not work, this does not mean it is not provable - -- the prover is REALLY incomplete *) + -- the prover is REALLY incomplete *) (* if debug then - begin - (* recompute the proofs *) - match witness_list_tags prover cnf_ff' with - | None -> failwith "abstraction is wrong" - | Some res -> () - end ; *) + begin + (* recompute the proofs *) + match witness_list_tags prover cnf_ff' with + | None -> failwith "abstraction is wrong" + | Some res -> () + end ; *) let res' = compact_proofs cnf_ff res cnf_ff' in let ff', res', ids = (ff', res', Mc.ids_of_formula ff') in let res' = dump_list spec.proof_typ spec.dump_proof res' in @@ -2009,7 +2009,7 @@ let micromega_genr prover tac = let goal_vars = List.map (fun (_, i) -> List.nth env (i - 1)) vars in let arith_args = goal_props @ goal_vars in let kill_arith = Tacticals.New.tclTHEN tac_arith tac in - (* Tacticals.New.tclTHEN + (* Tacticals.New.tclTHEN (Tactics.keep []) (Tactics.tclABSTRACT None*) Tacticals.New.tclTHENS diff --git a/plugins/micromega/csdpcert.ml b/plugins/micromega/csdpcert.ml index 90dd81adf4..a636fb0bdf 100644 --- a/plugins/micromega/csdpcert.ml +++ b/plugins/micromega/csdpcert.ml @@ -14,7 +14,7 @@ (* *) (************************************************************************) -open Num +open NumCompat open Sos open Sos_types open Sos_lib @@ -96,7 +96,7 @@ let real_nonlinear_prover d l = | Axiom_lt i -> poly_mul p y | Axiom_eq i -> poly_mul (poly_pow p 2) y | _ -> failwith "monoids") - m (poly_const (Int 1)) + m (poly_const Q.one) , List.map snd m )) (sets_of_list neq) in @@ -127,7 +127,7 @@ let real_nonlinear_prover d l = match List.map (function Axiom_eq i -> i | _ -> failwith "error") neq with - | [] -> Rational_lt (Int 1) + | [] -> Rational_lt Q.one | l -> Monoid l in List.fold_right (fun x y -> Product (x, y)) lt sq @@ -146,7 +146,7 @@ let real_nonlinear_prover d l = let pure_sos l = let l = List.map (fun (e, o) -> (Mc.denorm e, o)) l in (* If there is no strict inequality, - I should nonetheless be able to try something - over Z > is equivalent to -1 >= *) + I should nonetheless be able to try something - over Z > is equivalent to -1 >= *) try let l = List.combine l (CList.interval 0 (List.length l - 1)) in let lt, i = @@ -162,11 +162,11 @@ let pure_sos l = , List.fold_right (fun (c, p) rst -> Sum (Product (Rational_lt c, Square (term_of_poly p)), rst)) - polys (Rational_lt (Int 0)) ) + polys (Rational_lt Q.zero) ) in let proof = Sum (Axiom_lt i, pos) in - (* let s,proof' = scale_certificate proof in - let cert = snd (cert_of_pos proof') in *) + (* let s,proof' = scale_certificate proof in + let cert = snd (cert_of_pos proof') in *) S (Some proof) with (* | Sos.CsdpNotFound -> F "Sos.CsdpNotFound" *) | any -> @@ -184,8 +184,8 @@ let main () = try let (prover, poly) = (input_value stdin : provername * micromega_polys) in let cert = run_prover prover poly in - (* Printf.fprintf chan "%a -> %a" print_list_term poly output_csdp_certificate cert ; - close_out chan ; *) + (* Printf.fprintf chan "%a -> %a" print_list_term poly output_csdp_certificate cert ; + close_out chan ; *) output_value stdout (cert : csdp_certificate); flush stdout; Marshal.to_channel chan (cert : csdp_certificate) []; diff --git a/plugins/micromega/itv.ml b/plugins/micromega/itv.ml index 214edb46ba..74a9657038 100644 --- a/plugins/micromega/itv.ml +++ b/plugins/micromega/itv.ml @@ -8,12 +8,13 @@ (* * (see LICENSE file for the text of the license) *) (************************************************************************) -(** Intervals (extracted from mfourier.ml) *) +open NumCompat +open Q.Notations -open Num +(** Intervals (extracted from mfourier.ml) *) (** The type of intervals is *) -type interval = num option * num option +type interval = Q.t option * Q.t option (** None models the absence of bound i.e. infinity As a result, - None , None -> \]-oo,+oo\[ @@ -26,11 +27,11 @@ type interval = num option * num option let pp o (n1, n2) = ( match n1 with | None -> output_string o "]-oo" - | Some n -> Printf.fprintf o "[%s" (string_of_num n) ); + | Some n -> Printf.fprintf o "[%s" (Q.to_string n) ); output_string o ","; match n2 with | None -> output_string o "+oo[" - | Some n -> Printf.fprintf o "%s]" (string_of_num n) + | Some n -> Printf.fprintf o "%s]" (Q.to_string n) (** if then interval [itv] is empty, [norm_itv itv] returns [None] otherwise, it returns [Some itv] *) @@ -51,11 +52,11 @@ let inter i1 i2 = | None, Some _ -> o2 | Some n1, Some n2 -> Some (f n1 n2) in - norm_itv (inter max_num l1 l2, inter min_num r1 r2) + norm_itv (inter Q.max l1 l2, inter Q.min r1 r2) let range = function | None, _ | _, None -> None - | Some i, Some j -> Some (floor_num j -/ ceiling_num i +/ Int 1) + | Some i, Some j -> Some (Q.floor j -/ Q.ceiling i +/ Q.one) let smaller_itv i1 i2 = match (range i1, range i2) with diff --git a/plugins/micromega/itv.mli b/plugins/micromega/itv.mli index c7164f2c98..0dec639353 100644 --- a/plugins/micromega/itv.mli +++ b/plugins/micromega/itv.mli @@ -7,13 +7,13 @@ (* * GNU Lesser General Public License Version 2.1 *) (* * (see LICENSE file for the text of the license) *) (************************************************************************) -open Num +open NumCompat -type interval = num option * num option +type interval = Q.t option * Q.t option val pp : out_channel -> interval -> unit val inter : interval -> interval -> interval option -val range : interval -> num option +val range : interval -> Q.t option val smaller_itv : interval -> interval -> bool -val in_bound : interval -> num -> bool +val in_bound : interval -> Q.t -> bool val norm_itv : interval -> interval option diff --git a/plugins/micromega/mfourier.ml b/plugins/micromega/mfourier.ml index da75137185..838dab8ec8 100644 --- a/plugins/micromega/mfourier.ml +++ b/plugins/micromega/mfourier.ml @@ -8,8 +8,9 @@ (* * (see LICENSE file for the text of the license) *) (************************************************************************) +open NumCompat +open Q.Notations open Util -open Num open Polynomial open Vect @@ -61,11 +62,11 @@ let pp_cstr o (vect, bnd) = let l, r = bnd in ( match l with | None -> () - | Some n -> Printf.fprintf o "%s <= " (string_of_num n) ); + | Some n -> Printf.fprintf o "%s <= " (Q.to_string n) ); Vect.pp o vect; match r with | None -> output_string o "\n" - | Some n -> Printf.fprintf o "<=%s\n" (string_of_num n) + | Some n -> Printf.fprintf o "<=%s\n" (Q.to_string n) let pp_system o sys = System.iter (fun vect ibnd -> pp_cstr o (vect, !ibnd.bound)) sys @@ -121,12 +122,12 @@ let normalise_cstr vect cinfo = | None -> Contradiction | Some (l, r) -> ( match Vect.choose vect with - | None -> if Itv.in_bound (l, r) (Int 0) then Redundant else Contradiction + | None -> if Itv.in_bound (l, r) Q.zero then Redundant else Contradiction | Some (_, n, _) -> Cstr ( Vect.div n vect , let divn x = x // n in - if Int.equal (sign_num n) 1 then + if Int.equal (Q.sign n) 1 then {cinfo with bound = (Option.map divn l, Option.map divn r)} else { cinfo with @@ -139,7 +140,7 @@ let normalise_cstr vect cinfo = let count v = Vect.fold (fun (n, p) _ vl -> - let sg = sign_num vl in + let sg = Q.sign vl in assert (sg <> 0); if Int.equal sg 1 then (n, p + 1) else (n + 1, p)) (0, 0) v @@ -181,20 +182,20 @@ let system_list sys = System.fold (fun k bi l -> (k, !bi) :: l) s [] (** [add (v1,c1) (v2,c2) ] - precondition: (c1 <>/ Int 0 && c2 <>/ Int 0) + precondition: (c1 <>/ Q.zero && c2 <>/ Q.zero) @return a pair [(v,ln)] such that [v] is the sum of vector [v1] divided by [c1] and vector [v2] divided by [c2] Note that the resulting vector is not normalised. *) let add (v1, c1) (v2, c2) = - assert (c1 <>/ Int 0 && c2 <>/ Int 0); - let res = mul_add (Int 1 // c1) v1 (Int 1 // c2) v2 in + assert (c1 <>/ Q.zero && c2 <>/ Q.zero); + let res = mul_add (Q.one // c1) v1 (Q.one // c2) v2 in (res, count res) let add (v1, c1) (v2, c2) = let res = add (v1, c1) (v2, c2) in - (* Printf.printf "add(%a,%s,%a,%s) -> %a\n" pp_vect v1 (string_of_num c1) pp_vect v2 (string_of_num c2) pp_vect (fst res) ;*) + (* Printf.printf "add(%a,%s,%a,%s) -> %a\n" pp_vect v1 (Q.to_string c1) pp_vect v2 (Q.to_string c2) pp_vect (fst res) ;*) res (** To perform Fourier elimination, constraints are categorised depending on the sign of the variable to eliminate. *) @@ -207,11 +208,11 @@ let add (v1, c1) (v2, c2) = *) let split x (vect : vector) info (l, m, r) = - match get x vect with - | Int 0 -> + let vl = get x vect in + if Q.zero =/ vl then (* The constraint does not mention [x], store it in m *) (l, (vect, info) :: m, r) - | vl -> + else (* otherwise *) let cons_bound lst bd = match bd with @@ -219,7 +220,7 @@ let split x (vect : vector) info (l, m, r) = | Some bnd -> (vl, vect, {info with bound = (Some bnd, None)}) :: lst in let lb, rb = info.bound in - if Int.equal (sign_num vl) 1 then (cons_bound l lb, m, cons_bound r rb) + if Int.equal (Q.sign vl) 1 then (cons_bound l lb, m, cons_bound r rb) else (* sign_num vl = -1 *) (cons_bound l rb, m, cons_bound r lb) @@ -239,8 +240,8 @@ let project vr sys = let {neg = n1; pos = p1; bound = bound1; prf = prf1} = info1 and {neg = n2; pos = p2; bound = bound2; prf = prf2} = info2 in let bnd1 = Option.get (fst bound1) and bnd2 = Option.get (fst bound2) in - let bound = (bnd1 // v1) +/ (bnd2 // minus_num v2) in - let vres, (n, p) = add (vect1, v1) (vect2, minus_num v2) in + let bound = (bnd1 // v1) +/ (bnd2 // Q.neg v2) in + let vres, (n, p) = add (vect1, v1) (vect2, Q.neg v2) in ( vres , { neg = n ; pos = p @@ -270,11 +271,11 @@ let project vr sys = *) let project_using_eq vr c vect bound prf (vect', info') = - match get vr vect' with - | Int 0 -> (vect', info') - | c2 -> - let c1 = if c2 >=/ Int 0 then minus_num c else c in - let c2 = abs_num c2 in + let c2 = get vr vect' in + if Q.zero =/ c2 then (vect', info') + else + let c1 = if c2 >=/ Q.zero then Q.neg c else c in + let c2 = Q.abs c2 in let vres, (n, p) = add (vect, c1) (vect', c2) in let cst = bound // c1 in let bndres = @@ -315,14 +316,14 @@ let eval_vect map vect = let val_v = IMap.find v map in (sum +/ (val_v */ vl), rst) with Not_found -> (sum, Vect.set v vl rst)) - (Int 0, Vect.null) vect + (Q.zero, Vect.null) vect (** [restrict_bound n sum itv] returns the interval of [x] given that (fst itv) <= x * n + sum <= (snd itv) *) let restrict_bound n sum (itv : interval) = let f x = (x -/ sum) // n in let l, r = itv in - match sign_num n with + match Q.sign n with | 0 -> if in_bound itv sum then (None, None) (* redundant *) else failwith "SystemContradiction" @@ -339,7 +340,7 @@ let bound_of_variable map v sys = match inter bnd (restrict_bound vl sum !iref.bound) with | None -> Printf.fprintf stdout "bound_of_variable: eval_vecr %a = %s,%a\n" - Vect.pp vect (Num.string_of_num sum) Vect.pp rst; + Vect.pp vect (Q.to_string sum) Vect.pp rst; Printf.fprintf stdout "current interval: %a\n" Itv.pp !iref.bound; failwith "bound_of_variable: impossible" | Some itv -> itv) @@ -348,12 +349,12 @@ let bound_of_variable map v sys = (** [pick_small_value bnd] picks a value being closed to zero within the interval *) let pick_small_value bnd = match bnd with - | None, None -> Int 0 - | None, Some i -> if Int 0 <=/ floor_num i then Int 0 else floor_num i - | Some i, None -> if i <=/ Int 0 then Int 0 else ceiling_num i + | None, None -> Q.zero + | None, Some i -> if Q.zero <=/ Q.floor i then Q.zero else Q.floor i + | Some i, None -> if i <=/ Q.zero then Q.zero else Q.ceiling i | Some i, Some j -> - if i <=/ Int 0 && Int 0 <=/ j then Int 0 - else if ceiling_num i <=/ floor_num j then ceiling_num i (* why not *) + if i <=/ Q.zero && Q.zero <=/ j then Q.zero + else if Q.ceiling i <=/ Q.floor j then Q.ceiling i (* why not *) else i (** [solution s1 sys_l = Some(sn,\[(vn-1,sn-1);...; (v1,s1)\]\@sys_l)] @@ -373,8 +374,8 @@ let solve_sys black_v choose_eq choose_variable sys sys_l = fst (List.find (fun ((v, _, _, _), _) -> v <> black_v) eqs) in if debug then ( - Printf.printf "\nE %a = %s variable %i\n" Vect.pp vect - (string_of_num cst) v; + Printf.printf "\nE %a = %s variable %i\n" Vect.pp vect (Q.to_string cst) + v; flush stdout ); let sys' = elim_var_using_eq v vect cst ln sys in solve_sys sys' ((v, sys) :: sys_l) @@ -422,7 +423,7 @@ module EstimateElimVar = struct | Some bnd -> (info.neg + info.pos) :: lst in let lb, rb = info.bound in - if Int.equal (sign_num vl) 1 then + if Int.equal (Q.sign vl) 1 then xpart rl ((rl1, info) :: ltl) (cons_bound n lb) z (cons_bound p rb) else @@ -568,7 +569,7 @@ module Fourier = struct (* We add a dummy (fresh) variable for vector *) let fresh = List.fold_left (fun fr c -> max fr (Vect.fresh c.coeffs)) 0 l in let cstr = - {coeffs = Vect.set fresh (Int (-1)) vect; op = Eq; cst = Int 0} + {coeffs = Vect.set fresh Q.neg_one vect; op = Eq; cst = Q.zero} in match solve fresh choose_equality_var choose_variable (cstr :: l) with | Inr prf -> None (* This is an unsatisfiability proof *) @@ -619,28 +620,27 @@ module Proof = struct let pivot v (p1, c1) (p2, c2) = let {coeffs = v1; op = op1; cst = n1} = c1 and {coeffs = v2; op = op2; cst = n2} = c2 in - match (Vect.get v v1, Vect.get v v2) with - | Int 0, _ | _, Int 0 -> None - | a, b -> - if Int.equal (sign_num a * sign_num b) (-1) then - Some - ( add (p1, abs_num a) (p2, abs_num b) - , { coeffs = add (v1, abs_num a) (v2, abs_num b) - ; op = add_op op1 op2 - ; cst = (n1 // abs_num a) +/ (n2 // abs_num b) } ) - else if op1 == Eq then - Some - ( add (p1, minus_num (a // b)) (p2, Int 1) - , { coeffs = add (v1, minus_num (a // b)) (v2, Int 1) - ; op = add_op op1 op2 - ; cst = (n1 // minus_num (a // b)) +/ (n2 // Int 1) } ) - else if op2 == Eq then - Some - ( add (p2, minus_num (b // a)) (p1, Int 1) - , { coeffs = add (v2, minus_num (b // a)) (v1, Int 1) - ; op = add_op op1 op2 - ; cst = (n2 // minus_num (b // a)) +/ (n1 // Int 1) } ) - else None + let a, b = (Vect.get v v1, Vect.get v v2) in + if Q.zero =/ a || Q.zero =/ b then None + else if Int.equal (Q.sign a * Q.sign b) (-1) then + Some + ( add (p1, Q.abs a) (p2, Q.abs b) + , { coeffs = add (v1, Q.abs a) (v2, Q.abs b) + ; op = add_op op1 op2 + ; cst = (n1 // Q.abs a) +/ (n2 // Q.abs b) } ) + else if op1 == Eq then + Some + ( add (p1, Q.neg (a // b)) (p2, Q.one) + , { coeffs = add (v1, Q.neg (a // b)) (v2, Q.one) + ; op = add_op op1 op2 + ; cst = (n1 // Q.neg (a // b)) +/ (n2 // Q.one) } ) + else if op2 == Eq then + Some + ( add (p2, Q.neg (b // a)) (p1, Q.one) + , { coeffs = add (v2, Q.neg (b // a)) (v1, Q.one) + ; op = add_op op1 op2 + ; cst = (n2 // Q.neg (b // a)) +/ (n1 // Q.one) } ) + else None (* op2 could be Eq ... this might happen *) @@ -656,7 +656,7 @@ module Proof = struct | Cstr (v, info) -> Inl ((prf, cstr, v, info) :: acc) )) (Inl []) l - type oproof = (vector * cstr * num) option + type oproof = (vector * cstr * Q.t) option let merge_proof (oleft : oproof) (prf, cstr, v, info) (oright : oproof) = let l, r = info.bound in @@ -679,7 +679,7 @@ module Proof = struct (* There is a contradiction - it should show up by scaling up the vectors - any pivot should do*) match Vect.choose cstrr.coeffs with | None -> - Inr (add (prfl, Int 1) (prfr, Int 1), cstrr) (* this is wrong *) + Inr (add (prfl, Q.one) (prfr, Q.one), cstrr) (* this is wrong *) | Some (v, _, _) -> ( match pivot v (prfl, cstrl) (prfr, cstrr) with | None -> failwith "merge_proof : pivot is not possible" @@ -687,12 +687,12 @@ module Proof = struct let mk_proof hyps prf = (* I am keeping list - I might have a proof for the left bound and a proof for the right bound. - If I perform aggressive elimination of redundancies, I expect the list to be of length at most 2. - For each proof list, all the vectors should be of the form a.v for different constants a. - *) + If I perform aggressive elimination of redundancies, I expect the list to be of length at most 2. + For each proof list, all the vectors should be of the form a.v for different constants a. + *) let rec mk_proof prf = match prf with - | Assum i -> [(Vect.set i (Int 1) Vect.null, List.nth hyps i)] + | Assum i -> [(Vect.set i Q.one Vect.null, List.nth hyps i)] | Elim (v, prf1, prf2) -> let prfsl = mk_proof prf1 and prfsr = mk_proof prf2 in (* I take only the pairs for which the elimination is meaningful *) diff --git a/plugins/micromega/micromega_plugin.mlpack b/plugins/micromega/micromega_plugin.mlpack index e3aa0dab7d..2630e883c9 100644 --- a/plugins/micromega/micromega_plugin.mlpack +++ b/plugins/micromega/micromega_plugin.mlpack @@ -1,4 +1,5 @@ Micromega +NumCompat Mutils Itv Vect diff --git a/plugins/micromega/mutils.ml b/plugins/micromega/mutils.ml index 51f0328e4b..2e054a21c2 100644 --- a/plugins/micromega/mutils.ml +++ b/plugins/micromega/mutils.ml @@ -19,6 +19,9 @@ (* *) (************************************************************************) +open NumCompat +module Z_ = NumCompat.Z + module Int = struct type t = int @@ -159,24 +162,6 @@ let saturate_bin (type a) (module Set : Set.S with type elt = a) let s0 = List.fold_left (fun acc e -> Set.add e acc) Set.empty l in Set.elements (Set.diff (iterate Set.empty s0) s0) -open Num -open Big_int - -let ppcm x y = - let g = gcd_big_int x y in - let x' = div_big_int x g in - let y' = div_big_int y g in - mult_big_int g (mult_big_int x' y') - -let denominator = function - | Int _ | Big_int _ -> unit_big_int - | Ratio r -> Ratio.denominator_ratio r - -let numerator = function - | Ratio r -> Ratio.numerator_ratio r - | Int i -> Big_int.big_int_of_int i - | Big_int i -> i - let iterate_until_stable f x = let rec iter x = match f x with None -> x | Some x' -> iter x' in iter x @@ -207,24 +192,23 @@ module CoqToCaml = struct (* Swap left-right ? *) match i with XH -> 1 | XI i -> 1 + (2 * index i) | XO i -> 2 * index i - open Big_int - let rec positive_big_int p = match p with - | XH -> unit_big_int - | XI p -> add_int_big_int 1 (mult_int_big_int 2 (positive_big_int p)) - | XO p -> mult_int_big_int 2 (positive_big_int p) + | XH -> Z_.one + | XI p -> Z_.add Z_.one (Z_.mul Z_.two (positive_big_int p)) + | XO p -> Z_.mul Z_.two (positive_big_int p) let z_big_int x = match x with - | Z0 -> zero_big_int + | Z0 -> Z_.zero | Zpos p -> positive_big_int p - | Zneg p -> minus_big_int (positive_big_int p) + | Zneg p -> Z_.neg (positive_big_int p) let z x = match x with Z0 -> 0 | Zpos p -> index p | Zneg p -> -index p let q_to_num {qnum = x; qden = y} = - Big_int (z_big_int x) // Big_int (z_big_int (Zpos y)) + let open Q.Notations in + Q.of_bigint (z_big_int x) // Q.of_bigint (z_big_int (Zpos y)) end (** @@ -259,27 +243,24 @@ module CamlToCoq = struct (* this should be -1 *) Zneg (positive (-x)) - open Big_int - let positive_big_int n = - let two = big_int_of_int 2 in let rec _pos n = - if eq_big_int n unit_big_int then XH + if Z_.equal n Z_.one then XH else - let q, m = quomod_big_int n two in - if eq_big_int unit_big_int m then XI (_pos q) else XO (_pos q) + let q, m = Z_.quomod n Z_.two in + if Z_.equal Z_.one m then XI (_pos q) else XO (_pos q) in _pos n let bigint x = - match sign_big_int x with + match Z_.sign x with | 0 -> Z0 | 1 -> Zpos (positive_big_int x) - | _ -> Zneg (positive_big_int (minus_big_int x)) + | _ -> Zneg (positive_big_int (Z_.neg x)) let q n = - { Micromega.qnum = bigint (numerator n) - ; Micromega.qden = positive_big_int (denominator n) } + { Micromega.qnum = bigint (Q.num n) + ; Micromega.qden = positive_big_int (Q.den n) } end (** diff --git a/plugins/micromega/mutils.mli b/plugins/micromega/mutils.mli index 9badddc255..a03b03ed8e 100644 --- a/plugins/micromega/mutils.mli +++ b/plugins/micromega/mutils.mli @@ -8,6 +8,8 @@ (* * (see LICENSE file for the text of the license) *) (************************************************************************) +open NumCompat + module Int : sig type t = int @@ -28,9 +30,6 @@ module IMap : sig (** [from k m] returns the submap of [m] with keys greater or equal k *) end -val numerator : Num.num -> Big_int.big_int -val denominator : Num.num -> Big_int.big_int - module Cmp : sig val compare_list : ('a -> 'b -> int) -> 'a list -> 'b list -> int val compare_lexical : (unit -> int) list -> int @@ -53,19 +52,19 @@ val pp_list : module CamlToCoq : sig val positive : int -> Micromega.positive - val bigint : Big_int.big_int -> Micromega.z + val bigint : Z.t -> Micromega.z val n : int -> Micromega.n val nat : int -> Micromega.nat - val q : Num.num -> Micromega.q + val q : Q.t -> Micromega.q val index : int -> Micromega.positive val z : int -> Micromega.z - val positive_big_int : Big_int.big_int -> Micromega.positive + val positive_big_int : Z.t -> Micromega.positive end module CoqToCaml : sig - val z_big_int : Micromega.z -> Big_int.big_int + val z_big_int : Micromega.z -> Z.t val z : Micromega.z -> int - val q_to_num : Micromega.q -> Num.num + val q_to_num : Micromega.q -> Q.t val positive : Micromega.positive -> int val n : Micromega.n -> int val nat : Micromega.nat -> int @@ -96,7 +95,6 @@ module Hash : sig val hash_elt : ('a -> int) -> int -> 'a -> int end -val ppcm : Big_int.big_int -> Big_int.big_int -> Big_int.big_int val all_pairs : ('a -> 'a -> 'b) -> 'a list -> 'b list val try_any : (('a -> 'b option) * 'c) list -> 'a -> 'b option val is_sublist : ('a -> 'b -> bool) -> 'a list -> 'b list -> bool diff --git a/plugins/micromega/numCompat.ml b/plugins/micromega/numCompat.ml new file mode 100644 index 0000000000..82993cd730 --- /dev/null +++ b/plugins/micromega/numCompat.ml @@ -0,0 +1,174 @@ +(************************************************************************) +(* * The Coq Proof Assistant / The Coq Development Team *) +(* v * INRIA, CNRS and contributors - Copyright 1999-2019 *) +(* <O___,, * (see CREDITS file for the list of authors) *) +(* \VV/ **************************************************************) +(* // * This file is distributed under the terms of the *) +(* * GNU Lesser General Public License Version 2.1 *) +(* * (see LICENSE file for the text of the license) *) +(************************************************************************) + +module type ZArith = sig + type t + + val zero : t + val one : t + val two : t + val add : t -> t -> t + val sub : t -> t -> t + val mul : t -> t -> t + val div : t -> t -> t + val neg : t -> t + val sign : t -> int + val equal : t -> t -> bool + val compare : t -> t -> int + val power_int : t -> int -> t + val quomod : t -> t -> t * t + val ppcm : t -> t -> t + val gcd : t -> t -> t + val lcm : t -> t -> t + val to_string : t -> string +end + +module Z = struct + type t = Big_int.big_int + + open Big_int + + let zero = zero_big_int + let one = unit_big_int + let two = big_int_of_int 2 + let add = Big_int.add_big_int + let sub = Big_int.sub_big_int + let mul = Big_int.mult_big_int + let div = Big_int.div_big_int + let neg = Big_int.minus_big_int + let sign = Big_int.sign_big_int + let equal = eq_big_int + let compare = compare_big_int + let power_int = power_big_int_positive_int + let quomod = quomod_big_int + + let ppcm x y = + let g = gcd_big_int x y in + let x' = div_big_int x g in + let y' = div_big_int y g in + mult_big_int g (mult_big_int x' y') + + let gcd = gcd_big_int + + let lcm x y = + if eq_big_int x zero && eq_big_int y zero then zero + else abs_big_int (div_big_int (mult_big_int x y) (gcd x y)) + + let to_string = string_of_big_int +end + +module type QArith = sig + module Z : ZArith + + type t + + val of_int : int -> t + val zero : t + val one : t + val two : t + val ten : t + val neg_one : t + + module Notations : sig + val ( // ) : t -> t -> t + val ( +/ ) : t -> t -> t + val ( -/ ) : t -> t -> t + val ( */ ) : t -> t -> t + val ( =/ ) : t -> t -> bool + val ( <>/ ) : t -> t -> bool + val ( >/ ) : t -> t -> bool + val ( >=/ ) : t -> t -> bool + val ( </ ) : t -> t -> bool + val ( <=/ ) : t -> t -> bool + end + + val compare : t -> t -> int + val make : Z.t -> Z.t -> t + val den : t -> Z.t + val num : t -> Z.t + val of_bigint : Z.t -> t + val to_bigint : t -> Z.t + val neg : t -> t + + (* val inv : t -> t *) + val max : t -> t -> t + val min : t -> t -> t + val sign : t -> int + val abs : t -> t + val mod_ : t -> t -> t + val floor : t -> t + + (* val floorZ : t -> Z.t *) + val ceiling : t -> t + val round : t -> t + val pow2 : int -> t + val pow10 : int -> t + val power : int -> t -> t + val to_string : t -> string + val of_string : string -> t + val to_float : t -> float +end + +module Q : QArith with module Z = Z = struct + module Z = Z + + type t = Num.num + + open Num + + let of_int x = Int x + let zero = Int 0 + let one = Int 1 + let two = Int 2 + let ten = Int 10 + let neg_one = Int (-1) + + module Notations = struct + let ( // ) = div_num + let ( +/ ) = add_num + let ( -/ ) = sub_num + let ( */ ) = mult_num + let ( =/ ) = eq_num + let ( <>/ ) = ( <>/ ) + let ( >/ ) = ( >/ ) + let ( >=/ ) = ( >=/ ) + let ( </ ) = ( </ ) + let ( <=/ ) = ( <=/ ) + end + + let compare = compare_num + let make x y = Big_int x // Big_int y + + let numdom r = + let r' = Ratio.normalize_ratio (ratio_of_num r) in + (Ratio.numerator_ratio r', Ratio.denominator_ratio r') + + let num x = numdom x |> fst + let den x = numdom x |> snd + let of_bigint x = Big_int x + let to_bigint = big_int_of_num + let neg = minus_num + + (* let inv = *) + let max = max_num + let min = min_num + let sign = sign_num + let abs = abs_num + let mod_ = mod_num + let floor = floor_num + let ceiling = ceiling_num + let round = round_num + let pow2 n = power_num two (Int n) + let pow10 n = power_num ten (Int n) + let power x = power_num (Int x) + let to_string = string_of_num + let of_string = num_of_string + let to_float = float_of_num +end diff --git a/plugins/micromega/numCompat.mli b/plugins/micromega/numCompat.mli new file mode 100644 index 0000000000..183285e259 --- /dev/null +++ b/plugins/micromega/numCompat.mli @@ -0,0 +1,85 @@ +(************************************************************************) +(* * The Coq Proof Assistant / The Coq Development Team *) +(* v * INRIA, CNRS and contributors - Copyright 1999-2019 *) +(* <O___,, * (see CREDITS file for the list of authors) *) +(* \VV/ **************************************************************) +(* // * This file is distributed under the terms of the *) +(* * GNU Lesser General Public License Version 2.1 *) +(* * (see LICENSE file for the text of the license) *) +(************************************************************************) + +module type ZArith = sig + type t + + val zero : t + val one : t + val two : t + val add : t -> t -> t + val sub : t -> t -> t + val mul : t -> t -> t + val div : t -> t -> t + val neg : t -> t + val sign : t -> int + val equal : t -> t -> bool + val compare : t -> t -> int + val power_int : t -> int -> t + val quomod : t -> t -> t * t + val ppcm : t -> t -> t + val gcd : t -> t -> t + val lcm : t -> t -> t + val to_string : t -> string +end + +module type QArith = sig + module Z : ZArith + + type t + + val of_int : int -> t + val zero : t + val one : t + val two : t + val ten : t + val neg_one : t + + module Notations : sig + val ( // ) : t -> t -> t + val ( +/ ) : t -> t -> t + val ( -/ ) : t -> t -> t + val ( */ ) : t -> t -> t + val ( =/ ) : t -> t -> bool + val ( <>/ ) : t -> t -> bool + val ( >/ ) : t -> t -> bool + val ( >=/ ) : t -> t -> bool + val ( </ ) : t -> t -> bool + val ( <=/ ) : t -> t -> bool + end + + val compare : t -> t -> int + val make : Z.t -> Z.t -> t + val den : t -> Z.t + val num : t -> Z.t + val of_bigint : Z.t -> t + val to_bigint : t -> Z.t + val neg : t -> t + + (* val inv : t -> t *) + + val max : t -> t -> t + val min : t -> t -> t + val sign : t -> int + val abs : t -> t + val mod_ : t -> t -> t + val floor : t -> t + val ceiling : t -> t + val round : t -> t + val pow2 : int -> t + val pow10 : int -> t + val power : int -> t -> t + val to_string : t -> string + val of_string : string -> t + val to_float : t -> float +end + +module Z : ZArith +module Q : QArith with module Z = Z diff --git a/plugins/micromega/persistent_cache.ml b/plugins/micromega/persistent_cache.ml index d5b28cb03e..4777b5e231 100644 --- a/plugins/micromega/persistent_cache.ml +++ b/plugins/micromega/persistent_cache.ml @@ -82,9 +82,9 @@ module PHashtable (Key : HashedType) : PHashtable with type key = Key.t = struct with Unix.Unix_error (_, _, _) -> () (* Here, this is really bad news -- - there is a pending lock which could cause a deadlock. - Should it be an anomaly or produce a warning ? - *); + there is a pending lock which could cause a deadlock. + Should it be an anomaly or produce a warning ? + *); ignore (lseek fd pos SEEK_SET) (* We make the assumption that an acquired lock can always be released *) diff --git a/plugins/micromega/polynomial.ml b/plugins/micromega/polynomial.ml index f83b36d847..68aa739a6f 100644 --- a/plugins/micromega/polynomial.ml +++ b/plugins/micromega/polynomial.ml @@ -14,7 +14,8 @@ (* *) (************************************************************************) -open Num +open NumCompat +open Q.Notations open Mutils module Mc = Micromega @@ -23,8 +24,8 @@ let max_nb_cstr = ref max_int type var = int let debug = false -let ( <+> ) = add_num -let ( <*> ) = mult_num +let ( <+> ) = ( +/ ) +let ( <*> ) = ( */ ) module Monomial : sig type t @@ -153,13 +154,11 @@ end let pp_mon o (m, i) = if Monomial.is_const m then - if eq_num (Int 0) i then () else Printf.fprintf o "%s" (string_of_num i) - else - match i with - | Int 1 -> Monomial.pp o m - | Int -1 -> Printf.fprintf o "-%a" Monomial.pp m - | Int 0 -> () - | _ -> Printf.fprintf o "%s*%a" (string_of_num i) Monomial.pp m + if Q.zero =/ i then () else Printf.fprintf o "%s" (Q.to_string i) + else if Q.one =/ i then Monomial.pp o m + else if Q.neg_one =/ i then Printf.fprintf o "-%a" Monomial.pp m + else if Q.zero =/ i then () + else Printf.fprintf o "%s*%a" (Q.to_string i) Monomial.pp m module Poly : (* A polynomial is a map of monomials *) (* @@ -171,51 +170,51 @@ sig type t val pp : out_channel -> t -> unit - val get : Monomial.t -> t -> num + val get : Monomial.t -> t -> Q.t val variable : var -> t - val add : Monomial.t -> num -> t -> t - val constant : num -> t + val add : Monomial.t -> Q.t -> t -> t + val constant : Q.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 fold : (Monomial.t -> Q.t -> 'a -> 'a) -> t -> 'a -> 'a val factorise : var -> t -> t * t end = struct (*normalisation bug : 0*x ... *) module P = Map.Make (Monomial) open P - type t = num P.t + type t = Q.t P.t let pp o p = P.iter (fun mn i -> Printf.fprintf o "%a + " pp_mon (mn, i)) 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 + let get : Monomial.t -> t -> Q.t = + fun mn p -> try find mn p with Not_found -> Q.zero (* The polynomial 1.x *) - let variable : var -> t = fun x -> add (Monomial.var x) (Int 1) empty + let variable : var -> t = fun x -> add (Monomial.var x) Q.one empty (*The constant polynomial *) - let constant : num -> t = fun c -> add Monomial.const c empty + let constant : Q.t -> t = fun c -> add Monomial.const c empty (* The addition of a monomial *) - let add : Monomial.t -> num -> t -> t = + let add : Monomial.t -> Q.t -> t -> t = fun mn v p -> - if sign_num v = 0 then p + if Q.sign v = 0 then p else let vl = get mn p <+> v in - if sign_num vl = 0 then remove mn p else add mn vl p + if Q.sign vl = 0 then remove mn p else 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 = + let mult : Monomial.t -> Q.t -> t -> t = fun mn v p -> - if sign_num v = 0 then constant (Int 0) + if Q.sign v = 0 then constant Q.zero else fold (fun mn' v' res -> P.add (Monomial.prod mn mn') (v <*> v') res) @@ -227,7 +226,7 @@ end = struct 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 uminus : t -> t = fun p -> map (fun v -> Q.neg v) p let fold = P.fold let factorise x p = @@ -240,12 +239,12 @@ end = struct let mx = Monomial.prod m1 (Monomial.exp x (i - 1)) in (add mx v px, cx)) p - (constant (Int 0), constant (Int 0)) + (constant Q.zero, constant Q.zero) end type vector = Vect.t -type cstr = {coeffs : vector; op : op; cst : num} +type cstr = {coeffs : vector; op : op; cst : Q.t} and op = Eq | Ge | Gt @@ -256,8 +255,7 @@ let eval_op = function Eq -> ( =/ ) | Ge -> ( >=/ ) | Gt -> ( >/ ) let string_of_op = function Eq -> "=" | Ge -> ">=" | Gt -> ">" let output_cstr o {coeffs; op; cst} = - Printf.fprintf o "%a %s %s" Vect.pp coeffs (string_of_op op) - (string_of_num cst) + Printf.fprintf o "%a %s %s" Vect.pp coeffs (string_of_op op) (Q.to_string cst) let opMult o1 o2 = match (o1, o2) with Eq, _ | _, Eq -> Eq | Ge, _ | _, Ge -> Ge | Gt, Gt -> Gt @@ -308,11 +306,11 @@ module LinPoly = struct let _ = register Monomial.const end - let var v = Vect.set (MonT.register (Monomial.var v)) (Int 1) Vect.null + let var v = Vect.set (MonT.register (Monomial.var v)) Q.one Vect.null let of_monomial m = let v = MonT.register m in - Vect.set v (Int 1) Vect.null + Vect.set v Q.one Vect.null let linpol_of_pol p = Poly.fold @@ -324,7 +322,7 @@ module LinPoly = struct let pol_of_linpol v = Vect.fold (fun p vr n -> Poly.add (MonT.retrieve vr) n p) - (Poly.constant (Int 0)) v + (Poly.constant Q.zero) v let coq_poly_of_linpol cst p = let pol_of_mon m = @@ -332,13 +330,13 @@ module LinPoly = struct (fun x v p -> Mc.PEmul (Mc.PEpow (Mc.PEX (CamlToCoq.positive x), CamlToCoq.n v), p)) m - (Mc.PEc (cst (Int 1))) + (Mc.PEc (cst Q.one)) in Vect.fold (fun acc x v -> let mn = MonT.retrieve x in Mc.PEadd (Mc.PEmul (Mc.PEc (cst v), pol_of_mon mn), acc)) - (Mc.PEc (cst (Int 0))) + (Mc.PEc (cst Q.zero)) p let pp_var o vr = @@ -346,7 +344,7 @@ module LinPoly = struct with Not_found -> Printf.fprintf o "v%i" vr let pp o p = Vect.pp_gen pp_var o p - let constant c = if sign_num c = 0 then Vect.null else Vect.set 0 c Vect.null + let constant c = if Q.sign c = 0 then Vect.null else Vect.set 0 c Vect.null let is_linear p = Vect.for_all @@ -357,7 +355,7 @@ module LinPoly = struct let is_variable p = let (x, v), r = Vect.decomp_fst p in - if Vect.is_null r && v >/ Int 0 then Monomial.get_var (MonT.retrieve x) + if Vect.is_null r && v >/ Q.zero then Monomial.get_var (MonT.retrieve x) else None let factorise x p = @@ -431,17 +429,15 @@ module LinPoly = struct end module ProofFormat = struct - open Big_int - type prf_rule = | Annot of string * prf_rule | Hyp of int | Def of int - | Cst of Num.num + | Cst of Q.t | Zero | Square of Vect.t | MulC of Vect.t * prf_rule - | Gcd of Big_int.big_int * prf_rule + | Gcd of Z.t * prf_rule | MulPrf of prf_rule * prf_rule | AddPrf of prf_rule * prf_rule | CutPrf of prf_rule @@ -458,7 +454,7 @@ module ProofFormat = struct | Annot (s, p) -> Printf.fprintf o "(%a)@%s" output_prf_rule p s | Hyp i -> Printf.fprintf o "Hyp %i" i | Def i -> Printf.fprintf o "Def %i" i - | Cst c -> Printf.fprintf o "Cst %s" (string_of_num c) + | Cst c -> Printf.fprintf o "Cst %s" (Q.to_string c) | Zero -> Printf.fprintf o "Zero" | Square s -> Printf.fprintf o "(%a)^2" Poly.pp (LinPoly.pol_of_linpol s) | MulC (p, pr) -> @@ -469,8 +465,7 @@ module ProofFormat = struct | AddPrf (p1, p2) -> Printf.fprintf o "%a + %a" output_prf_rule p1 output_prf_rule p2 | CutPrf p -> Printf.fprintf o "[%a]" output_prf_rule p - | Gcd (c, p) -> - Printf.fprintf o "(%a)/%s" output_prf_rule p (string_of_big_int c) + | Gcd (c, p) -> Printf.fprintf o "(%a)/%s" output_prf_rule p (Z.to_string c) let rec output_proof o = function | Done -> Printf.fprintf o "." @@ -485,11 +480,11 @@ module ProofFormat = struct let rec pr_size = function | Annot (_, p) -> pr_size p - | Zero | Square _ -> Int 0 - | Hyp _ -> Int 1 - | Def _ -> Int 1 + | Zero | Square _ -> Q.zero + | Hyp _ -> Q.one + | Def _ -> Q.one | Cst n -> n - | Gcd (i, p) -> pr_size p // Big_int i + | Gcd (i, p) -> pr_size p // Q.of_bigint i | MulPrf (p1, p2) | AddPrf (p1, p2) -> pr_size p1 +/ pr_size p2 | CutPrf p -> pr_size p | MulC (v, p) -> pr_size p @@ -601,12 +596,12 @@ module ProofFormat = struct (id, ExProof (i, j, k, x, z, t, prf)) | Enum (i, p1, v, p2, pl) -> (* Why do I have top-level cuts ? *) - (* let p1 = implicit_cut p1 in - let p2 = implicit_cut p2 in - let (ids,prfs) = List.split (List.map (normalise_proof id) pl) in - (List.fold_left max 0 ids , - Enum(i,p1,v,p2,prfs)) - *) + (* let p1 = implicit_cut p1 in + let p2 = implicit_cut p2 in + let (ids,prfs) = List.split (List.map (normalise_proof id) pl) in + (List.fold_left max 0 ids , + Enum(i,p1,v,p2,prfs)) + *) let bds1, id, p1' = pr_rule_def_cut id (implicit_cut p1) in let bds2, id, p2' = pr_rule_def_cut id (implicit_cut p2) in let ids, prfs = List.split (List.map (normalise_proof id) pl) in @@ -649,13 +644,13 @@ module ProofFormat = struct if s1 = s2 then compare p1 p2 else String.compare s1 s2 | Hyp i, Hyp j -> Int.compare i j | Def i, Def j -> Int.compare i j - | Cst n, Cst m -> Num.compare_num n m + | Cst n, Cst m -> Q.compare n m | Zero, Zero -> 0 | Square v1, Square v2 -> Vect.compare v1 v2 | MulC (v1, p1), MulC (v2, p2) -> cmp_pair Vect.compare compare (v1, p1) (v2, p2) | Gcd (b1, p1), Gcd (b2, p2) -> - cmp_pair Big_int.compare_big_int compare (b1, p1) (b2, p2) + cmp_pair Z.compare compare (b1, p1) (b2, p2) | MulPrf (p1, q1), MulPrf (p2, q2) -> cmp_pair compare compare (p1, q1) (p2, q2) | AddPrf (p1, q1), MulPrf (p2, q2) -> @@ -672,11 +667,11 @@ module ProofFormat = struct | Annot (s, p) -> Annot (s, mul_cst_proof c p) | MulC (v, p') -> MulC (Vect.mul c v, p') | _ -> ( - match sign_num c with + match Q.sign c with | 0 -> Zero (* This is likely to be a bug *) | -1 -> MulC (LinPoly.constant c, p) (* [p] should represent an equality *) - | 1 -> if eq_num (Int 1) c then p else MulPrf (Cst c, p) + | 1 -> if Q.one =/ c then p else MulPrf (Cst c, p) | _ -> assert false ) let sMulC v p = @@ -698,7 +693,7 @@ module ProofFormat = struct match p with | Annot (s, p) -> dev_prf_rule p | Hyp _ | Def _ | Cst _ | Zero | Square _ -> - PrfRuleMap.singleton p (LinPoly.constant (Int 1)) + PrfRuleMap.singleton p (LinPoly.constant Q.one) | MulC (v, p) -> PrfRuleMap.map (fun v1 -> LinPoly.product v v1) (dev_prf_rule p) | AddPrf (p1, p2) -> @@ -716,9 +711,9 @@ module ProofFormat = struct let p2'' = prf_rule_of_map p2' in match p1'' with | Cst c -> PrfRuleMap.map (fun v1 -> Vect.mul c v1) p2' - | _ -> - PrfRuleMap.singleton (MulPrf (p1'', p2'')) (LinPoly.constant (Int 1)) ) - | _ -> PrfRuleMap.singleton p (LinPoly.constant (Int 1)) + | _ -> PrfRuleMap.singleton (MulPrf (p1'', p2'')) (LinPoly.constant Q.one) + ) + | _ -> PrfRuleMap.singleton p (LinPoly.constant Q.one) let simplify_prf_rule p = prf_rule_of_map (dev_prf_rule p) @@ -766,7 +761,7 @@ module ProofFormat = struct xid_of_hyp 0 l end - let cmpl_prf_rule norm (cst : num -> 'a) env prf = + let cmpl_prf_rule norm (cst : Q.t -> 'a) env prf = let rec cmpl = function | Annot (s, p) -> cmpl p | Hyp i | Def i -> Mc.PsatzIn (CamlToCoq.nat (Env.id_of_hyp i env)) @@ -783,7 +778,7 @@ module ProofFormat = struct cmpl prf let cmpl_prf_rule_z env r = - cmpl_prf_rule Mc.normZ (fun x -> CamlToCoq.bigint (numerator x)) env r + cmpl_prf_rule Mc.normZ (fun x -> CamlToCoq.bigint (Q.num x)) env r let rec cmpl_proof env = function | Done -> Mc.DoneProof @@ -810,7 +805,7 @@ module ProofFormat = struct | Hyp i | Def i -> env i | Cst n -> ( ( Vect.set 0 n Vect.null - , match Num.compare_num n (Int 0) with + , match Q.compare n Q.zero with | 0 -> Ge | 1 -> Gt | _ -> failwith "eval_prf_rule : negative constant" ) ) @@ -826,7 +821,7 @@ module ProofFormat = struct failwith "eval_prf_rule : not an equality" ) | Gcd (g, p) -> let v, op = eval_prf_rule env p in - (Vect.div (Big_int g) v, op) + (Vect.div (Q.of_bigint g) v, op) | MulPrf (p1, p2) -> let v1, o1 = eval_prf_rule env p1 in let v2, o2 = eval_prf_rule env p2 in @@ -839,7 +834,7 @@ module ProofFormat = struct let is_unsat (p, o) = let c, r = Vect.decomp_cst p in - if Vect.is_null r then not (eval_op o c (Int 0)) else false + if Vect.is_null r then not (eval_op o c Q.zero) else false let rec eval_proof env p = match p with @@ -882,7 +877,7 @@ module WithProof = struct let zero = ((Vect.null, Eq), ProofFormat.Zero) let const n = ((LinPoly.constant n, Ge), ProofFormat.Cst n) - let of_cstr (c, prf) = ((Vect.set 0 (Num.minus_num c.cst) c.coeffs, c.op), prf) + let of_cstr (c, prf) = ((Vect.set 0 (Q.neg c.cst) c.coeffs, c.op), prf) let product : t -> t -> t = fun ((p1, o1), prf1) ((p2, o2), prf2) -> @@ -897,7 +892,7 @@ module WithProof = struct | Eq -> ((LinPoly.product p p1, o1), ProofFormat.sMulC p prf1) | Gt | Ge -> let n, r = Vect.decomp_cst p in - if Vect.is_null r && n >/ Int 0 then + if Vect.is_null r && n >/ Q.zero then ((LinPoly.product p p1, o1), ProofFormat.mul_cst_proof n prf1) else ( if debug then @@ -908,34 +903,31 @@ module WithProof = struct let cutting_plane ((p, o), prf) = let c, p' = Vect.decomp_cst p in let g = Vect.gcd p' in - if - Big_int.eq_big_int Big_int.unit_big_int g - || c =/ Int 0 - || not (Big_int.eq_big_int (denominator c) Big_int.unit_big_int) - then None (* Nothing to do *) + if Z.equal Z.one g || c =/ Q.zero || not (Z.equal (Q.den c) Z.one) then None + (* Nothing to do *) else - let c1 = c // Big_int g in - let c1' = Num.floor_num c1 in + let c1 = c // Q.of_bigint g in + let c1' = Q.floor c1 in if c1 =/ c1' then None else match o with | Eq -> - Some ((Vect.set 0 (Int (-1)) Vect.null, Eq), ProofFormat.Gcd (g, prf)) + Some ((Vect.set 0 Q.neg_one Vect.null, Eq), ProofFormat.Gcd (g, prf)) | Gt -> failwith "cutting_plane ignore strict constraints" | Ge -> (* This is a non-trivial common divisor *) Some - ( (Vect.set 0 c1' (Vect.div (Big_int g) p), o) + ( (Vect.set 0 c1' (Vect.div (Q.of_bigint g) p), o) , ProofFormat.Gcd (g, prf) ) let construct_sign p = let c, p' = Vect.decomp_cst p in if Vect.is_null p' then Some - ( match sign_num c with + ( match Q.sign c with | 0 -> (true, Eq, ProofFormat.Zero) | 1 -> (true, Gt, ProofFormat.Cst c) - | _ (*-1*) -> (false, Gt, ProofFormat.Cst (minus_num c)) ) + | _ (*-1*) -> (false, Gt, ProofFormat.Cst (Q.neg c)) ) else None let get_sign l p = @@ -1007,7 +999,7 @@ module WithProof = struct | Some (c, p) -> Some (c, ProofFormat.simplify_prf_rule p) let is_substitution strict ((p, o), prf) = - let pred v = if strict then v =/ Int 1 || v =/ Int (-1) else true in + let pred v = if strict then v =/ Q.one || v =/ Q.neg_one else true in match o with Eq -> LinPoly.search_linear pred p | _ -> None let subst1 sys0 = @@ -1048,14 +1040,14 @@ module WithProof = struct , Some {cst = c2; var = v2; coeff = c2'} ) -> ( let good_coeff b o = match o with - | Eq -> Some (minus_num b) - | _ -> if b <=/ Int 0 then Some (minus_num b) else None + | Eq -> Some (Q.neg b) + | _ -> if b <=/ Q.zero then Some (Q.neg b) else None in match (good_coeff c1 o2, good_coeff c2 o1) with | None, _ | _, None -> None | Some c1, Some c2 -> let ext_mult c w = - if c =/ Int 0 then zero else mult (LinPoly.constant c) w + if c =/ Q.zero then zero else mult (LinPoly.constant c) w in Some (addition diff --git a/plugins/micromega/polynomial.mli b/plugins/micromega/polynomial.mli index 797ff5827d..357a2b10e1 100644 --- a/plugins/micromega/polynomial.mli +++ b/plugins/micromega/polynomial.mli @@ -9,6 +9,7 @@ (************************************************************************) open Mutils +open NumCompat module Mc = Micromega val max_nb_cstr : int ref @@ -81,7 +82,7 @@ module Poly : sig type t - val constant : Num.num -> t + val constant : Q.t -> t (** [constant c] @return the constant polynomial c *) @@ -101,24 +102,24 @@ module Poly : sig (** [uminus p] @return the polynomial -p i.e product by -1 *) - val get : Monomial.t -> t -> Num.num + val get : Monomial.t -> t -> Q.t (** [get mi p] @return the coefficient ai of the monomial mi. *) - val fold : (Monomial.t -> Num.num -> 'a -> 'a) -> t -> 'a -> 'a + val fold : (Monomial.t -> Q.t -> 'a -> 'a) -> t -> 'a -> 'a (** [fold f p a] folds f over the monomials of p with non-zero coefficient *) - val add : Monomial.t -> Num.num -> t -> t + val add : Monomial.t -> Q.t -> t -> t (** [add m n p] @return the polynomial n*m + p *) end -type cstr = {coeffs : Vect.t; op : op; cst : Num.num} +type cstr = {coeffs : Vect.t; op : op; cst : Q.t} (* Representation of linear constraints *) and op = Eq | Ge | Gt -val eval_op : op -> Num.num -> Num.num -> bool +val eval_op : op -> Q.t -> Q.t -> bool (*val opMult : op -> op -> op*) @@ -172,7 +173,7 @@ module LinPoly : sig @return 1.y where y is the variable index of the monomial x^1. *) - val coq_poly_of_linpol : (Num.num -> 'a) -> t -> 'a Mc.pExpr + val coq_poly_of_linpol : (Q.t -> 'a) -> t -> 'a Mc.pExpr (** [coq_poly_of_linpol c p] @param p is a multi-variate polynomial. @param c maps a rational to a Coq polynomial coefficient. @@ -206,7 +207,7 @@ module LinPoly : sig @return true if the polynomial is linear in x i.e can be written c*x+r where c is a constant and r is independent from x *) - val constant : Num.num -> t + val constant : Q.t -> t (** [constant c] @return the constant polynomial c *) @@ -216,9 +217,9 @@ module LinPoly : sig p is linear in x i.e x does not occur in b and a is a constant such that [pred a] *) - val search_linear : (Num.num -> bool) -> t -> var option + val search_linear : (Q.t -> bool) -> t -> var option - val search_all_linear : (Num.num -> bool) -> t -> var list + val search_all_linear : (Q.t -> bool) -> t -> var list (** [search_all_linear pred p] @return all the variables x such p = a.x + b such that p is linear in x i.e x does not occur in b and @@ -270,11 +271,11 @@ module ProofFormat : sig | Annot of string * prf_rule | Hyp of int | Def of int - | Cst of Num.num + | Cst of Q.t | Zero | Square of Vect.t | MulC of Vect.t * prf_rule - | Gcd of Big_int.big_int * prf_rule + | Gcd of Z.t * prf_rule | MulPrf of prf_rule * prf_rule | AddPrf of prf_rule * prf_rule | CutPrf of prf_rule @@ -287,20 +288,20 @@ module ProofFormat : sig (* x = z - t, z >= 0, t >= 0 *) - val pr_size : prf_rule -> Num.num + val pr_size : prf_rule -> Q.t val pr_rule_max_id : prf_rule -> int val proof_max_id : proof -> int val normalise_proof : int -> proof -> int * proof val output_prf_rule : out_channel -> prf_rule -> unit val output_proof : out_channel -> proof -> unit val add_proof : prf_rule -> prf_rule -> prf_rule - val mul_cst_proof : Num.num -> prf_rule -> prf_rule + val mul_cst_proof : Q.t -> prf_rule -> prf_rule val mul_proof : prf_rule -> prf_rule -> prf_rule val compile_proof : int list -> proof -> Micromega.zArithProof val cmpl_prf_rule : ('a Micromega.pExpr -> 'a Micromega.pol) - -> (Num.num -> 'a) + -> (Q.t -> 'a) -> int list -> prf_rule -> 'a Micromega.psatz @@ -332,7 +333,7 @@ module WithProof : sig val zero : t (** [zero] represents the tautology (0=0) *) - val const : Num.num -> t + val const : Q.t -> t (** [const n] represents the tautology (n>=0) *) val product : t -> t -> t diff --git a/plugins/micromega/simplex.ml b/plugins/micromega/simplex.ml index 54976221bc..15ab03964e 100644 --- a/plugins/micromega/simplex.ml +++ b/plugins/micromega/simplex.ml @@ -8,10 +8,9 @@ (* * (see LICENSE file for the text of the license) *) (************************************************************************) +open NumCompat +open Q.Notations open Polynomial -open Num - -(*open Util*) open Mutils type ('a, 'b) sum = Inl of 'a | Inr of 'b @@ -118,7 +117,7 @@ let output_vars o m = let unfeasible (rst : Restricted.t) tbl = Restricted.fold rst - (fun k v m -> if Vect.get_cst v >=/ Int 0 then m else IMap.add k () m) + (fun k v m -> if Vect.get_cst v >=/ Q.zero then m else IMap.add k () m) tbl IMap.empty let is_feasible rst tb = IMap.is_empty (unfeasible rst tb) @@ -138,7 +137,7 @@ let is_feasible rst tb = IMap.is_empty (unfeasible rst tb) let is_maximised_vect rst v = Vect.for_all (fun xi ai -> - if ai >/ Int 0 then false else Restricted.is_restricted xi rst) + if ai >/ Q.zero then false else Restricted.is_restricted xi rst) v (** [is_maximised rst v] @@ -161,11 +160,11 @@ let is_maximised rst v = *) type result = - | Max of num (** Maximum is reached *) + | Max of Q.t (** Maximum is reached *) | Ubnd of var (** Problem is unbounded *) | Feas (** Problem is feasible *) -type pivot = Done of result | Pivot of int * int * num +type pivot = Done of result | Pivot of int * int * Q.t type simplex = Opt of tableau * result (** For a row, x = ao.xo+...+ai.xi @@ -180,7 +179,7 @@ let rec find_pivot_column (rst : Restricted.t) (r : Vect.t) = match Vect.choose r with | None -> failwith "find_pivot_column" | Some (xi, ai, r') -> - if ai </ Int 0 then + if ai </ Q.zero then if Restricted.is_restricted xi rst then find_pivot_column rst r' (* ai.xi cannot be improved *) else (xi, -1) (* r is not restricted, sign of ai does not matter *) @@ -207,9 +206,9 @@ let find_pivot_row rst tbl j sgn = Restricted.fold rst (fun i' v res -> let aij = Vect.get j v in - if Int sgn */ aij </ Int 0 then + if Q.of_int sgn */ aij </ Q.zero then (* This would improve *) - let score' = Num.abs_num (Vect.get_cst v // aij) in + let score' = Q.abs (Vect.get_cst v // aij) in min_score res (i', score') else res) tbl None @@ -246,10 +245,10 @@ let find_pivot vr (rst : Restricted.t) tbl = let solve_column (c : var) (r : var) (e : Vect.t) : Vect.t = let a = Vect.get c e in - if a =/ Int 0 then failwith "Cannot solve column" + if a =/ Q.zero then failwith "Cannot solve column" else - let a' = Int (-1) // a in - Vect.mul a' (Vect.set r (Int (-1)) (Vect.set c (Int 0) e)) + let a' = Q.neg_one // a in + Vect.mul a' (Vect.set r Q.neg_one (Vect.set c Q.zero e)) (** [pivot_row r c e] @param c is such that c = e @@ -258,7 +257,7 @@ let solve_column (c : var) (r : var) (e : Vect.t) : Vect.t = let pivot_row (row : Vect.t) (c : var) (e : Vect.t) : Vect.t = let g = Vect.get c row in - if g =/ Int 0 then row else Vect.mul_add g e (Int 1) (Vect.set c (Int 0) row) + if g =/ Q.zero then row else Vect.mul_add g e Q.one (Vect.set c Q.zero row) let pivot_with (m : tableau) (v : var) (p : Vect.t) = IMap.map (fun (r : Vect.t) -> pivot_row r v p) m @@ -270,7 +269,7 @@ let pivot (m : tableau) (r : var) (c : var) = IMap.add c piv (pivot_with (IMap.remove r m) c piv) let adapt_unbounded vr x rst tbl = - if Vect.get_cst (IMap.find vr tbl) >=/ Int 0 then tbl else pivot tbl vr x + if Vect.get_cst (IMap.find vr tbl) >=/ Q.zero then tbl else pivot tbl vr x module BaseSet = Set.Make (struct type t = iset @@ -295,7 +294,7 @@ let simplex opt vr rst tbl = output_tableau stdout tbl; Printf.fprintf stdout "Error for variables %a\n" output_vars m end; - if (not opt) && Vect.get_cst (IMap.find vr tbl) >=/ Int 0 then + if (not opt) && Vect.get_cst (IMap.find vr tbl) >=/ Q.zero then Opt (tbl, Feas) else match find_pivot vr rst tbl with @@ -308,7 +307,7 @@ let simplex opt vr rst tbl = | Feas -> raise (Invalid_argument "find_pivot") ) | Pivot (i, j, s) -> if debug then begin - Printf.fprintf stdout "Find pivot for x%i(%s)\n" vr (string_of_num s); + Printf.fprintf stdout "Find pivot for x%i(%s)\n" vr (Q.to_string s); Printf.fprintf stdout "Leaving variable x%i\n" i; Printf.fprintf stdout "Entering variable x%i\n" j end; @@ -359,14 +358,13 @@ let push_real (opt : bool) (nw : var) (v : Vect.t) (rst : Restricted.t) | Feas -> Sat (t', None) | Max n -> if debug then begin - Printf.printf "The objective is maximised %s\n" (string_of_num n); + Printf.printf "The objective is maximised %s\n" (Q.to_string n); Printf.printf "%a = %a\n" LinPoly.pp_var nw pp_row (IMap.find nw t') end; - if n >=/ Int 0 then Sat (t', None) + if n >=/ Q.zero then Sat (t', None) else let v' = safe_find "push_real" nw t' in - Unsat - (Vect.set nw (Int 1) (Vect.set 0 (Int 0) (Vect.mul (Int (-1)) v'))) ) + Unsat (Vect.set nw Q.one (Vect.set 0 Q.zero (Vect.mul Q.neg_one v'))) ) open Mutils (** One complication is that equalities needs some pre-processing. @@ -381,7 +379,7 @@ let make_certificate vm l = (Vect.fold (fun acc x n -> let x', b = IMap.find x vm in - Vect.set x' (if b then n else Num.minus_num n) acc) + Vect.set x' (if b then n else Q.neg n) acc) Vect.null l) (** [eliminate_equalities vr0 l] @@ -397,11 +395,11 @@ let eliminate_equalities (vr0 : var) (l : Polynomial.cstr list) = | c :: l -> ( match c.op with | Ge -> - let v = Vect.set 0 (minus_num c.cst) c.coeffs in + let v = Vect.set 0 (Q.neg c.cst) c.coeffs in elim (idx + 1) (vr + 1) (IMap.add vr (idx, true) vm) l ((vr, v) :: acc) | Eq -> - let v1 = Vect.set 0 (minus_num c.cst) c.coeffs in - let v2 = Vect.mul (Int (-1)) v1 in + let v1 = Vect.set 0 (Q.neg c.cst) c.coeffs in + let v2 = Vect.mul Q.neg_one v1 in let vm = IMap.add vr (idx, true) (IMap.add (vr + 1) (idx, false) vm) in elim (idx + 1) (vr + 2) vm l ((vr, v1) :: (vr + 1, v2) :: acc) | Gt -> raise Strict ) @@ -419,7 +417,7 @@ let find_full_solution rst tbl = IMap.fold (fun vr v res -> Vect.set vr (Vect.get_cst v) res) tbl Vect.null let choose_conflict (sol : Vect.t) (l : (var * Vect.t) list) = - let esol = Vect.set 0 (Int 1) sol in + let esol = Vect.set 0 Q.one sol in let rec most_violating l e (x, v) rst = match l with | [] -> Some ((x, v), rst) @@ -476,7 +474,7 @@ let optimise obj l = let _, vm, l' = eliminate_equalities (vr0 + 1) l in let bound pos res = match res with - | Opt (_, Max n) -> Some (if pos then n else minus_num n) + | Opt (_, Max n) -> Some (if pos then n else Q.neg n) | Opt (_, Ubnd _) -> None | Opt (_, Feas) -> None in @@ -501,9 +499,7 @@ let make_farkas_certificate (env : WithProof.t IMap.t) vm v = begin try let x', b = IMap.find x vm in - mul_cst_proof - (if b then n else Num.minus_num n) - (snd (IMap.find x' env)) + mul_cst_proof (if b then n else Q.neg n) (snd (IMap.find x' env)) with Not_found -> (* This is an introduced hypothesis *) mul_cst_proof n (snd (IMap.find x env)) @@ -517,7 +513,7 @@ let make_farkas_proof (env : WithProof.t IMap.t) vm v = begin try let x', b = IMap.find x vm in - let n = if b then n else Num.minus_num n in + let n = if b then n else Q.neg n in let prf = IMap.find x' env in WithProof.mult (Vect.cst n) prf with Not_found -> @@ -526,7 +522,7 @@ let make_farkas_proof (env : WithProof.t IMap.t) vm v = end) WithProof.zero v -let frac_num n = n -/ Num.floor_num n +let frac_num n = n -/ Q.floor n type ('a, 'b) hitkind = | Forget @@ -538,38 +534,38 @@ type ('a, 'b) hitkind = let cut env rmin sol vm (rst : Restricted.t) tbl (x, v) = let n, r = Vect.decomp_cst v in let fn = frac_num n in - if fn =/ Int 0 then Forget (* The solution is integral *) + if fn =/ Q.zero then Forget (* The solution is integral *) else (* The cut construction is from: Letchford and Lodi. Strengthening Chvatal-Gomory cuts and Gomory fractional cuts. We implement the classic Proposition 2 from the "known results" - *) + *) (* Proposition 3 requires all the variables to be restricted and is therefore not always applicable. *) (* let ccoeff_prop1 v = frac_num v in - let ccoeff_prop3 v = - (* mixed integer cut *) - let fv = frac_num v in - Num.min_num fv (fn */ (Int 1 -/ fv) // (Int 1 -/ fn)) - in - let ccoeff_prop3 = - if Restricted.is_restricted x rst then ("Prop3", ccoeff_prop3) - else ("Prop1", ccoeff_prop1) - in *) - let n0_5 = Int 1 // Int 2 in + let ccoeff_prop3 v = + (* mixed integer cut *) + let fv = frac_num v in + Num.min_num fv (fn */ (Q.one -/ fv) // (Q.one -/ fn)) + in + let ccoeff_prop3 = + if Restricted.is_restricted x rst then ("Prop3", ccoeff_prop3) + else ("Prop1", ccoeff_prop1) + in *) + let n0_5 = Q.one // Q.two in (* If the fractional part [fn] is small, we construct the t-cut. If the fractional part [fn] is big, we construct the t-cut of the negated row. (This is only a cut if all the fractional variables are restricted.) - *) + *) let ccoeff_prop2 = let tmin = if fn </ n0_5 then (* t-cut *) - Num.ceiling_num (n0_5 // fn) + Q.ceiling (n0_5 // fn) else (* multiply by -1 & t-cut *) - minus_num (Num.ceiling_num (n0_5 // (Int 1 -/ fn))) + Q.neg (Q.ceiling (n0_5 // (Q.one -/ fn))) in ("Prop2", fun v -> frac_num (v */ tmin)) in @@ -651,7 +647,7 @@ let eliminate_variable (bounded, vr, env, tbl) x = let tv = var_of_vect t in (* x = z - t *) let xdef = Vect.add z (Vect.uminus t) in - let xp = ((Vect.set x (Int 1) (Vect.uminus xdef), Eq), Def vr) in + let xp = ((Vect.set x Q.one (Vect.uminus xdef), Eq), Def vr) in let zp = ((z, Ge), Def zv) in let tp = ((t, Ge), Def tv) in (* Pivot the current tableau using xdef *) @@ -662,11 +658,8 @@ let eliminate_variable (bounded, vr, env, tbl) x = (fun lp -> let (v, o), p = lp in let ai = Vect.get x v in - if ai =/ Int 0 then lp - else - WithProof.addition - (WithProof.mult (Vect.cst (Num.minus_num ai)) xp) - lp) + if ai =/ Q.zero then lp + else WithProof.addition (WithProof.mult (Vect.cst (Q.neg ai)) xp) lp) env in (* Add the variables to the environment *) diff --git a/plugins/micromega/simplex.mli b/plugins/micromega/simplex.mli index ff672edafd..8edea2d4b2 100644 --- a/plugins/micromega/simplex.mli +++ b/plugins/micromega/simplex.mli @@ -7,6 +7,8 @@ (* * GNU Lesser General Public License Version 2.1 *) (* * (see LICENSE file for the text of the license) *) (************************************************************************) + +open NumCompat open Polynomial (** Profiling *) @@ -23,7 +25,7 @@ val get_profile_info : unit -> profile_info (** Simplex interface *) -val optimise : Vect.t -> cstr list -> (Num.num option * Num.num option) option +val optimise : Vect.t -> cstr list -> (Q.t option * Q.t option) option val find_point : cstr list -> Vect.t option val find_unsat_certificate : cstr list -> Vect.t option diff --git a/plugins/micromega/sos.ml b/plugins/micromega/sos.ml index 772ed7a8c5..2b04bb80e2 100644 --- a/plugins/micromega/sos.ml +++ b/plugins/micromega/sos.ml @@ -9,7 +9,9 @@ (* ========================================================================= *) (* Nonlinear universal reals procedure using SOS decomposition. *) (* ========================================================================= *) -open Num + +open NumCompat +open Q.Notations open Sos_types open Sos_lib @@ -27,19 +29,19 @@ exception Sanity 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 + if Q.abs y </ Q.one // Q.ten then normalize (Q.ten */ y) - 1 + else if Q.abs y >=/ Q.one then normalize (y // Q.ten) + 1 else 0 in fun d x -> - if x =/ Int 0 then "0.0" + if x =/ Q.zero then "0.0" else - let y = abs_num x in + let y = Q.abs 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 (List.tl (explode (string_of_num k))) + let z = (Q.pow10 (-e) */ y) +/ Q.one in + let k = Q.round (Q.pow10 d */ z) in + (if x </ Q.zero then "-0." else "0.") + ^ implode (List.tl (explode (Q.to_string k))) ^ if e = 0 then "" else "e" ^ string_of_int e (* ------------------------------------------------------------------------- *) @@ -55,22 +57,22 @@ 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 vector = int * (int, Q.t) func +type matrix = (int * int) * (int * int, Q.t) func type monomial = (vname, int) func -type poly = (monomial, num) func +type poly = (monomial, Q.t) func (* ------------------------------------------------------------------------- *) (* Assignment avoiding zeros. *) (* ------------------------------------------------------------------------- *) -let ( |--> ) x y a = if y =/ Int 0 then a else (x |-> y) a +let ( |--> ) x y a = if y =/ Q.zero then a else (x |-> y) a (* ------------------------------------------------------------------------- *) (* This can be generic. *) (* ------------------------------------------------------------------------- *) -let element (d, v) i = tryapplyd v i (Int 0) +let element (d, v) i = tryapplyd v i Q.zero 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 @@ -82,12 +84,12 @@ 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 + if c =/ Q.zero then vector_0 n else ((n, List.fold_right (fun k -> k |-> c) (1 -- n) undefined) : vector) 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)) + if c =/ Q.zero then vector_0 n else (n, mapf (fun x -> c */ x) (snd v)) let vector_of_list l = let n = List.length l in @@ -102,15 +104,15 @@ let dimensions (m : matrix) = fst m let matrix_cmul c (m : matrix) = let i, j = dimensions m in - if c =/ Int 0 then matrix_0 (i, j) + if c =/ Q.zero 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_neg (m : matrix) = ((dimensions m, mapf Q.neg (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) + else ((d1, combine ( +/ ) (fun x -> x =/ Q.zero) (snd m1) (snd m2)) : matrix) let row k (m : matrix) = let i, j = dimensions m in @@ -150,21 +152,21 @@ let monomial_variables m = dom m (* ------------------------------------------------------------------------- *) 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_var x = (monomial_var x |=> Q.one : poly) +let poly_const c = if c =/ Q.zero 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 + if c =/ Q.zero then poly_0 else mapf (fun x -> c */ x) p -let poly_neg (p : poly) = (mapf minus_num p : poly) +let poly_neg (p : poly) = (mapf Q.neg p : poly) let poly_add (p1 : poly) (p2 : poly) = - (combine ( +/ ) (fun x -> x =/ Int 0) p1 p2 : poly) + (combine ( +/ ) (fun x -> x =/ Q.zero) 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 + if c =/ Q.zero 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 @@ -174,7 +176,7 @@ let poly_mul (p1 : poly) (p2 : poly) = let poly_square p = poly_mul p p let rec poly_pow p k = - if k = 0 then poly_const (Int 1) + if k = 0 then poly_const Q.one else if k = 1 then p else let q = poly_square (poly_pow p (k / 2)) in @@ -228,9 +230,9 @@ let string_of_monomial m = String.concat "*" 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 + if m = monomial_1 then Q.to_string c + else if c =/ Q.one then string_of_monomial m + else Q.to_string c ^ "*" ^ string_of_monomial m let string_of_poly (p : poly) = if p = poly_0 then "<<0>>" @@ -241,7 +243,7 @@ let string_of_poly (p : poly) = let s = List.fold_left (fun a (m, c) -> - if c </ Int 0 then a ^ " - " ^ string_of_cmonomial (minus_num c, m) + if c </ Q.zero then a ^ " - " ^ string_of_cmonomial (Q.neg c, m) else a ^ " + " ^ string_of_cmonomial (c, m)) "" cms in @@ -338,21 +340,19 @@ let token s = let decimal = let ( || ) = parser_or in let numeral = some isnum in - let decimalint = atleast 1 numeral >> o Num.num_of_string implode in + let decimalint = atleast 1 numeral >> o Q.of_string implode in let decimalfrac = atleast 1 numeral - >> fun s -> Num.num_of_string (implode s) // pow10 (List.length s) + >> fun s -> Q.of_string (implode s) // Q.pow10 (List.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 signed prs = a "-" ++ prs >> o Q.neg 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 + >> function h, [x] -> h */ Q.power 10 x | h, _ -> h let mkparser p s = let x, rst = p (explode s) in @@ -469,19 +469,19 @@ let run_csdp dbg obj mats = let scale_then = let common_denominator amat acc = - foldl (fun a m c -> lcm_num (denominator c) a) acc amat + foldl (fun a m c -> Z.lcm (Q.den c) a) acc amat and maximal_element amat acc = - foldl (fun maxa m c -> max_num maxa (abs_num c)) acc amat + foldl (fun maxa m c -> Q.max maxa (Q.abs c)) acc amat in fun solver obj mats -> - let cd1 = List.fold_right common_denominator mats (Int 1) - and cd2 = common_denominator (snd obj) (Int 1) in + let cd1 = Q.of_bigint @@ List.fold_right common_denominator mats Z.one + and cd2 = Q.of_bigint @@ common_denominator (snd obj) Z.one in let mats' = List.map (mapf (fun x -> cd1 */ x)) mats and obj' = vector_cmul cd2 obj in - let max1 = List.fold_right 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 max1 = List.fold_right maximal_element mats' Q.zero + and max2 = maximal_element (snd obj') Q.zero in + let scal1 = Q.pow2 (20 - int_of_float (log (Q.to_float max1) /. log 2.0)) + and scal2 = Q.pow2 (20 - int_of_float (log (Q.to_float max2) /. log 2.0)) in let mats'' = List.map (mapf (fun x -> x */ scal1)) mats' and obj'' = vector_cmul scal2 obj' in solver obj'' mats'' @@ -490,7 +490,7 @@ let scale_then = (* Round a vector to "nice" rationals. *) (* ------------------------------------------------------------------------- *) -let nice_rational n x = round_num (n */ x) // n +let nice_rational n x = Q.round (n */ x) // n let nice_vector n = mapa (nice_rational n) (* ------------------------------------------------------------------------- *) @@ -501,7 +501,7 @@ let nice_vector n = mapa (nice_rational n) let linear_program_basic a = let m, n = dimensions a in let mats = List.map (fun j -> diagonal (column j a)) (1 -- n) - and obj = vector_const (Int 1) m in + and obj = vector_const Q.one m in let rv, res = run_csdp false obj mats in if rv = 1 || rv = 2 then false else if rv = 0 then true @@ -521,8 +521,8 @@ let in_convex_hull pts pt = 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) ) + (fun pts j -> itern 1 pts (fun x i -> (i, j) |-> Q.of_int x)) + (iter (1, n) (fun i -> (v + i, i + 1) |-> Q.one) undefined) ) in linear_program_basic mat @@ -544,12 +544,14 @@ let minimal_convex_hull = (* 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_cmul c eq = + if c =/ Q.zero then Empty else mapf (fun d -> c */ d) eq + +let equation_add eq1 eq2 = combine ( +/ ) (fun x -> x =/ Q.zero) 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 + foldl (fun a v c -> a +/ (value v */ c)) Q.zero eq (* ------------------------------------------------------------------------- *) (* Eliminate all variables, in an essentially arbitrary order. *) @@ -574,11 +576,11 @@ let eliminate_all_equations one = 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 eq' = equation_cmul (Q.neg_one // 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) + let b = tryapplyd e v Q.zero in + if b =/ Q.zero then e + else equation_add e (equation_cmul (Q.neg b // a) eq) in eliminate ((v |-> eq') (mapf elim dun)) (List.map elim oeqs) in @@ -631,8 +633,8 @@ let diag 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 a11 </ Q.zero then failwith "diagonalize: not PSD" + else if a11 =/ Q.zero then if is_zero (row i m) then diagonalize (i + 1) m else failwith "diagonalize: not PSD" else @@ -659,21 +661,23 @@ let diag m = (* ------------------------------------------------------------------------- *) let deration d = - if d = [] then (Int 0, d) + if d = [] then (Q.zero, 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) + Q.make + (foldl (fun a i c -> Z.lcm a (Q.den c)) Z.one (snd l)) + (foldl (fun a i c -> Z.gcd a (Q.num c)) Z.zero (snd l)) in (c // (a */ a), mapa (fun x -> a */ x) l) in let d' = List.map adj d in let a = - List.fold_right (o lcm_num (o denominator fst)) d' (Int 1) - // List.fold_right (o gcd_num (o numerator fst)) d' (Int 0) + Q.make + (List.fold_right (o Z.lcm (o Q.den fst)) d' Z.one) + (List.fold_right (o Z.gcd (o Q.num fst)) d' Z.zero) in - (Int 1 // a, List.map (fun (c, l) -> (a */ c, l)) d') + (Q.one // a, List.map (fun (c, l) -> (a */ c, l)) d') (* ------------------------------------------------------------------------- *) (* Enumeration of monomials with given multidegree bound. *) @@ -702,11 +706,11 @@ let rec enumerate_monomials d vars = (* ------------------------------------------------------------------------- *) let rec enumerate_products d pols = - if d = 0 then [(poly_const num_1, Rational_lt num_1)] + if d = 0 then [(poly_const Q.one, Rational_lt Q.one)] else if d < 0 then [] else match pols with - | [] -> [(poly_const num_1, Rational_lt num_1)] + | [] -> [(poly_const Q.one, Rational_lt Q.one)] | (p, b) :: ps -> let e = multidegree p in if e = 0 then enumerate_products d ps @@ -736,7 +740,7 @@ let epoly_pmul p q acc = (* ------------------------------------------------------------------------- *) let epoly_of_poly p = - foldl (fun a m c -> (m |-> ((0, 0, 0) |=> minus_num c)) a) undefined p + foldl (fun a m c -> (m |-> ((0, 0, 0) |=> Q.neg c)) a) undefined p (* ------------------------------------------------------------------------- *) (* String for block diagonal matrix numbered k. *) @@ -796,7 +800,7 @@ let csdp nblocks blocksizes obj mats = if rv = 1 || rv = 2 then failwith "csdp: Problem is infeasible" else if rv = 3 then () (*Format.print_string "csdp warning: Reduced accuracy"; - Format.print_newline() *) + Format.print_newline() *) else if rv <> 0 then failwith ("csdp: error " ^ string_of_int rv) else (); res @@ -805,12 +809,12 @@ let csdp nblocks blocksizes obj mats = (* 3D versions of matrix operations to consider blocks separately. *) (* ------------------------------------------------------------------------- *) -let bmatrix_add = combine ( +/ ) (fun x -> x =/ Int 0) +let bmatrix_add = combine ( +/ ) (fun x -> x =/ Q.zero) let bmatrix_cmul c bm = - if c =/ Int 0 then undefined else mapf (fun x -> c */ x) bm + if c =/ Q.zero then undefined else mapf (fun x -> c */ x) bm -let bmatrix_neg = bmatrix_cmul (Int (-1)) +let bmatrix_neg = bmatrix_cmul Q.neg_one (* ------------------------------------------------------------------------- *) (* Smash a block matrix into components. *) @@ -839,7 +843,7 @@ let real_positivnullstellensatz_general linf d eqs leqs pol = in let monoid = if linf then - (poly_const num_1, Rational_lt num_1) + (poly_const Q.one, Rational_lt Q.one) :: List.filter (fun (p, c) -> multidegree p <= d) leqs else enumerate_products d leqs in @@ -850,7 +854,7 @@ let real_positivnullstellensatz_general linf d eqs leqs pol = let nons = List.combine mons (1 -- List.length mons) in ( mons , List.fold_right - (fun (m, n) -> m |-> ((-k, -n, n) |=> Int 1)) + (fun (m, n) -> m |-> ((-k, -n, n) |=> Q.one)) nons undefined ) in let mk_sqmultiplier k (p, c) = @@ -865,7 +869,7 @@ let real_positivnullstellensatz_general linf d eqs leqs pol = 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 c = if n1 = n2 then Q.one else Q.two in let e = tryapplyd a m undefined in (m |-> equation_add ((k, n1, n2) |=> c) e) a) nons) @@ -889,14 +893,14 @@ let real_positivnullstellensatz_general linf d eqs leqs pol = 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 = List.fold_right (fun v -> v |-> (v |=> Int 1)) pvs assig in + let allassig = List.fold_right (fun v -> v |-> (v |=> Q.one)) 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)) + let c = tryapplyd ass v Q.zero in + if c =/ Q.zero then m else ((b, j, i) |-> c) (((b, i, j) |-> c) m)) undefined allassig in let diagents = @@ -907,7 +911,7 @@ let real_positivnullstellensatz_general linf d eqs leqs pol = let mats = List.map mk_matrix qvars and obj = ( List.length pvs - , itern 1 pvs (fun v i -> i |--> tryapplyd diagents v (Int 0)) undefined ) + , itern 1 pvs (fun v i -> i |--> tryapplyd diagents v Q.zero) undefined ) in let raw_vec = if pvs = [] then vector_0 0 @@ -915,7 +919,7 @@ let real_positivnullstellensatz_general linf d eqs leqs pol = in let find_rounding d = if !debugging then ( - Format.print_string ("Trying rounding with limit " ^ string_of_num d); + Format.print_string ("Trying rounding with limit " ^ Q.to_string d); Format.print_newline () ) else (); let vec = nice_vector d raw_vec in @@ -930,16 +934,16 @@ let real_positivnullstellensatz_general linf d eqs leqs pol = (vec, List.map diag allmats) in let vec, ratdias = - if pvs = [] then find_rounding num_1 + if pvs = [] then find_rounding Q.one else tryfind find_rounding - (List.map Num.num_of_int (1 -- 31) @ List.map pow2 (5 -- 66)) + (List.map Q.of_int (1 -- 31) @ List.map Q.pow2 (5 -- 66)) in let newassigs = List.fold_right (fun k -> List.nth pvs (k - 1) |-> element vec k) (1 -- dim vec) - ((0, 0, 0) |=> Int (-1)) + ((0, 0, 0) |=> Q.neg_one) in let finalassigs = foldl (fun a v e -> (v |-> equation_eval newassigs e) a) newassigs allassig @@ -1017,7 +1021,7 @@ let monomial_order = let term_of_varpow x k = if k = 1 then Var x else Pow (Var x, k) let term_of_monomial m = - if m = monomial_1 then Const num_1 + if m = monomial_1 then Const Q.one else let m' = dest_monomial m in let vps = List.fold_right (fun (x, k) a -> term_of_varpow x k :: a) m' [] in @@ -1025,7 +1029,7 @@ let term_of_monomial m = let term_of_cmonomial (m, c) = if m = monomial_1 then Const c - else if c =/ num_1 then term_of_monomial m + else if c =/ Q.one then term_of_monomial m else Mul (Const c, term_of_monomial m) let term_of_poly p = @@ -1114,8 +1118,8 @@ let csdp obj mats = let rv, res = run_csdp !debugging obj mats in if rv = 1 || rv = 2 then failwith "csdp: Problem is infeasible" else if rv = 3 then () - (* (Format.print_string "csdp warning: Reduced accuracy"; - Format.print_newline()) *) + (* (Format.print_string "csdp warning: Reduced accuracy"; + Format.print_newline()) *) else if rv <> 0 then failwith ("csdp: error " ^ string_of_int rv) else (); res @@ -1162,7 +1166,7 @@ let sumofsquares_general_symmetry tool pol = match cls with | [] -> raise Sanity | [h] -> acc - | h :: t -> List.map (fun k -> (k |-> Int (-1)) (h |=> Int 1)) t @ acc + | h :: t -> List.map (fun k -> (k |-> Q.neg_one) (h |=> Q.one)) t @ acc in List.fold_right mk_eq eqvcls [] in @@ -1176,13 +1180,13 @@ let sumofsquares_general_symmetry tool pol = 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 + let c = if n1 = n2 then Q.one else Q.two 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 = List.fold_right (fun v -> v |-> (v |=> Int 1)) pvs assig in + let allassig = List.fold_right (fun v -> v |-> (v |=> Q.one)) pvs assig in let qvars = (0, 0) :: pvs in let diagents = end_itlist equation_add (List.map (fun i -> apply allassig (i, i)) (1 -- n)) @@ -1191,20 +1195,20 @@ let sumofsquares_general_symmetry tool pol = ( ( (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)) + let c = tryapplyd ass v Q.zero in + if c =/ Q.zero then m else ((j, i) |-> c) (((i, j) |-> c) m)) undefined allassig ) : matrix ) in let mats = List.map mk_matrix qvars and obj = ( List.length pvs - , itern 1 pvs (fun v i -> i |--> tryapplyd diagents v (Int 0)) undefined ) + , itern 1 pvs (fun v i -> i |--> tryapplyd diagents v Q.zero) 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_string ("Trying rounding with limit " ^ Q.to_string d); Format.print_newline () ) else (); let vec = nice_vector d raw_vec in @@ -1223,7 +1227,7 @@ let sumofsquares_general_symmetry tool pol = deration (diag mat) else tryfind find_rounding - (List.map Num.num_of_int (1 -- 31) @ List.map pow2 (5 -- 66)) + (List.map Q.of_int (1 -- 31) @ List.map Q.pow2 (5 -- 66)) in let poly_of_lin (d, v) = (d, foldl (fun a i c -> (List.nth lpps (i - 1) |-> c) a) undefined (snd v)) diff --git a/plugins/micromega/sos.mli b/plugins/micromega/sos.mli index ac75bd37f0..8a461b4c20 100644 --- a/plugins/micromega/sos.mli +++ b/plugins/micromega/sos.mli @@ -8,6 +8,7 @@ (* * (see LICENSE file for the text of the license) *) (************************************************************************) +open NumCompat open Sos_types type poly @@ -16,13 +17,10 @@ val poly_isconst : poly -> bool val poly_neg : poly -> poly val poly_mul : poly -> poly -> poly val poly_pow : poly -> int -> poly -val poly_const : Num.num -> poly +val poly_const : Q.t -> poly val poly_of_term : term -> poly val term_of_poly : poly -> term - -val term_of_sos : - positivstellensatz * (Num.num * poly) list -> positivstellensatz - +val term_of_sos : positivstellensatz * (Q.t * poly) list -> positivstellensatz val string_of_poly : poly -> string val real_positivnullstellensatz_general : @@ -31,6 +29,6 @@ val real_positivnullstellensatz_general : -> poly list -> (poly * positivstellensatz) list -> poly - -> poly list * (positivstellensatz * (Num.num * poly) list) list + -> poly list * (positivstellensatz * (Q.t * poly) list) list -val sumofsquares : poly -> Num.num * (Num.num * poly) list +val sumofsquares : poly -> Q.t * (Q.t * poly) list diff --git a/plugins/micromega/sos_lib.ml b/plugins/micromega/sos_lib.ml index 51221aa6b9..99c552e379 100644 --- a/plugins/micromega/sos_lib.ml +++ b/plugins/micromega/sos_lib.ml @@ -7,8 +7,6 @@ (* - Frédéric Besson (fbesson@irisa.fr) is using it to feed micromega *) (* ========================================================================= *) -open Num - (* ------------------------------------------------------------------------- *) (* Comparisons that are reflexive on NaN and also short-circuiting. *) (* ------------------------------------------------------------------------- *) @@ -28,32 +26,6 @@ let ( >? ) x y = cmp x y > 0 let o f g x = f (g x) (* ------------------------------------------------------------------------- *) -(* Some useful functions on "num" type. *) -(* ------------------------------------------------------------------------- *) - -let num_0 = Int 0 -and num_1 = Int 1 -and num_2 = Int 2 -and num_10 = Int 10 - -let pow2 n = power_num num_2 (Int n) -let pow10 n = power_num num_10 (Int n) - -let numdom r = - let r' = Ratio.normalize_ratio (ratio_of_num r) in - ( num_of_big_int (Ratio.numerator_ratio r') - , num_of_big_int (Ratio.denominator_ratio r') ) - -let numerator = o fst numdom -and denominator = o snd numdom - -let gcd_num n1 n2 = - num_of_big_int (Big_int.gcd_big_int (big_int_of_num n1) (big_int_of_num n2)) - -let lcm_num x y = - if x =/ num_0 && y =/ num_0 then num_0 else abs_num (x */ y // gcd_num x y) - -(* ------------------------------------------------------------------------- *) (* Various versions of list iteration. *) (* ------------------------------------------------------------------------- *) @@ -518,8 +490,8 @@ let deepen_until limit f n = let rec d_until f n = try (* if !debugging - then (print_string "Searching with depth limit "; - print_int n; print_newline()) ;*) + then (print_string "Searching with depth limit "; + print_int n; print_newline()) ;*) f n with Failure x -> (*if !debugging then (Printf.printf "solver error : %s\n" x) ; *) diff --git a/plugins/micromega/sos_lib.mli b/plugins/micromega/sos_lib.mli index 2bbcbf336b..7795808e12 100644 --- a/plugins/micromega/sos_lib.mli +++ b/plugins/micromega/sos_lib.mli @@ -9,9 +9,6 @@ (************************************************************************) val o : ('a -> 'b) -> ('c -> 'a) -> 'c -> 'b -val num_1 : Num.num -val pow10 : int -> Num.num -val pow2 : int -> Num.num val implode : string list -> string val explode : string -> string list val funpow : int -> ('a -> 'a) -> 'a -> 'a @@ -50,10 +47,6 @@ val sort : ('a -> 'a -> bool) -> 'a list -> 'a list val setify : 'a list -> 'a list val increasing : ('a -> 'b) -> 'a -> 'a -> bool val allpairs : ('a -> 'b -> 'c) -> 'a list -> 'b list -> 'c list -val gcd_num : Num.num -> Num.num -> Num.num -val lcm_num : Num.num -> Num.num -> Num.num -val numerator : Num.num -> Num.num -val denominator : Num.num -> Num.num val end_itlist : ('a -> 'a -> 'a) -> 'a list -> 'a val ( >> ) : ('a -> 'b * 'c) -> ('b -> 'd) -> 'a -> 'd * 'c val ( ++ ) : ('a -> 'b * 'c) -> ('c -> 'd * 'e) -> 'a -> ('b * 'd) * 'e diff --git a/plugins/micromega/sos_types.ml b/plugins/micromega/sos_types.ml index 988024968b..62699d8362 100644 --- a/plugins/micromega/sos_types.ml +++ b/plugins/micromega/sos_types.ml @@ -9,13 +9,13 @@ (************************************************************************) (* The type of positivstellensatz -- used to communicate with sos *) -open Num - type vname = string +open NumCompat + type term = | Zero - | Const of Num.num + | Const of Q.t | Var of vname | Opp of term | Add of (term * term) @@ -26,7 +26,7 @@ type term = let rec output_term o t = match t with | Zero -> output_string o "0" - | Const n -> output_string o (string_of_num n) + | Const n -> output_string o (Q.to_string n) | Var n -> Printf.fprintf o "v%s" n | Opp t -> Printf.fprintf o "- (%a)" output_term t | Add (t1, t2) -> Printf.fprintf o "(%a)+(%a)" output_term t1 output_term t2 @@ -42,9 +42,9 @@ type positivstellensatz = | Axiom_eq of int | Axiom_le of int | Axiom_lt of int - | Rational_eq of num - | Rational_le of num - | Rational_lt of num + | Rational_eq of Q.t + | Rational_le of Q.t + | Rational_lt of Q.t | Square of term | Monoid of int list | Eqmul of term * positivstellensatz @@ -55,9 +55,9 @@ let rec output_psatz o = function | Axiom_eq i -> Printf.fprintf o "Aeq(%i)" i | Axiom_le i -> Printf.fprintf o "Ale(%i)" i | Axiom_lt i -> Printf.fprintf o "Alt(%i)" i - | Rational_eq n -> Printf.fprintf o "eq(%s)" (string_of_num n) - | Rational_le n -> Printf.fprintf o "le(%s)" (string_of_num n) - | Rational_lt n -> Printf.fprintf o "lt(%s)" (string_of_num n) + | Rational_eq n -> Printf.fprintf o "eq(%s)" (Q.to_string n) + | Rational_le n -> Printf.fprintf o "le(%s)" (Q.to_string n) + | Rational_lt n -> Printf.fprintf o "lt(%s)" (Q.to_string n) | Square t -> Printf.fprintf o "(%a)^2" output_term t | Monoid l -> Printf.fprintf o "monoid" | Eqmul (t, ps) -> Printf.fprintf o "%a * %a" output_term t output_psatz ps diff --git a/plugins/micromega/sos_types.mli b/plugins/micromega/sos_types.mli index ca9a43b1d0..a0b9157880 100644 --- a/plugins/micromega/sos_types.mli +++ b/plugins/micromega/sos_types.mli @@ -8,13 +8,15 @@ (* * (see LICENSE file for the text of the license) *) (************************************************************************) +open NumCompat + (* The type of positivstellensatz -- used to communicate with sos *) type vname = string type term = | Zero - | Const of Num.num + | Const of Q.t | Var of vname | Opp of term | Add of (term * term) @@ -28,9 +30,9 @@ type positivstellensatz = | Axiom_eq of int | Axiom_le of int | Axiom_lt of int - | Rational_eq of Num.num - | Rational_le of Num.num - | Rational_lt of Num.num + | Rational_eq of Q.t + | Rational_le of Q.t + | Rational_lt of Q.t | Square of term | Monoid of int list | Eqmul of term * positivstellensatz diff --git a/plugins/micromega/vect.ml b/plugins/micromega/vect.ml index f53a7b42c9..198430295b 100644 --- a/plugins/micromega/vect.ml +++ b/plugins/micromega/vect.ml @@ -8,7 +8,8 @@ (* * (see LICENSE file for the text of the license) *) (************************************************************************) -open Num +open NumCompat +open Q.Notations open Mutils type var = int @@ -18,7 +19,7 @@ type var = int - values are all non-zero *) -type t = (var * num) list +type t = (var * Q.t) list type vector = t (** [equal v1 v2 = true] if the vectors are syntactically equal. *) @@ -33,32 +34,30 @@ let rec equal v1 v2 = let hash v = let rec hash i = function | [] -> i - | (vr, vl) :: l -> hash (i + Hashtbl.hash (vr, float_of_num vl)) l + | (vr, vl) :: l -> hash (i + Hashtbl.hash (vr, Q.to_float vl)) l in Hashtbl.hash (hash 0 v) let null = [] -let is_null v = match v with [] | [(0, Int 0)] -> true | _ -> false + +let is_null v = + match v with [] -> true | [(0, x)] when Q.zero =/ x -> true | _ -> false let pp_var_num pp_var o (v, n) = if Int.equal v 0 then - if eq_num (Int 0) n then () else Printf.fprintf o "%s" (string_of_num n) - else - match n with - | Int 1 -> pp_var o v - | Int -1 -> Printf.fprintf o "-%a" pp_var v - | Int 0 -> () - | _ -> Printf.fprintf o "%s*%a" (string_of_num n) pp_var v + if Q.zero =/ n then () else Printf.fprintf o "%s" (Q.to_string n) + else if Q.one =/ n then pp_var o v + else if Q.neg_one =/ n then Printf.fprintf o "-%a" pp_var v + else if Q.zero =/ n then () + else Printf.fprintf o "%s*%a" (Q.to_string n) pp_var v let pp_var_num_smt pp_var o (v, n) = if Int.equal v 0 then - if eq_num (Int 0) n then () else Printf.fprintf o "%s" (string_of_num n) - else - match n with - | Int 1 -> pp_var o v - | Int -1 -> Printf.fprintf o "(- %a)" pp_var v - | Int 0 -> () - | _ -> Printf.fprintf o "(* %s %a)" (string_of_num n) pp_var v + if Q.zero =/ n then () else Printf.fprintf o "%s" (Q.to_string n) + else if Q.one =/ n then pp_var o v + else if Q.neg_one =/ n then Printf.fprintf o "(- %a)" pp_var v + else if Q.zero =/ n then () + else Printf.fprintf o "(* %s %a)" (Q.to_string n) pp_var v let rec pp_gen pp_var o v = match v with @@ -75,36 +74,34 @@ let pp_smt o v = in Printf.fprintf o "(+ %a)" list v -let from_list (l : num list) = +let from_list (l : Q.t list) = let rec xfrom_list i l = match l with | [] -> [] | e :: l -> - if e <>/ Int 0 then (i, e) :: xfrom_list (i + 1) l + if e <>/ Q.zero then (i, e) :: xfrom_list (i + 1) l else xfrom_list (i + 1) l in xfrom_list 0 l -let zero_num = Int 0 - let to_list m = let rec xto_list i l = match l with | [] -> [] | (x, v) :: l' -> - if i = x then v :: xto_list (i + 1) l' else zero_num :: xto_list (i + 1) l + if i = x then v :: xto_list (i + 1) l' else Q.zero :: xto_list (i + 1) l in xto_list 0 m -let cons i v rst = if v =/ Int 0 then rst else (i, v) :: rst +let cons i v rst = if v =/ Q.zero then rst else (i, v) :: rst let rec update i f t = match t with - | [] -> cons i (f zero_num) [] + | [] -> cons i (f Q.zero) [] | (k, v) :: l -> ( match Int.compare i k with | 0 -> cons k (f v) l - | -1 -> cons i (f zero_num) t + | -1 -> cons i (f Q.zero) t | 1 -> (k, v) :: update i f l | _ -> failwith "compare_num" ) @@ -118,18 +115,17 @@ let rec set i n t = | 1 -> (k, v) :: set i n l | _ -> failwith "compare_num" ) -let cst n = if n =/ Int 0 then [] else [(0, n)] +let cst n = if n =/ Q.zero then [] else [(0, n)] let mul z t = - match z with - | Int 0 -> [] - | Int 1 -> t - | _ -> List.map (fun (i, n) -> (i, mult_num z n)) t + if z =/ Q.zero then [] + else if z =/ Q.one then t + else List.map (fun (i, n) -> (i, z */ n)) t let div z t = - if z <>/ Int 1 then List.map (fun (x, nx) -> (x, nx // z)) t else t + if z <>/ Q.one then List.map (fun (x, nx) -> (x, nx // z)) t else t -let uminus t = List.map (fun (i, n) -> (i, minus_num n)) t +let uminus t = List.map (fun (i, n) -> (i, Q.neg n)) t let rec add (ve1 : t) (ve2 : t) = match (ve1, ve2) with @@ -137,12 +133,12 @@ let rec add (ve1 : t) (ve2 : t) = | (v1, c1) :: l1, (v2, c2) :: l2 -> let cmp = Int.compare v1 v2 in if cmp == 0 then - let s = add_num c1 c2 in - if eq_num (Int 0) s then add l1 l2 else (v1, s) :: add l1 l2 + let s = c1 +/ c2 in + if Q.zero =/ s then add l1 l2 else (v1, s) :: add l1 l2 else if cmp < 0 then (v1, c1) :: add l1 ve2 else (v2, c2) :: add l2 ve1 -let rec xmul_add (n1 : num) (ve1 : t) (n2 : num) (ve2 : t) = +let rec xmul_add (n1 : Q.t) (ve1 : t) (n2 : Q.t) (ve2 : t) = match (ve1, ve2) with | [], _ -> mul n2 ve2 | _, [] -> mul n1 ve1 @@ -150,19 +146,19 @@ let rec xmul_add (n1 : num) (ve1 : t) (n2 : num) (ve2 : t) = let cmp = Int.compare v1 v2 in if cmp == 0 then let s = (n1 */ c1) +/ (n2 */ c2) in - if eq_num (Int 0) s then xmul_add n1 l1 n2 l2 + if Q.zero =/ s then xmul_add n1 l1 n2 l2 else (v1, s) :: xmul_add n1 l1 n2 l2 else if cmp < 0 then (v1, n1 */ c1) :: xmul_add n1 l1 n2 ve2 else (v2, n2 */ c2) :: xmul_add n1 ve1 n2 l2 let mul_add n1 ve1 n2 ve2 = - if n1 =/ Int 1 && n2 =/ Int 1 then add ve1 ve2 else xmul_add n1 ve1 n2 ve2 + if n1 =/ Q.one && n2 =/ Q.one then add ve1 ve2 else xmul_add n1 ve1 n2 ve2 let compare : t -> t -> int = Mutils.Cmp.compare_list (fun x y -> Mutils.Cmp.compare_lexical [ (fun () -> Int.compare (fst x) (fst y)) - ; (fun () -> compare_num (snd x) (snd y)) ]) + ; (fun () -> Q.compare (snd x) (snd y)) ]) (** [tail v vect] returns - [None] if [v] is not a variable of the vector [vect] @@ -181,28 +177,28 @@ let rec tail (v : var) (vect : t) = (* Hopeless *) -let get v vect = match tail v vect with None -> Int 0 | Some (vl, _) -> vl +let get v vect = match tail v vect with None -> Q.zero | Some (vl, _) -> vl let is_constant v = match v with [] | [(0, _)] -> true | _ -> false -let get_cst vect = match vect with (0, v) :: _ -> v | _ -> Int 0 +let get_cst vect = match vect with (0, v) :: _ -> v | _ -> Q.zero let choose v = match v with [] -> None | (vr, vl) :: rst -> Some (vr, vl, rst) let rec fresh v = match v with [] -> 1 | [(v, _)] -> v + 1 | _ :: v -> fresh v let variables v = List.fold_left (fun acc (x, _) -> ISet.add x acc) ISet.empty v -let decomp_cst v = match v with (0, vl) :: v -> (vl, v) | _ -> (Int 0, v) +let decomp_cst v = match v with (0, vl) :: v -> (vl, v) | _ -> (Q.zero, v) let rec decomp_at i v = match v with - | [] -> (Int 0, null) + | [] -> (Q.zero, null) | (vr, vl) :: r -> - if i = vr then (vl, r) else if i < vr then (Int 0, v) else decomp_at i r + if i = vr then (vl, r) else if i < vr then (Q.zero, v) else decomp_at i r -let decomp_fst v = match v with [] -> ((0, Int 0), []) | x :: v -> (x, v) +let decomp_fst v = match v with [] -> ((0, Q.zero), []) | x :: v -> (x, v) let rec subst (vr : int) (e : t) (v : t) = match v with | [] -> [] | (x, n) :: v' -> ( match Int.compare vr x with - | 0 -> mul_add n e (Int 1) v' + | 0 -> mul_add n e Q.one v' | -1 -> v | 1 -> add [(x, n)] (subst vr e v') | _ -> assert false ) @@ -227,25 +223,23 @@ let for_all p l = List.for_all (fun (v, n) -> p v n) l let decr_var i v = List.map (fun (v, n) -> (v - i, n)) v let incr_var i v = List.map (fun (v, n) -> (v + i, n)) v -open Big_int - let gcd v = let res = fold (fun c _ n -> - assert (Int.equal (compare_big_int (denominator n) unit_big_int) 0); - gcd_big_int c (numerator n)) - zero_big_int v + assert (Int.equal (Z.compare (Q.den n) Z.one) 0); + Z.gcd c (Q.num n)) + Z.zero v in - if Int.equal (compare_big_int res zero_big_int) 0 then unit_big_int else res + if Int.equal (Z.compare res Z.zero) 0 then Z.one else res let normalise v = - let ppcm = fold (fun c _ n -> ppcm c (denominator n)) unit_big_int v in + let ppcm = fold (fun c _ n -> Z.ppcm c (Q.den n)) Z.one v in let gcd = - let gcd = fold (fun c _ n -> gcd_big_int c (numerator n)) zero_big_int v in - if Int.equal (compare_big_int gcd zero_big_int) 0 then unit_big_int else gcd + let gcd = fold (fun c _ n -> Z.gcd c (Q.num n)) Z.zero v in + if Int.equal (Z.compare gcd Z.zero) 0 then Z.one else gcd in - List.map (fun (x, v) -> (x, v */ Big_int ppcm // Big_int gcd)) v + List.map (fun (x, v) -> (x, v */ Q.of_bigint ppcm // Q.of_bigint gcd)) v let rec exists2 p vect1 vect2 = match (vect1, vect2) with @@ -265,7 +259,7 @@ let dotproduct v1 v2 = else if x1 < x2 then dot acc v1' v2 else dot acc v1 v2' in - dot (Int 0) v1 v2 + dot Q.zero v1 v2 let map f v = List.map (fun (x, v) -> f x v) v @@ -276,18 +270,18 @@ let abs_min_elt v = Some (List.fold_left (fun (v1, vl1) (v2, vl2) -> - if abs_num vl1 </ abs_num vl2 then (v1, vl1) else (v2, vl2)) + if Q.abs vl1 </ Q.abs vl2 then (v1, vl1) else (v2, vl2)) (v, vl) r) let partition p = List.partition (fun (vr, vl) -> p vr vl) -let mkvar x = set x (Int 1) null +let mkvar x = set x Q.one null module Bound = struct - type t = {cst : num; var : var; coeff : num} + type t = {cst : Q.t; var : var; coeff : Q.t} let of_vect (v : vector) = match v with - | [(x, v)] -> if x = 0 then None else Some {cst = Int 0; var = x; coeff = v} + | [(x, v)] -> if x = 0 then None else Some {cst = Q.zero; var = x; coeff = v} | [(0, v); (x, v')] -> Some {cst = v; var = x; coeff = v'} | _ -> None end diff --git a/plugins/micromega/vect.mli b/plugins/micromega/vect.mli index 4b814cbb82..56c8ce87dd 100644 --- a/plugins/micromega/vect.mli +++ b/plugins/micromega/vect.mli @@ -8,7 +8,7 @@ (* * (see LICENSE file for the text of the license) *) (************************************************************************) -open Num +open NumCompat open Mutils type var = int @@ -50,18 +50,18 @@ val pp_smt : out_channel -> t -> unit val variables : t -> ISet.t (** [variables v] returns the set of variables with non-zero coefficients *) -val get_cst : t -> num +val get_cst : t -> Q.t (** [get_cst v] returns c i.e. the coefficient of the variable zero *) -val decomp_cst : t -> num * t +val decomp_cst : t -> Q.t * t (** [decomp_cst v] returns the pair (c,a1.x1+...+an.xn) *) -val decomp_at : int -> t -> num * t +val decomp_at : int -> t -> Q.t * t (** [decomp_cst v] returns the pair (ai, ai+1.xi+...+an.xn) *) -val decomp_fst : t -> (var * num) * t +val decomp_fst : t -> (var * Q.t) * t -val cst : num -> t +val cst : Q.t -> t (** [cst c] returns the vector v=c+0.x1+...+0.xn *) val is_constant : t -> bool @@ -74,33 +74,33 @@ val null : t val is_null : t -> bool (** [is_null v] returns whether [v] is the [null] vector i.e [equal v null] *) -val get : var -> t -> num +val get : var -> t -> Q.t (** [get xi v] returns the coefficient ai of the variable [xi]. [get] is also defined for the variable 0 *) -val set : var -> num -> t -> t +val set : var -> Q.t -> t -> t (** [set xi ai' v] returns the vector c+a1.x1+...ai'.xi+...+an.xn i.e. the coefficient of the variable xi is set to ai' *) val mkvar : var -> t (** [mkvar xi] returns 1.xi *) -val update : var -> (num -> num) -> t -> t +val update : var -> (Q.t -> Q.t) -> t -> t (** [update xi f v] returns c+a1.x1+...+f(ai).xi+...+an.xn *) val fresh : t -> int (** [fresh v] return the fresh variable with index 1+ max (variables v) *) -val choose : t -> (var * num * t) option +val choose : t -> (var * Q.t * t) option (** [choose v] decomposes a vector [v] depending on whether it is [null] or not. @return None if v is [null] @return Some(x,n,r) where v = r + n.x x is the smallest variable with non-zero coefficient n <> 0. *) -val from_list : num list -> t +val from_list : Q.t list -> t (** [from_list l] returns the vector c+a1.x1...an.xn from the list of coefficient [l=c;a1;...;an] *) -val to_list : t -> num list +val to_list : t -> Q.t list (** [to_list v] returns the list of all coefficient of the vector v i.e. [c;a1;...;an] The list representation is (obviously) not sparsed and therefore certain ai may be 0 *) @@ -114,7 +114,7 @@ val incr_var : int -> t -> t (** [incr_var i v] increments the variables of the vector [v] by the amount [i]. *) -val gcd : t -> Big_int.big_int +val gcd : t -> Z.t (** [gcd v] returns gcd(num(c),num(a1),...,num(an)) where num extracts the numerator of a rational value. *) @@ -130,17 +130,17 @@ val add : t -> t -> t @return c1+c1'+ (a1+a1').x1 + ... + (an+an').xn *) -val mul : num -> t -> t +val mul : Q.t -> t -> t (** [mul a v] is vector multiplication of vector [v] by a scalar [a]. @return a.v = a.c+a.a1.x1+...+a.an.xn *) -val mul_add : num -> t -> num -> t -> t +val mul_add : Q.t -> t -> Q.t -> t -> t (** [mul_add c1 v1 c2 v2] returns the linear combination c1.v1+c2.v2 *) val subst : int -> t -> t -> t (** [subst x v v'] replaces x by v in vector v' *) -val div : num -> t -> t +val div : Q.t -> t -> t (** [div c1 v1] returns the mutiplication by the inverse of c1 i.e (1/c1).v1 *) val uminus : t -> t @@ -148,36 +148,36 @@ val uminus : t -> t (** {1 Iterators} *) -val fold : ('acc -> var -> num -> 'acc) -> 'acc -> t -> 'acc +val fold : ('acc -> var -> Q.t -> 'acc) -> 'acc -> t -> 'acc (** [fold f acc v] returns f (f (f acc 0 c ) x1 a1 ) ... xn an *) -val fold_error : ('acc -> var -> num -> 'acc option) -> 'acc -> t -> 'acc option +val fold_error : ('acc -> var -> Q.t -> 'acc option) -> 'acc -> t -> 'acc option (** [fold_error f acc v] is the same as [fold (fun acc x i -> match acc with None -> None | Some acc' -> f acc' x i) (Some acc) v] but with early exit... *) -val find : (var -> num -> 'c option) -> t -> 'c option +val find : (var -> Q.t -> 'c option) -> t -> 'c option (** [find f v] returns the first [f xi ai] such that [f xi ai <> None]. If no such xi ai exists, it returns None *) -val for_all : (var -> num -> bool) -> t -> bool +val for_all : (var -> Q.t -> bool) -> t -> bool (** [for_all p v] returns /\_{i>=0} (f xi ai) *) -val exists2 : (num -> num -> bool) -> t -> t -> (var * num * num) option +val exists2 : (Q.t -> Q.t -> bool) -> t -> t -> (var * Q.t * Q.t) option (** [exists2 p v v'] returns Some(xi,ai,ai') if p(xi,ai,ai') holds and ai,ai' <> 0. It returns None if no such pair of coefficient exists. *) -val dotproduct : t -> t -> num +val dotproduct : t -> t -> Q.t (** [dotproduct v1 v2] is the dot product of v1 and v2. *) -val map : (var -> num -> 'a) -> t -> 'a list -val abs_min_elt : t -> (var * num) option -val partition : (var -> num -> bool) -> t -> t * t +val map : (var -> Q.t -> 'a) -> t -> 'a list +val abs_min_elt : t -> (var * Q.t) option +val partition : (var -> Q.t -> bool) -> t -> t * t module Bound : sig - type t = {cst : num; var : var; coeff : num} + type t = {cst : Q.t; var : var; coeff : Q.t} (** represents a0 + ai.xi *) val of_vect : vector -> t option diff --git a/plugins/micromega/zify.ml b/plugins/micromega/zify.ml index b3b627be14..dd8ea2c5ba 100644 --- a/plugins/micromega/zify.ml +++ b/plugins/micromega/zify.ml @@ -147,7 +147,7 @@ type classify_op = | OpConv (* e.g. Pos.ge == \x.y. Z.ge (Z.pos x) (Z.pos y) \x.y. Z.pos (Pos.add x y) == \x.y. Z.add (Z.pos x) (Z.pos y) Z.succ == (\x.x + 1) - *) + *) | OpOther (*let pp_classify_op = function @@ -1043,8 +1043,8 @@ let rec trans_expr env evd e = app_binop evd e binop a.(n - 2) prf1 a.(n - 1) prf2 | d -> mkvar evd inj e with Not_found -> - (* Feedback.msg_debug - Pp.(str "Not found " ++ Termops.Internal.debug_print_constr e); *) + (* Feedback.msg_debug + Pp.(str "Not found " ++ Termops.Internal.debug_print_constr e); *) mkvar evd inj e let trans_expr env evd e = |
