feat: Rewrite the Intervals module entirely

The existing Intervals module suffers from multiple drawbacks. It is
undocumented, uses a questionable implementation where explanations
associated with each internal bounds have unclear semantics (did I
mention it is not documented), and it has been the source of many
soundness bug in the past due to the way it uses exceptions to indicate
emptiness. This makes it hard to extend the module; for instance OCamlPro#1058
was delayed because we were not quite sure whether the implementation of
functions related to bit-vectors were correct.

This patch is a reimplementation of the Intervals module from scratch.
The new implementation keeps the representation used by the `intersect`
function in the existing implementation (and in a way can be thought of
as resolving the TODO in the existing implementation suggesting to
generalize that type). It is thoroughly documented, and is split between
a "core" module that provides safe functions to deal with explanations,
and specialized implementations for common functions (addition,
multiplication, etc.) using the "core" interface such that reasoning
about the implementation of addition etc. does not require thinking
about explanations at all. This makes it easier to extend the module
with new specialized functions.

The implementation is done through a (small) family of functors,
allowing separate implementations for real and integer intervals that
prevent accidentally mixing them up.

For the time being, the old interface is re-implemented on top of the
new interface (except where implementation details leaked too much) so
as to keep the changes mostly confined to the `Intervals` module.
Migrating to the new interface (and especially abandoning the use of
exceptions) will be done in a second step.

The patch is relatively big, but can't realistically be split into
smaller parts without ending up in intermediate states full of dead
code. I suggest reviewers first take a look at the documentation of the
`OrderedType`, `Interval` and `Union` signatures in the `Intervals_intf`
module (note that this includes some LaTeX, and might be easier to
read using `odoc` -- I tried to make sure the `odoc`-generated
documentation was usable). This should provide a good understanding of
the "core" functionality of the new implementation. The rest of the
review can be split into parts that should be mostly independent:

 - Implementation of the `OrderedType` interface for `Z.t` and `Q.t` in
   `ZEuclideanType` and `QAlgebraicType`;

 - Implementation of the concrete functions for addition,
   multiplication, etc. in the `Intervals` module;

 - Re-implementation of the old API in the `Intervals.Legacy` module
   (and minor related changes in other modules, notably
   `IntervalCalculus`);

 - Implementation of the "core" functionality in the `Intervals_core`
   module.

src/lib/dune

@@ -49,9 +49,10 @@
     ; reasoners
     Ac Arith Arrays_rel Bitv Ccx Shostak Relation Enum Enum_rel
     Fun_sat Fun_sat_frontend Inequalities Bitv_rel Th_util Adt Adt_rel
-    Instances IntervalCalculus Intervals Ite_rel Matching Matching_types
-    Polynome Records Records_rel Satml_frontend_hybrid Satml_frontend Satml
-    Sat_solver Sat_solver_sig Sig Sig_rel Theory Uf Use Rel_utils Bitlist
+    Instances IntervalCalculus Intervals_intf Intervals_core Intervals
+    Ite_rel Matching Matching_types Polynome Records Records_rel
+    Satml_frontend_hybrid Satml_frontend Satml Sat_solver Sat_solver_sig
+    Sig Sig_rel Theory Uf Use Rel_utils Bitlist
     ; structures
     Commands Errors Explanation Fpa_rounding
     Parsed Profiling Satml_types Symbols

src/lib/reasoners/inequalities.ml

@@ -121,10 +121,10 @@ module Container : Container_SIG = struct
       let ple0 = P.sub p1 p2 in
       match P.to_list ple0 with
       | ([], ctt) when is_le && Q.sign ctt > 0->
-        raise (Intervals.NotConsistent expl)
+        raise (Intervals.Legacy.NotConsistent expl)
 
       | ([], ctt) when not is_le  && Q.sign ctt >= 0 ->
-        raise (Intervals.NotConsistent expl)
+        raise (Intervals.Legacy.NotConsistent expl)
 
       | _ ->
         let p,c,d = P.normal_form ple0 in (* ple0 = (p + c) * d, and d > 0 *)

src/lib/reasoners/intervalCalculus.ml

@@ -39,7 +39,7 @@ module Q = Numbers.Q
 module L = Xliteral
 
 module Sy = Symbols
-module I = Intervals
+module I = Intervals.Legacy
 
 open Inequalities
 
