Using F#+Z3 for Inferring Congruence Equations

Using F# and Z3 for Inferring Congruence Equations

September 2008

The current version of the program can be found here

Introduction

Recently Andy King and Harald Sondergaard proposed a new method for synthesizing invariants of integer programs with bit-level semantics. The method is based on abstract interpretation. The main idea is to abstract relations between input and output bit vectors as a set of congruence equations modulo 2^w, where w is the length of the bit vectors. The computation of the abstract fixed point is based on congruence analysis and Boolean reasoning, implemented using SAT solvers.

King and Sondergaard compute invariants as follows: a program is divided in basic blocks. Each basic block is associated with a so-called summary. A summary represents an invariant for the given block of the program. Summaries are computed using a SAT solver: given a command and a summary of the input block, the algorithm incrementally computes a summary of a output block. Using this technique King and Sondergaard were able to synthesize invariants that no other analysis was able to derive.

In the following we show that one can use an SMT solver to get a more concise implementation of this abstract interpretation using directly bit vectors. We illustrate this idea on the same example that was used in the original paper - Wegner's fast bit counting algorithm. The algorithm takes a bit vector x as input a and returns the number c of 1 bits in x.

	y = x;
        c = 0;
	while (y != 0) {
	  y &= y - 1;
          c++;
        }

We give a top-down description of our solution. First we describe how we represent programs in Z3. A program P is as a tuple <n₁, n₂, E, m>, where n₁ denotes the number of program locations or basic blocks in P, n₂ denotes the number of variables in P, and set E represents the set of control-flow edges. The integer m is a modulus for later representations of congruence equations. A control-flow edge is modeled as a tuple <l₁, l₂, c>, where l₁ is the entry location, l₂ the exit location and c a command. A command consists of a guard and a sequence of assignments of terms to program variables.

Here is the complete specification of a program in Z3:

type location = int
type var = int
type term = | Var of var
            | Add of term * term
            | Sub of term * term
            | BvAnd of term * term
            | True
            | False
            | Num of int
            | Eq of term * term
            | Not of term
type guard      = term
type assignment = var * term
type command    = guard * assignment list
type edge       = { l_pre : location; l_post : location; command : command }

// location 0 is the initial location.
type program = { num_locations : int; num_vars : int; edges : edge list; m : int }

The following program graph represents Wegner's fast bit counting algorithm. It contains three locations and three variables (x, y and z).

The declaration of this program in Z3 looks as follows:

let sample_program =
  let varX = 0 in
  let varY = 1 in
  let varC = 2 in
  { num_locations = 3;
    num_vars = 3;
    edges =
    [
     {l_pre = 0; l_post = 1; command = (True, [(varY, Var varX); (varC, Num 0); (varX, Var varX)]) };
     {l_pre = 1; l_post = 1; command = (Not(Eq(Var varY, Num 0)), [(varY, BvAnd(Var varY, Sub(Var varY, Num 1))); (varC, Add(Var varC, Num 1)); (varX, Var varX)]) };
     {l_pre = 1; l_post = 2; command = (Eq(Var varY, Num 0), [(varX, Var varX); (varY, Var varY); (varC, Var varC)]) }
   ];
    m = 256;
  }

As mentioned previously, the abstract program associates with each program location a linear congruence. Naively speaking, the set of solutions of such a linear congruence describes the set of values that variables can have at the given program location in the concrete program. For describing the current state of the fixed point iteration of the abstract interpretation, we introduce the type abstract_interpretation_state. It contains the program and a sequence of congruences, one for each program location. It further contains a "todo" list which contains the list of locations that still need to be processed in order to reach the fixed point. The type todo_list also contains a Boolean array ("active") that keeps a flag for each location, indicating whether the location was changed during the fixed point computation. The precise definition of the type abstract_interpretation_state is as follows:

type todo_list = class
    val mutable list : int list;
    val active : bool array;
    new (num_items) = { list = []; active = Array.create num_items false }

    member this.add item =
      if not this.active.[item] then
        (this.active.[item] <- true; this.list <- item::this.list)

    member this.pop =
      match this.list with
      | [] -> None
      | item::rest -> (this.list <- rest; this.active.[item] <- false; Some item)

    member this.print =
      List.iter (printf "%d ") this.list; printf "\n"
  end

