diff --git a/libpoly/accessors.h b/libpoly/accessors.h new file mode 100644 index 0000000000000000000000000000000000000000..6fabad84456ae123eda99259fffe3d40a3e9ee17 --- /dev/null +++ b/libpoly/accessors.h @@ -0,0 +1,10 @@ +static size_t get_dyadic_interval_a_open(const lp_dyadic_interval_t * p){ + return p->a_open; +} + +static size_t get_dyadic_interval_b_open(const lp_dyadic_interval_t * p){ + return p->b_open; +} +static size_t get_dyadic_interval_is_point(const lp_dyadic_interval_t * p){ + return p->is_point; +} diff --git a/libpoly/dune b/libpoly/dune index ff84f0ec1dd6d8abaf2953947e863797a72ba088..63dc44780aecf1741dff29910f740409fde08109 100644 --- a/libpoly/dune +++ b/libpoly/dune @@ -43,12 +43,13 @@ (flags -w -9-27) (ctypes (external_library_name poly) - (build_flags_resolver (vendored (c_flags "-I." "-Ilibpoly" :standard) (c_library_flags "-lpoly -lgmp"))) + (build_flags_resolver (vendored (c_flags "-I." "-Ilibpoly" "-Werror=implicit-function-declaration" :standard) (c_library_flags "-lpoly -lgmp"))) (deps poly.h libpoly.a zarith.h poly.h algebraic_number.h upolynomial.h rational.h interval.h upolynomial.h output_language.h dyadic_interval.h dyadic_rational.h integer.h version.h value.h + accessors.h ) - (headers (include "poly.h" "algebraic_number.h" "string.h" "zarith.h" "upolynomial.h")) + (headers (include "poly.h" "algebraic_number.h" "string.h" "zarith.h" "upolynomial.h" "accessors.h")) (type_description (instance Type) (functor Type_description)) diff --git a/libpoly/function_description.ml b/libpoly/function_description.ml index 51e292267e03ef1d4c1c91a4529b6aa589d44d2b..033d21315f81388fcedad63b0db2053e6d57ad4d 100644 --- a/libpoly/function_description.ml +++ b/libpoly/function_description.ml @@ -16,113 +16,113 @@ module Functions (F : Ctypes.FOREIGN) = struct foreign "lp_trace_disable" (ocaml_string @-> returning void) let lp_algebraic_number_sgn = - foreign "lp_algebraic_number_sgn" (ptr lp_algebraic_number @-> returning int) + foreign "lp_algebraic_number_sgn" (ptr Algebraic_number.t @-> returning int) let lp_algebraic_number_cmp = foreign "lp_algebraic_number_cmp" - (ptr lp_algebraic_number @-> ptr lp_algebraic_number @-> returning int) + (ptr Algebraic_number.t @-> ptr Algebraic_number.t @-> returning int) let lp_algebraic_number_cmp_integer = foreign "lp_algebraic_number_cmp_integer" - (ptr lp_algebraic_number @-> MPZ.t @-> returning int) + (ptr Algebraic_number.t @-> MPZ.t @-> returning int) let lp_algebraic_number_cmp_rational = foreign "lp_algebraic_number_cmp_rational" - (ptr lp_algebraic_number @-> MPQ.t @-> returning int) + (ptr Algebraic_number.t @-> MPQ.t @-> returning int) let lp_algebraic_number_to_double = foreign "lp_algebraic_number_to_double" - (ptr lp_algebraic_number @-> returning double) + (ptr Algebraic_number.t @-> returning double) let lp_algebraic_number_to_rational = foreign "lp_algebraic_number_to_rational" - (ptr lp_algebraic_number @-> MPQ.t @-> returning void) + (ptr Algebraic_number.t @-> MPQ.t @-> returning void) let lp_algebraic_number_construct_zero = foreign "lp_algebraic_number_construct_zero" - (ptr lp_algebraic_number @-> returning void) + (ptr Algebraic_number.t @-> returning void) let lp_algebraic_number_construct_one = foreign "lp_algebraic_number_construct_one" - (ptr lp_algebraic_number @-> returning void) + (ptr Algebraic_number.t @-> returning void) let lp_algebraic_number_destruct = foreign "lp_algebraic_number_destruct" - (ptr lp_algebraic_number @-> returning void) + (ptr Algebraic_number.t @-> returning void) let lp_algebraic_number_construct_from_integer = foreign "lp_algebraic_number_construct_from_integer" - (ptr lp_algebraic_number @-> MPZ.t @-> returning void) + (ptr Algebraic_number.t @-> MPZ.t @-> returning void) let lp_algebraic_number_construct_from_rational = foreign "lp_algebraic_number_construct_from_rational" - (ptr lp_algebraic_number @-> MPQ.t @-> returning void) + (ptr Algebraic_number.t @-> MPQ.t @-> returning void) let lp_algebraic_number_add = foreign "lp_algebraic_number_add" - (ptr lp_algebraic_number @-> ptr lp_algebraic_number - @-> ptr lp_algebraic_number @-> returning void) + (ptr Algebraic_number.t @-> ptr Algebraic_number.t + @-> ptr Algebraic_number.t @-> returning void) let lp_algebraic_number_sub = foreign "lp_algebraic_number_sub" - (ptr lp_algebraic_number @-> ptr lp_algebraic_number - @-> ptr lp_algebraic_number @-> returning void) + (ptr Algebraic_number.t @-> ptr Algebraic_number.t + @-> ptr Algebraic_number.t @-> returning void) let lp_algebraic_number_mul = foreign "lp_algebraic_number_mul" - (ptr lp_algebraic_number @-> ptr lp_algebraic_number - @-> ptr lp_algebraic_number @-> returning void) + (ptr Algebraic_number.t @-> ptr Algebraic_number.t + @-> ptr Algebraic_number.t @-> returning void) let lp_algebraic_number_div = foreign "lp_algebraic_number_div" - (ptr lp_algebraic_number @-> ptr lp_algebraic_number - @-> ptr lp_algebraic_number @-> returning void) + (ptr Algebraic_number.t @-> ptr Algebraic_number.t + @-> ptr Algebraic_number.t @-> returning void) let lp_algebraic_number_neg = foreign "lp_algebraic_number_neg" - (ptr lp_algebraic_number @-> ptr lp_algebraic_number @-> returning void) + (ptr Algebraic_number.t @-> ptr Algebraic_number.t @-> returning void) let lp_algebraic_number_inv = foreign "lp_algebraic_number_inv" - (ptr lp_algebraic_number @-> ptr lp_algebraic_number @-> returning void) + (ptr Algebraic_number.t @-> ptr Algebraic_number.t @-> returning void) let lp_algebraic_number_pow = foreign "lp_algebraic_number_pow" - (ptr lp_algebraic_number @-> ptr lp_algebraic_number @-> uint + (ptr Algebraic_number.t @-> ptr Algebraic_number.t @-> uint @-> returning void) let lp_algebraic_number_positive_root = foreign "lp_algebraic_number_positive_root" - (ptr lp_algebraic_number @-> ptr lp_algebraic_number @-> uint + (ptr Algebraic_number.t @-> ptr Algebraic_number.t @-> uint @-> returning void) let lp_algebraic_number_is_rational = foreign "lp_algebraic_number_is_rational" - (ptr lp_algebraic_number @-> returning bool) + (ptr Algebraic_number.t @-> returning bool) let lp_algebraic_number_is_integer = foreign "lp_algebraic_number_is_integer" - (ptr lp_algebraic_number @-> returning bool) + (ptr Algebraic_number.t @-> returning bool) let lp_algebraic_number_ceiling = foreign "lp_algebraic_number_ceiling" - (ptr lp_algebraic_number @-> MPZ.t @-> returning void) + (ptr Algebraic_number.t @-> MPZ.t @-> returning void) let lp_algebraic_number_floor = foreign "lp_algebraic_number_floor" - (ptr lp_algebraic_number @-> MPZ.t @-> returning void) + (ptr Algebraic_number.t @-> MPZ.t @-> returning void) let lp_algebraic_number_to_string = foreign "lp_algebraic_number_to_string" - (ptr lp_algebraic_number @-> returning (ptr char)) + (ptr Algebraic_number.t @-> returning (ptr char)) let lp_algebraic_number_swap = foreign "lp_algebraic_number_swap" - (ptr lp_algebraic_number @-> ptr lp_algebraic_number @-> returning void) + (ptr Algebraic_number.t @-> ptr Algebraic_number.t @-> returning void) let lp_algebraic_number_hash_approx = foreign "lp_algebraic_number_hash_approx" - (ptr lp_algebraic_number @-> uint @-> returning size_t) + (ptr Algebraic_number.t @-> uint @-> returning size_t) let lp_upolynomial_construct = foreign "lp_upolynomial_construct" @@ -133,7 +133,7 @@ module Functions (F : Ctypes.FOREIGN) = struct let lp_upolynomial_roots_isolate = foreign "lp_upolynomial_roots_isolate" - (ptr lp_upolynomial @-> ptr lp_algebraic_number @-> ptr size_t + (ptr lp_upolynomial @-> ptr Algebraic_number.t @-> ptr size_t @-> returning void) let lp_upolynomial_to_string = @@ -147,6 +147,30 @@ module Functions (F : Ctypes.FOREIGN) = struct let lp_upolynomial_delete = foreign "lp_upolynomial_delete" (ptr lp_upolynomial @-> returning void) + let lp_upolynomial_destruct_copy = + foreign "lp_upolynomial_construct_copy" + (ptr lp_upolynomial @-> returning (ptr lp_upolynomial)) + + let get_dyadic_interval_a_open = + foreign "get_dyadic_interval_a_open" + (ptr Dyadic_interval.t @-> returning bool) + + let get_dyadic_interval_b_open = + foreign "get_dyadic_interval_b_open" + (ptr Dyadic_interval.t @-> returning bool) + + let get_dyadic_interval_is_point = + foreign "get_dyadic_interval_is_point" + (ptr Dyadic_interval.t @-> returning bool) + + let lp_dyadic_rational_get_num = + foreign "lp_dyadic_rational_get_num" + (ptr lp_dyadic_rational @-> MPZ.t @-> returning void) + + let lp_dyadic_rational_get_den = + foreign "lp_dyadic_rational_get_den" + (ptr lp_dyadic_rational @-> MPZ.t @-> returning void) + let lp_Z = foreign_value "lp_Z" (ptr lp_int_ring_t) let ml_z_mpz_init_set_z = diff --git a/libpoly/ocaml_poly.ml b/libpoly/ocaml_poly.ml index 51c49af0b016040c1625071aa200d2503115c864..bcb83f00ae237deee980377e343e0bd211b8cfb2 100644 --- a/libpoly/ocaml_poly.ml +++ b/libpoly/ocaml_poly.ml @@ -5,13 +5,19 @@ module PolyUtils = struct let trace_disable tag = C.Function.lp_trace_disable (ocaml_string_start tag) end +let finalise f p = + let typ = reference_type p in + let r = raw_address_of_ptr (to_voidp p) in + Gc.finalise_last (fun () -> f (from_voidp typ (ptr_of_raw_address r))) p; + p + module A = struct - type t = C.Type.lp_algebraic_number structure ptr + type t = C.Type.Algebraic_number.t structure ptr let mk () = allocate_n ~finalise:(fun p -> C.Function.lp_algebraic_number_destruct p) - C.Type.lp_algebraic_number ~count:1 + C.Type.Algebraic_number.t ~count:1 let zero () = let a = mk () in @@ -147,6 +153,37 @@ module A = struct let le a b = compare a b <= 0 let lt a b = compare a b < 0 let pp fmt a = Format.fprintf fmt "%s" (to_string a) + + type bound = Strict | Large + + type upolynomial_interval = + | IsPoint of Q.t + | U of C.Type.lp_upolynomial structure ptr * bound * Q.t * Q.t * bound + + let dyadic_to_q d = + let num = + return_z (fun mpz -> C.Function.lp_dyadic_rational_get_num d mpz) + in + let den = + return_z (fun mpz -> C.Function.lp_dyadic_rational_get_den d mpz) + in + Q.make num den + + let get_upolynomial_interval a = + let i = a |-> C.Type.Algebraic_number.i in + if C.Function.get_dyadic_interval_is_point i then + IsPoint (dyadic_to_q (i |-> C.Type.Dyadic_interval.a)) + else + let u = + finalise C.Function.lp_upolynomial_delete + !@(a |-> C.Type.Algebraic_number.f) + in + let a_open = C.Function.get_dyadic_interval_a_open i in + let b_open = C.Function.get_dyadic_interval_b_open i in + let to_bound open_ = if open_ then Strict else Large in + let a = dyadic_to_q (i |-> C.Type.Dyadic_interval.a) in + let b = dyadic_to_q (i |-> C.Type.Dyadic_interval.b) in + U (u, to_bound a_open, a, b, to_bound b_open) end module UPolynomial = struct @@ -192,7 +229,7 @@ module UPolynomial = struct let a = allocate_n ~finalise:(fun p -> C.Function.lp_algebraic_number_destruct p) - C.Type.lp_algebraic_number ~count:(degree u) + C.Type.Algebraic_number.t ~count:(degree u) in let n = allocate size_t Unsigned.Size_t.zero in C.Function.lp_upolynomial_roots_isolate u a n; diff --git a/libpoly/ocaml_poly.mli b/libpoly/ocaml_poly.mli index 89d5984195bdc5f8ab2a14b317b080f2880037ed..69a4e53d723a38736f52467e71224c6ab44e58e5 100644 --- a/libpoly/ocaml_poly.mli +++ b/libpoly/ocaml_poly.mli @@ -1,5 +1,5 @@ (** Algebraic number *) -module A : sig +module rec A : sig type t (** Algebraic number *) @@ -49,10 +49,20 @@ module A : sig val positive_pow : t -> Q.t -> t (** the numerator and denominator should fit in an int *) + + type bound = Strict | Large + + type upolynomial_interval = + | IsPoint of Q.t + | U of UPolynomial.t * bound * Q.t * Q.t * bound + + val get_upolynomial_interval : t -> upolynomial_interval + (** The polynomial which is used to define the algebraic number. The + polynomial is not necessarily unique *) end (** Polynomial in one variable *) -module UPolynomial : sig +and UPolynomial : sig type t val to_string : t -> string diff --git a/libpoly/type_description.ml b/libpoly/type_description.ml index 52dc6abbe2c35421daf02ada1ab05514977d4305..b86cb99a30bacc382929e7177db2fcaa7e1c88db 100644 --- a/libpoly/type_description.ml +++ b/libpoly/type_description.ml @@ -17,16 +17,40 @@ module Types (F : Ctypes.TYPE) = struct seal s; typedef s (name ^ "_t") - type lp_algebraic_number - - let lp_algebraic_number : lp_algebraic_number structure typ = - libpoly_typedef_structure_sized "lp_algebraic_number" - type lp_upolynomial let lp_upolynomial : lp_upolynomial structure typ = libpoly_typedef_structure_nosize "lp_upolynomial" + type lp_dyadic_rational + + let lp_dyadic_rational : lp_dyadic_rational structure typ = + libpoly_typedef_structure_nosize "lp_dyadic_rational" + + module Dyadic_interval = struct + type t + + let str : t structure typ = structure "lp_dyadic_interval_struct" + let t = typedef str "lp_dyadic_interval_t" + + (* the other fields are bit-fields *) + let a = field str "a" lp_dyadic_rational + let b = field str "b" lp_dyadic_rational + let () = seal str + end + + module Algebraic_number = struct + type t + + let str : t structure typ = structure "lp_algebraic_number_struct" + let t = typedef str "lp_algebraic_number_t" + let f = field str "f" (ptr lp_upolynomial) + let i = field str "I" Dyadic_interval.t + let sgn_at_a = field str "sgn_at_a" int + let sgn_at_b = field str "sgn_at_b" int + let () = seal str + end + type lp_int_ring_t let lp_int_ring_t : lp_int_ring_t structure typ = @@ -41,13 +65,9 @@ module Types (F : Ctypes.TYPE) = struct typedef (structure "__mpz_struct") "__mpz_struct" let mp_alloc = field _struct "_mp_alloc" int - let mp_size = field _struct "_mp_size" int - let mp_limb_t = field _struct "_mp_d" (ptr mp_limb_t) - let () = seal _struct - let t = ptr _struct end @@ -58,11 +78,8 @@ module Types (F : Ctypes.TYPE) = struct typedef (structure "__mpq_struct") "__mpq_struct" let mp_num = field _struct "_mp_num" MPZ._struct - let mp_den = field _struct "_mp_den" MPZ._struct - let () = seal _struct - let t = ptr _struct end diff --git a/src_colibri2/stdlib/std.ml b/src_colibri2/stdlib/std.ml index 5d9e610f10b55c5f3afc7395db4d3f09aa99f88b..92b2a16e2ddb0a470b2a819518586eb4a0eeab5a 100644 --- a/src_colibri2/stdlib/std.ml +++ b/src_colibri2/stdlib/std.ml @@ -166,7 +166,12 @@ module A = struct | Q a, A b -> -A.compare_q b a | A a, A b -> A.compare a b - let equal a b = compare a b = 0 + let equal a b = + match (a, b) with + | Q a, Q b -> Q.equal a b + | A a, Q b -> (not (Z.equal b.den Z.one)) && A.compare_q a b = 0 + | Q a, A b -> (not (Z.equal a.den Z.one)) && A.compare_q b a = 0 + | A a, A b -> A.compare a b = 0 let pp fmt = function | Q q -> Q.pp fmt q @@ -222,10 +227,11 @@ module A = struct let combine2 fq fa cv a b = match (a, b) with Q a, Q b -> Q (fq a b) | _ -> cv (fa (to_a a) (to_a b)) + (* todo special case for one, zero, ... *) let div = combine2 Q.div A.div normalize let add = combine2 Q.add A.add normalize let sub = combine2 Q.sub A.sub normalize - let mul = combine2 Q.mul A.add normalize + let mul = combine2 Q.mul A.mul normalize let ( + ) = add let ( - ) = sub let ( ~- ) = neg diff --git a/src_colibri2/stdlib/std.mli b/src_colibri2/stdlib/std.mli index 07e0c10974bd02a9a078d768fa7804af9541372f..e30dea97e47f2dfe0899faa8ab2850b8a3aaddfc 100644 --- a/src_colibri2/stdlib/std.mli +++ b/src_colibri2/stdlib/std.mli @@ -80,8 +80,12 @@ module Q : sig val shrink : t QCheck.Shrink.t end +(** Algebraic number with use of rational when possible *) module A : sig - type t = Q of Q.t | A of Ocaml_poly.A.t + (** {!A} allows to determine certainly if something is an integer, but is not + complete for rational. So Q is always used for integer and as best effort + for rational *) + type t = private Q of Q.t | A of Ocaml_poly.A.t val zero : t val one : t diff --git a/src_colibri2/tests/solve/models/dune.inc b/src_colibri2/tests/solve/models/dune.inc index fe6ebbbff1d6408d67bd7f36313dac537d9eb36c..5a9664c49d62b2a956a56851624bcfb829f2cdb0 100644 --- a/src_colibri2/tests/solve/models/dune.inc +++ b/src_colibri2/tests/solve/models/dune.inc @@ -4,3 +4,5 @@ (rule (alias runtest) (action (diff function.smt2.oracle function.smt2.res))) (rule (action (with-stdout-to get_value.smt2.res (run %{bin:colibri2} --size=50M --time=30s --max-steps 3500 --check-status colibri2 %{dep:get_value.smt2})))) (rule (alias runtest) (action (diff get_value.smt2.oracle get_value.smt2.res))) +(rule (action (with-stdout-to sqrt.smt2.res (run %{bin:colibri2} --size=50M --time=30s --max-steps 3500 --check-status colibri2 %{dep:sqrt.smt2})))) +(rule (alias runtest) (action (diff sqrt.smt2.oracle sqrt.smt2.res))) diff --git a/src_colibri2/tests/solve/models/sqrt.smt2 b/src_colibri2/tests/solve/models/sqrt.smt2 new file mode 100644 index 0000000000000000000000000000000000000000..1256c0abfef3a95f279ce8adff2be321f7d4fc42 --- /dev/null +++ b/src_colibri2/tests/solve/models/sqrt.smt2 @@ -0,0 +1,9 @@ +(set-logic ALL) +(set-info :status-colibri2 sat) + +(declare-fun b () Real) + +(assert (= (* b b) 2.0)) + +(check-sat) +(get-model) diff --git a/src_colibri2/tests/solve/models/sqrt.smt2.oracle b/src_colibri2/tests/solve/models/sqrt.smt2.oracle new file mode 100644 index 0000000000000000000000000000000000000000..1f8f0bdfa0102d566ee74eb95e4f0b9e438e4d83 --- /dev/null +++ b/src_colibri2/tests/solve/models/sqrt.smt2.oracle @@ -0,0 +1,2 @@ +sat +((b (epsilon (x) (and (= (1*x^2 + (-2)) 0) (< 181/128 x) (< x 1449/1024))))) diff --git a/src_colibri2/theories/LRA/dom_product.ml b/src_colibri2/theories/LRA/dom_product.ml index bb2bfced681703ebd413ff953c62dabfcd8eb40d..c0a816d02352e9addb770de56360103b47753c65 100644 --- a/src_colibri2/theories/LRA/dom_product.ml +++ b/src_colibri2/theories/LRA/dom_product.ml @@ -100,12 +100,13 @@ end = Pivot.WithUnsolved (struct include Product - let pp fmt p = Fmt.pf fmt "|%a|" pp p let of_one_node _ n = Product.monome n Q.one let subst p n q = - let p, c = subst p n q in - if Q.is_zero c then None else Some p + let p', c = subst p n q in + Debug.dprintf10 debug "subst: %a %a %a -> %a %a" Product.pp p Node.pp n + Product.pp q Product.pp p' Q.pp c; + if Q.is_zero c then None else Some p' let set d cl ~old_ ~new_ = SolveSign.iter_eqs d cl ~f:(fun sign -> @@ -120,7 +121,7 @@ end = Pivot.WithUnsolved (struct let q = f cl in Product.x_m_yc acc q c in - Product.fold add Product.one p + Product.fold add (Product.cst p.cst) p type info = data Node.M.t * data Node.M.t * data Node.M.t * Product.t @@ -140,20 +141,28 @@ end = Pivot.WithUnsolved (struct let solve (_, notz1, unk1, p1) (_, notz2, unk2, p2) : _ Pivot.solve_with_unsolved = let exception Solved of Node.t * Product.t in - try - (if - Node.M.is_empty notz1 && Node.M.is_empty unk1 && Node.M.is_empty notz2 - && Node.M.is_empty unk2 - then - (* all (strictly) positive: can use root *) - let p = Product.div p1 p2 in - match Product.extract p with - | One -> () - | Var (q, n, p) -> - let p = Product.power_cst (Q.inv (Q.neg q)) p in - raise (Solved (n, p)) - else if Node.M.is_empty unk1 then - (* side 1 with all different from zero *) + let all_zero (p : Product.t) : _ Pivot.solve_with_unsolved = + if Node.M.is_empty p.poly then + if A.is_zero p.cst then AlreadyEqual else Contradiction + else Subst (Node.M.map (fun _ -> Product.zero) p.poly) + in + if Product.is_zero p1 then all_zero p2 + else if Product.is_zero p2 then all_zero p1 + else if + Node.M.is_empty notz1 && Node.M.is_empty unk1 && Node.M.is_empty notz2 + && Node.M.is_empty unk2 + then + (* all (strictly) positive: can use root *) + let p = Product.div p1 p2 in + match Product.extract p with + | One -> AlreadyEqual + | Cst _ -> Contradiction + | Var (q, n, p) -> + let p = Product.power_cst (Q.inv (Q.neg q)) p in + Subst (Node.M.singleton n p) + else if Node.M.is_empty unk1 then + (* side 1 with all different from zero *) + try Node.M.iter (fun n q1 -> let q2 = Node.M.find_def Q.zero n p2.Product.poly in @@ -162,14 +171,16 @@ end = Pivot.WithUnsolved (struct assert (Q.equal (Node.M.find n p.poly) Q.minus_one); let p = Product.remove n p in raise (Solved (n, p)))) - p1.Product.poly - else if Node.M.is_num_elt 1 unk1 then - let n, q = Node.M.choose unk1 in - if Q.equal q Q.one && not (Node.M.mem n unk2) then - let p = Product.div p2 (Product.remove n p1) in - raise (Solved (n, p))); - Unsolved - with Solved (n, p) -> Subst (n, p) + p1.Product.poly; + Unsolved + with Solved (n, p) -> Subst (Node.M.singleton n p) + else if Node.M.is_num_elt 1 unk1 then + let n, q = Node.M.choose unk1 in + if Q.equal q Q.one && not (Node.M.mem n unk2) then + let p = Product.div p2 (Product.remove n p1) in + Subst (Node.M.singleton n p) + else Unsolved + else Unsolved end) and SolveSign : sig @@ -216,12 +227,12 @@ end = Pivot.WithUnsolved (struct | Plus_one -> AlreadyEqual | Minus_one -> Contradiction | Zero -> Contradiction - | Var (n, p) -> Subst (n, p) + | Var (n, p) -> Subst (Node.M.singleton n p) | NoNonZero -> Unsolved) | OneNonZero n -> if not (Node.M.mem n (Sign_product.nodes p2)) then let p = Sign_product.mul p2 (Sign_product.remove n p1) in - Subst (n, p) + Subst (Node.M.singleton n p) else Unsolved | Other -> Unsolved end) @@ -252,21 +263,22 @@ let factorize res a coef b d _ = let abs = Product.div abs common_abs in let sign = Sign_product.mul sign common_sign in match (Product.classify abs, Sign_product.classify sign) with - | ONE, PLUS_ONE -> (Q.add v cst, acc) - | ONE, MINUS_ONE -> (Q.sub v cst, acc) - | NODE n, NODE n' when Egraph.is_equal d n n' -> - (cst, (n, A.of_q v) :: acc) + | CST cst', PLUS_ONE -> (A.add (A.mul cst' v) cst, acc) + | CST cst', MINUS_ONE -> (A.sub (A.mul cst' v) cst, acc) + | NODE (cst', n), NODE (cst'', n') when Egraph.is_equal d n n' + -> + (cst, (n, A.mul (Sign_product.mul_cst cst'' cst') v) :: acc) | _ -> let n = node_of_product d abs sign in Egraph.register d n; - (cst, (n, A.of_q v) :: acc)) - (Q.zero, []) - [ (pa, sa, Q.one); (pb, sb, coef) ] + (cst, (n, v) :: acc)) + (A.zero, []) + [ (pa, sa, A.one); (pb, sb, coef) ] in - let p = Polynome.of_list (A.of_q cst) l in + let p = Polynome.of_list cst l in let n = Dom_polynome.node_of_polynome d p in let pp_coef fmt coef = - if Q.equal Q.one coef then Fmt.pf fmt "+" else Fmt.pf fmt "-" + if A.equal A.one coef then Fmt.pf fmt "+" else Fmt.pf fmt "-" in Debug.dprintf10 debug "[Factorize] %a = %a %a %a into %a" Node.pp res Node.pp a pp_coef coef Node.pp b Polynome.pp p; @@ -319,12 +331,12 @@ let converter d (f : Ground.t) = let a, b = IArray.extract2_exn args in reg a; reg b; - propagate_a_cb d ~res ~a ~b ~c:Q.one + propagate_a_cb d ~res ~a ~b ~c:A.one | { app = { builtin = Expr.Sub }; tyargs = []; args; _ } -> let a, b = IArray.extract2_exn args in reg a; reg b; - propagate_a_cb d ~res ~a ~b ~c:Q.minus_one + propagate_a_cb d ~res ~a ~b ~c:A.minus_one | { app = { builtin = Expr.Abs }; tyargs = []; args; _ } -> let a = IArray.extract1_exn args in reg a; @@ -342,29 +354,51 @@ let init env = match Dom_polynome.get_repr d n with | None -> () | Some p -> - if A.equal p.cst A.zero && Node.M.is_num_elt 1 p.poly then + if Node.M.is_empty p.poly then ( + (* p = cst *) + let p1 = Product.cst (A.abs p.cst) in + Debug.dprintf4 debug "[Polynome->Product] propagate %a is %a" + Node.pp n Polynome.pp p; + SolveAbs.assume_equality d n p1; + SolveSign.assume_equality d n (Sign_product.cst p.cst)) + else if A.equal p.cst A.zero && Node.M.is_num_elt 1 p.poly then ( + (* p = q * cl^1 *) let cl, q = Node.M.choose p.poly in - if A.equal q A.one then ( - let p1 = Product.of_list A.one [ (cl, Q.one) ] in - Debug.dprintf4 debug "[Polynome->Product] propagate %a is %a" - Node.pp n Node.pp cl; - SolveAbs.assume_equality d n p1; - SolveSign.assume_equality d n (Sign_product.of_one_node cl))); + let p1 = Product.of_list (A.abs q) [ (cl, Q.one) ] in + Debug.dprintf4 debug "[Polynome->Product] propagate %a is %a" + Node.pp n Polynome.pp p; + SolveAbs.assume_equality d n p1; + let sp = Sign_product.(mul (of_one_node cl) (cst q)) in + SolveSign.assume_equality d n sp)); let product_polynome d n = (* todo more efficient *) SolveAbs.iter_eqs d n ~f:(fun p -> - match Product.is_one_node p with - | None -> () - | Some cl -> + (* Look for the two following cases: *) + match Product.classify p with + | CST cst when A.is_zero cst -> Egraph.merge d n RealValue.zero + | CST cst -> + (* - |n| = cst /\ sgn(n)= ±1 *) + SolveSign.iter_eqs d n ~f:(fun p -> + match Sign_product.classify p with + | PLUS_ONE -> RealValue.set d n cst + | MINUS_ONE -> RealValue.set d n (A.neg cst) + | ZERO -> + (* |n| <> 0 and sgn(n) = 0 *) + Egraph.contradiction d + | NODE _ | PRODUCT -> ()) + | NODE (cst, cl) -> + (* |n| = cst*|cl| /\ sgn(n) = ±sgn(cl) *) SolveSign.iter_eqs d n ~f:(fun p -> - match Sign_product.is_one_node p with - | Some cl' when Egraph.is_equal d cl cl' -> - let p1 = Polynome.of_list A.zero [ (cl, A.one) ] in + match Sign_product.classify p with + | NODE (cst', cl') when Egraph.is_equal d cl cl' -> + let cst = Sign_product.mul_cst cst' cst in + let p1 = Polynome.of_list A.zero [ (cl, cst) ] in Debug.dprintf4 debug "[Product->Polynome] propagate %a is %a" Node.pp n Node.pp cl; Dom_polynome.assume_equality d n p1 - | _ -> ())) + | _ -> ()) + | PRODUCT -> ()) in SolveAbs.attach_eqs_change env product_polynome; SolveSign.attach_eqs_change env product_polynome; diff --git a/src_colibri2/theories/LRA/dom_product.mli b/src_colibri2/theories/LRA/dom_product.mli index ae452721d39969d4bd169a9f685993489b7c5f70..44d89ad1e6462b2d53dc7029ca79a2c4f54b04a0 100644 --- a/src_colibri2/theories/LRA/dom_product.mli +++ b/src_colibri2/theories/LRA/dom_product.mli @@ -33,4 +33,4 @@ end val init : Egraph.wt -> unit val propagate_a_cb : - _ Egraph.t -> res:Node.t -> a:Node.t -> c:Q.t -> b:Node.t -> unit + _ Egraph.t -> res:Node.t -> a:Node.t -> c:A.t -> b:Node.t -> unit diff --git a/src_colibri2/theories/LRA/fourier.ml b/src_colibri2/theories/LRA/fourier.ml index afc23ecee50033504762700ddaa4eb8102d88194..9f4d8605fc9db07dd421753fb465a1c8a5b3f303 100644 --- a/src_colibri2/theories/LRA/fourier.ml +++ b/src_colibri2/theories/LRA/fourier.ml @@ -81,6 +81,7 @@ let divide d (p : Polynome.t) = let cst, l = List.fold_left (fun (cst, acc) (abs, sign, v) -> + (* v * abs * sign *) let abs = Product.div abs common in let pos, sign = Sign_product.extract_cst sign in let v = @@ -88,9 +89,9 @@ let divide d (p : Polynome.t) = in match (Product.classify abs, Sign_product.classify sign) with | _, MINUS_ONE -> assert false - | ONE, PLUS_ONE -> (A.add v cst, acc) - | NODE n, NODE n' when Egraph.is_equal d n n' -> - (cst, (n, v) :: acc) + | CST cst', PLUS_ONE -> (A.add (A.mul cst' v) cst, acc) + | NODE (cst', n), NODE (cst'', n') when Egraph.is_equal d n n' -> + (cst, (n, A.mul (Sign_product.mul_cst cst'' cst') v) :: acc) | _ -> let n = Dom_product.node_of_product d abs sign in (cst, (n, v) :: acc)) diff --git a/src_colibri2/theories/LRA/pivot.ml b/src_colibri2/theories/LRA/pivot.ml index afcc44c928273ee55c431555732c25b4899d01b4..006ab309406a6487eba94f756ffaf53c44879b21 100644 --- a/src_colibri2/theories/LRA/pivot.ml +++ b/src_colibri2/theories/LRA/pivot.ml @@ -25,7 +25,7 @@ type 'a solve_with_unsolved = | AlreadyEqual | Contradiction | Unsolved - | Subst of Node.t * 'a + | Subst of 'a Node.M.t module WithUnsolved (P : sig type t @@ -35,11 +35,8 @@ module WithUnsolved (P : sig include Colibri2_popop_lib.Popop_stdlib.Datatype with type t := t val of_one_node : _ Egraph.t -> Node.t -> t - val is_one_node : t -> Node.t option - val subst : t -> Node.t -> t -> t option - val normalize : t -> f:(Node.t -> t) -> t type data @@ -54,15 +51,11 @@ module WithUnsolved (P : sig Egraph.wt -> (Egraph.rt -> Node.t -> Events.enqueue) -> unit val solve : info -> info -> t solve_with_unsolved - val set : Egraph.wt -> Node.t -> old_:t -> new_:t -> unit end) : sig val assume_equality : Egraph.wt -> Node.t -> P.t -> unit - val init : Egraph.wt -> unit - val get_repr : _ Egraph.t -> Node.t -> P.t option - val iter_eqs : _ Egraph.t -> Node.t -> f:(P.t -> unit) -> unit val attach_repr_change : @@ -132,6 +125,12 @@ end = struct | None -> P.of_one_node d cl | Some p -> p.repr) + let norm_dom d cl = function + | None -> + let r = P.of_one_node d cl in + { repr = r; eqs = P.S.empty } + | Some p -> p + module Th = struct let merged v1 v2 = Base.phys_equal v1 v2 @@ -193,33 +192,34 @@ end = struct Egraph.set_dom d dom cl2 p) and merge_one_new_eq d cl eq = + let eq = norm_product d eq in let po = Egraph.get_dom d dom cl in - match po with - | None -> - add_used_product d cl (P.nodes eq); - set_poly d cl { repr = eq; eqs = P.S.empty } [ (eq, eq) ] - | Some p -> ( - if (not (P.S.mem eq p.eqs)) && not (P.equal eq p.repr) then - match - solve d - (Base.List.cartesian_product - (part d (p.repr :: P.S.elements p.eqs)) - (part d [ eq ])) - with - | `Solved -> - (* The domains have been substituted, and possibly recursively *) - let eq = norm_product d eq in - merge_one_new_eq d cl eq - | `Not_solved -> - (* nothing to solve *) - let repr = p.repr in - let eqs = p.eqs |> P.S.add eq |> P.S.remove repr in - let p = { repr; eqs } in - add_used_product d cl (P.nodes eq); - set_poly d cl p [ (eq, eq) ]) + if Option.is_some po || Node.M.mem cl (P.nodes eq) then ( + let p = norm_dom d cl po in + if (not (P.S.mem eq p.eqs)) && not (P.equal eq p.repr) then + match + solve d + (Base.List.cartesian_product + (part d (p.repr :: P.S.elements p.eqs)) + (part d [ eq ])) + with + | `Solved -> + (* The domains have been substituted, and possibly recursively *) + merge_one_new_eq d cl eq + | `Not_solved -> + (* nothing to solve *) + let repr = p.repr in + let eqs = p.eqs |> P.S.add eq |> P.S.remove repr in + let p = { repr; eqs } in + add_used_product d cl (P.nodes eq); + set_poly d cl p [ (eq, eq) ]) + else ( + add_used_product d cl (P.nodes eq); + set_poly d cl { repr = eq; eqs = P.S.empty } [ (eq, eq) ]) and subst d cl eq = - Debug.dprintf4 debug "[Pivot] subst %a with %a" Node.pp cl P.pp eq; + Debug.dprintf5 debug "[Pivot:%s] subst %a with %a" P.name Node.pp cl P.pp + eq; let po = Egraph.get_dom d dom cl in match po with | None -> @@ -228,11 +228,7 @@ end = struct set_poly d cl p [ (eq, eq) ] | Some p -> assert (Option.equal Node.equal (P.is_one_node p.repr) (Some cl)); - subst_doms d cl eq; - let eq = norm_product d eq in - assert ( - P.equal eq (Base.Option.value_exn (Egraph.get_dom d dom cl)).repr); - merge_one_new_eq d cl eq + subst_doms d cl eq and subst_doms d cl p = let b = @@ -260,7 +256,7 @@ end = struct let new_cl, acc, chg = fold (Node.M.empty, P.S.empty, []) q.repr in - let repr = P.S.choose acc in + let repr = P.S.choose acc (* there is only one in acc *) in let new_cl, acc, chg = P.S.fold_left fold (new_cl, acc, chg) q.eqs in @@ -273,21 +269,23 @@ end = struct and part d l = List.map (fun p -> P.info d p) l and solve d l = - let exception Solved of Node.t * P.t in + let exception Solved of P.t Node.M.t in let criteria i1 i2 = let aux i1 i2 = match P.solve i1 i2 with | AlreadyEqual -> () | Contradiction -> Egraph.contradiction d | Unsolved -> () - | Subst (n, p) -> raise (Solved (n, p)) + | Subst m -> raise (Solved m) in aux i1 i2; aux i2 i1 in match List.iter (fun (a, b) -> criteria a b) l with - | exception Solved (n, p) -> + | exception Solved m -> + let n, p = Node.M.choose m in subst d n p; + Node.M.iter (fun n p -> merge_one_new_eq d n p) (Node.M.remove n m); `Solved | () -> `Not_solved @@ -328,7 +326,6 @@ end = struct let assume_equality d n (p : P.t) = Debug.dprintf5 debug "[Pivot %s] assume %a = %a" P.name Node.pp n P.pp p; let n = Egraph.find d n in - let p = norm_product d p in Th.merge_one_new_eq d n p let reshape d cl ~(f : P.t -> P.t option) = @@ -393,7 +390,6 @@ end = struct end let () = Events.register (module ChangeInfo) - let init d = ChangeInfo.init d let attach_eqs_change d ?node f = @@ -414,25 +410,18 @@ module Total (P : sig include Colibri2_popop_lib.Popop_stdlib.Datatype with type t := t val of_one_node : Node.t -> t - val is_one_node : t -> Node.t option - val subst : t -> Node.t -> t -> t - val normalize : t -> f:(Node.t -> t) -> t type data val nodes : t -> data Node.M.t - val solve : t -> t -> t solve_total - val set : Egraph.wt -> Node.t -> old_:t option -> new_:t -> unit end) : sig val assume_equality : Egraph.wt -> Node.t -> P.t -> unit - val init : Egraph.wt -> unit - val get_repr : _ Egraph.t -> Node.t -> P.t option val attach_repr_change : diff --git a/src_colibri2/theories/LRA/pivot.mli b/src_colibri2/theories/LRA/pivot.mli index 1112944c7f41b68a5ebcfcd04bec5d092270119c..df1009a39fc6994f571a76becb77d61720d6ceaa 100644 --- a/src_colibri2/theories/LRA/pivot.mli +++ b/src_colibri2/theories/LRA/pivot.mli @@ -24,7 +24,7 @@ type 'a solve_with_unsolved = | AlreadyEqual | Contradiction | Unsolved - | Subst of Node.t * 'a + | Subst of 'a Node.M.t module WithUnsolved (P : sig type t @@ -79,7 +79,6 @@ end) : sig (** Initialize the data-structure needed. *) val get_repr : _ Egraph.t -> Node.t -> P.t option - val iter_eqs : _ Egraph.t -> Node.t -> f:(P.t -> unit) -> unit val attach_repr_change : diff --git a/src_colibri2/theories/LRA/product.ml b/src_colibri2/theories/LRA/product.ml index 2139dd1a235600ee11297ccbee197ade713cc68a..318601dae4fccdcf855e5003dce70ec306a79112 100644 --- a/src_colibri2/theories/LRA/product.ml +++ b/src_colibri2/theories/LRA/product.ml @@ -40,13 +40,10 @@ module T = struct else Fmt.pf fmt "^%a" Q.pp q let pp fmt v = - (* let print_not_one fmt q = - * if not(Q.equal q Q.one) then Q.pp fmt q - * in *) - if Node.M.is_empty v.poly then Q.pp fmt Q.one + let print_not_one fmt q = if not (A.equal q A.one) then A.pp fmt q in + if Node.M.is_empty v.poly then A.pp fmt v.cst else - Fmt.pf fmt "@[%a@]" - (* print_not_one v.cst *) + Fmt.pf fmt "@[%a|%a|@]" print_not_one v.cst Fmt.(iter_bindings Node.M.iter (pair ~sep:nop Node.pp pp_exposant)) v.poly end @@ -59,8 +56,10 @@ include Popop_stdlib.MkDatatype (T) let invariant p = A.sign p.cst > 0 && Node.M.for_all (fun _ q -> not (Q.equal q Q.zero)) p.poly -(* let cst q = assert (not (Q.equal q Q.zero)); {cst = q; poly = Node.M.empty} - * (\* let zero = cst Q.zero *\) *) +let cst q = { cst = q; poly = Node.M.empty } +let zero = cst A.zero +let is_zero p = A.is_zero p.cst +let is_not_zero p = A.is_not_zero p.cst (** constructor *) let one = { cst = A.one; poly = Node.M.empty } @@ -69,11 +68,10 @@ let one = { cst = A.one; poly = Node.M.empty } * (\* let is_zero p = Q.equal p.cst Q.zero && Node.M.is_empty p.poly *\) * let is_one p = Q.equal p.cst Q.one && Node.M.is_empty p.poly *) -type extract = One | Var of Q.t * Node.t * t +type extract = One | Cst of A.t | Var of Q.t * Node.t * t let extract p = - if Node.M.is_empty p.poly then (* if Q.equal p.cst Q.one then *) One - (* else Cst p.cst *) + if Node.M.is_empty p.poly then if A.equal p.cst A.one then One else Cst p.cst else let x, q = Shuffle.chooseb Node.M.choose Node.M.choose_rnd p.poly in let p' = { p with poly = Node.M.remove x p.poly } in @@ -81,13 +79,13 @@ let extract p = let remove n p = { p with poly = Node.M.remove n p.poly } -type kind = ONE | NODE of Node.t | PRODUCT +type kind = CST of A.t | NODE of A.t * Node.t | PRODUCT let classify p = - if Node.M.is_empty p.poly then ONE + if Node.M.is_empty p.poly then CST p.cst else if Node.M.is_num_elt 1 p.poly then let node, k = Node.M.choose p.poly in - if Q.equal Q.one k then NODE node else PRODUCT + if Q.equal Q.one k then NODE (p.cst, node) else PRODUCT else PRODUCT let monome x c = @@ -95,32 +93,34 @@ let monome x c = let is_one_node p = (* cst = 0 and one empty monome *) - if (* Q.equal Q.zero p.cst && *) Node.M.is_num_elt 1 p.poly then + if A.equal A.one p.cst && Node.M.is_num_elt 1 p.poly then let node, k = Node.M.choose p.poly in if Q.equal Q.one k then Some node else None else None let mult_cst c p1 = - assert (A.is_not_zero c); - { p1 with cst = A.mul p1.cst c } + if A.is_zero c then zero else { p1 with cst = A.mul p1.cst c } let power_cst c p1 = - if Q.equal Q.zero c then one + if A.is_zero p1.cst then zero + else if Q.equal Q.zero c then one else if Q.equal Q.one c then p1 else let poly_mult c m = Node.M.map (fun c1 -> Q.mul c c1) m in - if Q.equal Q.zero c then one - else { cst = A.positive_pow p1.cst c; poly = poly_mult c p1.poly } + { cst = A.positive_pow p1.cst c; poly = poly_mult c p1.poly } (** Warning Node.M.union can be used only for defining an operation [op] that verifies [op 0 p = p] and [op p 0 = p] *) let mul p1 p2 = - let poly_add m1 m2 = - Node.M.union (fun _ c1 c2 -> Q.none_zero (Q.add c1 c2)) m1 m2 - in - { cst = A.mul p1.cst p2.cst; poly = poly_add p1.poly p2.poly } + if A.is_zero p1.cst || A.is_zero p2.cst then zero + else + let poly_add m1 m2 = + Node.M.union (fun _ c1 c2 -> Q.none_zero (Q.add c1 c2)) m1 m2 + in + { cst = A.mul p1.cst p2.cst; poly = poly_add p1.poly p2.poly } let div p1 p2 = + assert (is_not_zero p2); let poly_sub m1 m2 = Node.M.union_merge (fun _ c1 c2 -> @@ -131,18 +131,24 @@ let div p1 p2 = in { cst = A.div p1.cst p2.cst; poly = poly_sub p1.poly p2.poly } +(** x * y^c *) let x_m_yc p1 p2 c = - assert (not (Q.equal c Q.zero)); - let f a b = Q.add a (Q.mul c b) in - let poly m1 m2 = - Node.M.union_merge - (fun _ c1 c2 -> - match c1 with - | None -> Some (Q.mul c c2) - | Some c1 -> Q.none_zero (f c1 c2)) - m1 m2 - in - { cst = A.mul p1.cst (A.positive_pow p2.cst c); poly = poly p1.poly p2.poly } + if Q.equal c Q.zero then p1 + else if Q.equal c Q.one then mul p1 p2 + else + let f a b = Q.add a (Q.mul c b) in + let poly m1 m2 = + Node.M.union_merge + (fun _ c1 c2 -> + match c1 with + | None -> Some (Q.mul c c2) + | Some c1 -> Q.none_zero (f c1 c2)) + m1 m2 + in + { + cst = A.mul p1.cst (A.positive_pow p2.cst c); + poly = poly p1.poly p2.poly; + } let xc_m_yc p1 c1 p2 c2 = let c1_iszero = Q.equal c1 Q.zero in @@ -181,7 +187,13 @@ let subst_node p x y = let subst p x s = let poly, q = Node.M.find_remove x p.poly in - match q with None -> (p, Q.zero) | Some q -> (x_m_yc { p with poly } s q, q) + match q with + | None -> (p, Q.zero) + | Some q -> + if is_zero s then ( + assert (Q.sign q > 0); + (zero, q)) + else (x_m_yc { p with poly } s q, q) let fold f acc p = Node.M.fold_left f acc p.poly let iter f p = Node.M.iter f p.poly @@ -200,7 +212,7 @@ let of_map cst m = let common pa pb = { - cst = A.min pa.cst pb.cst; + cst = A.one; poly = Node.M.inter (fun _ a b -> Q.none_zero (Q.min (Q.floor a) (Q.floor b))) diff --git a/src_colibri2/theories/LRA/product.mli b/src_colibri2/theories/LRA/product.mli index 07d094a9cc38e0ad2673c42fe44ea6cf2bf889fa..c75c8602b55543ac16e997c5dfa05cb900844583 100644 --- a/src_colibri2/theories/LRA/product.mli +++ b/src_colibri2/theories/LRA/product.mli @@ -21,25 +21,29 @@ open Colibri2_popop_lib open Colibri2_core -(** Polynome *) +(** Product of absolute values *) type t = private { cst : A.t; poly : Q.t Node.M.t } include Popop_stdlib.Datatype with type t := t val invariant : t -> bool +val cst : A.t -> t +val zero : t val one : t val monome : Node.t -> Q.t -> t val is_one_node : t -> Node.t option +val is_zero : t -> bool val remove : Node.t -> t -> t type extract = | One (** p = 1 *) + | Cst of A.t (** p = cst <> 1 *) | Var of Q.t * Node.t * t (** p = x^q * p' *) val extract : t -> extract -type kind = ONE | NODE of Node.t | PRODUCT +type kind = CST of A.t | NODE of A.t * Node.t | PRODUCT val classify : t -> kind val power_cst : Q.t -> t -> t diff --git a/src_colibri2/theories/LRA/realValue.ml b/src_colibri2/theories/LRA/realValue.ml index 5b4b7bedc712842b3e092749f8f1183b555578e8..1ac0b9dcb6b85fee36707fbfb2346f76df5066de 100644 --- a/src_colibri2/theories/LRA/realValue.ml +++ b/src_colibri2/theories/LRA/realValue.ml @@ -160,12 +160,24 @@ let () = | Q v -> if Z.equal v.den Z.one then Z.pp_print fmt v.num else Fmt.pf fmt "(/ %a %a)" Z.pp_print v.num Z.pp_print v.den - | A v -> Fmt.pf fmt "(epsilon x. %s)" (Ocaml_poly.A.to_string v)) + | A v -> ( + match Ocaml_poly.A.get_upolynomial_interval v with + | Ocaml_poly.A.IsPoint _ -> + assert false (* absurd: by invariant of Std.A *) + | Ocaml_poly.A.U (u, ba, a, b, bb) -> + let pp_bound fmt = function + | Ocaml_poly.A.Strict -> Fmt.pf fmt "<" + | Large -> Fmt.pf fmt "≤" + in + Fmt.pf fmt "(epsilon (x) (and (= (%s) 0) (%a %a x) (%a x %a)))" + (Ocaml_poly.UPolynomial.to_string u) + pp_bound ba Q.pp a pp_bound bb Q.pp b)) let cst' c = index ~basename:(Format.asprintf "%aR" A.pp c) c let cst c = node (cst' c) let zero = cst A.zero let one = cst A.one +let set d n a = Egraph.set_value d n (nodevalue (index a)) let positive_int_sequence = let open Base.Sequence.Generator in @@ -185,7 +197,7 @@ let int_sequence = Base.Sequence.shift_right int_sequence_without_zero Z.zero let q_sequence d = let open Interp.SeqLim in let+ num = of_seq d int_sequence and* den = of_seq d positive_int_sequence in - (** algebraic non rational are not generated *) + (* algebraic non rational are not generated *) A.of_q (Q.make num den) let init_ty d = @@ -500,7 +512,7 @@ module Check = struct !<(A.of_bigint (bitv n a)) | { app = { builtin = Builtin.Int2BV n }; tyargs = []; args; _ } -> let a = IArray.extract1_exn args in - extract n (A.to_z (!>a)) + extract n (A.to_z !>a) | _ -> `None let init d = diff --git a/src_colibri2/theories/LRA/realValue.mli b/src_colibri2/theories/LRA/realValue.mli index 092aab0c570a44b9bd31e29d1000474801c94b4a..f35a0ea38479568a6b6da3ba9211a27085013a58 100644 --- a/src_colibri2/theories/LRA/realValue.mli +++ b/src_colibri2/theories/LRA/realValue.mli @@ -28,6 +28,7 @@ end include Value.S with type s = A.t +val set : Egraph.wt -> Node.t -> A.t -> unit val cst' : A.t -> t val cst : A.t -> Node.t val zero : Node.t diff --git a/src_colibri2/theories/LRA/sign_product.ml b/src_colibri2/theories/LRA/sign_product.ml index 7ef561b8a77da8297dfa1eddf7775571daa39705..72fd2aae7385de7855919efb5b6969cadfdfa00c 100644 --- a/src_colibri2/theories/LRA/sign_product.ml +++ b/src_colibri2/theories/LRA/sign_product.ml @@ -33,7 +33,7 @@ module T = struct match q with | OddNonZero -> () | Even -> Fmt.string fmt "²" - | Odd -> Fmt.string fmt "³" + | Odd -> Fmt.string fmt "¹" let pp fmt = function | Zero -> Fmt.pf fmt "0" @@ -63,11 +63,13 @@ let is_one_node = function else None let one = NZ { cst = Pos; poly = Node.M.empty } - let minus_one = NZ { cst = Neg; poly = Node.M.empty } - let zero = Zero +let cst a = + let sgn = A.sign a in + if sgn = 0 then zero else if sgn > 0 then one else minus_one + let mul_cst c1 c2 = match (c1, c2) with Pos, Pos | Neg, Neg -> Pos | Neg, Pos | Pos, Neg -> Neg @@ -180,14 +182,16 @@ let common pa pb = | Pos, Pos | Neg, Pos | Pos, Neg -> Pos); } -type kind = PLUS_ONE | MINUS_ONE | ZERO | NODE of Node.t | PRODUCT +let mul_cst (pos : cst) a = match pos with Pos -> a | Neg -> A.neg a + +type kind = PLUS_ONE | MINUS_ONE | ZERO | NODE of (cst * Node.t) | PRODUCT let classify = function | T.Zero -> ZERO | T.NZ p -> if Node.M.is_empty p.poly then match p.cst with Pos -> PLUS_ONE | Neg -> MINUS_ONE - else if T.equal_cst Pos p.cst && Node.M.is_num_elt 1 p.poly then + else if Node.M.is_num_elt 1 p.poly then let node, v = Node.M.choose p.poly in - match v with Even -> PRODUCT | Odd | OddNonZero -> NODE node + match v with Even -> PRODUCT | Odd | OddNonZero -> NODE (p.cst, node) else PRODUCT diff --git a/src_colibri2/theories/LRA/sign_product.mli b/src_colibri2/theories/LRA/sign_product.mli index 9f3515e62e833a6581f52fa31615b331241049f7..16f878b1a4349cb65dd1557cc0b44b081fbaca60 100644 --- a/src_colibri2/theories/LRA/sign_product.mli +++ b/src_colibri2/theories/LRA/sign_product.mli @@ -18,12 +18,16 @@ (* for more details (enclosed in the file licenses/LGPLv2.1). *) (*************************************************************************) -(** Product of node sign with one constant 1 or -1 *) +(** Represents product of sign of nodes *) type t include Colibri2_popop_lib.Popop_stdlib.Datatype with type t := t +type cst = Pos | Neg + +val cst : A.t -> t + val of_one_node : Node.t -> t (** Build a value that represent one node *) @@ -55,19 +59,17 @@ val extract : t -> extracted type extract_cst = Pos | Neg | Zero val extract_cst : t -> extract_cst * t - val remove : Node.t -> t -> t val common : t -> t -> t (** Return the common part *) val one : t - val minus_one : t - val zero : t +val mul_cst : cst -> A.t -> A.t -type kind = PLUS_ONE | MINUS_ONE | ZERO | NODE of Node.t | PRODUCT +type kind = PLUS_ONE | MINUS_ONE | ZERO | NODE of (cst * Node.t) | PRODUCT val classify : t -> kind diff --git a/src_colibri2/theories/LRA/simplex.ml b/src_colibri2/theories/LRA/simplex.ml index 56edb11b5cfa9078b59d10cb173b86ec3148d1b5..cd582e8806071b0be468844e2e99c3fb623a01b8 100644 --- a/src_colibri2/theories/LRA/simplex.ml +++ b/src_colibri2/theories/LRA/simplex.ml @@ -526,7 +526,7 @@ let converter d (f : Ground.t) = Egraph.register d n; Dom_interval.Propagate.sub d ~a ~b ~r:n; Dom_polynome.assume_equality d n p; - Dom_product.propagate_a_cb d ~a ~b ~c:Q.minus_one ~res:n); + Dom_product.propagate_a_cb d ~a ~b ~c:A.minus_one ~res:n); Dom_interval.Propagate.cmp d (build builtin) ~r:(Ground.node f) ~a:n ~b:RealValue.zero; n diff --git a/src_colibri2/theories/LRA/stages/compare_stage/compare_stage.ml b/src_colibri2/theories/LRA/stages/compare_stage/compare_stage.ml index 10eab50b40384b19a6db3af11b2f2782ea5b4170..a12b4cdb103242069f2b5ebcd954a7d0cc8a31e2 100644 --- a/src_colibri2/theories/LRA/stages/compare_stage/compare_stage.ml +++ b/src_colibri2/theories/LRA/stages/compare_stage/compare_stage.ml @@ -24,7 +24,7 @@ open Colibri2_stdlib.Std module type S = sig include Interval_sig.S - type split_heuristic = Singleton of Q.t | Splitted of t * t | NotSplitted + type split_heuristic = Singleton of A.t | Splitted of t * t | NotSplitted val split_heuristic : t -> split_heuristic end @@ -81,13 +81,13 @@ module Cmp (S1 : S) (S2 : S) : S = struct let mem q u = let b1 = S1.mem q u.s1 in let b2 = S2.mem q u.s2 in - print_if2 ((not b1) && b2) Q.pp q pp u "mem"; + print_if2 ((not b1) && b2) A.pp q pp u "mem"; b2 let absent q u = let b1 = S1.absent q u.s1 in let b2 = S2.absent q u.s2 in - print_if2 (b1 && not b2) Q.pp q pp u "absent"; + print_if2 (b1 && not b2) A.pp q pp u "absent"; b2 type is_comparable = S2.is_comparable = @@ -114,17 +114,11 @@ module Cmp (S1 : S) (S2 : S) : S = struct b2 let reals = mk S1.reals S2.reals - let integers = mk S1.integers S2.integers - let zero = mk S1.zero S2.zero - let mult_cst q s = { s1 = S1.mult_cst q s.s1; s2 = S2.mult_cst q s.s2 } - let add_cst q s = { s1 = S1.add_cst q s.s1; s2 = S2.add_cst q s.s2 } - let add s u = { s1 = S1.add s.s1 u.s1; s2 = S2.add s.s2 u.s2 } - let minus s u = { s1 = S1.minus s.s1 u.s1; s2 = S2.minus s.s2 u.s2 } let inter s u = @@ -143,39 +137,28 @@ module Cmp (S1 : S) (S2 : S) : S = struct match (b1, b2) with | _, None -> None | None, Some _ -> - print_if2 true pp s Q.pp u "except"; + print_if2 true pp s A.pp u "except"; None | Some s1, Some s2 -> Some { s1; s2 } let union s u = { s1 = S1.union s.s1 u.s1; s2 = S2.union s.s2 u.s2 } - let singleton q = { s1 = S1.singleton q; s2 = S2.singleton q } let from_convexe_hull b1 = { s1 = S1.from_convexe_hull b1; s2 = S2.from_convexe_hull b1 } - type split_heuristic = Singleton of Q.t | Splitted of t * t | NotSplitted + type split_heuristic = Singleton of A.t | Splitted of t * t | NotSplitted let split_heuristic _ = assert false - let gt q = { s1 = S1.gt q; s2 = S2.gt q } - let lt q = { s1 = S1.lt q; s2 = S2.lt q } - let le q = { s1 = S1.le q; s2 = S2.le q } - let ge q = { s1 = S1.ge q; s2 = S2.ge q } - let gt' q = { s1 = S1.gt' q.s1; s2 = S2.gt' q.s2 } - let lt' q = { s1 = S1.lt' q.s1; s2 = S2.lt' q.s2 } - let le' q = { s1 = S1.le' q.s1; s2 = S2.le' q.s2 } - let ge' q = { s1 = S1.ge' q.s1; s2 = S2.ge' q.s2 } - let choose s = S2.choose s.s2 - let invariant s = S1.invariant s.s1 && S2.invariant s.s2 let get_convexe_hull s = diff --git a/src_colibri2/theories/LRA/stages/stage1/interval.ml b/src_colibri2/theories/LRA/stages/stage1/interval.ml index fc23775b03b81c2a1cd827e02128f145b3667a67..dba566cc0f2aeb34548b8e267b4bb1a0b621d07d 100644 --- a/src_colibri2/theories/LRA/stages/stage1/interval.ml +++ b/src_colibri2/theories/LRA/stages/stage1/interval.ml @@ -21,17 +21,37 @@ open Colibri2_popop_lib open Colibri2_stdlib.Std +let compare_ord a b : Colibrics_lib.Ord.t = + let c = A.compare a b in + Int.(if c = 0 then Eq else if c <= 0 then Lt else Gt) + module Convexe = struct - include Colibrics_lib.Interval.Convexe + module Convexe = Colibrics_lib.Interval.Convexe.Make (struct + include A + + let infix_eq = equal + let infix_lseq = le + let infix_ls = lt + let infix_pl = add + let infix_mn = sub + let infix_as = mul + let infix_sl = div + let compare = compare_ord + let of_int = A.of_bigint + let floor a = A.to_z (A.floor a) + let ceil a = A.to_z (A.ceil a) + end) + + include Convexe include struct type bound = Colibrics_lib.Interval.Bound.t = Strict | Large [@@deriving eq, ord, hash] - type t = Colibrics_lib.Interval.Convexe.t = - | Finite of Q.t * bound * Q.t * bound - | InfFinite of Q.t * bound - | FiniteInf of Q.t * bound + type t = Convexe.t = + | Finite of A.t * bound * A.t * bound + | InfFinite of A.t * bound + | FiniteInf of A.t * bound | Inf [@@deriving eq, ord, hash] end @@ -47,25 +67,21 @@ module Convexe = struct in match t with | Finite (q1, b1, q2, b2) -> - Format.fprintf fmt "%a%a;%a%a" print_bound_left b1 Q.pp q1 Q.pp q2 + Format.fprintf fmt "%a%a;%a%a" print_bound_left b1 A.pp q1 A.pp q2 print_bound_right b2 | InfFinite (q1, b1) -> - Format.fprintf fmt "]-∞;%a%a" Q.pp q1 print_bound_right b1 + Format.fprintf fmt "]-∞;%a%a" A.pp q1 print_bound_right b1 | FiniteInf (q1, b1) -> - Format.fprintf fmt "%a%a;+∞[" print_bound_left b1 Q.pp q1 + Format.fprintf fmt "%a%a;+∞[" print_bound_left b1 A.pp q1 | Inf -> Format.fprintf fmt "]-∞;+∞[" include Popop_stdlib.MkDatatype (struct type nonrec t = t let equal = equal - let compare = compare - let hash = hash - let hash_fold_t = hash_fold_t - let pp = pp end) @@ -79,7 +95,7 @@ module Convexe = struct | Strict -> Strict | Large -> Large - let get_convexe_hull : t -> (Q.t * Bound.t) option * (Q.t * Bound.t) option = + let get_convexe_hull : t -> (A.t * Bound.t) option * (A.t * Bound.t) option = function | Finite (q1, b1, q2, b2) -> (Some (q1, of_bound b1), Some (q2, of_bound b2)) | InfFinite (q, b) -> (None, Some (q, of_bound b)) @@ -96,26 +112,26 @@ module Convexe = struct (* let nb_incr = 100 * let z_nb_incr = Z.of_int nb_incr * let choose_rnd rnd = function - * | {minv;maxv} when Q.equal Q.minus_inf minv -> - * Q.make (Z.sub (Z.of_int (rnd nb_incr)) (Q.to_bigint maxv)) Z.one - * | {minv;maxv} when Q.equal Q.inf maxv -> - * Q.make (Z.add (Z.of_int (rnd nb_incr)) (Q.to_bigint minv)) Z.one - * | {minv;maxv} when Q.equal minv maxv -> minv + * | {minv;maxv} when A.equal A.minus_inf minv -> + * A.make (Z.sub (Z.of_int (rnd nb_incr)) (A.to_bigint maxv)) Z.one + * | {minv;maxv} when A.equal A.inf maxv -> + * A.make (Z.add (Z.of_int (rnd nb_incr)) (A.to_bigint minv)) Z.one + * | {minv;maxv} when A.equal minv maxv -> minv * | {minv;maxv} -> - * let d = Q.sub maxv minv in + * let d = A.sub maxv minv in * let i = 1 + rnd (nb_incr - 2) in - * let x = Q.make (Z.of_int i) (Z.of_int 100) in - * let q = Q.add minv (Q.mul x d) in - * assert (Q.lt minv q); - * assert (Q.lt q maxv); + * let x = A.make (Z.of_int i) (Z.of_int 100) in + * let q = A.add minv (A.mul x d) in + * assert (A.lt minv q); + * assert (A.lt q maxv); * q * * let get_convexe_hull e = * let mk v b = - * match Q.classify v with - * | Q.ZERO | Q.NZERO -> Some (v,b) - * | Q.INF | Q.MINF -> None - * | Q.UNDEF -> assert false in + * match A.classify v with + * | A.ZERO | A.NZERO -> Some (v,b) + * | A.INF | A.MINF -> None + * | A.UNDEF -> assert false in * mk e.minv e.minb, mk e.maxv e.maxb * * let is_Large = function @@ -124,18 +140,18 @@ module Convexe = struct * * let inter_with_integer t = * let t = { - * minb = if Q.equal Q.minus_inf t.minv then Large else Strict; + * minb = if A.equal A.minus_inf t.minv then Large else Strict; * minv = - * if is_Large t.minb && Q.is_integer t.minv - * then Q.add t.minv Q.one - * else Q.ceil t.minv; - * maxb = if Q.equal Q.inf t.maxv then Large else Strict; + * if is_Large t.minb && A.is_integer t.minv + * then A.add t.minv A.one + * else A.ceil t.minv; + * maxb = if A.equal A.inf t.maxv then Large else Strict; * maxv = - * if is_Large t.maxb && Q.is_integer t.maxv - * then Q.sub t.maxv Q.one - * else Q.floor t.maxv; + * if is_Large t.maxb && A.is_integer t.maxv + * then A.sub t.maxv A.one + * else A.floor t.maxv; * } in - * if Q.lt t.maxv t.minv then None else begin + * if A.lt t.maxv t.minv then None else begin * assert (invariant t); * Some t * end *) @@ -144,15 +160,15 @@ end module ConvexeWithDivider = struct (** modulo = 0 -> no modulo *) - type t = { c: Convexe.t; divider: Q.t } + type t = { c: Convexe.t; divider: A.t } [@@ deriving eq, ord, hash] let pp fmt t = Format.fprintf fmt "%a%t" Convexe.pp t.c (fun fmt -> - if not (Q.equal t.divider Q.zero) - then Format.fprintf fmt "%%%a" Q.pp t.divider + if not (A.equal t.divider A.zero) + then Format.fprintf fmt "%%%a" A.pp t.divider ) @@ -164,41 +180,41 @@ module ConvexeWithDivider = struct let invariant e = Convexe.invariant e.c && - match Q.classify e.divider with - | Q.ZERO -> true - | Q.NZERO -> + match A.classify e.divider with + | A.ZERO -> true + | A.NZERO -> let large_finite_bound b v inf = match b with - | Strict -> Q.equal v inf + | Strict -> A.equal v inf | Large -> true in - large_finite_bound e.c.minb e.c.minv Q.minus_inf && - large_finite_bound e.c.maxb e.c.maxv Q.inf - | Q.INF | Q.MINF | Q.UNDEF -> false + large_finite_bound e.c.minb e.c.minv A.minus_inf && + large_finite_bound e.c.maxb e.c.maxv A.inf + | A.INF | A.MINF | A.UNDEF -> false let tighten_bounds e = - if Q.equal e.divider Q.zero then Some e + if A.equal e.divider A.zero then Some e else let round ~add ~divisible_up_to ~minus_inf ~minv ~minb = - if Q.equal minv minus_inf + if A.equal minv minus_inf then minb, minv else let minv' = divisible_up_to minv e.divider in Large, begin match minb with - | Strict when Q.equal minv' minv -> add minv e.divider + | Strict when A.equal minv' minv -> add minv e.divider | _ -> minv' end in let minb,minv = - round ~add:Q.add + round ~add:A.add ~divisible_up_to:Colibrics_lib.QUtils.divisible_up_to - ~minus_inf:Q.minus_inf ~minv:e.c.Convexe.minv ~minb:e.c.Convexe.minb + ~minus_inf:A.minus_inf ~minv:e.c.Convexe.minv ~minb:e.c.Convexe.minb in let maxb,maxv = - round ~add:Q.sub ~divisible_up_to:Colibrics_lib.QUtils.divisible_down_to - ~minus_inf:Q.inf ~minv:e.c.Convexe.maxv ~minb:e.c.Convexe.maxb + round ~add:A.sub ~divisible_up_to:Colibrics_lib.QUtils.divisible_down_to + ~minus_inf:A.inf ~minv:e.c.Convexe.maxv ~minb:e.c.Convexe.maxb in - if Q.lt maxv minv + if A.lt maxv minv then None else Some { c = { minv; minb; maxv; maxb }; divider = e.divider } @@ -222,14 +238,14 @@ module ConvexeWithDivider = struct let is_distinct e1 e2 = Convexe.is_distinct e1.c e2.c let divisible q d = - Q.equal d Q.zero || + A.equal d A.zero || Colibrics_lib.QUtils.divisible q d let is_included e1 e2 = Convexe.is_included e1.c e2.c && begin - Q.equal e2.divider Q.zero || - (not (Q.equal e1.divider Q.zero) && + A.equal e2.divider A.zero || + (not (A.equal e1.divider A.zero) && divisible e1.divider e2.divider ) end @@ -239,12 +255,12 @@ module ConvexeWithDivider = struct end let mult_cst q e = - assert (Q.is_real q); + assert (A.is_real q); { c = Convexe.mult_cst q e.c; - divider = Q.abs (Colibrics_lib.QUtils.mult_cst_divisible e.divider q) } + divider = A.abs (Colibrics_lib.QUtils.mult_cst_divisible e.divider q) } let add_cst q e = - if Q.equal q Q.zero then e + if A.equal q A.zero then e else let divider = Colibrics_lib.QUtils.union_divisible q e.divider in { c = Convexe.add_cst q e.c; divider } @@ -257,7 +273,7 @@ module ConvexeWithDivider = struct { c = Convexe.minus e1.c e2.c; divider = Colibrics_lib.QUtils.union_divisible e1.divider e2.divider } - let div_unknown c = { c; divider = Q.zero } + let div_unknown c = { c; divider = A.zero } let gt q = div_unknown (Convexe.gt q) @@ -286,11 +302,11 @@ module ConvexeWithDivider = struct let reals = div_unknown Convexe.reals let () = assert (invariant reals) - let is_reals e = Convexe.is_reals e.c && Q.equal e.divider Q.zero + let is_reals e = Convexe.is_reals e.c && A.equal e.divider A.zero let choose e = let q = Convexe.choose e.c in - if Q.equal e.divider Q.zero + if A.equal e.divider A.zero then q else Colibrics_lib.QUtils.divisible_down_to q e.divider @@ -325,25 +341,25 @@ end module ConvexeWithExceptions = struct - type t = {con: Convexe.t; exc: Q.S.t} + type t = {con: Convexe.t; exc: A.S.t} [@@deriving eq, ord] let pp fmt x = - if Q.S.is_empty x.exc then Convexe.pp fmt x.con + if A.S.is_empty x.exc then Convexe.pp fmt x.con else - let aux fmt m = Q.S.elements m - |> Format.(list ~sep:(const char ',') Q.pp) fmt + let aux fmt m = A.S.elements m + |> Format.(list ~sep:(const char ',') A.pp) fmt in Format.fprintf fmt "%a \ {%a}" Convexe.pp x.con aux x.exc let invariant x = Convexe.invariant x.con && - ( Q.S.is_empty x.exc || - ( Q.lt x.con.Convexe.minv (Q.S.min_elt x.exc) && - Q.lt (Q.S.max_elt x.exc) x.con.Convexe.maxv )) + ( A.S.is_empty x.exc || + ( A.lt x.con.Convexe.minv (A.S.min_elt x.exc) && + A.lt (A.S.max_elt x.exc) x.con.Convexe.maxv )) let hash e = 53 * (Convexe.hash e.con) + - 31 * (Q.S.fold_left (fun acc x -> 3*acc + Q.hash x) 5 e.exc) + 31 * (A.S.fold_left (fun acc x -> 3*acc + A.hash x) 5 e.exc) include Popop_stdlib.MkDatatype(struct type nonrec t = t @@ -351,100 +367,100 @@ module ConvexeWithExceptions = struct let hash = hash let pp = pp end) - let of_con con = {con;exc=Q.S.empty} + let of_con con = {con;exc=A.S.empty} let reals = of_con Convexe.reals let from_convexe f x = of_con (f x) let singleton = from_convexe Convexe.singleton - let zero = singleton Q.zero + let zero = singleton A.zero let le = from_convexe Convexe.le let lt = from_convexe Convexe.lt let ge = from_convexe Convexe.ge let gt = from_convexe Convexe.gt let is_singleton e = Convexe.is_singleton e.con - let is_reals e = Convexe.is_reals e.con && Q.S.is_empty e.exc - let mem x e = Convexe.mem x e.con && not (Q.S.mem x e.exc) + let is_reals e = Convexe.is_reals e.con && A.S.is_empty e.exc + let mem x e = Convexe.mem x e.con && not (A.S.mem x e.exc) let except e x = match Convexe.except e.con x with | None -> None | Some con -> if Convexe.mem x con - then Some {con; exc = Q.S.add x e.exc} + then Some {con; exc = A.S.add x e.exc} else Some {e with con} let union e1 e2 = {con=Convexe.union e1.con e2.con; - exc=Q.S.inter e1.exc e2.exc} + exc=A.S.inter e1.exc e2.exc} let normalize {con;exc} = match Convexe.is_singleton con with - | Some s -> if Q.S.mem s exc then None else Some (of_con con) + | Some s -> if A.S.mem s exc then None else Some (of_con con) | None -> - let _,has_min,exc = Q.S.split con.Convexe.minv exc in - let _,has_max,exc = Q.S.split con.Convexe.maxv exc in + let _,has_min,exc = A.S.split con.Convexe.minv exc in + let _,has_max,exc = A.S.split con.Convexe.maxv exc in Some {exc; con = {con with minb = if has_min then Strict else con.minb; maxb = if has_max then Strict else con.maxb }} let inter e1 e2 = match Convexe.inter e1.con e2.con with | None -> None | Some con -> - normalize { con; exc = Q.S.union e1.exc e2.exc } + normalize { con; exc = A.S.union e1.exc e2.exc } let is_comparable e1 e2 = Convexe.is_comparable e1.con e2.con (* || - * not (Q.S.equal e1.exc e2.exc) *) + * not (A.S.equal e1.exc e2.exc) *) let is_distinct e1 e2 = Convexe.is_distinct e1.con e2.con (* || - * not (Q.S.equal e1.exc e2.exc) *) + * not (A.S.equal e1.exc e2.exc) *) let is_included e1 e2 = Convexe.is_included e1.con e2.con || - Q.S.subset e1.exc e2.exc + A.S.subset e1.exc e2.exc let add e1 e2 = match Convexe.is_singleton e1.con, e1, Convexe.is_singleton e2.con, e2 with - | Some s1, _, Some s2, _ -> singleton (Q.add s1 s2) + | Some s1, _, Some s2, _ -> singleton (A.add s1 s2) | None, e, Some s, _ | Some s, _, None, e -> - {con=Convexe.add_cst s e.con; exc = Q.S.translate (Q.add s) e.exc} + {con=Convexe.add_cst s e.con; exc = A.S.translate (A.add s) e.exc} | _ -> of_con (Convexe.add e1.con e2.con) let add_cst s e = - {con = Convexe.add_cst s e.con; exc = Q.S.translate (Q.add s) e.exc} + {con = Convexe.add_cst s e.con; exc = A.S.translate (A.add s) e.exc} let minus e1 e2 = match Convexe.is_singleton e1.con, e1, Convexe.is_singleton e2.con, e2 with | Some s1, _, Some s2, _ -> - singleton (Q.sub s1 s2) + singleton (A.sub s1 s2) | None, e, Some s, _ -> - { con = Convexe.add_cst (Q.neg s) e.con; - exc = Q.S.translate (fun x -> Q.sub x s) e.exc } + { con = Convexe.add_cst (A.neg s) e.con; + exc = A.S.translate (fun x -> A.sub x s) e.exc } | Some s, _, None, e -> { con = Convexe.add_cst s e.con; exc = - Q.S.fold_left (fun acc exc -> Q.S.add (Q.sub s exc) acc) - Q.S.empty e.exc } + A.S.fold_left (fun acc exc -> A.S.add (A.sub s exc) acc) + A.S.empty e.exc } | _ -> of_con (Convexe.add e1.con e2.con) let mult_cst s e = {con = Convexe.mult_cst s e.con; exc = - Q.S.fold_left (fun acc exc -> Q.S.add (Q.mul s exc) acc) - Q.S.empty e.exc } + A.S.fold_left (fun acc exc -> A.S.add (A.mul s exc) acc) + A.S.empty e.exc } let choose e = - let con = if Q.S.is_empty e.exc then e.con + let con = if A.S.is_empty e.exc then e.con else (** by the invariant the intersection must succeed *) Opt.get_exn Impossible - (Convexe.inter e.con (Convexe.lt (Q.S.min_elt e.exc))) in + (Convexe.inter e.con (Convexe.lt (A.S.min_elt e.exc))) in Convexe.choose con let choose_rnd rnd e = - let con = if Q.S.is_empty e.exc then e.con + let con = if A.S.is_empty e.exc then e.con else (** by the invariant the intersection must succeed *) Opt.get_exn Impossible - (Convexe.inter e.con (Convexe.lt (Q.S.min_elt e.exc))) in + (Convexe.inter e.con (Convexe.lt (A.S.min_elt e.exc))) in Convexe.choose_rnd rnd con let get_convexe_hull e = Convexe.get_convexe_hull e.con @@ -461,15 +477,15 @@ module Union = struct | l -> let rec aux minb' minv' (l:t) = match minb', l with - | Large, [] when Q.equal minv' Q.inf -> false + | Large, [] when A.equal minv' A.inf -> false | _ , [] -> true - | _ , {minb=Large ; minv}::_ when Q.compare minv minv' <= 0 -> false - | Large , {minb=Strict; minv}::_ when Q.compare minv minv' <= 0 -> false - | Strict, {minb=Strict; minv}::_ when Q.compare minv minv' < 0 -> false + | _ , {minb=Large ; minv}::_ when A.compare minv minv' <= 0 -> false + | Large , {minb=Strict; minv}::_ when A.compare minv minv' <= 0 -> false + | Strict, {minb=Strict; minv}::_ when A.compare minv minv' < 0 -> false | _, ({maxv;maxb} as e)::l -> Convexe.invariant e && aux maxb maxv l in - aux Strict Q.minus_inf l + aux Strict A.minus_inf l let pp fmt l = Format.list ~sep:Format.(const string "∪") Convexe.pp fmt l @@ -494,7 +510,7 @@ module Union = struct let mem x e = List.exists (fun e -> Convexe.mem x e) e let lower_min_max e1 l1 t1 e2 l2 t2 = - let c = Q.compare e1.Convexe.minv e2.Convexe.minv in + let c = A.compare e1.Convexe.minv e2.Convexe.minv in match e1.Convexe.minb, e2.Convexe.minb with | Strict, Large when c = 0 -> e2,l2,t2,e1,l1,t1 @@ -510,7 +526,7 @@ module Union = struct let emin,lmin,_,emax,lmax,tmax = lower_min_max e1 l1 t1 e2 l2 t2 in (** look for an intersection *) - let c = Q.compare emin.maxv emax.minv in + let c = A.compare emin.maxv emax.minv in match emin.maxb, emax.minb with (** no intersection *) | Strict, Strict @@ -520,7 +536,7 @@ module Union = struct | _ -> (** intersection *) (** look for inclusion *) - let c = Q.compare emax.maxv emin.maxv in + let c = A.compare emax.maxv emin.maxv in match emax.maxb, emin.maxb with (** max included in min *) | Strict, Strict | Large, Large | Strict, Large @@ -545,7 +561,7 @@ module Union = struct (** order by the minimum *) let emin,lmin,tmin,emax,lmax,tmax = lower_min_max e1 l1 t1 e2 l2 t2 in (** look for an intersection *) - let c = Q.compare emin.maxv emax.minv in + let c = A.compare emin.maxv emax.minv in match emin.maxb, emax.minb with (** no intersection *) | Strict, Strict | Strict, Large | Large, Strict @@ -555,7 +571,7 @@ module Union = struct | _ -> (** intersection *) (** look for inclusion *) - let c = Q.compare emax.maxv emin.maxv in + let c = A.compare emax.maxv emin.maxv in match emax.maxb, emin.maxb with (** max included in min *) | Strict, Strict | Large, Large | Strict, Large @@ -594,7 +610,7 @@ module Union = struct assert (invariant t); t - let zero = singleton Q.zero + let zero = singleton A.zero let reals = [Convexe.reals] let is_reals = function @@ -633,9 +649,9 @@ module Union = struct let mult_neg q = List.rev_map (Convexe.mult_neg q) let mult_cst q t = - assert (Q.is_real q); - let c = Q.sign q in - if c = 0 then singleton Q.zero + assert (A.is_real q); + let c = A.sign q in + if c = 0 then singleton A.zero else if c > 0 then mult_pos q t else mult_neg q t @@ -648,13 +664,13 @@ module Union = struct let cons (t:Convexe.t) (l:t) = match t.maxb, l with | _,[] -> [t] - | Strict, ({minb=Strict} as e)::_ when Q.compare t.maxv e.minv <= 0 -> + | Strict, ({minb=Strict} as e)::_ when A.compare t.maxv e.minv <= 0 -> t::l - | _, e::_ when Q.compare t.maxv e.minv < 0 -> + | _, e::_ when A.compare t.maxv e.minv < 0 -> t::l | _, e::l -> - assert (Q.compare t.minv e.minv < 0); - assert (Q.compare t.maxv e.maxv < 0); + assert (A.compare t.minv e.minv < 0); + assert (A.compare t.maxv e.maxv < 0); {minb=t.minb; minv = t.maxv; maxv = e.maxv; maxb = e.maxb}::l let rec add_intemaxval t = function @@ -668,28 +684,28 @@ module Union = struct | None,_, None,_ -> List.fold_left (fun acc t -> union acc (add_intemaxval t t2)) [] t1 | Some q, _, None, t - | None, t, Some q, _ when Q.equal Q.zero q -> t + | None, t, Some q, _ when A.equal A.zero q -> t | Some q, _, None, t | None, t, Some q, _ -> add_cst q t - | Some q1,_, Some q2,_ -> singleton (Q.add q1 q2) + | Some q1,_, Some q2,_ -> singleton (A.add q1 q2) in assert ( invariant res ); res let minus t1 t2 = - add t1 (mult_neg Q.minus_one t2) + add t1 (mult_neg A.minus_one t2) (** TODO better heuristic *) let choose = function | [] -> assert false | {Convexe.minb=Large;minv}::_ -> minv | {maxb=Large;maxv}::_ -> maxv - | {minv;maxv}::_ when Q.equal Q.minus_inf minv -> - Q.make (Z.sub Z.one (Q.to_bigint maxv)) Z.one - | {minv;maxv}::_ when Q.equal Q.inf maxv -> - Q.make (Z.add Z.one (Q.to_bigint minv)) Z.one - | {minv;maxv}::_ -> Q.add maxv (Q.div_2exp (Q.sub minv maxv) 1) + | {minv;maxv}::_ when A.equal A.minus_inf minv -> + A.make (Z.sub Z.one (A.to_bigint maxv)) Z.one + | {minv;maxv}::_ when A.equal A.inf maxv -> + A.make (Z.add Z.one (A.to_bigint minv)) Z.one + | {minv;maxv}::_ -> A.add maxv (A.div_2exp (A.sub minv maxv) 1) let choose_rnd rnd l = Convexe.choose_rnd rnd (List.nth l (rnd (List.length l))) @@ -718,8 +734,8 @@ module ConvexeInfo (Info: sig val nothing: t end) = struct - type t = { minb : bound; minv: Q.t; mini: Info.t; - maxv: Q.t; maxb: bound; maxi: Info.t } + type t = { minb : bound; minv: A.t; mini: Info.t; + maxv: A.t; maxb: bound; maxi: Info.t } let get_info t = t.mini, t.maxi let to_convexe t = { Convexe.minb = t.minb; minv = t.minv; @@ -735,13 +751,13 @@ module ConvexeInfo (Info: sig let print_bound_right fmt = function | Large -> Format.fprintf fmt "]" | Strict -> Format.fprintf fmt "[" in - if equal_bound t.minb Large && equal_bound t.maxb Large && Q.equal t.minv t.maxv - then Format.fprintf fmt "{%a}" Q.pp t.minv + if equal_bound t.minb Large && equal_bound t.maxb Large && A.equal t.minv t.maxv + then Format.fprintf fmt "{%a}" A.pp t.minv else Format.fprintf fmt "%a%a;%a%a" print_bound_left t.minb - Q.pp t.minv - Q.pp t.maxv + A.pp t.minv + A.pp t.maxv print_bound_right t.maxb let compare e1 e2 = @@ -750,11 +766,11 @@ module ConvexeInfo (Info: sig | Large, Strict -> -1 | _ -> 0 in - let c = Q.compare e1.minv e2.minv in + let c = A.compare e1.minv e2.minv in if c = 0 then c else let c = compare_bound e1.minb e2.minb in if c = 0 then c else - let c = Q.compare e1.maxv e2.maxv in + let c = A.compare e1.maxv e2.maxv in if c = 0 then c else let c = compare_bound e1.maxb e2.maxb in if c = 0 then c else @@ -765,7 +781,7 @@ module ConvexeInfo (Info: sig let equal e1 e2 = Equal.physical e1.minb e2.minb && Equal.physical e1.maxb e2.maxb && - Q.equal e1.minv e2.minv && Q.equal e1.maxv e2.maxv && + A.equal e1.minv e2.minv && A.equal e1.maxv e2.maxv && Info.equal e1.mini e2.mini && Info.equal e1.maxi e2.maxi let hash e = @@ -784,15 +800,15 @@ module ConvexeInfo (Info: sig end) let invariant e = - not (Q.equal e.minv Q.inf) && - not (Q.equal e.maxv Q.minus_inf) && - let c = Q.compare e.minv e.maxv in + not (A.equal e.minv A.inf) && + not (A.equal e.maxv A.minus_inf) && + let c = A.compare e.minv e.maxv in if c = 0 then equal_bound e.minb Large && equal_bound e.maxb Large else c < 0 let is_singleton = function - | {minv;maxv} when Q.equal minv maxv -> Some minv + | {minv;maxv} when A.equal minv maxv -> Some minv | _ -> None let singleton ~min_info ?(max_info=min_info) q = @@ -801,25 +817,25 @@ module ConvexeInfo (Info: sig t let except e x = - let is_min = Q.equal e.minv x in - let is_max = Q.equal e.maxv x in + let is_min = A.equal e.minv x in + let is_max = A.equal e.maxv x in if is_min && is_max then None else if is_min then Some {e with minb=Strict} - else if Q.equal e.maxv x + else if A.equal e.maxv x then Some {e with maxb=Strict} else Some e let lower_min_max e1 e2 = - let c = Q.compare e1.minv e2.minv in + let c = A.compare e1.minv e2.minv in match e1.minb, e2.minb with | Strict, Large when c = 0 -> e2,e1, false | _ when c <= 0 -> e1,e2, true | _ -> e2,e1, false let bigger_min_max e1 e2 = - let c = Q.compare e1.maxv e2.maxv in + let c = A.compare e1.maxv e2.maxv in match e1.maxb, e2.maxb with | Strict, Large when c = 0 -> e1,e2 | _ when c >= 0 -> e2,e1 @@ -830,7 +846,7 @@ module ConvexeInfo (Info: sig (** order by the minimum *) let emin,emax,same = lower_min_max e1 e2 in (** emin.minv <= emax.minv *) (** look for inclusion *) - let c = Q.compare emin.maxv emax.minv in + let c = A.compare emin.maxv emax.minv in match emin.maxb, emax.minb with (** emin.minv <? e1 <? emin.maxv < emax.minv <? e2 <? emax.maxv *) | _ when c < 0 -> if same then `Lt else `Gt @@ -855,32 +871,32 @@ module ConvexeInfo (Info: sig let mem x e = (match e.minb with - | Strict -> Q.lt e.minv x - | Large -> Q.le e.minv x) + | Strict -> A.lt e.minv x + | Large -> A.le e.minv x) && (match e.maxb with - | Strict -> Q.lt x e.maxv - | Large -> Q.le x e.maxv) + | Strict -> A.lt x e.maxv + | Large -> A.le x e.maxv) let mult_pos q e = - {e with minv = Q.mul e.minv q; maxv = Q.mul e.maxv q} + {e with minv = A.mul e.minv q; maxv = A.mul e.maxv q} let mult_neg q e = { minb = e.maxb; maxb = e.minb; - minv = Q.mul e.maxv q; - maxv = Q.mul e.minv q; + minv = A.mul e.maxv q; + maxv = A.mul e.minv q; mini= e.mini; maxi = e.maxi } let mult_cst q e = - assert (Q.is_real q); - let c = Q.sign q in + assert (A.is_real q); + let c = A.sign q in if c = 0 then invalid_arg "mult_cst q = 0" else if c > 0 then mult_pos q e else mult_neg q e let add_cst q e = - {e with minv = Q.add e.minv q; maxv = Q.add e.maxv q} + {e with minv = A.add e.minv q; maxv = A.add e.maxv q} let mult_bound b1 b2 = match b1, b2 with @@ -889,29 +905,29 @@ module ConvexeInfo (Info: sig let add ~min_info ?(max_info=min_info) e1 e2 = let t = {minb = mult_bound e1.minb e2.minb; - minv = Q.add e1.minv e2.minv; - maxv = Q.add e1.maxv e2.maxv; + minv = A.add e1.minv e2.minv; + maxv = A.add e1.maxv e2.maxv; maxb = mult_bound e1.maxb e2.maxb; mini = min_info ; maxi = max_info} in assert (invariant t); t let minus ~min_info ?max_info e1 e2 = - add ~min_info ?max_info e1 (mult_neg Q.minus_one e2) + add ~min_info ?max_info e1 (mult_neg A.minus_one e2) let gt ~min_info q = - let t = {minb=Strict; minv = q; mini = min_info; maxv = Q.inf; maxb= Strict; maxi = Info.nothing} in + let t = {minb=Strict; minv = q; mini = min_info; maxv = A.inf; maxb= Strict; maxi = Info.nothing} in assert (invariant t); t let ge ~min_info q = - let t = {minb=Large; minv = q; mini = min_info; maxv = Q.inf; maxb= Strict; maxi = Info.nothing} in + let t = {minb=Large; minv = q; mini = min_info; maxv = A.inf; maxb= Strict; maxi = Info.nothing} in assert (invariant t); t let lt ~max_info q = - let t = {minb=Strict; minv = Q.minus_inf; mini = Info.nothing; maxv = q; maxb= Strict; maxi = max_info} in + let t = {minb=Strict; minv = A.minus_inf; mini = Info.nothing; maxv = q; maxb= Strict; maxi = max_info} in assert (invariant t); t let le ~max_info q = - let t = {minb=Strict; minv = Q.minus_inf; mini = Info.nothing; maxv = q; maxb= Large; maxi = max_info} in + let t = {minb=Strict; minv = A.minus_inf; mini = Info.nothing; maxv = q; maxb= Large; maxi = max_info} in assert (invariant t); t let union e1 e2 = @@ -931,7 +947,7 @@ module ConvexeInfo (Info: sig in if compare_bounds_inf_sup min max > 0 then None - else if Q.equal minv maxv && equal_bound minb Large && equal_bound maxb Large + else if A.equal minv maxv && equal_bound minb Large && equal_bound maxb Large then Some (singleton ~min_info:mini ~max_info:maxi minv) else Some {minv;minb;mini;maxv;maxb;maxi} @@ -944,55 +960,55 @@ module ConvexeInfo (Info: sig if the two arguments are equals, return the second *) - let zero ?max_info ~min_info = singleton ~min_info ?max_info Q.zero - let reals = {minb=Strict; minv=Q.minus_inf; mini=Info.nothing; - maxb=Strict; maxv=Q.inf; maxi=Info.nothing} + let zero ?max_info ~min_info = singleton ~min_info ?max_info A.zero + let reals = {minb=Strict; minv=A.minus_inf; mini=Info.nothing; + maxb=Strict; maxv=A.inf; maxi=Info.nothing} let () = assert (invariant reals) let is_reals = function | {minb=Strict; minv; maxb=Strict; maxv} - when Q.equal minv Q.minus_inf && - Q.equal maxv Q.inf -> true + when A.equal minv A.minus_inf && + A.equal maxv A.inf -> true | _ -> false let choose = function | {minb=Large;minv} -> minv | {maxb=Large;maxv} -> maxv - | {minv;maxv} when Q.equal Q.minus_inf minv && Q.equal Q.inf maxv -> - Q.zero - | {minv;maxv} when Q.equal Q.minus_inf minv -> - Q.make (Z.sub (Q.to_bigint maxv) Z.one) Z.one - | {minv;maxv} when Q.equal Q.inf maxv -> - Q.make (Z.add (Q.to_bigint minv) Z.one) Z.one + | {minv;maxv} when A.equal A.minus_inf minv && A.equal A.inf maxv -> + A.zero + | {minv;maxv} when A.equal A.minus_inf minv -> + A.make (Z.sub (A.to_bigint maxv) Z.one) Z.one + | {minv;maxv} when A.equal A.inf maxv -> + A.make (Z.add (A.to_bigint minv) Z.one) Z.one | {minv;maxv} -> - let q = Q.make (Z.add Z.one (Q.to_bigint minv)) Z.one in - if Q.lt q maxv then q - else Q.add maxv (Q.div_2exp (Q.sub minv maxv) 1) + let q = A.make (Z.add Z.one (A.to_bigint minv)) Z.one in + if A.lt q maxv then q + else A.add maxv (A.div_2exp (A.sub minv maxv) 1) let nb_incr = 100 let z_nb_incr = Z.of_int nb_incr let choose_rnd rnd = function - | {minv;maxv} when Q.equal Q.minus_inf minv -> - Q.make (Z.sub (Z.of_int (rnd nb_incr)) (Q.to_bigint maxv)) Z.one - | {minv;maxv} when Q.equal Q.inf maxv -> - Q.make (Z.add (Z.of_int (rnd nb_incr)) (Q.to_bigint minv)) Z.one - | {minv;maxv} when Q.equal minv maxv -> minv + | {minv;maxv} when A.equal A.minus_inf minv -> + A.make (Z.sub (Z.of_int (rnd nb_incr)) (A.to_bigint maxv)) Z.one + | {minv;maxv} when A.equal A.inf maxv -> + A.make (Z.add (Z.of_int (rnd nb_incr)) (A.to_bigint minv)) Z.one + | {minv;maxv} when A.equal minv maxv -> minv | {minv;maxv} -> - let d = Q.sub maxv minv in + let d = A.sub maxv minv in let i = 1 + rnd (nb_incr - 2) in - let x = Q.make (Z.of_int i) (Z.of_int 100) in - let q = Q.add minv (Q.mul x d) in - assert (Q.lt minv q); - assert (Q.lt q maxv); + let x = A.make (Z.of_int i) (Z.of_int 100) in + let q = A.add minv (A.mul x d) in + assert (A.lt minv q); + assert (A.lt q maxv); q let get_convexe_hull e = let mk v b = - match Q.classify v with - | Q.ZERO | Q.NZERO -> Some (v,b) - | Q.INF | Q.MINF -> None - | Q.UNDEF -> assert false in + match A.classify v with + | A.ZERO | A.NZERO -> Some (v,b) + | A.INF | A.MINF -> None + | A.UNDEF -> assert false in mk e.minv e.minb, mk e.maxv e.maxb end diff --git a/src_colibri2/theories/LRA/stages/stage1/interval.mli b/src_colibri2/theories/LRA/stages/stage1/interval.mli index c97378b34cd76fff4bed710ff15d49c822cdcddd..63c6b6683a69adeb38fab1b60440b61a2f19683e 100644 --- a/src_colibri2/theories/LRA/stages/stage1/interval.mli +++ b/src_colibri2/theories/LRA/stages/stage1/interval.mli @@ -18,26 +18,22 @@ (* for more details (enclosed in the file licenses/LGPLv2.1). *) (*************************************************************************) -module Convexe: sig +module Convexe : sig include Interval_sig.S - type split_heuristic = - | Singleton of Q.t - | Splitted of t * t - | NotSplitted - - val split_heuristic: t -> split_heuristic + type split_heuristic = Singleton of A.t | Splitted of t * t | NotSplitted + val split_heuristic : t -> split_heuristic end (* module ConvexeWithDivider: sig include Interval_sig.S val split_heuristic: t -> - [ `Singleton of Q.t + [ `Singleton of A.t | `Splitted of t * t | `NotSplitted ] - val divided_by: Q.t -> t + val divided_by: A.t -> t end module ConvexeWithExceptions: Interval_sig.S @@ -59,23 +55,23 @@ module ConvexeInfo(Info: sig val is_distinct: t -> t -> bool val is_included: t -> t -> bool - val mult_cst: Q.t -> t -> t - val add_cst : Q.t -> t -> t + val mult_cst: A.t -> t -> t + val add_cst : A.t -> t -> t val add: min_info:Info.t -> ?max_info:Info.t -> t -> t -> t val minus: min_info:Info.t -> ?max_info:Info.t -> t -> t -> t - val mem: Q.t -> t -> bool + val mem: A.t -> t -> bool - (** from Q.t *) - val singleton: min_info:Info.t -> ?max_info:Info.t -> Q.t -> t - val is_singleton: t -> Q.t option + (** from A.t *) + val singleton: min_info:Info.t -> ?max_info:Info.t -> A.t -> t + val is_singleton: t -> A.t option - val except: t -> Q.t -> t option + val except: t -> A.t -> t option - val gt: min_info:Info.t -> Q.t -> t - val ge: min_info:Info.t -> Q.t -> t - val lt: max_info:Info.t -> Q.t -> t - val le: max_info:Info.t -> Q.t -> t + val gt: min_info:Info.t -> A.t -> t + val ge: min_info:Info.t -> A.t -> t + val lt: max_info:Info.t -> A.t -> t + val le: max_info:Info.t -> A.t -> t (** > q, >= q, < q, <= q *) val union: t -> t -> t @@ -92,15 +88,15 @@ module ConvexeInfo(Info: sig (** R *) val is_reals: t -> bool - val choose: t -> Q.t + val choose: t -> A.t (** Nothing smart in this choosing *) - val choose_rnd : (int -> int) -> t -> Q.t + val choose_rnd : (int -> int) -> t -> A.t (** choose an element randomly (but non-uniformly), the given function is the random generator *) - val get_convexe_hull: t -> (Q.t * bound) option * (Q.t * bound) option + val get_convexe_hull: t -> (A.t * bound) option * (A.t * bound) option end *)