@@ -134,7 +134,9 @@ module Sim_Wrap = struct
     | Sim.Core.Unbounded _ -> assert false
     | Sim.Core.Max _ -> assert false
     | Sim.Core.Sat _ -> ()
-    | Sim.Core.Unsat ex ->
+    | Sim.Core.Unsat ex as res ->
+      Format.eprintf "??? %a@."
+        (Sim.Core.print res) simplex;
       let ex = Lazy.force ex in
       if get_debug_fm () then
         Printer.print_dbg
@@ -867,7 +869,7 @@ let rec tighten_ac are_eq x env expl =
       let env =
         if is_alien x then
           (* identity *)
-          let u = I.root n u in
+          let u = I.coerce ty (I.root n (I.coerce Treal u)) in
           let (pu, use_px) =
             try MX.n_find x env.monomes (* we know that x is a monome *)
             with Not_found -> I.undefined ty, SX.empty
@@ -884,7 +886,7 @@ let rec tighten_ac are_eq x env expl =
         Symbols.equal h (Symbols.Op Symbols.Mult) && n > 2 ->
       let env =
         if is_alien x then
-          let u = I.root n u in
+          let u = I.coerce ty (I.root n (I.coerce Treal u)) in
           let pu, use_px =
             try MX.n_find x env.monomes (* we know that x is a monome *)
             with Not_found -> I.undefined ty, SX.empty
@@ -1389,25 +1391,23 @@ let add_disequality are_eq env eqs p expl =
     env, eqs
   | ([a, x], v) ->
     let b = Q.div (Q.minus v) a in
-    let i1 = I.point b ty expl in
     let i2, use2 =
       try
         MX.n_find x env.monomes
       with Not_found -> I.undefined ty, SX.empty
     in
-    let i = I.exclude i1 i2 in
+    let i = I.exclude ~ex:expl b i2 in
     let env = MX.n_add x (i,use2) i2 env in
     let env = tighten_non_lin are_eq x use2 env expl in
     env, find_eq eqs x i env
   | _ ->
     let p, c, _ = P.normal_form_pos p in
-    let i1 = I.point (Q.minus c) ty expl in
     let i2 =
       try
         MP.n_find p env.polynomes
       with Not_found -> I.undefined ty
     in
-    let i = I.exclude i1 i2 in
+    let i = I.exclude ~ex:expl (Q.minus c) i2 in
     let env =
       if I.is_strict_smaller i i2 then
         update_monomes_from_poly p i (MP.n_add p i i2 env)
@@ -1819,11 +1819,16 @@ let new_terms _ = SE.empty
 
 let case_split_union_of_intervals =
   let aux acc uf i z =
-    if Uf.is_normalized uf z then
-      match I.bounds_of i with
-      | [] -> assert false
-      | [_] -> ()
-      | (_,(v, ex))::_ -> acc := Some (z, v, ex); raise Exit
+    if Uf.is_normalized uf z then (
+      ignore @@
+      I.fold (fun prev { lb = _ ; ub } ->
+          match prev, ub with
+          | None, Open ub -> Some (L.LT, ub)
+          | None, Closed ub -> Some (L.LE, ub)
+          | Some bnd, _ -> acc := Some (z, bnd); raise Exit
+          | None, _ -> None
+        ) None i
+    )
   in fun env uf ->
     let cs = ref None in
     try
@@ -1833,38 +1838,20 @@ let case_split_union_of_intervals =
     with Exit ->
     match !cs with
     | None -> assert false
-    | Some(_,None, _) -> assert false
-    | Some(r1,Some (n, eps), _) ->
+    | Some (r1, (pred, n)) ->
       let ty = X.type_info r1 in
       let r2 = alien_of (P.create [] n ty) in
-      let pred =
-        if Q.is_zero eps then L.LE else (assert (Q.is_m_one eps); L.LT)
-      in
       [LR.mkv_builtin true pred [r1; r2], true,
        Th_util.CS (Th_util.Th_arith, Q.one)]
 
 
 (*****)
 