type abstract_interpretation_state =
  { program     : program;
    congruences : LinearCongruence.t array;  // one linear congruence per program location.
    todo        : todo_list;
  }

In order to explain how to construct the initial abstract interpretation state, first we define a set of linear congruences. It is modeled using a matrix A, a vector b and a modulus m: A x ≡_m b. We represent matrix A as an array of rows and each row is an array of integers. The condition wf checks whether the linear congruence is well formed, i.e. whether coefficients are in Z_m = {0, 1, ..., m-1} and whether the number of rows in A corresponds to the number of rows in b. The matrix A has n*k columns, where n is the length of a bit vector, and k is the number of program variables.

type linear_congruence = { A : int [] []; b : int []; m : int }

module LinearCongruence = struct
  type t = linear_congruence

  let wf_coefficient m coefficient  = 0 <= coefficient && coefficient < m

  let wf lc =
    let in_bounds coefficient = wf_coefficient lc.m coefficient in
    Array.length lc.A = Array.length lc.b &&
    Array.for_all (Array.for_all in_bounds) lc.A &&
    Array.for_all in_bounds lc.b

  let printf lc =
    let mid = Array.length lc.A/2 in
    Array.iteri
      (fun idx vec ->
        printf "[";
        Array.iter (printf " %d") vec;
        printf " ]";
        if idx = mid then
         printf " = [ %d ] mod %d\n" lc.b.[idx] lc.m
	else
	  printf "   [ %d ]\n" lc.b.[idx]
      )
      lc.A

end

Next, we need to define the supremum and infimum of the abstract lattice. The supremum is given by the linear congruence 0 x ≡ 0. Every vector x is in the set of solutions of this congruence, i.e. it denotes the universal set "True" which is the supremum of the concrete lattice. The matrix A consists of one row and the dimensions of that row is determined by the number of program variables. The function call "zeroVec dim" returns the vector [| 0; 0; ... ;0 |] of a dimension dim. For example, "zeroVec 3" produces a vector [| 0; 0; 0 |]. The infimum of the abstract lattice denoting the empty set "False" is constructed similarly.

let constructTrue dim m1 = { A=[| zeroVec dim |]; b=[| 0 |]; m=m1 }

let constructFalse dim m1 = { A=[| zeroVec dim |]; b=[| 1 |]; m=m1 }

In the initial state of the abstract interpretation every single value is admitted at location 0 (the initial location) and all other locations are associated with the empty set. For this reason we associate the congruence 0 x ≡ 0 with location 0 and to all other locations we associate the congruence 0 x ≡ 1. Furthermore, we add location 0 into "todo" list.


// mk_initial_congruence (program:program) -> int(location) -> LinearCongruence.t
let mk_initial_congruence (program:program) location =
  if location = 0 then
    constructTrue (program.num_vars * (log2_ceiling program.m)) program.m
  else
    constructFalse program.num_vars program.m


// mk_initial_abstract_interpretation_state -> program -> abstract_interpretation_state
let mk_initial_abstract_interpretation_state program =
  {
    program = program;
    congruences = Array.init program.num_locations (fun idx -> mk_initial_congruence program idx);
    todo = let todo = new todo_list(program.num_locations) in todo.add(0); todo
  }

Incremental Construction of Congruences

