diff --git a/src/cdomains/affineEquality/sparseImplementation/listMatrix.ml b/src/cdomains/affineEquality/sparseImplementation/listMatrix.ml index 29a319a6e4..75b557fb20 100644 --- a/src/cdomains/affineEquality/sparseImplementation/listMatrix.ml +++ b/src/cdomains/affineEquality/sparseImplementation/listMatrix.ml @@ -4,6 +4,8 @@ open RatOps open Batteries +module M = Messages + let timing_wrap = Vector.timing_wrap module type SparseMatrix = @@ -16,6 +18,10 @@ sig val rref_vec: t -> vec -> t Option.t val rref_matrix: t -> t -> t Option.t + + val linear_disjunct: t -> t -> t + (** [linear_disjunct m1 m2] returns a matrix that contains the linear disjunct of [m1] and [m2]. + The result is in rref. If [m1] and [m2] are not linearly disjunct, an exception is raised. *) end module type SparseMatrixFunctor = @@ -296,29 +302,34 @@ module ListMatrix: SparseMatrixFunctor = @param v A vector with number of entries equal to the number of columns of [v]. *) let rref_vec m v = - if is_empty m then (* In this case, v is normalized and returned *) - BatOption.map (fun (_, value) -> init_with_vec @@ div_row v value) (V.find_first_non_zero v) - else (* We try to normalize v and check if a contradiction arises. If not, we insert v at the appropriate place in m (depending on the pivot) *) - let pivot_positions = get_pivot_positions m in - (* filtered_pivots are only the pivots which have a non-zero entry in the corresponding column of v. Only those are relevant to subtract from v *) - let filtered_pivots = List.rev @@ fst @@ List.fold_left (fun (res, pivs_tail) (col_idx, value) -> - let pivs_tail = List.drop_while (fun (_, piv_col, _) -> piv_col < col_idx) pivs_tail in (* Skipping until possible match of both cols *) - match pivs_tail with - | [] -> (res, []) - | (row_idx, piv_col, row) :: ps when piv_col = col_idx -> ((row_idx, piv_col, row, value) :: res, ps) - | _ -> (res, pivs_tail) - ) ([], pivot_positions) (V.to_sparse_list v) in - let v_after_elim = List.fold_left (fun acc (row_idx, pivot_position, piv_row, v_at_piv) -> - sub_scaled_row acc piv_row v_at_piv - ) v filtered_pivots in - match V.find_first_non_zero v_after_elim with (* now we check for contradictions and finally insert v *) - | None -> Some m (* v is zero vector and was therefore already covered by m *) - | Some (idx, value) -> - if idx = (num_cols m - 1) then - None - else - let normalized_v = V.map_f_preserves_zero (fun x -> x /: value) v_after_elim in - Some (insert_v_according_to_piv m normalized_v idx pivot_positions) + let res = + if is_empty m then (* In this case, v is normalized and returned *) + BatOption.map (fun (_, value) -> init_with_vec @@ div_row v value) (V.find_first_non_zero v) + else (* We try to normalize v and check if a contradiction arises. If not, we insert v at the appropriate place in m (depending on the pivot) *) + let pivot_positions = get_pivot_positions m in + (* filtered_pivots are only the pivots which have a non-zero entry in the corresponding column of v. Only those are relevant to subtract from v *) + let filtered_pivots = List.rev @@ fst @@ List.fold_left (fun (res, pivs_tail) (col_idx, value) -> + let pivs_tail = List.drop_while (fun (_, piv_col, _) -> piv_col < col_idx) pivs_tail in (* Skipping until possible match of both cols *) + match pivs_tail with + | [] -> (res, []) + | (row_idx, piv_col, row) :: ps when piv_col = col_idx -> ((row_idx, piv_col, row, value) :: res, ps) + | _ -> (res, pivs_tail) + ) ([], pivot_positions) (V.to_sparse_list v) + in + let v_after_elim = List.fold_left (fun acc (row_idx, pivot_position, piv_row, v_at_piv) -> + sub_scaled_row acc piv_row v_at_piv + ) v filtered_pivots in + match V.find_first_non_zero v_after_elim with (* now we check for contradictions and finally insert v *) + | None -> Some m (* v is zero vector and was therefore already covered by m *) + | Some (idx, value) -> + if idx = (num_cols m - 1) then + None + else + let normalized_v = V.map_f_preserves_zero (fun x -> x /: value) v_after_elim in + Some (insert_v_according_to_piv m normalized_v idx pivot_positions) + in + if M.tracing then M.trace "rref_vec" "rref_vec: m:\n%s, v: %s => res:\n%s" (show m) (V.show v) (match res with None -> "None" | Some r -> show r); + res let rref_vec m v = timing_wrap "rref_vec" (rref_vec m) v @@ -372,4 +383,129 @@ module ListMatrix: SparseMatrixFunctor = let is_covered_by m1 m2 = timing_wrap "is_covered_by" (is_covered_by m1) m2 + (** Direct implementation of https://doi.org/10.1007/BF00268497 , chapter 5.2 Calculation of Linear Disjunction + also available at https://www-apr.lip6.fr/~mine/enseignement/mpri/attic/2014-2015/exos/karr.pdf + only difference is the implementation being optimized for sparse matrices in row representation + *) + let linear_disjunct m1 m2 = + let maxcols = num_cols m1 in + let inverse_termorder = fun x y -> y - x in + let rev_matrix = List.map (fun x -> V.of_sparse_list (V.length x) (List.rev @@ V.to_sparse_list x) ) in + let del_col m i = List.map (fun v -> V.tail_afterindex v i) m in + let safe_get_row m i = + try List.nth m i with + | Invalid_argument _ -> V.zero_vec (num_cols m) (* if row is empty, we return zero *) + in + let safe_remove_row m i = + try remove_row m i with (* remove_row can fail for sparse representations *) + | Invalid_argument _ -> m (* if row is empty, we return the original matrix *) + in + + let col_and_rc m colidx rowidx = + let col = get_col_upper_triangular m colidx in + let rc = try V.nth col rowidx with (* V.nth could be integrated into get_col for the last few bits of performance... *) + | Invalid_argument _ -> A.zero (* if col is empty, we return zero *) in + col, rc + in + + let push_col m colidx col = + List.mapi (fun idx row -> + match V.nth col idx with + | valu when A.equal A.zero valu -> row (* if the value is zero, we do not change the row *) + | valu -> V.push_first row colidx valu + | exception _ -> row + ) m + in + + let case_two a r col_b = + let a_r = get_row a r in + let res = map2i (fun i x y -> if i < r then + V.map2_f_preserves_zero (fun u j -> u +: y *: j) x a_r + else x) a col_b in + if M.tracing then M.trace "linear_disjunct_cases" "case_two: \na:\n%s, r:%d,\n col_b: %s, a_r: %s, => res:\n%s" (show a) r (V.show col_b) (V.show a_r) (show res); + res + in + + let case_three col1 col2 m1 m2 result ridx cidx = (* no new pivots at ridx/cidx *) + let sub_and_lastterm c1 c2 = (* return last element/idx pair that differs*) + let len = V.length c1 in + let c1 = V.to_sparse_list c1 in + let c2 = V.to_sparse_list c2 in + let rec sub_and_last_aux (acclist,acc) c1 c2 = + match c1, c2 with + | (i,_)::_,_ when i >= ridx -> (acclist,acc) (* we are done, no more entries in c1 that are relevant *) + | (i1, v1) :: xs1, (i2, v2) :: xs2 when i1 = i2 -> + let res = A.sub v1 v2 in + let acc = if A.equal res A.zero then acc else Some (i1, v1, v2) in + sub_and_last_aux ((i1,res)::acclist,acc) xs1 xs2 + | (i1, v1) :: xs1, (i2, v2) :: xs2 when i1 < i2 -> sub_and_last_aux ((i1,v1)::acclist,Some (i1,v1,A.zero)) xs1 ((i2, v2)::xs2) + | (i1, v1) :: xs1, (i2, v2) :: xs2 (* when i1 > i2 *)-> sub_and_last_aux ((i2,A.neg v2)::acclist,Some (i2,A.zero,v2)) ((i1, v1)::xs1) xs2 + | (i,v)::xs ,[] -> sub_and_last_aux ((i,v)::acclist,Some (i,v,A.zero)) xs [] + | [], (i,v)::xs -> sub_and_last_aux ((i,v)::acclist,Some (i,A.zero,v)) [] xs + | [], [] -> (acclist,acc) + in + let resl,rest = sub_and_last_aux ([],None) c1 c2 in + if M.tracing then M.trace "linear_disjunct_cases" "sub_and_last: ridx: %d c1: %s, c2: %s, resultlist: %s, result_pivot: %s" ridx (V.show col1) (V.show col2) (String.concat "," (List.map (fun (i,v) -> Printf.sprintf "(%d,%s)" i (A.to_string v)) resl)) (match rest with None -> "None" | Some (i,v1,v2) -> Printf.sprintf "(%d,%s,%s)" i (A.to_string v1) (A.to_string v2)); + V.of_sparse_list len (List.rev resl), rest + in + let coldiff,lastdiff = sub_and_lastterm col1 col2 in + match lastdiff with + | None -> + let sameinboth=get_col_upper_triangular m1 cidx in + if M.tracing then M.trace "linear_disjunct_cases" "case_three: no difference found, cidx: %d, ridx: %d, coldiff: %s, sameinboth: %s" cidx ridx (V.show coldiff) (V.show sameinboth); + (del_col m1 cidx, del_col m2 cidx, push_col result cidx sameinboth, ridx) (* No difference found -> (del_col m1 cidx, del_col m2 cidx, push hd to result, ridx)*) + | Some (idx,x,y) -> + let r1 = safe_get_row m1 idx in + let r2 = safe_get_row m2 idx in + let resrow = safe_get_row result idx in + let diff = x -: y in + let multiply_by_t termorder m t = + map2i (fun i x c -> if i <= ridx then + let beta = c /: diff in + V.map2_f_preserves_zero_helper (termorder) (fun u j -> u -: (beta *: j)) x t + else x) m coldiff + in + let transformed_res = multiply_by_t (inverse_termorder) result resrow in + let transformed_a = multiply_by_t (-) m1 r1 in + let alpha = get_col_upper_triangular transformed_a cidx in + let res = push_col transformed_res cidx alpha in + if M.tracing then M.trace "linear_disjunct_cases" "case_three: found difference at ridx: %d idx: %d, x: %s, y: %s, diff: %s, m1: \n%s, m2:\n%s, res:\n%s" + ridx idx (A.to_string x) (A.to_string y) (A.to_string diff) (show m1) (show m2) (show @@ rev_matrix res); + safe_remove_row (transformed_a) idx, safe_remove_row (multiply_by_t (-) m2 r2) idx, safe_remove_row (res) idx, ridx - 1 + in + + let rec lindisjunc_aux currentrowindex currentcolindex m1 m2 result = + if M.tracing then M.trace "linear_disjunct" "result so far: \n%s, currentrowindex: %d, currentcolindex: %d, m1: \n%s, m2:\n%s " + (show @@ rev_matrix result) currentrowindex currentcolindex (show m1) (show m2); + if currentcolindex >= maxcols then result + else + let col1, rc1 = col_and_rc m1 currentcolindex currentrowindex in + let col2, rc2 = col_and_rc m2 currentcolindex currentrowindex in + match Z.to_int @@ A.get_num rc1, Z.to_int @@ A.get_num rc2 with + | 1, 1 -> lindisjunc_aux + (currentrowindex + 1) (currentcolindex+1) + (del_col m1 currentrowindex) (del_col m2 currentrowindex) + (List.mapi (fun idx row -> if idx = currentrowindex then V.push_first row currentcolindex A.one else row) result) + | 1, 0 -> let beta = get_col_upper_triangular m2 currentcolindex in + if M.tracing then M.trace "linear_disjunct_cases" "case 1,0: currentrowindex: %d, currentcolindex: %d, m1: \n%s, m2:\n%s , beta %s" currentrowindex currentcolindex (show m1) (show m2) (V.show beta); + lindisjunc_aux + (currentrowindex) (currentcolindex+1) + (safe_remove_row (case_two m1 currentrowindex col2) currentrowindex) (safe_remove_row m2 currentrowindex) + (safe_remove_row (push_col result currentcolindex beta) currentrowindex) + | 0, 1 -> let beta = get_col_upper_triangular m1 currentcolindex in + if M.tracing then M.trace "linear_disjunct_cases" "case 0,1: currentrowindex: %d, currentcolindex: %d, m1: \n%s, m2:\n%s , beta %s" currentrowindex currentcolindex (show m1) (show m2) (V.show beta); + lindisjunc_aux + (currentrowindex) (currentcolindex+1) + (safe_remove_row m1 currentrowindex) (safe_remove_row (case_two m2 currentrowindex col1) currentrowindex) + (safe_remove_row (push_col result currentcolindex beta) currentrowindex) + | 0, 0 -> let m1 , m2, result, currentrowindex = case_three col1 col2 m1 m2 result currentrowindex currentcolindex in + lindisjunc_aux currentrowindex (currentcolindex+1) m1 m2 result (* we need to process m1, m2 and result *) + | a,b -> failwith ("matrix not in rref m1: " ^ (string_of_int a) ^ (string_of_int b)^(show m1) ^ " m2: " ^ (show m2)) + in + (* create a totally empty intial result, with dimensions rows x cols *) + let pseudoempty = BatList.make (max (num_rows m1) (num_rows m1)) (V.zero_vec (num_cols m1)) in + let res = rev_matrix @@ lindisjunc_aux 0 0 m1 m2 pseudoempty in + if M.tracing then M.tracel "linear_disjunct" "linear_disjunct between \n%s and \n%s =>\n%s" (show m1) (show m2) (show res); + res + end diff --git a/src/cdomains/affineEquality/sparseImplementation/sparseVector.ml b/src/cdomains/affineEquality/sparseImplementation/sparseVector.ml index 88879e21de..d1c6a8b859 100644 --- a/src/cdomains/affineEquality/sparseImplementation/sparseVector.ml +++ b/src/cdomains/affineEquality/sparseImplementation/sparseVector.ml @@ -1,13 +1,19 @@ open Vector open RatOps +module M = Messages + open Batteries module type SparseVector = sig include Vector + val push_first: t -> int -> num -> t + val is_zero_vec: t -> bool + val tail_afterindex: t -> int -> t + val insert_zero_at_indices: t -> (int * int) list -> int -> t val remove_at_indices: t -> int list -> t @@ -23,6 +29,8 @@ sig val map2_f_preserves_zero: (num -> num -> num) -> t -> t -> t + val map2_f_preserves_zero_helper: (int -> int -> int) -> (num -> num -> num) -> t -> t -> t + val find2i_f_false_at_zero: (num -> num -> bool) -> t -> t -> int val apply_with_c_f_preserves_zero: (num -> num -> num) -> num -> t -> t @@ -86,7 +94,8 @@ module SparseVector: SparseVectorFunctor = let show v = let rec sparse_list_str i l = if i >= v.len then "]" - else match l with + else + match l with | [] -> (A.to_string A.zero) ^" "^ (sparse_list_str (i + 1) l) | (idx, value) :: xs -> if i = idx then (A.to_string value) ^" "^ sparse_list_str (i + 1) xs @@ -126,6 +135,13 @@ module SparseVector: SparseVectorFunctor = | Some (idx, value) when idx = n -> value | _ -> A.zero + let push_first v n num = + if n >= v.len then raise (Invalid_argument "Index out of bounds") + else let res = + {v with entries = (n,num)::v.entries} in + if M.tracing then M.trace "push_first" "pushed %s at index %d, new length: %d, resulting in %s" (A.to_string num) n res.len (res.entries |> List.map (fun (i, x) -> Printf.sprintf "(%d, %s)" i (A.to_string x)) |> String.concat ", "); + res + (** [set_nth v n num] returns [v] where the [n]-th entry has been set to [num]. @raise Invalid_argument if [n] is out of bounds. @@ -156,6 +172,22 @@ module SparseVector: SparseVectorFunctor = in {entries = add_indices_helper v.entries indices 0 []; len = v.len + num_zeros} + (** + [tail_afterindex v n] returns the vector starting after the [n]-th entry, i.e. all entries with index > [n]. + @raise Invalid_argument if [n] is out of bounds. + *) + let tail_afterindex v n = + if n >= v.len then raise (Invalid_argument "Index out of bounds") + else + match v.entries with + | [] -> v (* If the vector is empty, return it as is *) + | (headidx,headval) :: _ -> + if M.tracing then M.trace "tail_afterindex" "headidx: %d, n: %d, v.len: %d" headidx n v.len; + if headidx > n then v + else + let entries = List.tl v.entries in + {entries; len = v.len } + (** [remove_nth v n] returns [v] where the [n]-th entry is removed, decreasing the length of the vector by one. @raise Invalid_argument if [n] is out of bounds @@ -257,13 +289,15 @@ module SparseVector: SparseVectorFunctor = {v with entries = entries'} (** - [map2_f_preserves_zero f v v'] returns the mapping of [v] and [v'] specified by [f]. + [map2_f_preserves_zero termorder f v v'] returns the mapping of [v] and [v'] specified by [f]. Note that [f] {b must} be such that [f 0 0 = 0]! + [termorder] is a function specifying, if the entries of [v] and [v'] are ordered in increasing or decreasing index order. + @raise Invalid_argument if [v] and [v'] have unequal lengths *) - let map2_f_preserves_zero f v v' = + let map2_f_preserves_zero_helper termorder f v v' = let f_rem_zero acc idx e1 e2 = let r = f e1 e2 in if r =: A.zero then acc else (idx, r) :: acc @@ -274,7 +308,7 @@ module SparseVector: SparseVectorFunctor = | [], (yidx, yval) :: ys -> aux (f_rem_zero acc yidx A.zero yval) [] ys | (xidx, xval) :: xs, [] -> aux (f_rem_zero acc xidx xval A.zero) xs [] | (xidx, xval) :: xs, (yidx, yval) :: ys -> - match xidx - yidx with + match termorder xidx yidx with | d when d < 0 -> aux (f_rem_zero acc xidx xval A.zero) xs v2 | d when d > 0 -> aux (f_rem_zero acc yidx A.zero yval) v1 ys | _ -> aux (f_rem_zero acc xidx xval yval) xs ys @@ -282,6 +316,17 @@ module SparseVector: SparseVectorFunctor = if v.len <> v'.len then raise (Invalid_argument "Unequal lengths") else {v with entries = List.rev (aux [] v.entries v'.entries)} + (** + [map2_f_preserves_zero f v v'] returns the mapping of [v] and [v'] specified by [f]. + + Note that [f] {b must} be such that [f 0 0 = 0]! + + The entries of [v] and [v'] are assumed to be ordered in increasing index order. + + @raise Invalid_argument if [v] and [v'] have unequal lengths + *) + let map2_f_preserves_zero f v v'= map2_f_preserves_zero_helper (-) f v v' + let map2_f_preserves_zero f v1 v2 = timing_wrap "map2_f_preserves_zero" (map2_f_preserves_zero f v1) v2 (** @@ -321,6 +366,5 @@ module SparseVector: SparseVectorFunctor = let rev v = let entries = List.rev_map (fun (idx, value) -> (v.len - 1 - idx, value)) v.entries in - {entries; len = v.len} - + {entries; len = v.len} end diff --git a/src/cdomains/apron/affineEqualityDomain.apron.ml b/src/cdomains/apron/affineEqualityDomain.apron.ml index 8d15a9d249..e556ab533b 100644 --- a/src/cdomains/apron/affineEqualityDomain.apron.ml +++ b/src/cdomains/apron/affineEqualityDomain.apron.ml @@ -170,7 +170,8 @@ struct Array.modify (( *:) lcm) row; let int_arr = Array.map Mpqf.get_num row in let div = Array.fold_left Z.gcd int_arr.(0) int_arr in - Array.modify (fun x -> Z.div x div) int_arr; + if not @@ Z.equal div Z.zero then + Array.modify (fun x -> Z.div x div) int_arr; int_arr in let vec_to_constraint arr env = @@ -282,47 +283,6 @@ struct res let join a b = - let rec lin_disjunc r s a b = - if s >= Matrix.num_cols a then a else - let case_two a r col_b = - let a_r = Matrix.get_row a r in - let a = Matrix.map2i (fun i x y -> if i < r then Vector.map2_f_preserves_zero (fun u j -> u +: y *: j) x a_r else x) a col_b; in - Matrix.remove_row a r - in - let case_three a b col_a col_b max = - let col_a, col_b = Vector.keep_vals col_a max, Vector.keep_vals col_b max in - let col_diff = Vector.map2_f_preserves_zero (-:) col_a col_b in - if Vector.is_zero_vec col_diff then - (a, b, max) - else - ( - let col_a = Vector.rev col_a in - let col_b = Vector.rev col_b in - let i = Vector.find2i_f_false_at_zero (<>:) col_a col_b in - let (x, y) = Vector.nth col_a i, Vector.nth col_b i in - let r, diff = Vector.length col_a - (i + 1), x -: y in - let a_r, b_r = Matrix.get_row a r, Matrix.get_row b r in - let multiply_by_t m t = - Matrix.map2i (fun i' x c -> if i' <= max then (let beta = c /: diff in Vector.map2_f_preserves_zero (fun u j -> u -: (beta *: j)) x t) else x) m col_diff - in - Matrix.remove_row (multiply_by_t a a_r) r, Matrix.remove_row (multiply_by_t b b_r) r, (max - 1) - ) - in - let col_a, col_b = Matrix.get_col_upper_triangular a s, Matrix.get_col_upper_triangular b s in - let nth_zero v i = match Vector.nth v i with - | exception Invalid_argument _ -> Mpqf.zero - | x -> x - in - let a_rs, b_rs = nth_zero col_a r, nth_zero col_b r in - if not (Z.equal (Mpqf.get_den a_rs) Z.one) || not (Z.equal (Mpqf.get_den b_rs) Z.one) then failwith "Matrix not in rref form" else - begin match Z.to_int @@ Mpqf.get_num a_rs, Z.to_int @@ Mpqf.get_num b_rs with - | 1, 1 -> lin_disjunc (r + 1) (s + 1) a b - | 1, 0 -> lin_disjunc r (s + 1) (case_two a r col_b) b - | 0, 1 -> lin_disjunc r (s + 1) a (case_two b r col_a) - | 0, 0 -> let new_a, new_b, new_r = case_three a b col_a col_b r in - lin_disjunc new_r (s + 1) new_a new_b - | _ -> failwith "Matrix not in rref form" end - in if is_bot a then b else if is_bot b then @@ -334,9 +294,9 @@ struct let sup_env = Environment.lce a.env b.env in let mod_x = dim_add (Environment.dimchange a.env sup_env) x in let mod_y = dim_add (Environment.dimchange b.env sup_env) y in - {d = Some (lin_disjunc 0 0 mod_x mod_y); env = sup_env} + {d = Some (Matrix.linear_disjunct mod_x mod_y); env = sup_env} | x, y when Matrix.equal x y -> {d = Some x; env = a.env} - | x, y -> {d = Some(lin_disjunc 0 0 x y); env = a.env} + | x, y -> {d = Some(Matrix.linear_disjunct x y); env = a.env} let join a b = timing_wrap "join" (join a) b