-let bnd_to_simplex_bound ((bnd, explanation) : I.bnd) : Sim.Core.bound option =
-  match bnd with
-  | None -> None
-  | Some (bnd, offset) ->
-    let bvalue =
-      if Q.equal offset Q.one
-      then Sim.Core.R2.lower bnd
-      else if Q.equal offset Q.m_one
-      then Sim.Core.R2.upper bnd
-      else if Q.equal offset Q.zero
-      then Sim.Core.R2.of_r bnd
-      else assert false (* alt-ergo style *)
-    in Some {bvalue; explanation}
-
 let int_constraints_from_map_intervals =
   let aux p xp i uf acc =
     if Uf.is_normalized uf xp && I.is_point i == None
        && P.type_info p == Ty.Tint
-    then (p, I.bounds_of i) :: acc
+    then (p, I.integer_hull i) :: acc
     else acc
   in
   fun env uf ->
@@ -1878,36 +1865,37 @@ let fm_simplex_unbounded_integers_encoding env uf =
   let simplex = Sim.Core.empty ~is_int:true ~check_invs:true in
   let int_ctx = int_constraints_from_map_intervals env uf in
   List.fold_left
-    (fun simplex (p, uints) ->
-       match uints with
-       | [] ->
-         Printer.print_err "Intervals already empty !!!";
-         assert false
-       | _::_::_ ->
-         Printer.print_err
-           "case-split over unions of intervals is needed !!!";
-         assert false
-
-       | [(lb, ex_lb), (ub, ex_ub)] ->
-         let l, c = P.to_list p in
-         assert (Q.is_zero c);
-         assert (lb == None || ub == None);
-         let min = bnd_to_simplex_bound (lb, ex_lb) in
-         let max = bnd_to_simplex_bound (ub, ex_ub) in
-         match l with
-         | [] -> assert false
-         | [c, x] ->
-           assert (Q.is_one c);
-           Sim.Assert.var simplex x ?min ?max |> fst
-         | _ ->
-           let l = List.rev_map (fun (c, x) -> x, c) (List.rev l) in
-           let xp = alien_of p in
-           let sim_p =
-             match Sim.Core.poly_of_slake simplex xp with
-             | Some res -> res
-             | None -> Sim.Core.P.from_list l
-           in
-           Sim.Assert.poly simplex sim_p xp ?min ?max |> fst
+    (fun simplex (p, (lb, ub)) ->
+       let l, c = P.to_list p in
+       assert (Q.is_zero c);
+       let min =
+         match lb with
+         | Some (lb, ex_lb) ->
+           let lb = Q.from_z lb in
+           Some Sim.Core.{ bvalue = R2.of_r lb ; explanation = ex_lb }
+         | None -> None
+       in
+       let max =
+         match ub with
+         | Some (ub, ex_ub) ->
+           let ub = Q.from_z ub in
+           Some Sim.Core.{ bvalue = R2.of_r ub ; explanation = ex_ub }
+         | None -> None
+       in
+       match l with
+       | [] -> assert false
+       | [c, x] ->
+         assert (Q.is_one c);
+         Sim.Assert.var simplex x ?min ?max |> fst
+       | _ ->
+         let l = List.rev_map (fun (c, x) -> x, c) (List.rev l) in
+         let xp = alien_of p in
+         let sim_p =
+           match Sim.Core.poly_of_slake simplex xp with
+           | Some res -> res
+           | None -> Sim.Core.P.from_list l
+         in
+         Sim.Assert.poly simplex sim_p xp ?min ?max |> fst
 
     ) simplex int_ctx
 
@@ -2006,15 +1994,15 @@ let middle_value env ~is_max ty p bound =
       begin
         try
           if is_max then
-            Intervals.new_borne_sup Ex.empty bound ~is_le:false i
+            I.new_borne_sup Ex.empty bound ~is_le:false i
           else
-            Intervals.new_borne_inf Ex.empty bound ~is_le:false i
-        with Intervals.NotConsistent _ -> assert false
+            I.new_borne_inf Ex.empty bound ~is_le:false i
+        with I.NotConsistent _ -> assert false
       end
     | Some i, None -> i
-    | None, _ -> Intervals.point Q.zero ty Ex.empty
+    | None, _ -> I.point Q.zero ty Ex.empty
   in
-  let q = Option.get (Intervals.pick ~is_max interval) in
+  let q = I.pick ~is_max interval in
   alien_of (P.create [] q ty)
 
 let optimizing_objective env uf Objective.Function.{ e; is_max; _ } =