We now briefly describe how the incremental construction of congruences works. For a given program location l, let F_l be a formula denoting a set of values that program variables can have in the location l. Let e=(l_pre,l_post,c) be a control-flow edge and let T_e be a formula constructed from the command c associated with the edge e that represents the transition relation of command c. We will explain how to construct these formula later. Given the formula F_{l_pre}(x) and the formula T_e(x, x') we iteratively construct a concruence F_{l_post}(x') that abstracts the successor states of F_{l_pre}(x) under command c. Initially we set F_{l_post}(x') to the current concruence associated with location l_post. We then iteratively check the satisfiability of H (x, x') ≡ F_{l_pre}(x) AND T_e(x, x') AND NOT F_{l_post}(x') by calling Z3. If H (x, x') is unsatisfiable, then the formula F_{l_post} covers successor states of F_{l_pre} under command c and we are done. However, if the formula H(x, x') is satisfiable, then there are some values (v,v') that make it satisfiable. The values v' describe states at location l_post that are not yet covered by formula F_{l_post}. We need to refine the concruence assoicated with l_post so that it encompasses the values v'.

For checking satisfiability of formula H (x, x') we use Z3 with the parameter for "MODEL" set to "true". It means that if the formula is satisfiable, Z3 will also return one of its models. We then extract the values v' from the model that are used for the refinement of the congruence associated with location l_post. This refinement is done using operations join and project. We incrementally continue with this procedure until we find the most precise abstraction.

The above procedure of propagating constraints to all locations is implemented in the following loop. It takes the abstract interpretation state and first checks the emptiness of the todo list. If the todo list is empty, the algorithms ends and each location contains a congruence that describes an invariant for that location in the concrete program. However, if there still exists some location in the todo_list, we call function "propagate_location location state" which takes the abstract interpretation state and the location from the todo list and propagates the constraints for all its outgoing control-flow edges:

let propagate_location location (state:abstract_interpretation_state) =
  List.iter (propagate_edge location state) state.program.edges

let rec propagate (state:abstract_interpretation_state) =
  match state.todo.pop with
  | None -> ()
  | Some location ->
      propagate_location location state;
      propagate state

We explain now how the basic ideas described at the beginning of this section are encoded in the function propagate_location. For this we use the function propagate_edge. Let fmla be a Z3 formula that corresponds to H (x, x') (later in the text we explain all the technical details about constructing Z3 formulas from edges and linear congruences). Using Z3 we try to construct a model for fmla. If a model is not found, then the formula is unsatisfiable and we stop with further abstraction refinements.

If a model is found, then we have to update the linear congruence that corresponds to the location l_post. Let lc2 be a linear congruence corresponding to location l_post. We update lc2 as follows: first, based on the returned model we construct a new linear congruence lc1. The set of solutions for the new linear congruence for l_post should be the union of the set of solutions for lc2 and of the set of solutions for lc1. This new linear congruence is derived by computing the join of lc1 and lc2. Note that the join does not simply correspond to the union of the set of solutions, because linear congruences are not closed under union. Basically, in order to compute the join we compute the set of generators for both congruences and then we take the union of those generators. We explain the details later.

Once the new lc2 is computed, we update it in the abstract interpretation state. We add the start and end point of the edge in the todo list since we do not know whether we already reached a fixed point or not.

The call "let ty = z3.MkBvType(uint32 (log2_ceiling state.congruences.[0].m))" creates a bit-vector type of size ⌈log ₂m⌉ since we use the theory of bit-vectors to handle modular arithmetic.

let propagate_edge location (state:abstract_interpretation_state) (edge:edge) =
  if location = edge.l_pre then
    let fmlLocation =  mk_formulaFromCongruence_bitwise 0 state.congruences.[location] in
    let ty = z3.MkBvType(uint32 (log2_ceiling state.program.m)) in
    let fmlEdge =  mk_formula z3 state.program ty edge in
    let fmlPost =  mk_formulaFromCongruence_bitwise state.program.num_vars state.congruences.[edge.l_post] in
    let fml = z3.MkAnd [| fmlLocation ; fmlEdge; z3.MkNot(fmlPost) |] in
    let model = ref (null : TypeSafeModel) in
    z3.Push();
    z3.AssertCnstr(fml);
    (match z3.CheckAndGetModel(model) with
    | LBool.True ->
        let model = !model in
          assert (model <> null);
          model.Display(Console.Out);
          let p_vars = get_primed_vars state.program in
          let lc1 = extract_linear_congruence_bitwise model p_vars state.program.m  in
          let lc2 = constructJoin lc1 state.congruences.[edge.l_post] in
            state.congruences.[edge.l_post] <- lc2;
            state.todo.add edge.l_post;
            state.todo.add edge.l_pre;
            model.Dispose();
    | LBool.False ->
    | _ ->
	assert false
    );
    z3.Pop()

We still need to explain how to handle primed variables. We use the following command to create a variable of the type ty in Z3:

z3.MkConst(z3.MkSymbol(some integer identifier), ty))

For a program P we use integers 0,...,P.num_vars - 1 as identifiers for unprimed program variables. For primed variables we use integers P.num_vars,...,P.num_vars + P.num_vars - 1. Calling function get_primed_vars on a program P returns the list of Z3 variables that corresponds to the list of primed variables of program P.

let get_primed_vars (program:program) =
  let m = program.m in
  let ty = z3.MkBvType(uint32 (log2_ceiling m)) in
  List.init program.num_vars (fun idx -> z3.MkConst(z3.MkSymbol(program.num_vars + idx), ty))

Constructing a Z3 Formula for the Transition Relation of a Control-flow Edge

Every control-flow edge in a program P contains a command c. This command consists of two parts: a guard and a sequence of variable assignments. Assuming that we are able to construct the Z3 formula for each part, the resulting edge formula would be a conjunction of both those formulas. Let mk_term be a function that takes a term t and returns the corresponding Z3 term. It can be easily defined using pattern matching. The guard that is associated with a command can be translated directely into a Z3 formula using function mk_term. The assigments are translated into equations between primed and unprimed variables.

// mk_term takes a term and returns z3 term
let rec mk_term (z3:TypeSafeContext) ty term =
  match term with
  | Var v -> z3.MkConst(z3.MkSymbol(v), ty)
  | Add (t1,t2) -> z3.MkBvAdd ((mk_term z3 ty t1), (mk_term z3 ty t2))
  | Sub (t1, t2) -> z3.MkBvSub ((mk_term z3 ty t1), (mk_term z3 ty t2))
  | BvAnd (t1, t2) -> z3.MkBvAnd((mk_term z3 ty t1), (mk_term z3 ty t2))
  | True -> z3.MkTrue()
  | False -> z3.MkFalse()
  | Num n -> z3.MkNumeral(n, ty)
  | Eq (t1, t2) -> z3.MkEq((mk_term z3 ty t1), (mk_term z3 ty t2))
  | Not t1 -> z3.MkNot (mk_term z3 ty t1)

let mk_formula (z3:TypeSafeContext) (program:program) ty (edge:edge) =
  let (guard, assignment) = edge.command in
  let mk_var v = z3.MkConst(z3.MkSymbol(v + program.num_vars),ty) in
  let assignment = List.map (fun (v,t) -> z3.MkEq(mk_var v, mk_term z3 ty t)) assignment in
  let assignment = Array.of_list assignment in
  z3.MkAnd(mk_term z3 ty guard, z3.MkAnd assignment)

Constructing a Z3 Formula from a Linear Congruence

In this section we explain how to construct a Z3 formula that corresponds to the congruence A x ≡_m b. The function mk_formulaFromCongruence_bitwise has two parameters: the linear congruence lc and the integer offset. The last parameter is used for handling primed and unprimed variables. Let (l1, l2, c) be a control-flow edge. If function mk_formulaFromCongruence_bitwise is invoked for constructing the formula corresponding to the congruence at location l1 then the offset has value zero and if it is invoked for location l2 then the resulting formula is supposed to range over primed variables, thus, the offset has value program.num_vars.

Let mk_eq be a function that takes the ith row of the matrix A and returns the Z3 formula that corresponds to A.[i] x ≡_m b.[i]. Note that this is a single linear congruence, not a set. Then, the Z3 formula that corresponds to the whole system A x ≡_m b is a conjunction of all such row formulas.

For a given row (r.[0]; ... ; r.[n-1]) and index idx the corresponding Z3 formula is given by the equality a = b.[idx], where a is the product of the sum of all r.[j] and the variable at position j.

Recall that the variable vector has dimension n*k, where n is the length of a bit vector and k is the number of program variables. For our sample program the variable vector looks as follows: (x₀, ..., x₇, y₀, ..., y₇, c₀,..., c₇). Let x be a variable vector in the linear congruence system A x ≡_m b. Then the Z3 term corresponding to x_j is constructed as follows:

using j, determine the program variable to which x_j belongs, let var be this variable and let idx be the index of x_j in var, i.e. x_j = var_idx)
construct the bitvector (Z3 term) v that corresponds to var (offset is used for distinguishing primed and unprimed variables)
using z3.MkBvExtract, extract idx-th bit from v (the newly derived bitvector has size 1)
using z3.MkBvZeroExt, extend this vector with zeros to the equivalent bit vector of standard size (n*k)

let mk_formulaFromCongruence_bitwise offset (lc: LinearCongruence.t) =
  assert (LinearCongruence.wf lc);
  let bv_size = log2_ceiling lc.m                      in
  assert ((getDimension1 lc.A) % bv_size = 0);
  let ty = z3.MkBvType(uint32 bv_size)                 in
  let mk_var (id:int) =
    let var = id / bv_size                             in
    let idx = id % bv_size                             in
    let v  = z3.MkConst(z3.MkSymbol(offset + var), ty) in
    z3.MkBvZeroExt(uint32 (bv_size-1), z3.MkBvExtract(uint32 idx, uint32 idx, v))
  in
  let add a b = match a with None -> b | Some a -> z3.MkBvAdd(a,b) in
  let getOption a = match a with None -> z3.MkNumeral(0, ty) | Some a -> a in
  let mk_eq idx row =
    let mutable a = None in
    for i = 0 to Array.length row - 1 do
      match row.[i] with
      | 0 -> ()
      | 1 -> a <- Some (add a (mk_var i))
      | c -> a <- Some (add a (z3.MkBvMul(z3.MkNumeral(c, ty), mk_var i)))
    done;
    z3.MkEq(getOption a, z3.MkNumeral(lc.b.[idx],ty))
  in
  z3.MkAnd (Array.mapi mk_eq lc.A)

Constructing a Linear Congruence from a Model

Given a model that is produced by Z3, we want to construct the linear congruence whose set of solutions corresponds to that model. The function extract_linear_congruence_bitwise takes as parameters the model, the list of variables in which we are interested, and the modulus m. Let model(x) be the numerical value that the model assigns to variable x. Then we construct the linear congruence from the model as follows: A is the identity matrix and vector b consists of model(x) values. The identity matrix of the dimension d is constructed by invoking function idMatrix d.

let ithPlace dim place value = Array.init dim (fun x -> if x <> place then 0 else value)
let diagMatrix dim v = Array.init dim (fun row -> ithPlace dim row v)
let idMatrix dim = diagMatrix dim 1

Vector b is constructed from the returned values as follows: let var be a variable that belongs to the list of input parameters of extract_linear_congruence_bitwise. The numerical value of var in the model can be determined with the following command:

model.GetNumeralValueInt(model.Eval(var))

The corresponding value is then passed to the function expand_val which takes an integer as parameter and creates a list that corresponds to the binary representation of that integer. Finally, vector b is constructed by concatenating all those lists of binary numbers.


let extract_linear_congruence_bitwise (model:TypeSafeModel) vars m =
  let bv_size = log2_ceiling m                                    in
  let dimension = List.length vars                                in
  let extract_val var = model.GetNumeralValueInt(model.Eval(var)) in
  let expand_val v = List.init bv_size (fun idx -> if 0 <> ((1 <<< idx) &&& v) then 1 else 0) in
  let vals = List.map (extract_val >> expand_val) vars            in
  let vals = List.concat vals                                     in
  assert (List.for_all (LinearCongruence.wf_coefficient m) vals);
  {
    A = idMatrix (dimension*bv_size);
    b = Array.of_list vals;
    m = m
   }

Constructing the Join of two Linear Congruences

To construct the join of two linear congruences, we will use results presented in the paper by King and Sondergaard. They mostly rely on the results shown in the paper by Müller-Olm and Seidl on analysis of modular arithmetic. Each congruence system Ax =_m b defines an affine subset of Z_m. This affine subset can be described through the set of its generators. Calculating the join of two affine subset corresponds to calculating the set that is generated by the union of their generators. Let lc be a linear congruence. Then True join lc = True and False join lc = lc.

let constructJoin lc1 lc2 =
  if isTrue lc1 || isTrue lc2 then
    constructTrue (getDimension1 lc1.A) lc1.m
  else if isFalse lc1 && isFalse lc2 then
    constructFalse (getDimension1 lc1.A) lc1.m
  else if isFalse lc1 then
    lc2
  else if isFalse lc2 then
    lc1
  else
    constructNonTrivialJoin lc1 lc2

Computing the join for non-trivial congruences is more involved. For this purpose we use Proposition 3 in the paper by King and Sondergaard from which an algorithm for computing the join can be derived. This algorithm looks as follows:

Given A1 x = b1, A2 x = b2 create the matrix:


   [ 1     1   0   0  0 1]
   [ -b1   0  A1   0  0 0]
   [ 0   -b2   0  A2  0 0]
   [ 0     0  -I  -I  I 0]

Triangulate that matrix
Find the coordinates for projection of this new triangulated matrix
Return the projection of the triangulated matrix onto the variable vector

let constructNonTrivialJoin lc1 lc2 =
  assert (lc1.m = lc2.m);
  let dim = getDimension1 lc1.A in
  let m1 = triangulate (createJointMatrix lc1 lc2) lc1.m in
  //printf "triangulated matrix\n";
  //printMatrix m1;
  let rawNum = findDimForProjection m1 dim in
  if (rawNum = -1) then
    constructTrue dim lc1.m
  else
    let (matJ, vecJ) = projection m1 rawNum dim in
    {
      A = matJ;
      b = vecJ;
      m = lc1.m
    }

We briefly explain each of those steps. The first step is a straight-forward encoding of the bigger matrix using the pattern shown above.

let (@@) a b = Array.append a b
let (@@@) m1 m2 =
  assert (getDimension0 m1 = getDimension0 m2);
  Array.map2 Array.append m1 m2

let createJointMatrix lc1 lc2 =
  assert (getDimension1 lc1.A >= getDimension1 lc2.A);
  let dim = getDimension1 lc1.A in
  let m = lc1.m in
  let (rows1, cols1) = getDimensions lc1.A in
  let (rows2, cols2) = getDimensions lc2.A in
  let b1 = transpose (negateVector m lc1.b) in
  let b2 = transpose (negateVector m lc2.b) in
  let zV dim = zeroVec dim in
  let zC dim = zeroCol dim in
  let dM dim v = diagMatrix dim v in
  let zM rows cols = zeroMatrix rows cols in
  let matrix =
    [| [| 1;       1 |]     @@  zV cols1      @@  zV cols2       @@  zV dim       @@ [| 1|] |]  @@
    (b1       @@@ zC rows1 @@@ lc1.A          @@@ zM rows1 cols2 @@@ zM rows1 dim @@@ zC rows1) @@
    (zC rows2 @@@ b2       @@@ zM rows2 cols1 @@@ lc2.A          @@@ zM rows2 dim @@@ zC rows2) @@
    (zC dim   @@@ zC dim   @@@ dM dim (m-1)   @@@ dM dim (m-1)   @@@ dM dim 1     @@@ zC dim)
  in
  matrix

To triangulate the matrix obtained in the first step we use the algorithm described in the paper by Müller-Olm and Seidl. To obtain the triangulated matrix we cannot just apply standard Gaussian eliminations, since there are zero divisors in Z_m. Therefore, we could not use the matrix library that exists in F#, because it does not support modular arithmetic.

The algorithm of Müller-Olm and Seidl takes modular arithmetic into account. The implementation corresponds closely to their algorithm. At the end we invoke the function updateMatrix which removes potential zero vectors from the matrix.

let leading vec = Array.find_index (fun x -> x <> 0) vec

let vecDiff vec vec1 = Array.map2 (fun x y -> x - y) vec vec1
				
let multScalVec num vec = Array.map (fun x -> num * x) vec

let normVec m vec = Array.map (fun x -> norm x m) vec

let incr1 (a, b) = (a+1, b)

let rec rank num =
  match num with
  | 0 -> (0, 0)
  | _ when num % 2 <> 0 -> (0, num)
  | _ -> incr1 (rank (num / 2))

let leadingIndexEntryRankNonInvert vec =
  let index = leading vec     in
  let entry = vec.[index]     in
  let (rnk, inv) = rank entry in
  (index, entry, rnk, inv)

let is_zero vec = Array.for_all (fun x -> x = 0) vec

let is_non_zero vec =  not (is_zero vec)

let addInSet set inx vec =
  if 0 <= inx && inx < Array.length set then
    set.[inx] <- vec


let coeff r r1 i = (power2 (r - r1)) * i

let updateMatrix mat = Array.filter is_non_zero mat

let rec triangulateSet (set : int [] []) vec m =
  if is_non_zero vec then
    let (inx, ent, rnk, inv) = leadingIndexEntryRankNonInvert vec in
    assert (0 <= inx && inx < Array.length set);
    if is_zero set.[inx] then
      addInSet set inx vec
    else
      let vec1 = set.[inx] in
      let (inx1, ent1, rnk1, inv1) = leadingIndexEntryRankNonInvert vec1 in
      if rnk1 <= rnk then
	let vec2 = normVec m (vecDiff (multScalVec inv1 vec) (multScalVec (coeff rnk rnk1 inv) vec1)) in
	triangulateSet set vec2 m
      else
	let vec2 = normVec m (vecDiff (multScalVec inv vec1) (multScalVec (coeff rnk1 rnk inv1) vec)) in
	(addInSet set inx vec);
        triangulateSet set vec2 m


let triangulate ( mat : int [] []) m =
  let num_rows = getDimension0 mat in
  let num_cols = getDimension1 mat in
  let set =  Array.init num_cols (fun _ -> [||]) in
  for i = 0 to num_rows - 1 do
    if is_non_zero mat.[i] then
      let l = leading mat.[i] in
      if set.[l] = [||] then
	addInSet set l mat.[i]
      else
	triangulateSet set mat.[i] m
  done;
  updateMatrix set

The last step is to project the triangulated matrix onto the variable vector. Since the matrix is in triangulated form, we use Proposition 2 from the paper by King and Sondergaard. First we calculate the dimension for the projection which is the top-most row j of the matrix A for which (A(j,1), . . . , A(j,i-1)) = 0 holds.

let findDimForProjection (mat: int [] []) dim =
  let n = Array.length mat.[0] in
  let mutable i = Array.length mat - 1 in
  let mutable j = 0 in
  while (0 <= i) do
    if (is_zero (Array.sub mat.[i] 0 (n - 1 - dim))) then i <- i - 1
    else begin
      j <- i; i <- -2
    end
  done;
  //
  // (j = 0, i = - 2 ) = return 1 if there are at least two rows
  //      if is_zero (Array.sub mat.[i] 0 (n - 1 - dim) holds, it is covered
  //      if is_zero (Array.sub mat.[i] 0 (n - 1 - dim) does not holds, it is degeneric case, return -1
  // (j = 0, i = - 1 ) = return 0
  // (j = k > 0, i = - 2 ) = return k + 1 if not last, if last return number of rows
  //
  if (i = -1) then
    0
  else if (j = (Array.length mat - 1)) then
    if (is_zero (Array.sub mat.[j] 0 (n - 1 - dim))) then
      j
    else -1
  else j + 1

Once we find the dimension for the projection, we just extract the submatrix from the triangulated matrix

let projection (mat: int [] []) m n =
  let lM = getDimension0 mat - 1 in
  let lN = getDimension1 mat - 1 in
  let mat1 = Array.init ( lM - m + 1) (fun _ -> [||]) in
  let b1 = Array.init ( lM - m  + 1) (fun _ -> 0 ) in
  for i = m to lM  do
    let vec = Array.sub mat.[i] (lN - n ) n in
    mat1.[i - m] <- vec
  done;
  for i = m to lM do
    b1.[i - m] <- mat.[i].[lN]
  done;
  (mat1, b1)

Acknowledgements

Nikolaj Bjørner, James Margetson, Daniel Balasubramanian, Utkarsh Upadhyay, Thomas Wies and Damien Zufferey who looked at previous versions of the code above and made helpful suggestions.