diff --git a/src/bin/common/parse_command.ml b/src/bin/common/parse_command.ml index 5e799b7df0..c4c7e55c46 100644 --- a/src/bin/common/parse_command.ml +++ b/src/bin/common/parse_command.ml @@ -158,6 +158,7 @@ module Debug = struct | Fpa | Gc | Interpretation + | Intervals | Matching | Sat | Split @@ -174,7 +175,7 @@ module Debug = struct let all = [ Debug; Ac; Adt; Arith; Arrays; Bitv; Sum; Ite; Cc; Combine; Constr; Explanation; Fm; Fpa; Gc; - Interpretation; Matching; Sat; Split; Triggers; + Interpretation; Intervals; Matching; Sat; Split; Triggers; Types; Typing; Uf; Unsat_core; Use; Warnings; Commands; Optimize ] @@ -196,6 +197,7 @@ module Debug = struct | Fpa -> "fpa" | Gc -> "gc" | Interpretation -> "interpretation" + | Intervals -> "intervals" | Matching -> "matching" | Sat -> "sat" | Split -> "split" @@ -228,6 +230,7 @@ module Debug = struct | Fpa -> Options.set_debug_fpa verbosity | Gc -> Options.set_debug_gc true | Interpretation -> Options.set_debug_interpretation true + | Intervals -> Options.set_debug_intervals true | Matching -> Options.set_debug_matching verbosity | Sat -> Options.set_debug_sat true | Split -> Options.set_debug_split true diff --git a/src/lib/dune b/src/lib/dune index 3391e92faf..922a6f2c9e 100644 --- a/src/lib/dune +++ b/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 diff --git a/src/lib/reasoners/inequalities.ml b/src/lib/reasoners/inequalities.ml index f4def0a2ed..d51a45e9c0 100644 --- a/src/lib/reasoners/inequalities.ml +++ b/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 *) diff --git a/src/lib/reasoners/intervalCalculus.ml b/src/lib/reasoners/intervalCalculus.ml index ac25097a0f..5636cd82f9 100644 --- a/src/lib/reasoners/intervalCalculus.ml +++ b/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 @@ -492,8 +492,8 @@ let generic_add x j use is_mon env = let p, c = P.separate_constant p0 in let p, c0, d = P.normal_form_pos p in assert (Q.sign d <> 0 && Q.sign c0 = 0); - let j = I.add j (I.point (Q.minus c) ty Explanation.empty) in - let j = I.scale (Q.inv d) j in + (* x = (p + c0) * d <-> p = x * 1/d - c0 / d *) + let j = I.affine_scale ~coef:(Q.inv d) ~const:(Q.div (Q.minus c) d) j in try MP.n_add p j (MP.n_find p env.polynomes) env with Not_found -> MP.n_add p j (I.undefined ty) env @@ -769,14 +769,14 @@ let rec init_monomes_of_poly are_eq env p use_p expl = update_monome are_eq expl use_p env x ) env (fst (P.to_list p)) -and init_alien are_eq expl p (normal_p, c, d) ty use_x env = +and init_alien are_eq expl p (normal_p, c, d) use_x env = let env = init_monomes_of_poly are_eq env p use_x expl in let i = intervals_from_monomes env p in let i = try + (* p = (normal_p + c) * d *) let old_i = MP.n_find normal_p env.polynomes in - let old_i = I.scale d - (I.add old_i (I.point c ty Explanation.empty)) in + let old_i = I.affine_scale ~const:(Q.mult c d) ~coef:d old_i in I.intersect i old_i with Not_found -> i in @@ -822,10 +822,10 @@ and update_monome are_eq expl use_x env x = let (pa', ca, da) as npa = P.normal_form_pos pa in let (pb', cb, db) as npb = P.normal_form_pos pb in let env, ia = - init_alien are_eq expl pa npa ty use_x env in + init_alien are_eq expl pa npa use_x env in let ia = I.add_explanation ia ea in (* take repr into account*) let env, ib = - init_alien are_eq expl pb npb ty use_x env in + init_alien are_eq expl pb npb use_x env in let ib = I.add_explanation ib eb in (* take repr into account*) let ia, ib = match cannot_be_equal_to_zero env pb ib with | Some (ex, _) when Q.equal ca cb @@ -867,7 +867,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 +884,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 @@ -920,19 +920,19 @@ and tighten_non_lin are_eq x use_x env expl = let update_monomes_from_poly p i env = let lp, _ = P.to_list p in - let ty = P.type_info p in List.fold_left (fun env (a,x) -> let np = P.remove x p in let (np,c,d) = P.normal_form_pos np in + (* a * x + (np + c) * d = 0 -> a * x = -d * np - d * c *) try let inp = MP.n_find np env.polynomes in let new_ix = I.scale - (Q.div Q.one a) + (Q.inv a) (I.add i - (I.scale (Q.minus d) - (I.add inp - (I.point c ty Explanation.empty)))) in + (I.affine_scale + ~coef:(Q.minus d) ~const:(Q.mult (Q.minus d) c) inp)) + in let old_ix, ux = MX.n_find x env.monomes in let ix = I.intersect old_ix new_ix in MX.n_add x (ix, ux) old_ix env @@ -1389,25 +1389,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 +1817,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 +1836,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 +1863,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 +1992,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; _ } = @@ -2261,8 +2247,7 @@ let domain_matching _lem_name tr sbt env uf optimized = let p = poly_of rr in let p', c', d = P.normal_form_pos p in let env, i' = best_interval_of optimized env p' in - let ic = I.point c' (P.type_info p') Explanation.empty in - let i = I.scale d (I.add i' ic) in + let i = I.affine_scale ~coef:d ~const:(Q.mult d c') i' in begin match I.match_interval lb ub i idoms with | None -> raise (Sem_match_fails env) | Some idoms -> idoms, maps_to, env, uf diff --git a/src/lib/reasoners/intervals.ml b/src/lib/reasoners/intervals.ml index 92385dc2dd..4418f0e3b2 100644 --- a/src/lib/reasoners/intervals.ml +++ b/src/lib/reasoners/intervals.ml @@ -1,7 +1,7 @@ (**************************************************************************) (* *) (* Alt-Ergo: The SMT Solver For Software Verification *) -(* Copyright (C) 2013-2023 --- OCamlPro SAS *) +(* Copyright (C) 2024 --- OCamlPro SAS *) (* *) (* This file is distributed under the terms of OCamlPro *) (* Non-Commercial Purpose License, version 1. *) @@ -19,1309 +19,1016 @@ (* *) (* CNRS - INRIA - Universite Paris Sud *) (* *) -(* Until 2013, some parts of this code were released under *) -(* the Apache Software License version 2.0. *) -(* *) (* --------------------------------------------------------------- *) (* *) (* More details can be found in the directory licenses/ *) (* *) (**************************************************************************) -module Ex = Explanation - -type borne = - | Strict of (Q.t * Ex.t) - | Large of (Q.t * Ex.t) - | Pinfty | Minfty - -type t = { - ints : (borne * borne) list; - is_int : bool; - expl: Ex.t -} - -exception NotConsistent of Ex.t -exception No_finite_bound - -(*BISECT-IGNORE-BEGIN*) -module Debug = struct - - let print_borne fmt = function - | Minfty -> Format.fprintf fmt "-inf" - | Pinfty -> Format.fprintf fmt "+inf" - | Strict (v, e) | Large (v, e) -> - Format.fprintf fmt "%s" (Q.to_string v); - if Options.(get_verbose () || get_unsat_core ()) then - Format.fprintf fmt " %a" Ex.print e - - let print_interval fmt (b1,b2) = - let c1, c2 = match b1, b2 with - | Large _, Large _ -> '[', ']' - | Large _, _ -> '[', '[' - | _, Large _ -> ']', ']' - | _, _ -> ']', '[' - in - Format.fprintf fmt "%c%a;%a%c" c1 print_borne b1 print_borne b2 c2 +open Intervals_intf - let print_list fmt = function - | [] -> Format.fprintf fmt "[empty]" - | e::l -> - print_interval fmt e; - List.iter (Format.fprintf fmt " U %a" print_interval) l +(* Let us pretend we are using [Logs]. *) +module Log = struct + type ('a, 'b) msgf = + (?header:string -> + ('a, Stdlib.Format.formatter, unit, 'b) Stdlib.format4 -> 'a) + -> 'b - let print fmt { ints; expl = e; _ } = - print_list fmt ints; - if Options.(get_verbose () || get_unsat_core ()) then - Format.fprintf fmt " %a" Ex.print e + type 'a log = ('a, unit) msgf -> unit + let module_name = "AltErgoLib.Intervals" + + let debug : 'a log = + fun k -> + if Options.get_debug_intervals () then + k (fun ?header:function_name fmt -> + let header = Option.is_some function_name in + Printer.print_dbg ~header ~module_name ?function_name fmt + ) end -(*BISECT-IGNORE-END*) - -let print = Debug.print -let pretty_print = Debug.print - - -let large_borne_of n e = Large (n, e) -let strict_borne_of n e = Strict (n, e) - -let undefined_int = {ints = [Minfty, Pinfty]; is_int = true ; expl = Ex.empty} -let undefined_real = {ints = [Minfty, Pinfty]; is_int = false; expl = Ex.empty} - -let undefined ty = match ty with - | Ty.Tint -> undefined_int - | Ty.Treal -> undefined_real - | _ -> assert false - -let is_undefined t = match t.ints with - | [Minfty, Pinfty] -> true - | _ -> false - -let point b ty e = { - ints = [Large (b, e), Large (b, e)]; - is_int = ty == Ty.Tint; - expl = e -} - -let is_point { ints = l; expl = e; _ } = - match l with - | [Large (v1, e1) , Large (v2, e2)] when Q.equal v1 v2 -> - Some (v1, Ex.union e (Ex.union e1 e2)) - | _ -> None - -let finite_size { ints = l; is_int = is_int; _ } = - if not is_int then None - else - try - let acc = ref [] in - List.iter - (fun (b1, b2) -> - match b1, b2 with - | Minfty, _ | _, Pinfty -> raise Exit - | Large (v1, _) , Large (v2, _) -> acc := (v1, v2) :: !acc - | _, _ -> assert false - )l; - let res = - List.fold_left - (fun n (v1, v2) -> Q.add n (Q.add (Q.sub v2 v1) Q.one)) Q.zero !acc - in - Some res - with Exit -> None -let borne_inf = function - | { ints = (Large (v, ex), _) :: _; _ } -> v, ex, true - | { ints = (Strict (v, ex), _) :: _; _ } -> v, ex, false - | _ -> raise No_finite_bound - -let only_borne_inf = function - | { ints = (inf, _) :: _ ; _ } as t -> { t with ints = [(inf, Pinfty)] } - | _ -> assert false - -let borne_sup { ints; _ } = - match List.rev ints with - | (_, Large (v, ex))::_ -> v, ex, true - | (_, Strict (v, ex))::_ -> v, ex, false - | _ -> raise No_finite_bound - -let only_borne_sup t = - let rec aux = function - | [] -> assert false - | [ (_, sup) ] -> { t with ints = [(Minfty, sup)] } - | _ :: tl -> aux tl - in aux t.ints - -let explain_borne = function - | Large (_, e) | Strict (_, e) -> e - | _ -> Ex.empty - -let add_expl_to_borne b e = - if Ex.is_empty e then b - else - match b with - | Large (n, e') -> Large (n, Ex.union e e') - | Strict (n, e') -> Strict (n, Ex.union e e') - | Pinfty | Minfty -> b - -let add_expl_zero i expl = - if Ex.is_empty expl then i - else - let res = - List.rev_map (fun x -> - match x with - | Large (c1,e1), Large (c2,e2) when Q.sign c1 = 0 && Q.sign c2 = 0 -> - Large (Q.zero, Ex.union e1 expl), Large (Q.zero, Ex.union e2 expl) - | _ -> x - ) i.ints +module Ring(C : Core)(RT : RingType) = struct + include C.Union(RT) + + type monotony = Inc | Dec + + (* Helper to keep the code of [mul] and [ediv] readable. *) + let map2_mon_to_set f p1 u1 p2 u2 = + let f_lb, f_ub = + match p1 with + | Inc -> (fun i1 -> f i1.lb), (fun i1 -> f i1.ub) + | Dec -> (fun i1 -> f i1.ub), (fun i1 -> f i1.lb) in - { i with ints = List.rev res } - -let int_of_borne_inf b = - match b with - | Minfty | Pinfty -> b - | Large (v, e) -> - Large ((if Numbers.Q.is_int v then v else Numbers.Q.ceiling v), e) - | Strict (v, e) -> - if Numbers.Q.is_int v then Large (Q.add v Q.one, e) + map_mon_to_set (match p2 with + | Inc -> fun i1 -> approx_map_inc_to_set (f_lb i1) (f_ub i1) u2 + | Dec -> fun i1 -> approx_map_dec_to_set (f_lb i1) (f_ub i1) u2 + ) u1 + + let trace1 fname u u' = + Log.debug (fun m -> + m "@[%s(%a) = %a@]" + fname (Fmt.box pp) u (Fmt.box pp) u' + ); + u' + + let trace2 fname u1 u2 u = + Log.debug (fun m -> + m "@[%s(%a, %a) = %a@]" + fname (Fmt.box pp) u1 (Fmt.box pp) u2 (Fmt.box pp) u + ); + u + + let add u1 u2 = + trace2 "add" u1 u2 @@ of_set_nonempty @@ + map2_mon_to_set RT.add Inc u1 Inc u2 + + let mul u1 u2 = + trace2 "mul" u1 u2 @@ of_set_nonempty @@ + trisection_map_to_set RT.zero u1 + (fun lt1 -> + trisection_map_to_set RT.zero u2 + (map2_mon_to_set RT.mul Dec lt1 Dec) + (fun () -> interval_set { lb = RT.zero ; ub = RT.zero }) + (map2_mon_to_set RT.mul Inc lt1 Dec)) + (fun _ -> interval_set { lb = RT.zero ; ub = RT.zero }) + (fun gt1 -> + trisection_map_to_set RT.zero u2 + (map2_mon_to_set RT.mul Dec gt1 Inc) + (fun () -> interval_set { lb = RT.zero ; ub = RT.zero }) + (map2_mon_to_set RT.mul Inc gt1 Inc)) + + let pow n u = + if n <= 0 then Fmt.invalid_arg "pow: nonpositive exponent %d" n; + if n mod 2 = 0 then + of_set_nonempty @@ + trisection_map_to_set RT.zero u + (map_dec_to_set (RT.pow n)) + (fun () -> interval_set { lb = RT.zero ; ub = RT.zero }) + (map_inc_to_set (RT.pow n)) else - let v' = Numbers.Q.ceiling v in - assert (Q.compare v' v > 0); - Large (v', e) - -let int_of_borne_sup b = - match b with - | Minfty | Pinfty -> b - | Large (v, e) -> - Large ((if Numbers.Q.is_int v then v else Numbers.Q.floor v), e) - | Strict (v, e) -> - if Numbers.Q.is_int v then Large (Q.sub v Q.one, e) + map_strict_inc (RT.pow n) u + + let neg u = trace1 "neg" u @@ map_strict_dec RT.neg u +end + +module Field(C : Core)(FT : FieldType) = struct + include Ring(C)(FT) + + let inv ?(inv0 = Interval.full) u = + (* NB: trisection is required even though we are decreasing on both side, + because we are not decreasing on ]-oo, +oo[. *) + trace1 "inv" u @@ of_set_nonempty @@ + trisection_map_to_set FT.zero u + (map_dec_to_set FT.inv) + (fun () -> interval_set inv0) + (map_dec_to_set FT.inv) + + let div ?inv0 x y = + mul x (inv ?inv0 y) +end + +module AlgebraicField(C : Core)(AT : AlgebraicType) = struct + include Field(C)(AT) + + let root n u = + if n <= 0 then Fmt.invalid_arg "root: nonpositive radical %d" n; + of_set_checked @@ + if n mod 2 = 0 then + trisection_map_to_set AT.zero u + (fun _ -> empty_set) + (fun () -> interval_set { lb = AT.zero ; ub = AT.zero }) + (fun gt -> + let pos = + approx_map_inc_to_set (AT.root_default n) (AT.root_excess n) gt + in + let neg_int i = { lb = AT.neg i.ub ; ub = AT.neg i.lb } in + union_set (map_set neg_int pos) pos) else - let v' = Numbers.Q.floor v in - assert (Q.compare v' v < 0); - Large (v', e) - -let int_bornes (l, u) = int_of_borne_inf l, int_of_borne_sup u - -let compare_bounds b1 ~is_low1 b2 ~is_low2 = - match b1, b2 with - | Minfty, Minfty | Pinfty, Pinfty -> 0 - | Minfty, _ | _, Pinfty -> -1 - | _, Minfty | Pinfty, _ -> 1 - - | Large (v1, _), Large (v2, _) -> Q.compare v1 v2 - - | Strict (v1, _), Strict (v2, _) -> - let c = Q.compare v1 v2 in - if c <> 0 then c - else if is_low1 == is_low2 then 0 (* bl_bl or bu_bu *) - else if is_low1 then 1 (* implies not is_low2 *) - else -1 (* implies not is_low1 and is_low2 *) - - | Strict (v1, _), Large (v2, _) -> - let c = Q.compare v1 v2 in - if c <> 0 then c else if is_low1 then 1 else -1 - - | Large (v1, _), Strict (v2, _) -> - let c = Q.compare v1 v2 in - if c <> 0 then c else if is_low2 then -1 else 1 - -let zero_endpoint b = - match b with - | Large (v, _) -> Numbers.Q.is_zero v - | _ -> false - -let min_of_lower_bounds b1 b2 = - if compare_bounds b1 ~is_low1:true b2 ~is_low2:true <= 0 then b1 else b2 - -let max_of_upper_bounds b1 b2 = - if compare_bounds b1 ~is_low1:false b2 ~is_low2:false >= 0 then b1 else b2 - -let zero_large = Large (Q.zero, Ex.empty) - -let low_borne_pos_strict b = - compare_bounds b ~is_low1:true zero_large ~is_low2:true > 0 - -let up_borne_pos_strict b = - compare_bounds b ~is_low1:false zero_large ~is_low2:false > 0 - -let low_borne_neg_strict b = - compare_bounds b ~is_low1:true zero_large ~is_low2:true < 0 - -let up_borne_neg_strict b = - compare_bounds b ~is_low1:false zero_large ~is_low2:false < 0 - - -(* should be removed: probably buggy when mixing lower and upper bounds *) -let pos_borne b = match b with - | Pinfty -> true - | Minfty -> false - | Strict (v, _) | Large (v, _) -> Q.sign v >= 0 - -(* should be removed: probably buggy when mixing lower and upper bounds *) -let neg_borne b = match b with - | Pinfty -> false - | Minfty -> true - | Strict (v, _) | Large (v, _) -> Q.sign v <= 0 - -(* TODO: generalize the use of this type and the function joint below - to other operations on intervals *) -type kind = - | Empty of Explanation.t - | Int of (borne * borne) - -let join l glob_ex = (* l should not be empty *) - let rec j_aux _todo _done = - match _todo, _done with - | [], [] -> assert false - | [], _ -> List.rev _done, None - | [Empty ex], [] -> [], Some ex - | (Int b) :: l, _ -> j_aux l (b :: _done) - | (Empty ex) :: l, _ -> - let _done = match _done with - | [] -> [] - | (low, up) :: r -> (low, add_expl_to_borne up ex) :: r - in - let _todo = match l with - | [] -> [] - | (Empty ex') :: r -> (Empty (Ex.union ex ex')) :: r - | (Int (low, up)) :: r -> (Int (add_expl_to_borne low ex, up)) :: r - in - j_aux _todo _done - in - match j_aux l [] with - | [], None -> assert false - | l , None -> l - | [], Some ex -> raise (NotConsistent (Ex.union ex glob_ex)); - | _ , Some _ -> assert false - - -let intersect = - let rec step is_int l1 l2 acc = - match l1, l2 with - | [], _ | _, [] -> List.rev acc - - | (_, up1) :: r1, (lo2, _) :: _ when - compare_bounds up1 ~is_low1:false lo2 ~is_low2:true < 0 -> - (* No overlap. (lo1, up1) is smaller *) - let nexpl = Ex.union (explain_borne up1) (explain_borne lo2) in - step is_int r1 l2 ((Empty nexpl) :: acc) - - | (lo1, _) :: _, (_, up2) :: r2 when - compare_bounds lo1 ~is_low1:true up2 ~is_low2:false > 0 -> - (* No overlap. (lo2, up2) is smaller *) - let nexpl = Ex.union (explain_borne up2) (explain_borne lo1) in - step is_int l1 r2 ((Empty nexpl) :: acc) - - | (lo1, up1)::r1, (lo2, up2)::r2 -> - let cll = compare_bounds lo1 ~is_low1:true lo2 ~is_low2:true in - let cuu = compare_bounds up1 ~is_low1:false up2 ~is_low2:false in - if cll <= 0 && cuu >= 0 then (* (lo1, up1) subsumes (lo2, up2) *) - step is_int l1 r2 ((Int (lo2,up2))::acc) (* ex of lo1 and up1 ? *) + approx_map_inc_to_set (AT.root_default n) (AT.root_excess n) u +end + +module EuclideanRing(C : Core)(ET : EuclideanType) = struct + include Ring(C)(ET) + + let ediv ?(div0 = Interval.full) u1 u2 = + trace2 "ediv" u1 u2 @@ of_set_nonempty @@ + trisection_map_to_set ET.zero u2 + (fun lt2 -> + trisection_map_to_set ET.zero u1 + (fun lt1 -> map2_mon_to_set ET.ediv Dec lt1 Inc lt2) + (fun () -> interval_set { lb = ET.zero ; ub = ET.zero }) + (fun gt1 -> map2_mon_to_set ET.ediv Dec gt1 Dec lt2)) + (fun () -> interval_set div0) + (fun gt2 -> + trisection_map_to_set ET.zero u1 + (fun lt1 -> map2_mon_to_set ET.ediv Inc lt1 Inc gt2) + (fun () -> interval_set { lb = ET.zero ; ub = ET.zero }) + (fun gt1 -> map2_mon_to_set ET.ediv Inc gt1 Dec gt2)) +end + +(* EuclideanType interface for integers: all intervals are closed, but we still + need to handle infinities. *) +module ZEuclideanType = struct + type finite = Z.t + + let pp_finite = Z.pp_print + + type t = Neg_infinite | Finite of finite | Pos_infinite + + let pp ppf = function + | Neg_infinite -> Format.fprintf ppf "-oo" + | Finite n -> pp_finite ppf n + | Pos_infinite -> Format.fprintf ppf "+oo" + + let equal x y = + match x, y with + | Neg_infinite, Neg_infinite -> true + | Neg_infinite, _ | _, Neg_infinite -> false + + | Pos_infinite, Pos_infinite -> true + | Pos_infinite, _ | _, Pos_infinite -> false + + | Finite x, Finite y -> Z.equal x y + + let compare x y = + match x, y with + | Neg_infinite, Neg_infinite -> 0 + | Neg_infinite, _ -> -1 + | _, Neg_infinite -> 1 + + | Finite x, Finite y -> Z.compare x y + + | Pos_infinite, Pos_infinite -> 0 + | Pos_infinite, _ -> 1 + | _, Pos_infinite -> -1 + + let minfty = Neg_infinite + + let finite n = Finite n + + let pinfty = Pos_infinite + + let value_opt = function Finite n -> Some n | _ -> None + + let succ = function + | Finite x -> Finite (Z.succ x) + | inf -> inf + + let pred = function + | Finite x -> Finite (Z.pred x) + | inf -> inf + + let view = function + | Neg_infinite | Pos_infinite -> Unbounded + | Finite x -> Closed x + + let zero = Finite Z.zero + + let sign = function + | Neg_infinite -> -1 + | Finite n -> Z.sign n + | Pos_infinite -> 1 + + let neg = function + | Neg_infinite -> Pos_infinite + | Finite n -> Finite (Z.neg n) + | Pos_infinite -> Neg_infinite + + let add x y = + match x, y with + | Pos_infinite, Neg_infinite | Neg_infinite, Pos_infinite -> + Fmt.invalid_arg "add: %a + %a is undefined" pp x pp y + + | Pos_infinite, _ | _, Pos_infinite -> Pos_infinite + | Neg_infinite, _ | _, Neg_infinite -> Neg_infinite + + | Finite x, Finite y -> Finite (Z.add x y) + + let mul x y = + match x, y with + | (Pos_infinite | Neg_infinite), (Pos_infinite | Neg_infinite) -> + if sign x = sign y then Pos_infinite else Neg_infinite + + | Finite x, Finite y -> Finite (Z.mul x y) + + | Finite fin, (Pos_infinite | Neg_infinite as inf) + | (Pos_infinite | Neg_infinite as inf), Finite fin -> + let s = Z.sign fin in + if s > 0 then inf + else if s < 0 then neg inf + else Finite Z.zero + + let pow n x = + if n <= 0 then Fmt.invalid_arg "pow: nonpositive exponent %d" n; + match x with + | Pos_infinite -> Pos_infinite + | Neg_infinite -> + if n mod 2 = 0 then Pos_infinite else Neg_infinite + | Finite x -> Finite (Z.pow x n) + + let ediv x y = + match x, y with + | Pos_infinite, Pos_infinite + | Neg_infinite, Neg_infinite -> + Pos_infinite + + | Neg_infinite, Pos_infinite + | Pos_infinite, Neg_infinite -> + Neg_infinite + + | Finite x, Finite y -> Finite (Z.ediv x y) + + | Pos_infinite, Finite x -> + if Z.sign x > 0 then Pos_infinite else Neg_infinite + + | Neg_infinite, Finite x -> + if Z.sign x < 0 then Pos_infinite else Neg_infinite + + | Finite x, Pos_infinite -> + if Z.sign x < 0 then Finite Z.minus_one else Finite Z.zero + + | Finite x, Neg_infinite -> + if Z.sign x < 0 then Finite Z.one else Finite Z.zero +end + +(* AlgebraicType interface for reals + + In order to deal with open intervals, we allow bounds of the form [x + eps] + and [x - eps] where [eps] is an infinitesimal. Effectively, [x + eps] is a + strict lower bound, and [x - eps] is a strict upper bound. *) +module QAlgebraicType = struct + type finite = Q.t + + let pp_finite = Q.pp_print + + type kind = Lt | Eq | Gt + (* We store bounds as pairs [(kind, v)] where [(Lt, v)] denotes [v - eps] or + the upper bound [x < v]] and [(Gt, v)] denotes [v + eps] or the lower bound + [x > v]. *) + + let pp_kind pp ppf (k, x) = + match k with + | Lt -> Format.fprintf ppf "%a-eps" pp x + | Eq -> pp ppf x + | Gt -> Format.fprintf ppf "%a+eps" pp x + + let equal_kind x y = + match x, y with + | Lt, Lt -> true + | Lt, _ | _, Lt -> false + + | Gt, Gt -> true + | Gt, _ | _, Gt -> false + + | Eq, Eq -> true + + let compare_kind x y = + match x, y with + | Lt, Lt -> 0 + | Lt, _ -> -1 + | _, Lt -> 1 + + | Eq, Eq -> 0 + + | Gt, Gt -> 0 + | Gt, _ -> 1 + | _, Gt -> -1 + + type t = Neg_infinite | Finite of kind * finite | Pos_infinite + + let of_bigint = function + | ZEuclideanType.Neg_infinite -> Neg_infinite + | Finite n -> Finite (Eq, Q.of_bigint n) + | Pos_infinite -> Pos_infinite + + let floor = function + | Neg_infinite -> ZEuclideanType.Neg_infinite + | Finite (k, x) -> + (* If [y < x] with [x] an integer, then [floor(y) <= x - 1] *) + if Z.equal (Q.den x) Z.one then + match k with + | Lt -> Finite (Z.pred @@ Q.num x) + | Eq | Gt -> Finite (Q.num x) else - if cll >= 0 && cuu <= 0 then (* (lo2, up2) subsumes (lo1, up1) *) - step is_int r1 l2 ((Int(lo1,up1))::acc) (* ex of lo2 and up2 ? *) + Finite (Z.fdiv (Q.num x) (Q.den x)) + | Pos_infinite -> Pos_infinite + + let ceiling = function + | Neg_infinite -> ZEuclideanType.Neg_infinite + | Finite (k, x) -> + (* If [x < y] with [x] an integer, then [x + 1 <= ceiling(y)] *) + if Z.equal (Q.den x) Z.one then + match k with + | Lt | Eq -> Finite (Q.num x) + | Gt -> Finite (Z.succ @@ Q.num x) else - if cll <= 0 && cuu <= 0 then (* lo1 <= lo2 <= up1 <= up2 *) - step is_int r1 l2 ((Int(lo2,up1))::acc) (* ex of lo1 and up2 ? *) + Finite (Z.cdiv (Q.num x) (Q.den x)) + | Pos_infinite -> Pos_infinite + + let pp ppf = function + | Neg_infinite -> Format.fprintf ppf "-oo" + | Finite (k, x) -> pp_kind pp_finite ppf (k, x) + | Pos_infinite -> Format.fprintf ppf "+oo" + + let equal x y = + match x, y with + | Neg_infinite, Neg_infinite -> true + | Neg_infinite, _ | _, Neg_infinite -> false + + | Pos_infinite, Pos_infinite -> true + | Pos_infinite, _ | _, Pos_infinite -> false + + | Finite (kx, x), Finite (ky, y) -> + Q.equal x y && equal_kind kx ky + + let compare x y = + match x, y with + | Neg_infinite, Neg_infinite -> 0 + | Neg_infinite, _ -> -1 + | _, Neg_infinite -> 1 + + | Finite (kx, x), Finite (ky, y) -> + let c = Q.compare x y in + if c <> 0 then c else + compare_kind kx ky + + | Pos_infinite, Pos_infinite -> 0 + | Pos_infinite, _ -> 1 + | _, Pos_infinite -> -1 + + let minfty = Neg_infinite + + let finite x = Finite (Eq, x) + + let pinfty = Pos_infinite + + let value_opt = function Finite (Eq, x) -> Some x | _ -> None + + let succ = function + | Finite (Lt, x) -> Finite (Eq, x) + | Finite (Eq, x) -> Finite (Gt, x) + | Finite (Gt, _) as x -> + Fmt.invalid_arg "succ: %a" pp x + | inf -> inf + + let pred = function + | Finite (Lt, _) as x -> + Fmt.invalid_arg "pred: %a" pp x + | Finite (Eq, x) -> Finite (Lt, x) + | Finite (Gt, x) -> Finite (Eq, x) + | inf -> inf + + let view = function + | Neg_infinite | Pos_infinite -> Unbounded + | Finite (Eq, x) -> Closed x + | Finite (_, x) -> Open x + + let zero = Finite (Eq, Q.zero) + + let sign_kind = function + | Lt -> -1 + | Eq -> 0 + | Gt -> 1 + + let kind_of_sign s = + if s < 0 then Lt else if s > 0 then Gt else Eq + + let neg_kind k = kind_of_sign (-sign_kind k) + + let neg = function + | Pos_infinite -> Neg_infinite + | Neg_infinite -> Pos_infinite + | Finite (k, x) -> Finite (neg_kind k, Q.neg x) + + let add x y = + match x, y with + | Pos_infinite, Pos_infinite + | Pos_infinite, Finite _ + | Finite _, Pos_infinite -> + Pos_infinite + + | Neg_infinite, Neg_infinite + | Neg_infinite, Finite _ + | Finite _, Neg_infinite -> + Neg_infinite + + | Finite (kx, qx), Finite (ky, qy) -> + (* Can never add lower bound to upper bound *) + let skx = sign_kind kx and sky = sign_kind ky in + if (skx * sky < 0) then + Fmt.invalid_arg "add: %a + %a" pp x pp y; + Finite (kind_of_sign (sign_kind kx + sign_kind ky), Q.add qx qy) + + | Pos_infinite, Neg_infinite + | Neg_infinite, Pos_infinite -> + Fmt.invalid_arg "add: %a + %a is undefined" pp x pp y + + let sign = function + | Neg_infinite -> -1 + | Pos_infinite -> 1 + | Finite (k, x) -> + let s = Q.sign x in + if s <> 0 then s else + sign_kind k + + let mul x y = + match x, y with + | Pos_infinite, Pos_infinite + | Neg_infinite, Neg_infinite -> + Pos_infinite + + | Pos_infinite, Neg_infinite + | Neg_infinite, Pos_infinite -> + Neg_infinite + + | (Pos_infinite | Neg_infinite as i), (Finite _ as f) + | (Finite _ as f), (Pos_infinite | Neg_infinite as i) -> + let s = sign f in + if s > 0 then i + else if s < 0 then neg i + else Finite (Eq, Q.zero) + + | Finite (kx, qx), Finite (ky, qy) -> + let sx = sign x and sy = sign y in + if sx = 0 || sy = 0 then + Finite (Eq, Q.zero) else - if cll >= 0 && cuu >= 0 then (* lo2 <= lo1 <= up2 <= up1 *) - step is_int l1 r2 (Int((lo1,up2))::acc) (* ex of lo2 and up1 ? *) - else assert false - in - fun - ({ints=l1; expl=e1; is_int=is_int} as uints1) - { ints = l2; expl = e2; _ } -> - (*l1 and l2 are supposed to be normalized *) - let expl = Ex.union e1 e2 in - let l = step is_int l1 l2 [] in - let res = { uints1 with ints = join l expl; expl } in - assert (res.ints != []); - res - -let new_borne_sup expl b ~is_le uints = - let aux b expl = - let b = (if is_le then large_borne_of else strict_borne_of) b expl in - if uints.is_int then int_of_borne_sup b else b - in - intersect - { ints = [Minfty, aux b expl]; - is_int = uints.is_int; - expl = Ex.empty } uints - -let new_borne_inf expl b ~is_le uints = - let aux b expl = - let b = (if is_le then large_borne_of else strict_borne_of) b expl in - if uints.is_int then int_of_borne_inf b else b - in - intersect - { ints = [aux b expl, Pinfty]; - is_int = uints.is_int; - expl = Ex.empty } uints - -type interval_class = - | P - | M - | N - | Z - -let class_of l u = - if zero_endpoint l && zero_endpoint u then Z - else if pos_borne l && pos_borne u then begin - assert (up_borne_pos_strict u); - P - end - else if neg_borne l && neg_borne u then begin - assert (low_borne_neg_strict l); - N - end - else begin - assert (low_borne_neg_strict l); - assert (up_borne_pos_strict u); - M - end - -let union_bornes is_int l = - let rec aux is_int l acc = - match l with - | [] -> acc - | [e] -> e::acc - | (l1, u1)::((l2, u2)::r as r2) -> - if compare_bounds u1 ~is_low1:false l2 ~is_low2:true < 0 then - match is_int, u1, l2 with - | true, Large(x,_), Large (y, _) when Q.equal (Q.sub y x) Q.one -> - aux is_int ((l1, u2)::r) acc - | _ -> - (* the only case where we put things in acc *) - aux is_int r2 ((l1, u1)::acc) + let skx = sign_kind kx and sky = sign_kind ky in + (* Check that [x] and [y] do not have opposite epsilon if they + have the same sign, or identical epsilon if they have opposite + signs: this matches the monotonicity of multiplication. *) + assert (skx = 0 || sky = 0 || sx * sy = skx * sky); + let k = kind_of_sign (skx * sy + sky * sx) in + Finite (k, Q.mul qx qy) + + let pow n x = + if n <= 0 then Fmt.invalid_arg "pow: nonpositive exponent %d" n; + match x with + | Pos_infinite -> Pos_infinite + | Neg_infinite -> + if n mod 2 = 0 then Pos_infinite else Neg_infinite + | Finite (k, qx) -> + let qpow = Q.make (Z.pow (Q.num qx) n) (Z.pow (Q.den qx) n) in + if n mod 2 = 0 && sign x < 0 then + Finite (neg_kind k, qpow) else - if compare_bounds u1 ~is_low1:false u2 ~is_low2:false > 0 then - aux is_int ((l1, u1)::r) acc + Finite (k, qpow) + + let inv = function + | Neg_infinite -> Finite (Lt, Q.zero) + | Pos_infinite -> Finite (Gt, Q.zero) + | Finite (k, x) -> + if Q.(x = zero) then + match k with + | Lt -> Neg_infinite + | Eq -> Fmt.invalid_arg "inv: 1/0 is undefined" + | Gt -> Pos_infinite else - aux is_int ((l1, u2)::r) acc - in - List.rev (aux is_int l []) - -let union_intervals uints = - let l = List.fast_sort (fun (l1, _) (l2, _) -> - compare_bounds l1 ~is_low1:true l2 ~is_low2:true) uints.ints in - {uints with ints = union_bornes uints.is_int l} - -let minus_borne = function - | Minfty -> Pinfty - | Pinfty -> Minfty - | Large (v, e) -> Large (Numbers.Q.minus v, e) - | Strict (v, e) -> Strict (Numbers.Q.minus v, e) - -let rev_normalize_int_bounds rl ex n = - let l = - List.rev_map - (fun b -> - let b = int_bornes b in - match b with - | Large (v, ex1), Large (w, ex2) when Q.compare w v < 0 -> - Empty (Ex.union ex1 ex2) - | Strict (v, ex1), Large (w, ex2) - | Large (v, ex1) , Strict (w, ex2) - | Strict (v, ex1), Strict (w, ex2) when Q.compare w v <= 0 -> - Empty (Ex.union ex1 ex2) - | _ -> Int b - ) rl - in - if Q.compare n Q.zero > 0 (* !!! this test should be checked *) then join l ex - else List.rev (join (List.rev l) ex) - - -let exclude = - let rec complement l prev acc = - match l with - | (b1,b2)::r -> - let bu = match b1 with - | Strict v -> Large v - | Large v -> Strict v - | _ -> b1 + Finite (neg_kind k, Q.inv x) + + let root_default_num n v = + if n = 2 then Numbers.Q.sqrt_default v else Numbers.Q.root_default v n + + let root_excess_num n v = + if n = 2 then Numbers.Q.sqrt_excess v else Numbers.Q.root_excess v n + + let root_default n = function + | Neg_infinite -> Neg_infinite + | Pos_infinite -> Pos_infinite + | Finite (_, x) -> + match + if Q.sign x >= 0 then root_default_num n x + else Option.map Q.neg (root_excess_num n (Q.neg x)) + with + | None -> Neg_infinite + | Some x -> Finite (Eq, x) + + let root_excess n = function + | Neg_infinite -> Neg_infinite + | Pos_infinite -> Pos_infinite + | Finite (_, x) -> + match + if Q.sign x >= 0 then root_excess_num n x + else Option.map Q.neg (root_default_num n (Q.neg x)) + with + | None -> Pos_infinite + | Some x -> Finite (Eq, x) + + let ( + ) = add + + let ( * ) = mul + + let (~$$) = of_bigint +end + +module Core = Intervals_core.Make(struct + include Explanation + + let pp = print + end) + +type 'a union = 'a Core.union + +module Real = AlgebraicField(Core)(QAlgebraicType) + +module Int = EuclideanRing(Core)(ZEuclideanType) + +module Legacy = struct + type t = Real of Real.t | Int of Int.t + + let pretty_print ppf = function + | Real u -> Real.pp ppf u + | Int u -> Int.pp ppf u + + let print = pretty_print + + let of_real u = Real u + + let of_int u = Int u + + exception NotConsistent of Explanation.t + + let check = function + | NonEmpty u -> u + | Empty ex -> raise @@ NotConsistent ex + + exception No_finite_bound + + let undefined ty = + match ty with + | Ty.Tint -> of_int @@ Int.of_interval @@ Int.Interval.full + | Treal -> of_real @@ Real.of_interval @@ Real.Interval.full + | _ -> invalid_arg "undefined" + + let point v ty ex = + match ty with + | Ty.Treal -> + of_real @@ Real.of_interval ~ex @@ Real.Interval.singleton v + | Tint when Z.equal (Q.den v) Z.one -> + of_int @@ Int.of_interval ~ex @@ Int.Interval.singleton @@ Q.num v + | _ -> Fmt.invalid_arg "point: %a" Q.pp_print v + + let doesnt_contain_0 u = + Option.map (fun ex -> (ex, [])) @@ + match u with + | Real u -> Real.( + match intersect u (of_interval @@ Interval.singleton Q.zero) with + | NonEmpty _ -> None + | Empty ex -> Some ex + ) + | Int u -> Int.( + match intersect u (of_interval @@ Interval.singleton Z.zero) with + | NonEmpty _ -> None + | Empty ex -> Some ex + ) + + let is_strict_smaller u1 u2 = + match u1, u2 with + | Real u1, Real u2 -> Real.subset ~strict:true u1 u2 + | Int u1, Int u2 -> Int.subset ~strict:true u1 u2 + | Int _, Real _ | Real _, Int _ -> invalid_arg "is_strict_smaller" + + let new_borne_sup ex b ~is_le u = + match u with + | Real u -> + let ub = if is_le then Closed b else Open b in + of_real @@ check @@ Real.intersect u @@ Real.of_bounds ~ex Unbounded ub + | Int u -> + let b = Z.fdiv (Q.num b) (Q.den b) in + let ub = if is_le then Closed b else Open b in + of_int @@ check @@ Int.intersect u @@ Int.of_bounds ~ex Unbounded ub + + let new_borne_inf ex b ~is_le u = + match u with + | Real u -> + let lb = if is_le then Closed b else Open b in + of_real @@ check @@ Real.intersect u @@ Real.of_bounds ~ex lb Unbounded + | Int u -> + let b = Z.cdiv (Q.num b) (Q.den b) in + let lb = if is_le then Closed b else Open b in + of_int @@ check @@ Int.intersect u @@ Int.of_bounds ~ex lb Unbounded + + let only_borne_sup= function + | Real u -> + let ub, ex = Real.upper_bound u in + of_real @@ Real.of_bounds ~ex Unbounded ub + | Int u -> + let ub, ex = Int.upper_bound u in + of_int @@ Int.of_bounds ~ex Unbounded ub + + let only_borne_inf = function + | Real u -> + let lb, ex = Real.lower_bound u in + of_real @@ Real.of_bounds ~ex lb Unbounded + | Int u -> + let lb, ex = Int.lower_bound u in + of_int @@ Int.of_bounds ~ex lb Unbounded + + let is_point = function + | Real u -> Real.value_opt u + | Int u -> ( + match Int.value_opt u with + | Some (v, ex) -> Some (Q.of_bigint v, ex) + | _ -> None + ) + + let intersect u1 u2 = + match u1, u2 with + | Real u1, Real u2 -> of_real @@ check @@ Real.intersect u1 u2 + | Int u1, Int u2 -> of_int @@ check @@ Int.intersect u1 u2 + | Int _, Real _ | Real _, Int _ -> invalid_arg "intersect" + + let exclude ?ex v u = + match u with + | Real u -> ( + let v = Real.Interval.singleton v in + match Real.intersect u @@ check @@ Real.of_complement ?ex v with + | NonEmpty u -> of_real u + | Empty ex -> raise @@ NotConsistent ex + ) + | Int u -> ( + if Z.equal (Q.den v) Z.one then + let v = Int.Interval.singleton @@ Q.num v in + match Int.intersect u @@ check @@ Int.of_complement ?ex v with + | NonEmpty u -> of_int u + | Empty ex -> raise @@ NotConsistent ex + else + Int u + ) + + let map2 fname freal fint u1 u2 = + match u1, u2 with + | Real u1, Real u2 -> of_real (freal u1 u2) + | Int u1, Int u2 -> of_int (fint u1 u2) + | _ -> + Fmt.invalid_arg "%s: cannot %s int and real" + fname fname + + let mult = map2 "mult" Real.mul Int.mul + + let power n = function + | Real u -> of_real @@ Real.pow n u + | Int u -> of_int @@ Int.pow n u + + let root n = function + | Real u -> of_real @@ check @@ Real.root n u + | Int _ -> invalid_arg "root: cannot take integral root" + + let add = map2 "add" Real.add Int.add + + let neg = function + | Real u -> of_real @@ Real.neg u + | Int u -> of_int @@ Int.neg u + + let affine_scale ~const ~coef u = + let coef, u = + let c = Q.sign coef in + if c < 0 then Q.neg coef, neg u + else if c > 0 then coef, u + else Fmt.invalid_arg "affine: scale cannot be zero" + in + let open QAlgebraicType in + let scale = finite coef in + let offset = finite const in + match u with + | Real u -> + of_real @@ Real.map_strict_inc (fun b -> b * scale + offset) u + | Int u -> + of_int @@ check @@ + Int.partial_map_inc + (fun lb -> ceiling @@ ~$$lb * scale + offset) + (fun ub -> floor @@ ~$$ub * scale + offset) + u + + let scale scale u = affine_scale ~const:Q.zero ~coef:scale u + + let finite_size u = + match u with + | Int u -> ( + try + Some ( + Q.of_bigint @@ + Int.fold (fun acc i -> + match Int.Interval.view i with + | { lb = Closed lb ; ub = Closed ub } -> Z.(acc + ub - lb + one) + | _ -> raise Exit + ) Z.zero u + ) + with Exit -> None + ) + | _ -> None + + let integer_hull = function + | Real _ -> invalid_arg "integer_hull" + | Int u -> + let lb, ex_lb = Int.lower_bound u in + let ub, ex_ub = Int.upper_bound u in + let lb = + match lb with + | Closed v -> Some (v, ex_lb) + | _ -> None in - let bl = match b2 with - | Strict v -> Large v - | Large v -> Strict v - | _ -> b2 + let ub = + match ub with + | Closed v -> Some (v, ex_ub) + | _ -> None in - if bu == Minfty then complement r bl acc - else complement r bl ((prev, bu)::acc) - | [] -> - List.rev (if prev == Pinfty then acc else (prev, Pinfty)::acc) - in - fun uints1 uints2 -> - let l_c = complement uints1.ints Minfty [] in - let l_c = - if uints2.is_int then List.rev (List.rev_map int_bornes l_c) else l_c - in - let uints1_c = union_intervals {uints1 with ints = l_c} in - intersect uints1_c uints2 - -let add_borne b1 b2 = - match b1,b2 with - | Minfty, Pinfty | Pinfty, Minfty -> assert false - | Minfty, _ | _, Minfty -> Minfty - | Pinfty, _ | _, Pinfty -> Pinfty - | Large (v1, e1), Large (v2, e2) -> - Large (Q.add v1 v2, Ex.union e1 e2) - | (Large (v1, e1) | Strict (v1, e1)), (Large (v2, e2) | Strict (v2, e2)) -> - Strict (Q.add v1 v2, Ex.union e1 e2) - -let translate c ((b1, b2) as b) = - if Q.(equal zero) c then b - else begin - let tmp = Large (c, Ex.empty) in - (add_borne b1 tmp, add_borne b2 tmp) - end - -let scale_interval_zero n (b1, b2) = - assert (Q.sign n = 0); - Large (Q.zero, explain_borne b1), Large (Q.zero, explain_borne b2) - -let scale_borne_non_zero n b = - assert (Q.sign n > 0); - match b with - | Pinfty | Minfty -> b - | Large (v, e) -> Large (Numbers.Q.mult n v, e) - | Strict (v, e) -> Strict (Numbers.Q.mult n v, e) - -let scale_interval_pos n (b1, b2) = - scale_borne_non_zero n b1, scale_borne_non_zero n b2 - -let scale_interval_neg n (b1, b2) = - minus_borne (scale_borne_non_zero (Numbers.Q.minus n) b2), - minus_borne (scale_borne_non_zero (Numbers.Q.minus n) b1) - - -let affine_scale ~const ~coef uints = - Options.tool_req 4 "TR-Arith-Axiomes scale"; - let sgn = Q.sign coef in - let aux = - if sgn = 0 then scale_interval_zero - else if sgn > 0 then scale_interval_pos - else scale_interval_neg - in - let rl = List.rev_map (fun bornes -> - translate const (aux coef bornes) - ) uints.ints in - let l = - if uints.is_int then rev_normalize_int_bounds rl uints.expl coef - else List.rev rl - in - let res = union_intervals { uints with ints = l } in - assert (res.ints != []); - res - -let scale coef uints = - affine_scale ~const:Q.zero ~coef uints - -let coerce ty uints = - match ty with - | Ty.Treal -> { uints with is_int = false; } - | Ty.Tint -> - if uints.is_int then uints - else scale Q.one { uints with is_int = true; } - | _ -> assert false - -let add_interval is_int l (b1,b2) = - List.fold_right - (fun (b1', b2') l -> - let l1 = ((add_borne b1 b1'),(add_borne b2 b2'))::l in - union_bornes is_int (l1) - ) l [] - -let add {ints = l1; is_int = is_int; expl = e1} { ints = l2; expl = e2; _ } = - let l = - List.fold_left - (fun l bs -> let i = add_interval is_int l1 bs in i@l) [] l2 - in - let res = union_intervals { ints = l ; is_int; expl = Ex.union e1 e2 } in - assert (res.ints != []); - res - -let sub i1 i2 = add i1 (scale Numbers.Q.m_one i2) - -let merge i1 i2 = - union_intervals - {ints = List.rev_append i1.ints i2.ints; is_int = i1.is_int; - expl = Explanation.union i1.expl i2.expl} - -let contains i q = - List.exists - (fun (b1, b2) -> - begin - match b1 with - | Minfty -> true - | Pinfty -> assert false - | Large(v, _) -> Q.compare v q <= 0 - | Strict(v, _) -> Q.compare v q < 0 - end - && begin - match b2 with - | Pinfty -> true - | Minfty -> assert false - | Large(v, _) -> Q.compare v q >= 0 - | Strict(v, _) -> Q.compare v q > 0 - end - ) i.ints - -let doesnt_contain_0 = - let rec explain_no_zero l = - match l with - | [] -> assert false - | (b1, b2) :: l -> - let pos_strict_b1 = low_borne_pos_strict b1 in - let pos_strict_b2 = up_borne_pos_strict b2 in - if pos_strict_b1 && pos_strict_b2 then - (* there is no negative values at all *) - Some (Ex.union (explain_borne b1) (explain_borne b2), []) - else - begin - (* we know l does not contain 0. So, these asserts should hold *) - assert (not pos_strict_b1); - assert (not pos_strict_b2); - assert (low_borne_neg_strict b1); - assert (up_borne_neg_strict b2); - match l with - | [] -> - (* there is no positive values at all *) - Some (Ex.union (explain_borne b1) (explain_borne b2), []) - | (c1,_)::_ -> - if low_borne_pos_strict c1 then - Some (Ex.union (explain_borne b2) (explain_borne c1), []) - else explain_no_zero l - end - in - fun int -> - if contains int Q.zero then None - else explain_no_zero int.ints - - -let is_positive { ints; expl; _ } = - match ints with - | [] -> assert false - | (lb,_)::_ -> if pos_borne lb then Some (expl, []) else None - -let root_default_num v n = - if n = 2 then Numbers.Q.sqrt_default v else Numbers.Q.root_default v n - -let root_exces_num v n = - if n = 2 then Numbers.Q.sqrt_excess v else Numbers.Q.root_excess v n - -(* should be removed and replaced with compare_bounds, with makes distinction - between lower and upper bounds *) -let compare_bornes b1 b2 = - match b1, b2 with - | Minfty, Minfty | Pinfty, Pinfty -> 0 - | Minfty, _ | _, Pinfty -> -1 - | Pinfty, _ | _, Minfty -> 1 - | Strict (v1, _), Strict (v2, _) | Large (v1, _), Large (v2, _) - | Strict (v1, _), Large (v2, _) | Large (v1, _), Strict (v2, _) -> - Q.compare v1 v2 - - -let is_strict_smaller = - let rec aux l1 l2 nb_eq sz_l1 sz_l2 = - match l1, l2 with - [], _ -> true, nb_eq, sz_l1, (sz_l2 + List.length l2) - | _, [] -> false, nb_eq, (sz_l1 + List.length l1), sz_l2 - | b1::r1, b2::r2 -> - let lo1, up1 = b1 in - let lo2, up2 = b2 in - let c_l1_l2 = compare_bounds lo1 ~is_low1:true lo2 ~is_low2:true in - let c_u1_u2 = compare_bounds up1 ~is_low1:false up2 ~is_low2:false in - let c_l1_u2 = compare_bounds lo1 ~is_low1:true up2 ~is_low2:false in - let c_u1_l2 = compare_bounds up1 ~is_low1:false lo2 ~is_low2:true in - if c_l1_l2 = 0 && c_u1_u2 = 0 then - aux r1 r2 (nb_eq + 1) (sz_l1 + 1) (sz_l2 + 1) - else - if c_l1_l2 >= 0 && c_u1_u2 <= 0 then (* without being equal *) - (* b1 \subset b2! look for inclusion of r1 in l2 *) - aux r1 l2 nb_eq (sz_l1 + 1) sz_l2 - else - if c_l1_u2 >= 0 then - (*ignore b2, and look for inclusion of l1 in r2*) - aux l1 r2 nb_eq sz_l1 (sz_l2 + 1) - else - if c_u1_l2 < 0 then - raise Exit(* b1 is not included in any part of l2*) - else - if c_l1_l2 <= 0 && c_u1_u2 >= 0 then (* without being equal *) - raise Exit (*no inclusion, we have b2 \subset b1 !! *) - else - if c_l1_l2 < 0 && c_u1_u2 < 0 || - c_l1_l2 > 0 && c_u1_u2 > 0 then - raise Exit (* intersection and differences are not empty *) - else - assert false - in - fun i1 i2 -> - try - let res, nb_eq, sz_l1, sz_l2 = aux i1.ints i2.ints 0 0 0 in - (* if res is true, we should check that i1 and i2 are not equal *) - res && (sz_l1 <> sz_l2 || nb_eq <> sz_l1) - with Exit -> false - - -let mult_borne b1 b2 = - match b1,b2 with - | Minfty, Pinfty | Pinfty, Minfty -> assert false - | Minfty, b | b, Minfty -> - if compare_bornes b (large_borne_of Q.zero Ex.empty) = 0 - then b - else if pos_borne b then Minfty - else Pinfty - | Pinfty, b | b, Pinfty -> - if compare_bornes b (large_borne_of Q.zero Ex.empty) = 0 - then b - else if pos_borne b then Pinfty - else Minfty - | Strict (_, e1), Large (v, e2) - | Large (v, e1), Strict (_, e2) when Numbers.Q.is_zero v -> - Large (Q.zero, Ex.union e1 e2) - - | Strict (v1, e1), Strict (v2, e2) | Strict (v1, e1), Large (v2, e2) - | Large (v1, e1), Strict (v2, e2) -> - Strict (Numbers.Q.mult v1 v2, Ex.union e1 e2) - | Large (v1, e1), Large (v2, e2) -> - Large (Numbers.Q.mult v1 v2, Ex.union e1 e2) - -let mult_borne_inf b1 b2 = - match b1,b2 with - | Minfty, Pinfty | Pinfty, Minfty -> Minfty - | _, _ -> mult_borne b1 b2 - -let mult_borne_sup b1 b2 = - match b1,b2 with - | Minfty, Pinfty | Pinfty, Minfty -> Pinfty - | _, _ -> mult_borne b1 b2 - - -let mult_bornes (a,b) (c,d) = - (* see ../extra/intervals_mult.png *) - (* put the explanation of both bounds for multiplication. Putting - just one of them is probably incorrect. When a bound is [0,0], put - the explanation of that bound instead of empty. - TODO: make a deeper study!!!*) - let ex_a_b = Ex.union (explain_borne a) (explain_borne b) in - let ex_c_d = Ex.union (explain_borne c) (explain_borne d) in - let all_ex = Ex.union ex_a_b ex_c_d in - match class_of a b, class_of c d with - | P, P -> - mult_borne_inf a c, mult_borne_sup b d, all_ex - | P, M -> - mult_borne_inf b c, mult_borne_sup b d, all_ex - | P, N -> - mult_borne_inf b c, mult_borne_sup a d, all_ex - | M, P -> - mult_borne_inf a d, mult_borne_sup b d, all_ex - | M, M -> - min_of_lower_bounds (mult_borne_inf a d) (mult_borne_inf b c), - max_of_upper_bounds (mult_borne_sup a c) (mult_borne_sup b d), - all_ex - | M, N -> - mult_borne_inf b c, mult_borne_sup a c, all_ex - | N, P -> - mult_borne_inf a d, mult_borne_sup b c, all_ex - | N, M -> - mult_borne_inf a d, mult_borne_sup a c, all_ex - | N, N -> - mult_borne_inf b d, mult_borne_sup a c, all_ex - | Z, (P | M | N | Z) -> (a, b, ex_a_b) - | (P | M | N ), Z -> (c, d, ex_c_d) - -let rec power_borne_inf p b = - match p with - | 1 -> b - | p -> mult_borne_inf b (power_borne_inf (p-1) b) - -let rec power_borne_sup p b = - match p with - | 1 -> b - | p -> mult_borne_sup b (power_borne_sup (p-1) b) - -let max_merge b1 b2 = - let ex = Ex.union (explain_borne b1) (explain_borne b2) in - let max = max_of_upper_bounds b1 b2 in - match max with - | Minfty | Pinfty -> max - | Large (v, _) -> Large (v, ex) - | Strict (v, _) -> Strict (v, ex) - -let power_bornes p (b1,b2) = - if neg_borne b1 && pos_borne b2 then - match p with - | 0 -> assert false - | p when p mod 2 = 0 -> - (* max_merge to have explanations !!! *) - let m = max_merge (power_borne_sup p b1) (power_borne_sup p b2) in - (Large (Q.zero, Ex.empty), m) - | _ -> (power_borne_inf p b1, power_borne_sup p b2) - else if pos_borne b1 && pos_borne b2 then - (power_borne_inf p b1, power_borne_sup p b2) - else if neg_borne b1 && neg_borne b2 then - match p with - | 0 -> assert false - | p when p mod 2 = 0 -> (power_borne_inf p b2, power_borne_sup p b1) - | _ -> (power_borne_inf p b1, power_borne_sup p b2) - else assert false - -let int_div_of_borne_inf min_f b = - match b with - | Minfty | Pinfty -> b - | Large (v, e) -> - Large ((if Numbers.Q.is_int v then v else (*Q.floor*) min_f v), e) - - | Strict (v, e) -> - (* this case really happens ?? *) - if Numbers.Q.is_int v then Large (Q.add v Q.one, e) - else - let v' = (*Q.floor*) min_f v in (* correct ? *) - assert (Q.compare v' v > 0); - Large (v', e) - -let int_div_of_borne_sup max_f b = - match b with - | Minfty | Pinfty -> b - | Large (v, e) -> - Large ((if Numbers.Q.is_int v then v else (*Q.floor*) max_f v), e) - - | Strict (v, e) -> - (* this case really happens ?? *) - if Numbers.Q.is_int v then Large (Q.sub v Q.one, e) - else - let v' = (*Q.floor*) max_f v in (* correct ? *) - assert (Q.compare v' v < 0); - Large (v', e) - -(* we use int_div_bornes for division of integer intervals instead of - int_bornes because the div op is Euclidean division is this case *) -let int_div_bornes min_f max_f l u = - int_div_of_borne_inf min_f l, int_div_of_borne_sup max_f u - -let mult u1 u2 = - Options.tool_req 4 "TR-Arith-Axiomes mult"; - let resl, expl = - List.fold_left - (fun (l', expl) b1 -> - List.fold_left - (fun (l, ex) b2 -> - let bl, bu, ex' = mult_bornes b1 b2 in - let bl = add_expl_to_borne bl ex' in - let bu = add_expl_to_borne bu ex' in - (bl, bu)::l, Ex.union ex ex') (l', expl) u2.ints) - ([], Ex.empty) u1.ints - in - let res = - union_intervals - { ints= resl; is_int = u1.is_int; - expl = Ex.union expl (Ex.union u1.expl u2.expl) } in - assert (res.ints != []); - res - - -let power n u = - Options.tool_req 4 "TR-Arith-Axiomes power"; - let l = List.map (power_bornes n) u.ints in - union_intervals { u with ints = l } - -let root_default_borne is_int x n = - match x with - | Pinfty -> Pinfty - | Minfty -> Minfty - | Large (v, e) | Strict (v, e) -> - let sign, s = - if Q.sign v >= 0 then (fun q -> q), root_default_num v n - else Numbers.Q.minus, root_exces_num (Numbers.Q.minus v) n - in - match s with - | None -> Minfty - | Some s -> - let s = sign s in - if is_int then - let cs = Numbers.Q.ceiling s in - let cs2 = Numbers.Q.power cs n in - if Q.compare v cs2 <= 0 then Large (cs, e) - else Large (Q.add cs Q.one, e) - else Large (s, e) - -let root_exces_borne is_int x n = - match x with - | Pinfty -> Pinfty - | Minfty -> Minfty - | Large (v, e) | Strict (v, e) -> - let sign, s = - if Q.sign v >= 0 then (fun d -> d), root_exces_num v n - else Numbers.Q.minus, root_default_num (Numbers.Q.minus v) n + (lb, ub) + + let map_bound f = function + | Unbounded -> Unbounded + | Open x -> Open (f x) + | Closed x -> Closed (f x) + + let lower_bound = function + | Real u -> Real.lower_bound u + | Int u -> + let lb, ex_lb = Int.lower_bound u in + map_bound Q.of_bigint lb, ex_lb + + let borne_inf u = + let lb, ex = lower_bound u in + match lb with + | Open b -> b, ex, false + | Closed b -> b, ex, true + | Unbounded -> raise No_finite_bound + + let upper_bound = function + | Real u -> Real.upper_bound u + | Int u -> + let ub, ex_ub = Int.upper_bound u in + map_bound Q.of_bigint ub, ex_ub + + let borne_sup u = + let ub, ex = upper_bound u in + match ub with + | Open b -> b, ex, false + | Closed b -> b, ex, true + | Unbounded -> raise No_finite_bound + + let div = map2 "div" Real.div Int.ediv + + let coerce ty u = + match ty, u with + | Ty.Treal, Real _ -> u + | Treal, Int u -> + of_real @@ Real.map_strict_inc QAlgebraicType.of_bigint u + | Tint, Int _ -> u + | Tint, Real u -> + of_int @@ check @@ + Int.partial_map_inc QAlgebraicType.ceiling QAlgebraicType.floor u + | _ -> Fmt.invalid_arg "coerce" + + let contains u v = + match u with + | Real u -> Real.( + match intersect u (of_interval @@ Interval.singleton v) with + | NonEmpty _ -> true + | Empty _ -> false + ) + | Int u -> Int.( + if Z.equal (Q.den v) Z.one then + let v = Q.num v in + match intersect u (of_interval @@ Interval.singleton v) with + | NonEmpty _ -> true + | Empty _ -> false + else + false + ) + + let add_explanation u ex = + match u with + | Real u -> Real (Real.add_explanation ~ex u) + | Int u -> Int (Int.add_explanation ~ex u) + + let equal u1 u2 = + match u1, u2 with + | Real u1, Real u2 -> Real.equal u1 u2 + | Int u1, Int u2 -> Int.equal u1 u2 + | _ -> invalid_arg "equal" + + let pick ~is_max u = + match u with + | Real u -> + let pick_interval lb ub = + match lb, ub with + | Unbounded, Unbounded -> Q.zero + + | _, Closed ub when is_max -> ub + + | Unbounded, Open ub when is_max -> Q.(ub - ~$1) + | (Open lb | Closed lb), Open ub when is_max -> + Q.(ub - min ~$1 ((ub - lb) / ~$2)) + + | (Open lb | Closed lb), Unbounded when is_max -> Q.(lb + ~$1) + + | Closed lb, _ -> lb + | Open lb, Unbounded -> Q.(lb + ~$1) + | Open lb, (Open ub | Closed ub) -> + Q.(lb + min ~$1 ((ub - lb) / ~$2)) + | Unbounded, (Open ub | Closed ub) -> Q.(ub - ~$1) + in + let lb, _ = Real.lower_bound u in + let ub, _ = Real.upper_bound u in + Real.fold (fun acc i -> + let i = Real.Interval.view i in + (if is_max then max else min) acc (pick_interval i.lb i.ub) + ) (pick_interval lb ub) u + | Int u -> + let pick_interval lb ub = + match lb, ub with + | Unbounded, Unbounded -> Z.zero + + | _, Closed ub when is_max -> ub + | Closed lb, Unbounded when is_max -> Z.succ lb + + | Closed lb, _ -> lb + | Unbounded, Closed ub -> Z.pred ub + + | Open _, _ | _, Open _ -> assert false + in + let lb, _ = Int.lower_bound u in + let ub, _ = Int.upper_bound u in + let n = + Int.fold (fun acc i -> + let i = Int.Interval.view i in + (if is_max then max else min) acc (pick_interval i.lb i.ub) + ) (pick_interval lb ub) u + in + Q.of_bigint n + + let fold f acc u = + match u with + | Real u -> + Real.fold (fun acc i -> f acc (Real.Interval.view i)) acc u + | Int u -> + Int.fold (fun acc i -> + f acc + (Real.Interval.view + { lb = QAlgebraicType.of_bigint i.lb + ; ub = QAlgebraicType.of_bigint i.ub }) + ) acc u + + type interval_matching = + ((Q.t * bool) option * (Q.t * bool) option * Ty.t) Var.Map.t + + module MV = Var.Map + module Sy = Symbols + + let consistent_bnds low up = + match low, up with + | Some (q_low, str_low), Some (q_up, str_up) -> + let c = Q.compare q_up q_low in + c > 0 || (c = 0 && not str_low && not str_up) + | _ -> true + + let new_up_bound idoms s ty q is_strict = + let old_low, old_up, ty = + try MV.find s idoms + with Not_found -> None,None,ty in - match s with - | None -> Pinfty - | Some s -> - let s = sign s in - if is_int then - let cs = Numbers.Q.floor s in - let cs2 = Numbers.Q.power cs n in - if Q.compare v cs2 >= 0 then Large (cs, e) - else Large (Q.sub cs Q.one, e) - else Large (s, e) - -let sqrt_interval is_int (l, ex) (b1,b2) = - let l1 = minus_borne (root_exces_borne is_int b2 2) in - let u1 = minus_borne (root_default_borne is_int b1 2) in - let l2 = root_default_borne is_int b1 2 in - let u2 = root_exces_borne is_int b2 2 in - let c1 = compare_bornes l1 u1 in - let c2 = compare_bornes l2 u2 in - if c1 > 0 then - if c2 > 0 then - l, Ex.union ex (Ex.union (explain_borne b1) (explain_borne b2)) - else (l2,u2)::l, ex - else - if c2 > 0 then (l1, u1)::l, ex - else (union_bornes is_int [(l1,u1); (l2, u2)]) @ l, ex - -let sqrt {ints = l; is_int = is_int; expl} = - Options.tool_req 4 "TR-Arith-Axiomes sqrt"; - let l, expl = List.fold_left (sqrt_interval is_int) ([], expl) l in - if l == [] then raise (NotConsistent expl); - let res = union_intervals { ints = l; is_int; expl} in - assert (res.ints != []); - res - -let root_interval is_int (b1,b2) n = - let u,l = (root_default_borne is_int b1 n, root_exces_borne is_int b2 n) in - if compare_bornes u l <= 0 then Int(u,l) - else Empty (Ex.union (explain_borne b1) (explain_borne b2)) - -let rec root n ({ints = l; is_int = is_int; expl } as u) = - Options.tool_req 4"TR-Arith-Axiomes root"; - if n mod 2 = 0 then root (n/2) (sqrt u) - else - let l = List.rev_map (fun bs -> root_interval is_int bs n) l in - let l = join (List.rev l) expl in - let res = union_intervals {u with ints = l; is_int = is_int} in - assert (res.ints != []); - res - -let inv_borne_inf b is_int ~other = - match b with - | Pinfty -> assert false - | Minfty -> - if is_int then Large (Q.zero, explain_borne other) - else Strict (Q.zero, explain_borne other) - | Strict (c, _) | Large (c, _) when Q.sign c = 0 -> Pinfty - | Strict (v, e) -> Strict (Q.div Q.one v, e) - | Large (v, e) -> Large (Q.div Q.one v, e) - -let inv_borne_sup b is_int ~other = - match b with - | Minfty -> assert false - | Pinfty -> - if is_int then Large (Q.zero, explain_borne other) - else Strict (Q.zero, explain_borne other) - | Strict (c, _) | Large (c, _) when Q.sign c = 0 -> Minfty - | Strict (v, e) -> Strict (Q.div Q.one v, e) - | Large (v, e) -> Large (Q.div Q.one v, e) - -let inv_bornes (l, u) is_int = - inv_borne_sup u is_int ~other:l, inv_borne_inf l is_int ~other:u - - -let inv ({ ints = l; is_int; _ } as u) = - match doesnt_contain_0 u with - | None -> { u with ints = [Minfty, Pinfty] } - | Some (ex, _) -> - let l' = - List.fold_left - (fun acc (l,u) -> - let l = add_expl_to_borne l ex in - let u = add_expl_to_borne u ex in - (inv_bornes (l, u) is_int) :: acc - ) [] l + let new_up = + match old_up with + | None -> Some (q, is_strict) + | Some (r, str) -> + let c = Q.compare r q in + if c < 0 then old_up + else if c > 0 then Some (q, is_strict) + else + if is_strict == str || is_strict then old_up + else Some (q, is_strict) in - assert (l' != []); - (* ! SHOULD NOT try to simplify here if is_int is true *) - union_intervals { u with ints = l' } - - -type sign_of_interval = Zero | Pos | Neg | Mixed - -let sign_of_interval { ints; _ } = - match ints, List.rev ints with - | [], _ | _, [] -> assert false - | (inf, _)::_, (_,sup)::_ -> - match inf, sup with - | Pinfty, _ | _, Minfty -> assert false - | Minfty, Pinfty -> Mixed - | Large(v,_), Large(v',_) -> - if Q.compare v Q.zero > 0 then Pos - else if Q.compare v' Q.zero < 0 then Neg - else if Numbers.Q.is_zero v && Numbers.Q.is_zero v' then Zero - else Mixed - - | (Strict(v,_) | Large(v,_)), (Strict(v',_) | Large(v',_)) -> - if Q.compare v Q.zero >= 0 then Pos - else if Q.compare v' Q.zero <= 0 then Neg - else Mixed - - | (Strict(v,_) | Large(v,_)), Pinfty -> - if Q.compare v Q.zero >= 0 then Pos - else Mixed - - | Minfty, (Strict(v',_) | Large(v',_)) -> - if Q.compare v' Q.zero <= 0 then Neg - else Mixed - - -let div i1 i2 = - Options.tool_req 4 "TR-Arith-Axiomes div"; - let inv_i2 = inv i2 in - if is_undefined inv_i2 then inv_i2 - else - let i1 = match doesnt_contain_0 i2 with - | Some (ex, _) -> add_expl_zero i1 ex - | None -> i1 + if new_up == old_up then idoms + else + if consistent_bnds old_low new_up then MV.add s (old_low, new_up, ty) idoms + else raise Exit + + let new_low_bound idoms s ty q is_strict = + let old_low, old_up, ty = + try MV.find s idoms + with Not_found -> None,None,ty in - let ({ ints = l; is_int = is_int; _ } as i) = mult i1 inv_i2 in - assert (l != []); - if is_int then - (* not just int_bornes because it's Euclidean division *) - let min_f, max_f = - match sign_of_interval i2 with - | Zero -> assert false (* inv_i2 is not undefined *) - | Pos -> Numbers.Q.floor, Numbers.Q.floor - | Neg -> Numbers.Q.ceiling, Numbers.Q.ceiling - | Mixed -> Numbers.Q.floor, Numbers.Q.ceiling - in - let rl = List.rev_map (fun (l,u) -> int_div_bornes min_f max_f l u) l in - union_intervals { i with ints = List.rev rl } - else { i with ints = l } - -let abs = - let zero_inf_r = - new_borne_inf Ex.empty Q.zero ~is_le:true (undefined Ty.Treal) - in - let zero_inf_i = - new_borne_inf Ex.empty Q.zero ~is_le:true (undefined Ty.Tint) - in - fun i -> - let xx = if i.is_int then zero_inf_i else zero_inf_r in - intersect (merge i (scale Numbers.Q.m_one i)) xx - -let mk_closed l u llarge ularge lexp uexp ty = - let lb = if llarge then Large(l, lexp) else Strict (l, lexp) in - let ub = if ularge then Large(u, uexp) else Strict (u, uexp) in - { ints = [lb, ub]; is_int = ty == Ty.Tint; expl = Ex.union lexp uexp } - -type bnd = (Numbers.Q.t * Numbers.Q.t) option * Explanation.t - -let bnd_of_borne b ex0 low = - match b with - | Pinfty when not low -> None, ex0 - | Minfty when low -> None, ex0 - | Pinfty | Minfty -> assert false - | Large (c, ex) -> Some (c, Q.zero), Ex.union ex0 ex - | Strict (c, ex) -> - Some (c, if low then Q.one else Numbers.Q.m_one), Ex.union ex0 ex - -let bounds_of env = - let ex = env.expl in - match env.ints with - | [] -> [(None, ex), (None, ex)] - | l -> - List.rev - (List.rev_map - (fun (b1, b2) -> bnd_of_borne b1 ex true, bnd_of_borne b2 ex false) l) - -let add_explanation i ex = - if Ex.is_empty ex then i - else - let rl = - List.rev_map (fun (l, u) -> - add_expl_to_borne l ex, add_expl_to_borne u ex) i.ints + let new_low = + match old_low with + | None -> Some (q, is_strict) + | Some (r, str) -> + let c = Q.compare r q in + if c > 0 then old_low + else if c < 0 then Some (q, is_strict) + else + if is_strict == str || is_strict then old_low + else Some (q, is_strict) in - {i with ints = List.rev rl; expl = Ex.union i.expl ex} - -let equal i1 i2 = - try - List.iter2 - (fun (b1,c1) (b2,c2) -> - if compare_bounds b1 ~is_low1:true b2 ~is_low2:true <> 0 || - compare_bounds c1 ~is_low1:false c2 ~is_low2:false <> 0 then - raise Exit - )i1.ints i2.ints; - true - with Exit | Invalid_argument _ -> false - -let min_bound { ints; _ } = - match ints with - | [] -> assert false - | (b, _) :: _ -> b - -let max_bound { ints; _} = - match List.rev ints with - | [] -> assert false - | (_, b) :: _ -> b - -type interval_matching = - ((Q.t * bool) option * (Q.t * bool) option * Ty.t) Var.Map.t - -module MV = Var.Map -module Sy = Symbols - -let consistent_bnds low up = - match low, up with - | Some (q_low, str_low), Some (q_up, str_up) -> - let c = Q.compare q_up q_low in - c > 0 || (c = 0 && not str_low && not str_up) - | _ -> true - -let new_up_bound idoms s ty q is_strict = - let old_low, old_up, ty = - try MV.find s idoms - with Not_found -> None,None,ty - in - let new_up = - match old_up with - | None -> Some (q, is_strict) - | Some (r, str) -> - let c = Q.compare r q in - if c < 0 then old_up - else if c > 0 then Some (q, is_strict) - else - if is_strict == str || is_strict then old_up - else Some (q, is_strict) - in - if new_up == old_up then idoms - else - if consistent_bnds old_low new_up then MV.add s (old_low, new_up, ty) idoms - else raise Exit - -let new_low_bound idoms s ty q is_strict = - let old_low, old_up, ty = - try MV.find s idoms - with Not_found -> None,None,ty - in - let new_low = - match old_low with - | None -> Some (q, is_strict) - | Some (r, str) -> - let c = Q.compare r q in - if c > 0 then old_low - else if c < 0 then Some (q, is_strict) - else - if is_strict == str || is_strict then old_low - else Some (q, is_strict) - in - if new_low == old_low then idoms - else - if consistent_bnds new_low old_up then MV.add s (new_low, old_up, ty) idoms - else raise Exit - -let new_var idoms s ty = - if MV.mem s idoms then idoms - else MV.add s (None, None, ty) idoms - -let match_interval_upper {Sy.sort; is_open; kind; is_lower} i imatch = - assert (not is_lower); - match kind, max_bound i with - | Sy.Unbounded, _ -> imatch (* ? var *) - | Sy.VarBnd _, Minfty -> assert false - | Sy.VarBnd s, Pinfty -> new_var imatch s sort - | Sy.VarBnd s, Strict (v, _) -> new_low_bound imatch s sort v false - | Sy.VarBnd s, Large (v, _) -> new_low_bound imatch s sort v is_open - - | Sy.ValBnd _, Minfty -> assert false - | Sy.ValBnd _, Pinfty -> raise Exit - | Sy.ValBnd vl, Strict (v, _) -> - let c = Q.compare v vl in - if c > 0 then raise Exit; - imatch - - | Sy.ValBnd vl, Large (v, _) -> - let c = Q.compare v vl in - if c > 0 || c = 0 && is_open then raise Exit; - imatch - - -let match_interval_lower {Sy.sort; is_open; kind; is_lower} i imatch = - assert (is_lower); - match kind, min_bound i with - | Sy.Unbounded, _ -> imatch (* ? var *) - | Sy.VarBnd _, Pinfty -> assert false - | Sy.VarBnd s, Minfty -> new_var imatch s sort - | Sy.VarBnd s, Strict (v, _) -> new_up_bound imatch s sort v false - | Sy.VarBnd s, Large (v, _) -> new_up_bound imatch s sort v is_open - - | Sy.ValBnd _, Minfty -> raise Exit - | Sy.ValBnd _, Pinfty -> assert false - | Sy.ValBnd vl, Strict (v, _) -> - let c = Q.compare v vl in - if c < 0 then raise Exit; - imatch - - | Sy.ValBnd vl, Large (v, _) -> - let c = Q.compare v vl in - if c < 0 || c = 0 && is_open then raise Exit; - imatch -let match_interval lb ub i accu = - try Some (match_interval_upper ub i (match_interval_lower lb i accu)) - with Exit -> None - -(* Assumes: the input set of intervals is normalized. *) -let pick ~is_max { ints; is_int; _ } = - let ints = if is_max then List.rev ints else ints in - match ints with - | [] -> None - | (Minfty, Pinfty) :: _ -> Some Q.zero - | (_, Large (q, _)) :: _ when is_max -> Some q - | (_, Strict(q, _)) :: _ when is_max && is_int -> - (* By normalization, an integer interval of the form |p, q) has to - contain at least one integer and thus [q-1] is an element of this - interval. *) - Some Q.(q - ~$1) - - | (Large (q, _), _) :: _ when not is_max -> Some q - | (Strict (q, _), _) :: _ when not is_max && is_int -> - (* By normalization, an integer interval of the form (q, p| has to - contain at least one integer and thus [q+1] is an element of this - interval. *) - Some Q.(q + ~$1) - - | (Minfty, (Strict (q, _) | Large (q, _))) :: _ -> Some Q.(q - ~$1) - | ((Strict (q, _) | Large (q, _)), Pinfty) :: _ -> Some Q.(q + ~$1) - | ((Strict (q1, _) | Large (q1, _)), (Strict (q2, _) | Large (q2, _))) :: _ -> - begin - assert (not is_int); - Some Q.((q1 + q2) / ~$2) - end - | (_, Minfty) :: _ | (Pinfty, _) :: _ -> - (* As the set of intervals is normalized, it cannot contain - empty intervals. *) - assert false - -(*****************) - -(* Some debug code for Intervals: commented by default - - let no_inclusion = - let not_included (b1, c1) (b2, c2) = - not ( - compare_bounds b1 ~is_low1:true b2 ~is_low2:true >= 0 && - compare_bounds c1 ~is_low1:false c2 ~is_low2:false <= 0 - ) - in - let b_inc_list d l = - List.iter (fun e -> - assert (not_included d e); - assert (not_included e d) - ) l - in - let rec aux todo = - match todo with - | [] -> assert false - | [e] -> () - | d::l -> b_inc_list d l; aux l - in - fun i -> - (*fprintf fmt "[no_inclusion] i = %a@." print i;*) - aux i.ints - - let not_mergeable = - let rec aux is_int l = - match l with - | [] -> assert false - | [e] -> () - | (_,a)::(((b,_)::_) as l) -> - begin - match a, b with - | Minfty, _ | _, Pinfty -> assert false (*should not happen*) - | Pinfty, _ | _, Minfty -> assert false - (*should not happen or not norm*) - | Large(q1,_) , Large(q2,_) -> - assert (Q.compare q1 q2 < 0); (* otherwise, we can merge *) - if is_int then - (* the gap between q1 and q2 should be > 1. Otherwise, - we can merge *) - assert (Q.compare (Q.sub q2 q1) Q.one > 0) - - | Strict(q1,_), Large(q2,_) -> - assert (not is_int); - assert (Q.compare q1 q2 < 0) (* otherwise, we can merge *) - - | Large(q1,_) , Strict(q2,_) -> - assert (not is_int); - assert (Q.compare q1 q2 < 0) (* otherwise, we can merge *) - - | Strict(q1,_) , Strict(q2,_) -> - assert (not is_int); - assert (Q.compare q1 q2 <= 0) (* otherwise, we can merge *) - end; - aux is_int l; - in - fun i -> - (*fprintf fmt "[no_mergeable] i = %a@." print i;*) - aux i.is_int i.ints - - let assert_is_normalized i = - not_mergeable i; - no_inclusion i; - i - - let exclude i j = - try - let k = exclude i j in - k |> assert_is_normalized - with Assert_failure _ -> assert false - - let intersect i j = - try - let k = intersect i j in - k |> assert_is_normalized - with Assert_failure _ -> assert false - - let mult i j = - try mult i j |> assert_is_normalized - with Assert_failure _ -> assert false - - let power i j = - try power i j |> assert_is_normalized - with Assert_failure _ -> assert false - - let sqrt i = - try sqrt i |> assert_is_normalized - with Assert_failure _ -> assert false - - let root n i = - try root n i |> assert_is_normalized - with Assert_failure _ -> assert false - - let add i j = - try - (*fprintf fmt "@.i = %a@." print i; - fprintf fmt "j = %a@." print j;*) - let k = add i j in - (*fprintf fmt "res = %a@." print k;*) - k |> assert_is_normalized - with Assert_failure _ -> assert false - - let scale q i = - try scale q i |> assert_is_normalized - with Assert_failure _ -> assert false - - let sub i j = - try sub i j |> assert_is_normalized - with Assert_failure _ -> assert false - - let merge i j = - try merge i j |> assert_is_normalized - with Assert_failure _ -> assert false - - let abs i = - try abs i |> assert_is_normalized - with Assert_failure _ -> assert false - - let div i j = - try div i j |> assert_is_normalized - with Assert_failure _ -> assert false -*) + if new_low == old_low then idoms + else + if consistent_bnds new_low old_up then MV.add s (new_low, old_up, ty) idoms + else raise Exit + + let new_var idoms s ty = + if MV.mem s idoms then idoms + else MV.add s (None, None, ty) idoms + + let match_interval_upper {Symbols.sort; is_open; kind; is_lower} u imatch = + assert (not is_lower); + match kind, fst (upper_bound u) with + | Unbounded, _ -> imatch (* ? var *) + | VarBnd s, Unbounded -> new_var imatch s sort + | VarBnd s, Open v -> new_low_bound imatch s sort v false + | VarBnd s, Closed v -> new_low_bound imatch s sort v is_open + + | ValBnd _, Unbounded -> raise Exit + | ValBnd vl, Open v -> + if Q.compare v vl > 0 then raise Exit; + imatch + | ValBnd vl, Closed v -> + let c = Q.compare v vl in + if c > 0 || c = 0 && is_open then raise Exit; + imatch + + let match_interval_lower {Symbols.sort; is_open; kind; is_lower} u imatch = + assert is_lower; + match kind, fst (lower_bound u) with + | Unbounded, _ -> imatch (* ? var *) + + | VarBnd s, Unbounded -> new_var imatch s sort + | VarBnd s, Open v -> new_up_bound imatch s sort v false + | VarBnd s, Closed v -> new_up_bound imatch s sort v is_open + + | ValBnd _, Unbounded -> raise Exit + | ValBnd vl, Open v -> + if Q.compare v vl < 0 then raise Exit; + imatch + | ValBnd vl, Closed v -> + let c = Q.compare v vl in + if c < 0 || c = 0 && is_open then raise Exit; + imatch + + let match_interval lb ub u acc = + try Some (match_interval_upper ub u (match_interval_lower lb u acc)) + with Exit -> None +end + +module DebugExplanations : Explanations with type t = string list = struct + (** This module implements the [Explanations] interface by storing + explanations as (sorted) lists of strings for debugging purposes. *) + + type t = string list + + let pp ppf ex = + Format.fprintf ppf "@[[%a]@]" + (Format.pp_print_list + ~pp_sep:(fun ppf () -> Format.fprintf ppf ",@ ") + Format.pp_print_string) ex + + let is_empty = function + | [] -> true + | _ -> false + + let empty = [] + + let union l1 l2 = + List.rev_append l1 l2 |> List.sort_uniq String.compare + + let compare l1 l2 = + List.compare String.compare l1 l2 +end diff --git a/src/lib/reasoners/intervals.mli b/src/lib/reasoners/intervals.mli index 9e88828d4c..f5a26faf4f 100644 --- a/src/lib/reasoners/intervals.mli +++ b/src/lib/reasoners/intervals.mli @@ -1,7 +1,7 @@ (**************************************************************************) (* *) (* Alt-Ergo: The SMT Solver For Software Verification *) -(* Copyright (C) 2013-2023 --- OCamlPro SAS *) +(* Copyright (C) 2024 --- OCamlPro SAS *) (* *) (* This file is distributed under the terms of OCamlPro *) (* Non-Commercial Purpose License, version 1. *) @@ -19,140 +19,137 @@ (* *) (* CNRS - INRIA - Universite Paris Sud *) (* *) -(* Until 2013, some parts of this code were released under *) -(* the Apache Software License version 2.0. *) -(* *) (* --------------------------------------------------------------- *) (* *) (* More details can be found in the directory licenses/ *) (* *) (**************************************************************************) -type t - -exception NotConsistent of Explanation.t -exception No_finite_bound +open Intervals_intf -val undefined : Ty.t -> t +(** This module provides implementations of union-of-intervals over reals and + integers. *) -val is_undefined : t -> bool +type 'a union +(** Polymorphic union type. This allows writing conversion functions between + integer and real unions. *) -val point : Numbers.Q.t -> Ty.t -> Explanation.t -> t +module Real : AlgebraicField + with type explanation := Explanation.t + and type value := Q.t + and type 'a union = 'a union +(** Union-of-intervals over real numbers. *) -val doesnt_contain_0 : t -> Th_util.answer +module Int : EuclideanRing + with type explanation := Explanation.t + and type value := Z.t + and type 'a union = 'a union +(** Union-of-intervals over integers. *) -val is_positive : t -> Th_util.answer +module Legacy : sig + (** The [Legacy] module reimplements (most of) the old legacy [Intervals] + module on top of the current implementation to ease the transition. *) -val is_strict_smaller : t -> t -> bool + type t -val new_borne_sup : Explanation.t -> Numbers.Q.t -> is_le : bool -> t -> t + exception NotConsistent of Explanation.t + exception No_finite_bound -val new_borne_inf : Explanation.t -> Numbers.Q.t -> is_le : bool -> t -> t + val undefined : Ty.t -> t -val only_borne_sup : t -> t -(** Keep only the upper bound of the interval, - setting the lower bound to minus infty. *) + val point : Numbers.Q.t -> Ty.t -> Explanation.t -> t -val only_borne_inf : t -> t -(** Keep only the lower bound of the interval, - setting the upper bound to plus infty. *) + val doesnt_contain_0 : t -> Th_util.answer -val is_point : t -> (Numbers.Q.t * Explanation.t) option + val is_strict_smaller : t -> t -> bool -val intersect : t -> t -> t + val new_borne_sup : Explanation.t -> Numbers.Q.t -> is_le : bool -> t -> t -val exclude : t -> t -> t + val new_borne_inf : Explanation.t -> Numbers.Q.t -> is_le : bool -> t -> t -val mult : t -> t -> t + val only_borne_sup : t -> t + (** Keep only the upper bound of the interval, + setting the lower bound to minus infty. *) -val power : int -> t -> t + val only_borne_inf : t -> t + (** Keep only the lower bound of the interval, + setting the upper bound to plus infty. *) -val sqrt : t -> t + val is_point : t -> (Numbers.Q.t * Explanation.t) option -val root : int -> t -> t + val intersect : t -> t -> t -val add : t -> t -> t + val exclude : ?ex:Explanation.t -> Q.t -> t -> t -val scale : Numbers.Q.t -> t -> t + val mult : t -> t -> t -val affine_scale : const:Numbers.Q.t -> coef:Numbers.Q.t -> t -> t -(** Perform an affine transformation on the given bounds. - Suposing input bounds (b1, b2), this will return - (const + coef * b1, const + coef * b2). - This function is useful to avoid the incorrect roundings that - can take place when scaling down an integer range. *) + val power : int -> t -> t -val sub : t -> t -> t + val root : int -> t -> t -val merge : t -> t -> t + val add : t -> t -> t -val abs : t -> t + val scale : Numbers.Q.t -> t -> t -val pretty_print : Format.formatter -> t -> unit + val affine_scale : const:Numbers.Q.t -> coef:Numbers.Q.t -> t -> t + (** Perform an affine transformation on the given bounds. + Suposing input bounds (b1, b2), this will return + (const + coef * b1, const + coef * b2). + This function is useful to avoid the incorrect roundings that + can take place when scaling down an integer range. *) -val print : Format.formatter -> t -> unit + val pretty_print : Format.formatter -> t -> unit -val finite_size : t -> Numbers.Q.t option + val print : Format.formatter -> t -> unit -val borne_inf : t -> Numbers.Q.t * Explanation.t * bool -(** bool is true when bound is large. Raise: No_finite_bound if no - finite lower bound *) + val finite_size : t -> Numbers.Q.t option -val borne_sup : t -> Numbers.Q.t * Explanation.t * bool -(** bool is true when bound is large. Raise: No_finite_bound if no - finite upper bound*) + val integer_hull : + t -> + (Numbers.Z.t * Explanation.t) option * + (Numbers.Z.t * Explanation.t) option -val div : t -> t -> t + val borne_inf : t -> Numbers.Q.t * Explanation.t * bool + (** bool is true when bound is large. Raise: No_finite_bound if no + finite lower bound *) -val coerce : Ty.t -> t -> t -(** Coerce an interval to the given type. The main use of that function is - to round a rational interval to an integer interval. This is particularly - useful to avoid roudning too many times when manipulating intervals that - at the end represent an integer interval, but whose intermediate state do - not need to represent integer intervals (e.g. computing the interval for - an integer polynome from the intervals of the monomes). *) + val borne_sup : t -> Numbers.Q.t * Explanation.t * bool + (** bool is true when bound is large. Raise: No_finite_bound if no + finite upper bound*) -val mk_closed : - Numbers.Q.t -> Numbers.Q.t -> bool -> bool -> - Explanation.t -> Explanation.t -> Ty.t -> t -(** - takes as argument in this order: - - a lower bound - - an upper bound - - a bool that says if the lower bound it is large (true) or strict - - a bool that says if the upper bound it is large (true) or strict - - an explanation of the lower bound - - an explanation of the upper bound - - a type Ty.t (Tint or Treal *) + val div : t -> t -> t -type bnd = (Numbers.Q.t * Numbers.Q.t) option * Explanation.t -(* - None <-> Infinity - - the first number is the real bound - - the second number if +1 (resp. -1) for strict lower (resp. upper) bound, - and 0 for large bounds -*) + val coerce : Ty.t -> t -> t + (** Coerce an interval to the given type. The main use of that function is + to round a rational interval to an integer interval. This is particularly + useful to avoid roudning too many times when manipulating intervals that + at the end represent an integer interval, but whose intermediate state do + not need to represent integer intervals (e.g. computing the interval for + an integer polynome from the intervals of the monomes). *) -val bounds_of : t -> (bnd * bnd) list + val contains : t -> Numbers.Q.t -> bool -val contains : t -> Numbers.Q.t -> bool + val add_explanation : t -> Explanation.t -> t -val add_explanation : t -> Explanation.t -> t + val equal : t -> t -> bool -val equal : t -> t -> bool + val pick : is_max:bool -> t -> Numbers.Q.t + (** [pick ~is_max t] returns an elements of the set of intervals [t]. If + [is_max] is [true], we pick the largest element of [t], if it exists. + We look for the smallest element if [is_max] is [false]. *) -val pick : is_max:bool -> t -> Numbers.Q.t option -(** [pick ~is_max t] returns an elements of the set of intervals [t]. If - [is_max] is [true], we pick the largest element of [t], if it exists. - We look for the smallest element if [is_max] is [false]. *) + val fold : + ('a -> Q.t bound interval -> 'a) -> 'a -> t -> 'a -type interval_matching = - ((Numbers.Q.t * bool) option * (Numbers.Q.t * bool) option * Ty.t) - Var.Map.t + type interval_matching = + ((Numbers.Q.t * bool) option * (Numbers.Q.t * bool) option * Ty.t) + Var.Map.t -(** matchs the given lower and upper bounds against the given interval, and - update the given accumulator with the constraints. Returns None if - the matching problem is inconsistent -*) -val match_interval: - Symbols.bound -> Symbols.bound -> t -> interval_matching -> - interval_matching option + (** matchs the given lower and upper bounds against the given interval, and + update the given accumulator with the constraints. Returns None if + the matching problem is inconsistent + *) + val match_interval: + Symbols.bound -> Symbols.bound -> t -> interval_matching -> + interval_matching option +end diff --git a/src/lib/reasoners/intervals_core.ml b/src/lib/reasoners/intervals_core.ml new file mode 100644 index 0000000000..9365044df6 --- /dev/null +++ b/src/lib/reasoners/intervals_core.ml @@ -0,0 +1,807 @@ +(**************************************************************************) +(* *) +(* Alt-Ergo: The SMT Solver For Software Verification *) +(* Copyright (C) 2024 --- OCamlPro SAS *) +(* *) +(* This file is distributed under the terms of OCamlPro *) +(* Non-Commercial Purpose License, version 1. *) +(* *) +(* As an exception, Alt-Ergo Club members at the Gold level can *) +(* use this file under the terms of the Apache Software License *) +(* version 2.0. *) +(* *) +(* --------------------------------------------------------------- *) +(* *) +(* The Alt-Ergo theorem prover *) +(* *) +(* Sylvain Conchon, Evelyne Contejean, Francois Bobot *) +(* Mohamed Iguernelala, Stephane Lescuyer, Alain Mebsout *) +(* *) +(* CNRS - INRIA - Universite Paris Sud *) +(* *) +(* --------------------------------------------------------------- *) +(* *) +(* More details can be found in the directory licenses/ *) +(* *) +(**************************************************************************) + +open Intervals_intf + +(** Pretty-printer for a bound when used as a lower bound. *) +let pp_lower_bound pp_a ppf lb = + match lb with + | Unbounded -> Format.fprintf ppf "]-oo" + | Open x -> Format.fprintf ppf "]%a" pp_a x + | Closed x -> Format.fprintf ppf "[%a" pp_a x + +(** Pretty-printer for a bound when used as an upper bound. *) +let pp_upper_bound pp_a ppf ub = + match ub with + | Unbounded -> Format.fprintf ppf "+oo[" + | Open x -> Format.fprintf ppf "%a[" pp_a x + | Closed x -> Format.fprintf ppf "%a]" pp_a x + +module EqLtLeNotations(OT : OrderedType) = struct + (** This module contains convenient redefinitions of [=], [<], [<=], [min] and + [max] for use with a specific [OrderedType]. *) + + let[@inline always] ( = ) x y = OT.equal x y + + let[@inline always] ( < ) x y = OT.compare x y < 0 + + let[@inline always] ( <= ) x y = OT.compare x y <= 0 + + let[@inline always] min x y = if x <= y then x else y + + let[@inline always] max x y = if y <= x then x else y +end + +module Interval(OT : OrderedType) = struct + type value = OT.finite + + type bnd = OT.t + + type t = OT.t interval + + let pp ppf i = + Format.fprintf ppf "@[%a;@ %a@]" + (pp_lower_bound OT.pp_finite) (OT.view i.lb) + (pp_upper_bound OT.pp_finite) (OT.view i.ub) + + open EqLtLeNotations(OT) + + let equal i1 i2 = + i1.lb = i2.lb && i1.ub = i2.ub + + let of_lower_bound = function + | Unbounded -> OT.minfty + | Closed b -> OT.finite b + | Open b -> OT.succ (OT.finite b) + + let of_upper_bound = function + | Unbounded -> OT.pinfty + | Closed b -> OT.finite b + | Open b -> OT.pred (OT.finite b) + + let of_bounds lb ub = + let lb = of_lower_bound lb in + let ub = of_upper_bound ub in + if ub < lb then invalid_arg "of_bounds" + else { lb ; ub } + + let view { lb ; ub } = { lb = OT.view lb ; ub = OT.view ub } + + let singleton v = of_bounds (Closed v) (Closed v) + + let full = of_bounds Unbounded Unbounded + + let value_opt { lb ; ub } = + if lb = ub then OT.value_opt lb else None +end + +module Make(Ex : Explanations) (*: Core with type explanation = Ex.t *)= struct + (* This module implements the core functionality of the union-of-intervals + representation. + + It is the only module that knows about the details of the internal + representation of the [union] type; all other modules only use the API + exposed through the [Core] signature. + + As such, its correction is critical to the correction of the implementation + as a whole; unfortunately, it is also where most of the complexity lies. + + The implementation is actually contained in the [Union] functor; the [Core] + module itself only contain (polymorphic) type definitions to be able to + convert values between [Union]s of different types. *) + + type explanation = Ex.t + + type 'a union = 'a * ('a interval, Ex.t) kind list * 'a + (* We represent intervals as a list of maximal disjoint allowed intervals with + precise bounds and explanations justifying that the (nonempty) gaps between + each consecutive allowed interval is forbidden. + + This induces a loss of + precision: for instance, if [0, 4] and [10, 20] are allowed and [6, 7] is + forbidden with justification [ex], this will "widen" the forbidden interval + to [4+eps, 10-eps] instead; if we find that the associated value must be + equal to [5], we will include [ex] in the conflict even though it is not + strictly necessary. On the other hand, there are diminishing returns: we + don't want to keep track too precisely of the forbidden intervals if they + are not going to occur in a conflict. + + Note that this is always correct, because it is equivalent to adding + explanations to forbidden intervals (if {m I} is impossible in any model + where {m e} holds, it is also impossible in any model where {b both} {m e} + and {m e'} hold). It also does not change the evaluation in the current + context: we keep track of the bounds of allowed intervals precisely. + + We also store global lower and upper bounds for the union that are always + valid. This is useful for the [trisection_map_to_set] function to + communicate information about the range of values to its arguments. *) + + module Union(OT : OrderedType) = struct + module Interval = Interval(OT) + + open EqLtLeNotations(OT) + + type nonrec 'a union = 'a union + + type bnd = OT.t + + type t = OT.t union + + let pp ppf (glb, u, gub) = + let rec loop ex = function + | [] -> + if Ex.is_empty ex then () + else Fmt.pf ppf "!%a (<= %a)" Ex.pp ex OT.pp gub + | Empty ex' :: u' -> loop (Ex.union ex ex') u' + | NonEmpty i :: u' -> + if not (Ex.is_empty ex) then Fmt.pf ppf "!%a U@ " Ex.pp ex; + Interval.pp ppf i; + if not (Lists.is_empty u') then Fmt.pf ppf " U@ "; + loop ex u' + in + begin match u with + | NonEmpty _ :: _ -> () + | _ -> Fmt.pf ppf "(%a <=)@ " OT.pp glb + end; + loop Ex.empty u + + let pp = Fmt.box pp + + (** {2 Creation} *) + + let of_interval ?(ex = Ex.empty) i = + if Ex.is_empty ex then + (i.lb, [ NonEmpty i ], i.ub) + else + let u = if i.ub = OT.pinfty then [] else [ Empty ex ] in + let u = NonEmpty i :: u in + let u = if i.lb = OT.minfty then u else Empty ex :: u in + (OT.minfty, u, OT.pinfty) + + let[@inline always] of_bounds ?ex lb ub = + of_interval ?ex @@ Interval.of_bounds lb ub + + let empty ex u = + if Ex.is_empty ex then u else Empty ex :: u + + let of_complement ?(ex = Ex.empty) i = + if i.ub = OT.pinfty then + if i.lb = OT.minfty then + Empty ex + else + NonEmpty ( + let u = + NonEmpty { lb = OT.minfty ; ub = OT.pred i.lb } :: empty ex [] + in + (OT.minfty, u, OT.pinfty) + ) + else + NonEmpty ( + let u = empty ex [NonEmpty { lb = OT.succ i.ub ; ub = OT.pinfty }] in + let u = + if i.lb = OT.minfty then u + else NonEmpty { lb = OT.minfty ; ub = OT.pred i.lb } :: u + in + (OT.minfty, u, OT.pinfty) + ) + + (** {2 Inspection} *) + + let lower_bound (_, u, _) = + let rec lower_bound ex = function + | [] -> invalid_arg "lower_bound: empty union" + | NonEmpty i :: _ -> OT.view i.lb, ex + | Empty ex' :: u -> lower_bound (Ex.union ex ex') u + in lower_bound Ex.empty u + + let upper_bound (_, u, _) = + let rec upper_bound ub ex = function + | [] -> OT.view ub, ex + | Empty ex' :: u -> upper_bound ub (Ex.union ex ex') u + | NonEmpty i :: u -> upper_bound i.ub Ex.empty u + in + let rec init = function + | [] -> invalid_arg "upper_bound: empty union" + | Empty _ :: u -> init u + | NonEmpty i :: u -> upper_bound i.ub Ex.empty u + in + init u + + let value_opt (_, u, _) = + let rec loop ex v = function + | [] -> Option.map (fun v -> (v, ex)) v + | Empty ex' :: u -> loop (Ex.union ex ex') v u + | NonEmpty i :: u -> + match v with + | None when i.lb = i.ub -> loop ex (OT.value_opt i.lb) u + | _ -> None + in loop Ex.empty None u + + (** {2 Iteration} *) + + let to_seq (_, u, _) = + List.to_seq u |> Seq.filter_map (function + | NonEmpty i -> Some i + | Empty _ -> None) + + let fold f acc (_, u, _) = + List.fold_left (fun acc -> function + | NonEmpty i -> f acc i + | Empty _ -> acc + ) acc u + + let iter f (_, u, _) = + List.iter (function + | NonEmpty i -> f i + | Empty _ -> () + ) u + + (** {2 Comparison} *) + + (* Note: [subset] and [equal] could be implemented outside of this module as + they do not deal with explanations. There are the exception to the rule + that only functions related to explanations are implemented in [Union], + because it just makes sense to have them here. *) + + let rec subset_seq ~strict i1 s1 i2 s2 = + if i2.ub < i1.lb then + (* i1 is entirely after i2; skip i2 *) + match Seq.uncons s2 with + | None -> false + | Some (i2', s2') -> subset_seq ~strict:false i1 s1 i2' s2' + else + i2.lb <= i1.lb && i1.ub <= i2.ub && + let strict = strict && i1.lb = i2.lb && i1.ub = i2.ub in + match Seq.uncons s1 with + | None -> not strict || Option.is_some (Seq.uncons s2) + | Some (i1', s1') -> + if strict then + match Seq.uncons s2 with + | None -> false + | Some (i2', s2') -> subset_seq ~strict:true i1' s1' i2' s2' + else + subset_seq ~strict:false i1' s1' i2 s2 + + let subset ?(strict = false) u1 u2 = + match Seq.uncons (to_seq u1), Seq.uncons (to_seq u2) with + | None, None -> not strict + | Some _, None -> false + | None, Some _ -> true + | Some (i1, s1), Some (i2, s2) -> + subset_seq ~strict i1 s1 i2 s2 + + let equal u1 u2 = + Seq.equal Interval.equal (to_seq u1) (to_seq u2) + + let checked ((glb, u', gub) as u) = + let rec loop ex = function + | [] -> Empty ex + | Empty ex' :: u' -> loop (Ex.union ex ex') u' + | NonEmpty _ :: _ -> NonEmpty u + in + if gub < glb then Empty Ex.empty else loop Ex.empty u' + + let nonempty = function + | NonEmpty u -> u + | Empty _ -> invalid_arg "nonempty" + + let checked_opt = function + | None -> Empty Ex.empty + | Some u -> checked u + + let add_explanation ~ex ((_, u, _) as u') = + if Ex.is_empty ex then u' else + let[@tail_mod_cons] rec loop ex' = function + | [] -> [ Empty ex' ] + | [ NonEmpty i ] as u when i.ub = OT.pinfty -> Empty ex' :: u + | NonEmpty i :: u -> Empty ex' :: NonEmpty i :: loop ex u + | Empty ex'' :: u -> loop (Ex.union ex' ex'') u + in + let u = + match u with + | NonEmpty i :: u when i.lb = OT.minfty -> + if i.ub = OT.pinfty then [ NonEmpty i ] else NonEmpty i :: loop ex u + | _ -> loop ex u + in + (OT.minfty, u, OT.pinfty) + + let min_ex ex1 ex2 = + if Stdlib.(Ex.compare ex1 ex2 <= 0) then ex1 else ex2 + + (* Computing the intersection follows a fairly natural, if a bit tedious, + process: we are simply iterating on the two unions at the same time and + computing the intersections as we go. + + Some care need to be taken of to make sure to use appropriate + explanations. When one of the unions is enough to forbid a given + forbidden interval, we only want to use that one (for instance, the + intersection of [-oo; 4] with explanation [ex1] and [-oo; 10] with + explanation [ex2] only needs explanation [ex1]). We only need to take the + union of the explanations when the explanation in one union is used to + skip over an interval that was allowed in the other union; for instance, + the intersection of + + [0; 1] U !ex1 U [2; 3] + + and + + !ex2 U [4; 10] + + must be + + !(union ex1 ex2) U [2; 3] + + We need [ex1] to forbid ]1; 2[ and [ex2] to forbid [0; 1] and [2; 4]. + + In addition, the current implementation respects the known global + lower/upper bounds, which adds a low amount of complexity (everything + before/after the global lower/upper bound can be removed without + justification). An alternative implementation could ignore the global + lower and upper bounds and always return [-oo, +oo] as an enclosing + interval; which might improve efficiency of further computations. *) + let intersect = + let[@tail_mod_cons] rec loop ex1 u1 ex2 u2 gub = + match u1, u2 with + | [], [] -> (empty [@tailcall false]) (min_ex ex1 ex2) [] + | Empty ex1' :: u1', _ -> loop (Ex.union ex1 ex1') u1' ex2 u2 gub + | _, Empty ex2' :: u2' -> loop ex1 u1 (Ex.union ex2 ex2') u2' gub + | NonEmpty i1 :: _, _ when gub < i1.lb -> loop ex1 [] ex2 u2 gub + | _, NonEmpty i2 :: _ when gub < i2.lb -> loop ex1 u1 ex2 [] gub + | [], NonEmpty _ :: _ -> (empty [@tailcall false]) ex1 [] + | NonEmpty _ :: _, [] -> (empty [@tailcall false]) ex2 [] + | NonEmpty i1 :: u1', NonEmpty i2 :: u2' -> + let i1, u1' = + if gub < i1.ub then { i1 with ub = gub }, [] else i1, u1' + and i2, u2' = + if gub < i2.ub then { i2 with ub = gub }, [] else i2, u2' + in + if i2.ub < i1.lb then loop ex1 u1 (Ex.union ex1 ex2) u2' gub + else if i1.ub < i2.lb then loop (Ex.union ex1 ex2) u1' ex2 u2 gub + else + let ub = min i1.ub i2.ub in + let lb, ex = + if i1.lb < i2.lb then i2.lb, ex2 + else if i2.lb < i1.lb then i1.lb, ex1 + else i1.lb, min_ex ex1 ex2 + in + if Ex.is_empty ex then cons_loop { lb ; ub } i1 u1' i2 u2' gub + else Empty ex :: cons_loop { lb ; ub } i1 u1' i2 u2' gub + and[@tail_mod_cons] cons_loop i i1 u1 i2 u2 gub = + NonEmpty i :: + let lb = OT.succ i.ub in + let u1' = if i.ub < i1.ub then NonEmpty { i1 with lb } :: u1 else u1 in + let u2' = if i.ub < i2.ub then NonEmpty { i2 with lb } :: u2 else u2 in + loop Ex.empty u1' Ex.empty u2' gub + in + fun (glb1, u1, gub1) (glb2, u2, gub2) -> + let glb = max glb1 glb2 and gub = min gub1 gub2 in + if gub < glb then Empty Ex.empty else + let rec gist_ge glb ex = function + | [] -> empty ex [] + | Empty ex' :: u' -> gist_ge glb (Ex.union ex ex') u' + | NonEmpty i :: u' when i.ub < glb -> gist_ge glb Ex.empty u' + | NonEmpty i :: u' as u -> + if glb < i.lb then empty ex u + else NonEmpty { i with lb = glb } :: u' + in + let u1 = gist_ge glb Ex.empty u1 and u2 = gist_ge glb Ex.empty u2 in + checked (glb, loop Ex.empty u1 Ex.empty u2 gub, gub) + + (** {2 Image by monotone functions} *) + + (* Note that [partial_map_inc] and [partial_map_dec] are currently + implemented in the next section using [approx_map_inc_to_set] and + [approx_map_dec_to_set] for simplicity. *) + + let map_strict_inc f (glb, u, gub) = + (* Use a [ref] here instead of returning a pair to allow [@tail_mod_cons] + annotation in [loop]. *) + let f_gub = ref None in + let[@tail_mod_cons] rec loop = function + | [] -> [] + | Empty ex :: u' -> Empty ex :: loop u' + | NonEmpty i :: u' -> + let f_i = { lb = f i.lb ; ub = f i.ub } in + if Lists.is_empty u' then f_gub := Some f_i.ub; + NonEmpty f_i :: loop u' + in + let f_u = loop u in + let f_glb = match f_u with NonEmpty i :: _ -> i.lb | _ -> f glb in + let f_gub = match !f_gub with Some gub -> gub | None -> f gub in + (f_glb, f_u, f_gub) + + let map_strict_dec f (glb, u, gub) = + let rec loop f_glb acc = function + | [] -> (f gub, acc, f_glb) + | Empty ex' :: u' -> loop f_glb (Empty ex' :: acc) u' + | NonEmpty i :: u' -> + let f_i = { lb = f i.ub ; ub = f i.lb } in + loop f_glb (NonEmpty f_i :: acc) u' + in + let init = function + | NonEmpty i :: u' -> + let f_i = { lb = f i.ub ; ub = f i.lb } in + loop f_i.ub [NonEmpty f_i] u' + | u -> loop (f glb) [] u + in init u + + (** {2 Image by arbitrary functions} *) + + type 'a basic = + | Allow of OT.t interval + | Forbid of OT.t interval * Ex.t + + (* The ['a set] type represents a loose union of [basic] intervals. There + are no restrictions on the underlying [basic] intervals (although they + should be valid, as per the {!Intervals} signature). + + Because we mostly use the [set] type as an intermediate accumulator when + building unions, we expose it as functions from list of [basic] intervals + to [basic] intervals. The underlying list can be recovered by calling the + [set] with an empty list as argument, and we will loosely identify the + [set] and its representation as a list. *) + + type set = OT.t basic list -> OT.t basic list + + let empty_set = fun acc -> acc + + let interval_set i acc = Allow i :: acc + + let union_set s1 s2 acc = s1 (s2 acc) + + let map_set f s acc = + (* Make sure to only call [f] on the intervals that come from [s]! *) + List.rev_append ( + List.rev_map (function + | Allow i -> Allow (f i) + | Forbid (i, ex) -> Forbid (f i, ex) + ) (s []) + ) acc + + let of_set = + (* [of_set] is used to convert an unnormalized [set] to a normalized + [union]. Since a [set] can be empty, but an [union] cannot, we return + an [option] here. + + We maintain a list of the [basic] sets, sorted by decreasing upper + bound (it would be more natural to sort by increasing lower bound, but + sorting by decreasing upper bound allows us to build the resulting + [union] backwards in a tail-recursive way). *) + + (* [prepend_forbid gub lb u] is used to add all the forbidden intervals + from [fs] to the initial segment before [u]. + + [u] must begin with an allowed segment that starts at [lb]; any + forbidden intervals that start at [lb] or after can be skipped. *) + let rec prepend_forbid gub lb u = function + | [] -> (lb, u, gub) + | (flb, _) :: fs when lb <= flb -> + prepend_forbid gub lb u fs + | (glb, ex) :: fs -> + let glb, ex = + List.fold_left (fun (glb, ex) (lb', ex') -> + if lb <= lb' then (glb, ex) else (min glb lb', Ex.union ex ex') + ) (glb, ex) fs + in + (glb, Empty ex :: u, gub) + in + (* Main loop. + + We assume the input list to be sorted by decreasing upper bound, and + the upper bound of intervals in the input list are currently [<= ub]. + + We have built the union [u], and are processing an allowed interval + { lb ; ub } to be added at the start of [u]. We don't add it yet, + because it may need to be merged with other allowed intervals we have + not processed yet: we have only processed intervals with an upper bound + greater than or equal to [ub], but there may be allowed intervals that + end inside [lb; ub] (or at [pred lb]). + + The list [fs] collects a set of forbidden half-intervals that must be + added before [lb; ub] (unless they get subsumed by merging). *) + let rec loop gub lb ub u fs = function + | [] -> + prepend_forbid gub lb (NonEmpty { lb ; ub } :: u) fs + | Forbid (i, _) :: s when lb <= i.lb -> + (* [i] is subsumed by [lb; ub] *) + loop gub lb ub u fs s + | Forbid (i, ex) :: s -> + (* [i] starts before [lb]; need to keep it *) + loop gub lb ub u ((i.lb, ex) :: fs) s + | Allow i :: s when lb <= i.ub || OT.succ i.ub = lb -> + (* No gap: merge intervals ([i.ub] is smaller by sorting)*) + loop gub (min i.lb lb) ub u fs s + | Allow i :: s -> + (* Gap with another allowed interval *) + assert (OT.succ i.ub < lb); + (* Need to filter the [flb] again because the forbidden intervals + might have been added before the merging that subsumed it. *) + let ex = + List.fold_left (fun ex (flb, ex') -> + if flb < lb then Ex.union ex ex' else ex + ) Ex.empty fs + in + let forbid = List.filter (fun (flb, _) -> flb < i.lb) fs in + let u = empty ex @@ NonEmpty { lb ; ub } :: u in + loop gub i.lb i.ub u forbid s + in + (* Look for the first allowed interval (which has the highest upper bound, + since the input list is sorted by decreasing order). + + Note that we assume that ties are broken by sorting allowed intervals + before forbidden intervals. *) + let rec init lb ex gub fs = function + | [] -> (lb, [ Empty ex ], gub) + | Allow i :: s -> + assert (i.ub < gub); + let fs = List.filter (fun (lb, _) -> lb < i.lb) fs in + loop gub i.lb i.ub [Empty ex] fs s + | Forbid (i, ex') :: s -> + (* NB: all allowed intervals have a strictly smaller upper bound due + to our tie-breaking. *) + init (min lb i.lb) (Ex.union ex ex') gub ((i.lb, ex') :: fs) s + in + fun s -> + let s = + List.filter + (function Allow _ -> true | Forbid (_, ex) -> not (Ex.is_empty ex)) + (s []) + in + (* Sort by decreasing upper bound, ties going to allowed intervals *) + let rev_cmp b1 b2 = + let (Allow i1 | Forbid (i1, _)) = b1 in + let (Allow i2 | Forbid (i2, _)) = b2 in + let c = OT.compare i2.ub i1.ub in + if c <> 0 then c else + match b1, b2 with + | Allow _, Forbid _ -> -1 + | Forbid _, Allow _ -> 1 + | Allow _, Allow _ | Forbid _, Forbid _ -> 0 + in + match List.fast_sort rev_cmp s with + | [] -> None + | Allow i :: s -> + Some (loop i.ub i.lb i.ub [] [] s) + | Forbid (i, ex) :: s -> + (* NB: all allowed intervals have a strictly smaller upper bound due + to our tie-breaking. *) + Some (init i.lb ex i.ub [(i.lb, ex)] s) + + let of_set_checked s = checked_opt (of_set s) + + let of_set_nonempty s = nonempty @@ checked_opt @@ of_set s + + let forbid_set ~ex s acc = + if Ex.is_empty ex then acc + else + List.rev_append + (List.rev_map + (function Allow i | Forbid (i, _) -> Forbid (i, ex)) (s [])) + acc + + (* Note that for forbidden intervals we call [f] with a slightly + larger interval that includes the allowed parts. *) + let map_to_set f (glb, u, gub) acc = + let rec loop acc lb ex = function + | [] -> forbid_set ~ex (f { lb ; ub = gub }) acc + | Empty ex' :: u -> loop acc lb (Ex.union ex ex') u + | NonEmpty i :: u -> + let acc = + if Ex.is_empty ex then acc + else forbid_set ~ex (f { i with lb }) acc + in + loop (f i acc) i.lb Ex.empty u + in loop acc glb Ex.empty u + + let union_hull i1 i2 = + { lb = min i1.lb i2.lb ; ub = max i1.ub i2.ub } + + let hull_set = function + | [] -> None + | (Allow i | Forbid (i, _)) :: s -> + Some ( + List.fold_left (fun { lb = glb ; ub = gub } b -> + let { lb ; ub } = match b with Allow i | Forbid (i, _) -> i in + { lb = min lb glb ; ub = max ub gub } + ) i s + ) + + (* Similar to [map_to_set] except that we are keeping track of the + global hull of each call to [f i] where [i] is allowed. + + Then we know that if [x] is between two consecutive intervals + [i] and [i'] (with [i.ub < i'.lb]), [f x] is included in the + hull of [union (f i) (f i')]. + + Note that we "overflow" from the strict gap between [i] and + [i'] but that is OK: the image of [f i] and [f i'] are allowed, + which will subsume the corresponding forbidden part from the + hulls when converting back to an union. *) + let map_mon_to_set f (glb, u, gub) acc = + let rec loop acc prev_hull lb ex = function + | [] -> forbid_set ~ex (f { lb ; ub = gub }) acc + | Empty ex' :: u -> loop acc prev_hull lb (Ex.union ex ex') u + | NonEmpty i :: u -> + let f_i = f i [] in + match hull_set f_i with + | None -> loop acc prev_hull lb ex u + | Some hull_f -> + let acc = + if Ex.is_empty ex then acc + else Forbid (union_hull prev_hull hull_f, ex) :: acc + in + let acc = List.rev_append f_i acc in + loop acc hull_f i.lb Ex.empty u + in + let rec init ex = function + | [] -> forbid_set ~ex (f { lb = glb ; ub = gub }) acc + | Empty ex' :: u -> init (Ex.union ex ex') u + | NonEmpty i :: u -> + let f_i = f i [] in + match hull_set f_i with + | None -> init ex u + | Some hull_f -> + let acc = + if Ex.is_empty ex then acc + else forbid_set ~ex (f { i with lb = glb }) acc + in + let acc = List.rev_append f_i acc in + loop acc hull_f i.lb Ex.empty u + in init Ex.empty u + + let approx_map_inc_to_set f_lb f_ub u acc = + map_mon_to_set (fun i -> + let lb = f_lb i.lb and ub = f_ub i.ub in + if OT.compare lb ub > 0 then empty_set + else interval_set { lb ; ub } + ) u acc + + let approx_map_dec_to_set f_lb f_ub u acc = + map_mon_to_set (fun i -> + let ub = f_ub i.lb and lb = f_lb i.ub in + if OT.compare lb ub > 0 then empty_set + else interval_set { lb ; ub } + ) u acc + + let map_inc_to_set f u = approx_map_inc_to_set f f u + + let map_dec_to_set f u = approx_map_dec_to_set f f u + + (* We could generalize [map_strict_inc] and [map_strict_dec] to implement + [partial_map_inc] and [partial_map_dec] but that would make the code more + complex due to needing to potentially merge consecutive intervals. + Instead, we centralize such logic in [of_set] (but still provide the + [partial_map_inc] and [partial_map_dec] interfaces to be able to switch + to a specialized implementation later). *) + + let partial_map_inc f_lb f_ub u = + of_set_checked (approx_map_inc_to_set f_lb f_ub u) + + let partial_map_dec f_lb f_ub u = + of_set_checked (approx_map_dec_to_set f_lb f_ub u) + + (* [uncons] and [cons] functions provide a convenient "view" of an [union]. + They allow uniform reasoning between allowed and forbidden intervals in + [trisection_map_to_set]. *) + + (** [uncons u] returns a triple [i, ex, u'] such that [i] is the first + interval of [u] and [u'] is the part of [u] strictly after [i.ub] if + any. + + The explanation [ex] is [None] if [i] is allowed, and [Some] with the + corresponding explanation if [i] is forbidden. *) + let uncons (glb, u, gub) = + let rec loop ex = function + | [] -> { lb = glb ; ub = gub }, Some ex, None + | Empty ex' :: u' -> loop (Ex.union ex ex') u' + | NonEmpty i :: u' as u -> + if Ex.is_empty ex then + match u' with + | [] -> (i, None, None) + | _ -> (i, None, Some (OT.succ i.ub, u', gub)) + else ({ lb = glb ; ub = OT.pred i.lb }, Some ex, Some (i.lb, u, gub)) + in loop Ex.empty u + + (** [cons i ex u] prepends the interval [i] (as an allowed interval if [ex] + is [None], and a forbidden interval if it is [Some]) to the union [u]. + + [u] must start after i, i.e. we must have [i.ub < (hull u).lb]. + + [cons] is an inverse of [uncons], but [uncons] is not an inverse of + [cons] because [cons] can merge consecutive intervals. *) + + let cons i ex u = + match u with + | None -> ( + match ex with + | None-> (i.lb, [ NonEmpty i ], i.ub) + | Some ex-> + if Ex.is_empty ex then (i.lb, [], i.ub) + else (i.lb, [ Empty ex ], i.ub) + ) + | Some (glb, u, gub) -> + assert (i.ub < glb); + match ex, u with + | None, NonEmpty i' :: u' when OT.succ i.ub = i'.lb -> + (i.lb, NonEmpty { lb = i.lb ; ub = i'.ub } :: u', gub) + | None, _ -> (i.lb, NonEmpty i :: u, gub) + | Some ex, Empty ex' :: u' -> + (i.lb, Empty (Ex.union ex ex') :: u', gub) + | Some ex, _ -> (i.lb, empty ex u, gub) + + + let trisection_map_to_set v ((_, _, gub) as u) f_lt f_eq f_gt acc = + if gub < v then + (* Optimization to avoid traversing the whole union *) + f_lt u acc + else + let rec trisection acc u = + let i, ex, u' = uncons u in + if i.ub < v then + (* The interval is strictly before the value *) + match u' with + | None -> + Some (cons i ex None), acc + | Some u' -> + let before, acc = trisection acc u' in + Some (cons i ex before), acc + else + (* v <= i.ub *) + let before = + (* [before] is the portion of [i] strictly before [v], if any *) + if v <= i.lb then None + else Some (cons { i with ub = OT.pred v } ex None) + in + let acc = + (* if [v] is inside [i], add [f_eq] to the mix with appropriate + explanation. *) + if v < i.lb || Option.fold ~none:false ~some:Ex.is_empty ex then + acc + else + match ex with + | None -> f_eq () acc + | Some ex -> forbid_set ~ex (f_eq ()) acc + in + let after = + (* [after] is the portion of [u] strictly after [v], if any *) + if v = i.ub then u' + else if v < i.lb then Some (cons i ex u') + else Some (cons { i with lb = OT.succ v } ex u') + in + let acc = + match after with + | Some after -> f_gt after acc + | None -> acc + in + before, acc + in + let before, acc = trisection acc u in + match before with + | Some before -> f_lt before acc + | None -> acc + end +end diff --git a/src/lib/reasoners/intervals_core.mli b/src/lib/reasoners/intervals_core.mli new file mode 100644 index 0000000000..0297d91490 --- /dev/null +++ b/src/lib/reasoners/intervals_core.mli @@ -0,0 +1,33 @@ +(**************************************************************************) +(* *) +(* Alt-Ergo: The SMT Solver For Software Verification *) +(* Copyright (C) 2024 --- OCamlPro SAS *) +(* *) +(* This file is distributed under the terms of OCamlPro *) +(* Non-Commercial Purpose License, version 1. *) +(* *) +(* As an exception, Alt-Ergo Club members at the Gold level can *) +(* use this file under the terms of the Apache Software License *) +(* version 2.0. *) +(* *) +(* --------------------------------------------------------------- *) +(* *) +(* The Alt-Ergo theorem prover *) +(* *) +(* Sylvain Conchon, Evelyne Contejean, Francois Bobot *) +(* Mohamed Iguernelala, Stephane Lescuyer, Alain Mebsout *) +(* *) +(* CNRS - INRIA - Universite Paris Sud *) +(* *) +(* --------------------------------------------------------------- *) +(* *) +(* More details can be found in the directory licenses/ *) +(* *) +(**************************************************************************) + +open Intervals_intf + +(** This module implements union of intervals with explanations. See the + {!Intervals_intf.Core} signature. *) + +module Make(Ex : Explanations) : Core with type explanation = Ex.t diff --git a/src/lib/reasoners/intervals_intf.ml b/src/lib/reasoners/intervals_intf.ml new file mode 100644 index 0000000000..f29da80648 --- /dev/null +++ b/src/lib/reasoners/intervals_intf.ml @@ -0,0 +1,735 @@ +(**************************************************************************) +(* *) +(* Alt-Ergo: The SMT Solver For Software Verification *) +(* Copyright (C) 2024 --- OCamlPro SAS *) +(* *) +(* This file is distributed under the terms of OCamlPro *) +(* Non-Commercial Purpose License, version 1. *) +(* *) +(* As an exception, Alt-Ergo Club members at the Gold level can *) +(* use this file under the terms of the Apache Software License *) +(* version 2.0. *) +(* *) +(* --------------------------------------------------------------- *) +(* *) +(* The Alt-Ergo theorem prover *) +(* *) +(* Sylvain Conchon, Evelyne Contejean, Francois Bobot *) +(* Mohamed Iguernelala, Stephane Lescuyer, Alain Mebsout *) +(* *) +(* CNRS - INRIA - Universite Paris Sud *) +(* *) +(* --------------------------------------------------------------- *) +(* *) +(* More details can be found in the directory licenses/ *) +(* *) +(**************************************************************************) + +(** Interface for union-of-interval modules. *) + +type 'a interval = { lb : 'a ; ub : 'a } +(** The type of {b closed} intervals over type ['a]. + + An interval [{ lb ; ub }] represents the closed interval {m [lb, ub]}, i.e. + a value [v] is in the interval iff [lb <= v <= ub]. *) + +type 'a bound = + | Open of 'a (** An open (strict) finite bound. *) + | Closed of 'a (** A closed (large) finite bound. *) + | Unbounded (** An infinite bound. *) +(** The type ['a bound] is used to create and inspect intervals; see + {!OrderedType.view} and {!Interval.of_bounds}. *) + +(** Type signature for explanations. This is mostly intended for Alt-Ergo's + built-in {!Explanation} module, but it can be convenient to use a + different module for debugging. *) +module type Explanations = sig + type t + (** The type of explanations. + + Explanation encode predicates that may hold (be [true]) or not (be + [false]) in a model. *) + + val pp : t Fmt.t + (** Pretty-printer for explanations. *) + + val is_empty : t -> bool + (** An explanation is empty if it is [true] in all models. *) + + val empty : t + (** The [empty] explanation. Must satisfy [is_empty empty = true]. *) + + val union : t -> t -> t + (** The union [union ex ex'] of two explanations [ex] and [ex'] is [true] in + any model where both [ex] and [ex'] are [true]. + + Note that [union ex empty = union empty ex = ex]. *) + + val compare : t -> t -> int + (** Arbitrary comparison function on explanations. In case multiple + explanations for a given fact are possible, the smaller explanation + according to this comparison function is preferred. + + It is recommended, but not required, that {!empty} be the lowest + explanation for [compare]. *) +end + +module type OrderedType = sig + (** An ordered type of finite values, extended with positive and negative + infinites as well as successor and predecessor functions to represent + strict lower and upper bounds. + + Negative and positive infinites are used to represent half-planes as + intervals; we require that all values are larger than the negative + infinite and smaller than the positive infinite. + + Values of the form {m x + \epsilon} and {m x - \epsilon} are used to + transform open bounds in the {!type-finite} type into closed bounds on the + extended type {!t}, with the equivalence: + + {math + y <= x - \epsilon \Leftrightarrow y < x + } + + and + + {math + x + \epsilon <= y \Leftrightarrow y < x + } + + This signature does not expose {m x + \epsilon} and {m x - \epsilon} + values directly; instead, we use {m \mathrm{succ}} and {m \mathrm{pred}} + functions which can be thought of as addition and subtraction of + {m \epsilon}, respectively. *) + + type finite + (** The type of finite values. *) + + val pp_finite : finite Fmt.t + (** Pretty-printer for finite values. *) + + type t + (** The type of extended values used to represent interval bounds. *) + + val pp : t Fmt.t + (** Pretty-printer for extended values. *) + + val equal : t -> t -> bool + (** Equality on extended values. Must be compatible with [compare]. *) + + val compare : t -> t -> int + (** The comparison function on extended values is an extension of the regular + order on finite values. *) + + val view : t -> finite bound + (** [view t] converts the internal representation to a [bound] for + examination. *) + + val minfty : t + (** [minfty] is an extended value {m -\infty} such that for any extended value + {m x}, {m -\infty \le x}. *) + + val finite : finite -> t + (** Finite values are included in [t]. We will identify a finite value in + [finite] and its representation in [t]. *) + + val pinfty : t + (** [pinfty] is an extended value {m +\infty} such that for any extended value + {m x}, {m x \le +\infty}. *) + + val value_opt : t -> finite option + (** [value_opt] is the partial inverse of [finite]. *) + + val succ : t -> t + (** Each element of the ordered type [t] has a successor. The successor + of an element is always greater than the element itself: + + {math + \forall x, x \le \mathrm{succ}(x) + } + + We say that {m x} is a {e finite upper bound} if it is {e strictly} + smaller than its successor, i.e. if {m x < \mathrm{succ}(x)}, and we + require that all {b finite} values are finite upper bounds (however there + may be finite upper bounds that are not finite values). + + [succ] must be the inverse of [pred] below. *) + + val pred : t -> t + (** Each element of the ordered type [t] has a predecessor. The predecessor + of an element is always smaller than the element itself: + + {math + \forall x, \mathrm{pred}(x) \le x + } + + We say that {m x} is a {e finite lower bound} if it is {e strictly} + greater than its predecessor, i.e. if {m \mathrm{pred}(x) < x}, and we + require that all {b finite} values are finite lower bounds (however there + may be finite lower bounds that are not finite values). + + [pred] must be the inverse of [succ] above. *) +end + +(** Intervals over an {!OrderedType}. *) +module type Interval = sig + type value + (** The type of values in the interval. *) + + type bnd + (** The type of the interval bounds. *) + + type t = bnd interval + (** The type of {b closed} intervals with bounds of type [bound]. + + Note that the {b values} of the interval may have a different type; for + instance, {m -\infty} and {m +\infty} are valid bounds for real intervals, + but are not themselves reals. + + We say that an extended value {m x} is a {b valid lower bound} if it + satisfies: + + {math + \forall y, y < x \Rightarrow \mathrm{pred}(x) < x + } + + Similarly, we say that an extended value {m x} is a {b valid upper bound} + if it satisfies: + + {math + \forall y, x < y \Rightarrow x < \mathrm{succ}(x) + } + + Remark that valid lower bounds are either {m -\infty} or finite lower + bounds, and valid upper bounds are either {m +\infty} or finite upper + bounds. + + We require that an interval [{ lb ; ub }] be such that [lb] is a valid + lower bound, [ub] is a valid upper bound, and [lb <= ub]. *) + + val pp : t Fmt.t + (** Pretty-printer for intervals. *) + + val equal : t -> t -> bool + (** Equality test for intervals. *) + + val of_bounds : value bound -> value bound -> t + (** Build an interval from a pair of lower and upper bounds. + + @raises Invalid_argument if the upper bound is smaller than the lower + bound. *) + + val view : t -> value bound interval + (** Returns a view of the interval using the [bound] type for convenient + examination. *) + + val singleton : value -> t + (** The interval [singleton v] contains the singleton {m \{ v \}}. + + This is an equivalent shortcut for [of_bounds (Closed v) (Closed v)]. *) + + val full : t + (** The full interval {m [-\infty, +\infty]} contains all values. + + This is an equivalent shortcut for [of_bounds Unbounded Unbounded]. *) + + val value_opt : t -> value option + (** If the interval is a singleton containing a single value, returns that + value; otherwise, returns [None]. *) +end + +type ('a, 'b) kind = + | NonEmpty of 'a (** A {b non-empty} value of type ['a]. *) + | Empty of 'b (** An empty value with additional justification. *) +(** The [kind] type is an equivalent to the [result] type from the standard + library with more appropriate constructor names. + + It is used to distinguish between unions of intervals that are empty in the + current context (with an explanation abstracting a set of model where the + union is also empty) and unions of intervals that have at least one value in + the current context. *) + +module type Union = sig + (** Signature for union-of-intervals implementations {b with explanations}. *) + + (** Given an explanation {m e} (see {!Intervals_intf.Explanations}), we say + that an interval {m I} is forbidden for a variable {m x} by {m e} if + {m M(x)} is never in {m I} when {m e} is [true], i.e. if the following + implication holds: + + {math \forall M, M(e) \Rightarrow M(x) \not\in I} + + Given a current context {m C} (a set of explanations) and an implicit + variable {m x}, we say that an interval {m I} is {e forbidden} if it is + forbidden by an union of explanations in {m C}, and that it is {e allowed} + otherwise. + + The [union] type below internally keeps track of both allowed and + forbidden intervals, where each forbidden interval is annotated by a + specific explanation (true in the current context) that forbids it, but + the functions in this module provide abstractions to deal with the unions + without having to worry about explanations (and forbidden intervals) + directly. + + In order to avoid having to worry about variables while reasoning about + intervals, we define a contextual evaluation function for unions (and set + of intervals); the evaluation of an allowed interval {m I} is always {m I} + in all contexts, and the evaluation of a forbidden interval {m I} is {m I} + in contexts where the associated explanation {b does not hold}, and the + empty set otherwise (in particular, it is the empty set in the current + context). + + Note that the evaluation of a forbidden interval is a subset of the + interval in any context. *) + + type explanation + (** The type of explanations. *) + + type value + (** The type of values. *) + + type bnd + (** The type of bounds. *) + + type 'a union + + type t = bnd union + (** The type of unions. *) + + module Interval : Interval + with type value = value + and type bnd = bnd + (** Intervals over the [value] type. *) + + val pp : t Fmt.t + (** Pretty-printer for unions. *) + + (** {1 Creation} + + The following functions are used to create unions. *) + + val of_interval : + ?ex:explanation -> bnd interval -> t + (** Build an union from an interval. + + [of_interval ?ex:None i] evaluates to [i] in all contexts. + + [of_interval ~ex i] evaluates to [i] in contexts where [ex] holds + (including the current context) and to the full set otherwise. *) + + val of_bounds : + ?ex:explanation -> value bound -> value bound -> t + (** [of_bounds ?ex lb ub] is a shortcut for + [of_interval ?ex @@ Interval.of_bounds lb ub] *) + + val of_complement : + ?ex:explanation -> bnd interval -> (t, explanation) kind + (** Build an union from the complement of an interval. + + The explanation, if provided, justifies that the value is {b not} in the + interval; the resulting union evaluates to the complement of the interval + in any model where the explanation holds, and to the full set otherwise. + + If the provided interval covers the full domain, the resulting union is + empty and an [Empty] kind is returned; otherwise, the resulting union is + [NonEmpty]. *) + + (** {1 Inspection} + + The following functions are used to inspect the global bounds of an + union of intervals. *) + + val lower_bound : t -> value bound * explanation + (** [lower_bound u] returns a pair [lb, ex] of a global lower bound and an + explanation that justifies that this lower bound is valid. *) + + val upper_bound : t -> value bound * explanation + (** [upper_bound u] returns a pair [lb, ex] of a global upper bound and an + explanation that justifies that this upper bound is valid. *) + + val value_opt : t -> (value * explanation) option + (** If [u] is a singleton {b in the current context}, [value_opt u] returns a + pair [(v, ex)] where [v] is the value of [u] in the current context and + [ex] justifies [v] being both a lower and upper bound. + + This is more efficient than checking if the values returned by + [lower_bound] and [upper_bound] are equal. *) + + (** {1 Iteration} + + The following functions are used to iterate over the maximal disjoint + intervals composing the union {b in the current context}. *) + + val to_seq : t -> bnd interval Seq.t + (** Convert an union of interval to a sequence of maximal disjoint and + non-adjacent intervals. + + {b Warning}: The sequence of intervals is only valid in the current + context. *) + + val iter : (bnd interval -> unit) -> t -> unit + (** [iter f u] calls [f] on each of the maximal disjoint and non-adjacent + intervals that make up the union [u] {b in the current context}. *) + + val fold : ('a -> bnd interval -> 'a) -> 'a -> t -> 'a + (** [fold f acc u] folds [f] over each of the maximal disjoint and + non-adjacent intervals that make up the union [u] {b in the current + context}. *) + + (** {1 Comparison} + + The following functions are used to compare intervals {b in the current + context}. *) + + val subset : ?strict:bool -> t -> t -> bool + (** [subset ?strict u1 u2] returns [true] if [u1] is a subset of [u2] {b in + the current context}. [subset] ignores all explanations. + + If set, the [strict] flag ([false] by default) checks for strict + inclusion. *) + + val equal : t -> t -> bool + (** [equal u1 u2] returns [true] if [u1] is equal to [u2] {b in the current + context}. [equal] ignores all explanations. *) + + (** {1 Set manipulation} *) + + val add_explanation : ex:explanation -> t -> t + (** [add_explanation ~ex u] adds the explanation [ex] to the union [u]. + + More precisely: [add_explanation ~ex u] is equivalent to [u] in models + where [ex] holds, and is [Interval.full] otherwise. *) + + val intersect : t -> t -> (t, explanation) kind + (** [intersect u1 u2] computes the intersection of [u1] and [u2]. + + The resulting intersection is valid in all models. *) + + (** {1 Image by monotone functions} + + The following functions can compute the image of an union by monotone + function on bounds. *) + + val map_strict_inc : ('a -> bnd) -> 'a union -> t + (** [map_strict_inc f u] computes the image of [u] by the {b strictly + increasing} function [f]. + + This function is very efficient and should be used when possible. *) + + val map_strict_dec : ('a -> bnd) -> 'a union -> t + (** [map_strict_dec f u] computes the image of [u] by the {b strictly + decreasing} function [f]. + + This function is very efficient and should be used when possible. *) + + val partial_map_inc : + ('a -> bnd) -> ('a -> bnd) -> 'a union -> (t, explanation) kind + (** [partial_map_inc f_lb f_ub u] computes the image of [u] by a partial + (weakly) increasing function [f]. + + [f] is represented by the pair [(f_lb, f_ub)] of functions that are such + that for any [x], [f_lb x = f x = f_ub x] if [x] is in the domain of + [f], and [f_ub x < f_lb x] otherwise. + + {b Warning}: The functions [f_lb] and [f_ub] must themselve be (weakly) + increasing. Moreover, the inequality [f_ub x <= f_lb x] must hold + everywhere; if [f_lb] and [f_ub] are imprecise approximations of [f], use + {!approx_map_inc_to_set} instead. *) + + val partial_map_dec : + ('a -> bnd) -> ('a -> bnd) -> 'a union -> (t, explanation) kind + (** [partial_map_dec f_lb f_ub u] computes the image of [u] by a partial + (weakly) decreasing function [f]. + + [f] is represented by the pair [(f_lb, f_ub)] of functions that are such + that for any [x], [f_lb x = f x = f_ub x] if [x] is in the domain of + [f], and [f_ub x < f_lb x] otherwise. + + {b Warning}: The functions [f_lb] and [f_ub] must themselve be (weakly) + decreasing. Moreover, the inequality [f_ub x <= f_lb x] must hold + everywhere; if [f_lb] and [f_ub] are imprecise approximations of [f], use + {!approx_map_dec_to_set} instead. *) + + (** {1 Image by arbitrary functions} + + When computing the image of an union by non-monotone functions and + multi-variate functions, work needs to be done to convert intermediate + results into a normalized form of disjoint intervals (because applying the + function to initially disjoint intervals can create arbitrary unions). + + In such cases, the following family of [to_set] functions should be used. + The use the [set] type to represent intermediate results, which can be + converted back to an union using the [of_set] family of functions as + needed. *) + + type set + (** The [set] type represents a (possibly empty) set of intervals. + + Unlike an [union], it is not normalized: it is an intermediate type not + intended to be manipulated directly, except in the process of building an + union. Use the [of_set] family of functions to convert it to an union as + needed. *) + + val empty_set : set + (** [empty_set] represents an empty set of intervals. *) + + val interval_set : bnd interval -> set + (** [interval_set i] represents the interval [i] as a [set]. *) + + val union_set : set -> set -> set + (** [union_set f1 f2] computes the union of the sets of intervals represented + by [f1] and [f2]. *) + + val map_set : (bnd interval -> bnd interval) -> set -> set + (** [map_set f s] applies the function [f] to all the intervals in [s]. *) + + val of_set_checked : set -> (t, explanation) kind + (** Converts a [set] to a nonempty union. Returns [Empty ex] if the set is + currently empty, and [NonEmpty u] otherwise. + + If a [NonEmpty] value is returned, the resulting union is guaranteed to be + nonempty in the current context. *) + + val of_set_nonempty : set -> t + (** Converts a [set] to an union, assuming that it is nonempty in the current + context. + + @raise Invalid_argument if the [set] represents an union that is empty in + the current context. *) + + val map_to_set : ('a interval -> set) -> 'a union -> set + (** [map_to_set f u] computes the image of [u] by [f]. + + There are no restrictions on [f] (except that it be isotone, i.e. it must + be compatible with inclusion), which means that [map_to_set] needs + to call [f] on currently allowed intervals, but also on some intervals + that are currently impossible but would be possible in other contexts + (depending on explanations). + + When possible, prefer using a more specialized variant of + [map_to_set] that use properties of the function [f] to avoid + certain calls to [f]. *) + + val map_mon_to_set : ('a interval -> set) -> 'a union -> set + (** [map_mon_to_set] is a variant of [map_to_set] when the + function [f] is monotone. + + More precisely, we require that for any pair of intervals [(i1, i2)] the + interval hull of [f i1] and [f i2] is included in the image of the + interval hull of [i1] and [i2] by [f]. + + It is often more convenient to use [approx_map_inc_to_set] or + [approx_map_dec_to_set] instead, to lift a function on values to a + function on unions without going through intervals. *) + + val approx_map_inc_to_set : ('a -> bnd) -> ('a -> bnd) -> 'a union -> set + (** [approx_map_inc_to_set f_lb f_ub u] is a variant of [map_to_set] that is + appropriate for (weakly) increasing functions. + + The (weakly) increasing function [f] is represented by a pair + [(f_lb, f_ub)] of functions such that for all [x] in the domain of [f], + [f_lb x <= f x <= f_ub x] (and [f_ub x < f_lb x] otherwise). + + Note that, unlike [partial_map_inc], this function can be used with + imprecise approximations: we allow [f_lb x < f x < f_ub x] (this is useful + for radicals). *) + + val approx_map_dec_to_set : ('a -> bnd) -> ('a -> bnd) -> 'a union -> set + (** [approx_map_inc_to_set f_lb f_ub u] is a variant of [map_to_set] that is + appropriate for (weakly) decreasing functions. + + The (weakly) decreasing function [f] is represented by a pair + [(f_lb, f_ub)] of functions such that for all [x] in the domain of [f], + [f_lb x <= f x <= f_ub x] (and [f_ub x < f_lb x] otherwise). + + Note that, unlike [partial_map_dec], this function can be used with + imprecise approximations: we allow [f_lb x < f x < f_ub x]. *) + + val map_inc_to_set : ('a -> bnd) -> 'a union -> set + (** [map_inc_to_set f u] is a variant of [approx_map_inc_to_set] when the + underlying function is total and precisely known. + + See also {!map_strict_inc}. + + Note that strict monotony is {b not} required for [map_inc_to_set]. *) + + val map_dec_to_set : ('a -> bnd) -> 'a union -> set + (** [map_dec_to_set f u] is a variant of [approx_map_dec_to_set] when the + underlying function is total and precisely known. + + See also {!map_strict_dec}. + + Note that strict monotony is {b not} required for [map_dec_to_set]. *) + + val trisection_map_to_set : + bnd -> t -> (t -> set) -> (unit -> set) -> (t -> set) -> set + (** [trisection_map_to_set v u f_lt f_eq f_gt] constructs an union of + intervals by combining the image of [f_lt] on the fragment of [u] that + only contains values strictly less than [v], the image of [f_gt] on the + fragment of [u] that only contains values strictly greater than [v], and + the image of [f_eq] if [v] is contained in [u]. + + It is helpful to build piecewise monotone functions such as + multiplication, where trisection around [0] can be used. *) +end + +(** Polymorphic {!module-type-Union} implementations. *) +module type Core = sig + type explanation + (** The type of explanations. *) + + type 'a union + (** Normalized non-empty union of intervals over ['a] values. *) + + module Union(OT : OrderedType) : Union + with type 'a union = 'a union + and type value := OT.finite + and type bnd = OT.t + and type explanation := explanation +end + +module type RingType = sig + (** An ordered ring is an ordered type with addition and multiplication that + follows ring laws. *) + + include OrderedType + + val add : t -> t -> t + (** [add] will only be called with values that are compatible with its + monotonicity: its argument can never be two bounds of different kinds + (upper or lower). *) + + val zero : t + (** Neutral element for addition. *) + + val neg : t -> t + (** Inverse for addition: for all [u], [add u (neg u) = zero]. *) + + val mul : t -> t -> t + (** [mul] will be only be called with values that are compatible with + its monotonicity: its arguments can be two bounds with the same sign + and kind (upper or lower), or two bounds with opposite signs and + opposite kinds (upper or lower), but never two bounds with the same + sign and opposite kinds or opposite signs and the same kind. + + It is recommended to program defensively and raise an assertion + failure if [mul] is ever called with two bounds of the same sign and + opposite kinds or opposite signs and the same kind. *) + + val pow : int -> t -> t + (** [pow n x] raises [x] to the [n]-th power. + + @raise Invalid_arg if [n] is nonpositive. *) +end + +(** Union-of-intervals over a {!RingType}. *) +module type Ring = sig + include Union (** @inline *) + + (** {1 Ring interface} *) + + val neg : t -> t + (** [neg u] evaluates to {m \{ -x \mid x \in S \}} when [u] evaluates to + {m S}. *) + + val add : t -> t -> t + (** [add u1 u2] evaluates to {m \{ x + y \mid x \in S_1, y \in S_2 \}} when + [u1] evaluates to {m S_1} and [u2] evaluates to {m S_2}. *) + + val mul : t -> t -> t + (** [mul u1 u2] evaluates to {m \{ x \times y \mid x \in S_1, y \in S_2 \}} + when [u1] evaluates to {m S_1} and [u2] evaluates to {m S_2}. *) + + val pow : int -> t -> t + (** [pow n u] evaluates to {m \{ x^n \mid x \in S \}} when [u] evaluates to + {m S}. + + @raise Invalid_argument if [n] is nonpositive. *) +end + +module type FieldType = sig + (** An ordered field type is an ordered ring type equipped with a + multiplicative inverse. *) + + include RingType + + val inv : t -> t + (** Inverse for multiplication: for all [u], [mul u (inv u) = one]. + + [inv zero] is undefined. *) +end + +(** Union-of-intervals over a {!FieldType}. *) +module type Field = sig + include Ring (** @inline *) + + (** {1 Field interface} *) + + val inv : ?inv0:bnd interval -> t -> t + (** [inv u] computes the image of [u] by {m x \mapsto x^{-1}}. + + [inv0] is used to define the inverse of [0] and defaults to the full + interval. *) + + val div : ?inv0:bnd interval -> t -> t -> t + (** [div u1 u2] computes {m \{ x \times y^{-1} \mid x \in S_1, y \in S_2\}} + when [u1] evaluates to {m S_1} and [u2] evaluates to {m S_2}. + + [inv0] is used when the divisor can be [0]. It defaults to the full + interval. *) +end + +module type AlgebraicType = sig + (** An ordered algebraic type is an ordered field with an (approximation of) + the n-th root. *) + + include FieldType + + val root_default : int -> t -> t + (** [root_default n x] is an under-approximation of the [n]-th root of [x]. *) + + val root_excess : int -> t -> t + (** [root_default n x] is an over-approximation of the [n]-th root of [x]. *) +end + +(** Union-of-intervals over an {!AlgebraicType}. *) +module type AlgebraicField = sig + include Field (** @inline *) + + (** {1 Algebraic operations} *) + + val root : int -> t -> (t, explanation) kind + (** [root n u] computes an over-approximation of the [n]-th root over [u]. + + Note that the resulting set can be empty because not all values + necessarily have a [n]-th root (for instance negative numbers have no real + square root). + + @raise Invalid_argument if [n] is nonpositive. *) +end + +module type EuclideanType = sig + (** An ordered euclidean type is an ordered ring equipped with an euclidean + division. *) + + include RingType + + val ediv : t -> t -> t + (** [ediv x y] returns the euclidean division of [x] by [y]. + + @raise Division_by_zero if [y] is [zero]. *) +end + +(** Union-of-intervals over an {!EuclideanType}. *) +module type EuclideanRing = sig + include Ring (** @inline *) + + (** {1 Euclidean division} *) + + val ediv : ?div0:bnd interval -> t -> t -> t + (** [ediv u1 u2] computes the image of euclidean division of a value in [u1] + by a value in [u2]. + + [div0] (defaults to the full interval) is used as the image of division by + [zero]. *) +end diff --git a/src/lib/util/options.ml b/src/lib/util/options.ml index c9d90a8210..31d8d41a77 100644 --- a/src/lib/util/options.ml +++ b/src/lib/util/options.ml @@ -141,6 +141,7 @@ let debug_combine = ref false let debug_constr = ref false let debug_explanations = ref false let debug_fm = ref false +let debug_intervals = ref false let debug_fpa = ref 0 let debug_gc = ref false let debug_interpretation = ref false @@ -171,6 +172,7 @@ let set_debug_combine b = debug_combine := b let set_debug_constr b = debug_constr := b let set_debug_explanations b = debug_explanations := b let set_debug_fm b = debug_fm := b +let set_debug_intervals b = debug_intervals := b let set_debug_fpa i = debug_fpa := i let set_debug_gc b = debug_gc := b let set_debug_interpretation b = debug_interpretation := b @@ -201,6 +203,7 @@ let get_debug_combine () = !debug_combine let get_debug_constr () = !debug_constr let get_debug_explanations () = !debug_explanations let get_debug_fm () = !debug_fm +let get_debug_intervals () = !debug_intervals let get_debug_fpa () = !debug_fpa let get_debug_gc () = !debug_gc let get_debug_interpretation () = !debug_interpretation diff --git a/src/lib/util/options.mli b/src/lib/util/options.mli index 7f6e0fbced..ee2dd0ed12 100644 --- a/src/lib/util/options.mli +++ b/src/lib/util/options.mli @@ -117,6 +117,9 @@ val set_debug_explanations : bool -> unit (** Set [debug_fm] accessible with {!val:get_debug_fm} *) val set_debug_fm : bool -> unit +(** Set [debug_intervals] accessible with {!val:get_debug_intervals} *) +val set_debug_intervals : bool -> unit + (** Set [debug_fpa] accessible with {!val:get_debug_fpa} Possible values are @@ -495,6 +498,9 @@ val get_debug_uf : unit -> bool (** Get the debugging flag of inequalities. *) val get_debug_fm : unit -> bool +(** Get the debugging flag of intervals. *) +val get_debug_intervals : unit -> bool + (** Get the debugging value of floating-point. *) val get_debug_fpa : unit -> int (** Default to [0]. *) diff --git a/src/lib/util/printer.ml b/src/lib/util/printer.ml index 76de673f8f..fe71d62144 100644 --- a/src/lib/util/printer.ml +++ b/src/lib/util/printer.ml @@ -251,7 +251,10 @@ let print_wrn ?(flushed=true) ?(header=(Options.get_output_with_headers ())) let print_dbg ?(flushed=true) ?(header=(Options.get_output_with_headers ())) ?(module_name="") ?(function_name="") s = let fmt = Options.Output.get_fmt_diagnostic () in - Format.fprintf fmt "@[%s" (pp_smt clean_dbg_print); + if header then + Format.fprintf fmt "@[%s" (pp_smt clean_dbg_print) + else + Format.fprintf fmt "@[%s" (pp_smt clean_dbg_print); if header then begin let fname = if String.equal function_name "" @@ -267,11 +270,11 @@ let print_dbg ?(flushed=true) ?(header=(Options.get_output_with_headers ())) force_new_line fmt; if Options.get_output_with_colors () then Format.fprintf fmt - "@{@{[Debug]%s%s@}@} @[" + "@{@{[Debug]%s%s@}@}@ @[" mname fname else Format.fprintf fmt - "[Debug]%s%s @[" mname fname + "[Debug]%s%s@ @[" mname fname end; if flushed || Options.get_output_with_forced_flush () then Format.kfprintf flush fmt s else Format.fprintf fmt s