Making use of F#'s math libraries together with Z3

A short note by Byron Cook

Recent work on F#'s math libraries, together with the latest release of Z3 make for a pretty powerful mixture. In particular I find it interesting that its so easy to combine F#'s polymorphic matrix code together with the power of Z3. I recently used F#'s new matrix syntax and the new Z3 release in order to re-implement the rank function synthesis engine used within TERMINATOR. The result turned out to be so concise that I thought it would be interesting to the larger F# community. I expect that, in the future, Don will probably pick up this example and use it as an F# sample. Thus, if you're looking for an up-to-date version of this example check the F# distribution.

At the high-level we're going to build a tool that takes in a mathematical relation represented as the conjunction of linear inequalities. As an example consider "x>0 and x' = x-1 and y'>y", which is a relation stating that the new value of x is always one less than the old value of x, that x is always positive, and that y goes up. We're out to automatically prove that this relation is well-founded, meaning that if you apply it pointwise to any infinite sequence of pairs (x0,y0),(x1,y1),......... that the relation will eventually not hold on a pair. See recent lecture notes (lecture 1, lecture 2, and lecture 3) for more information.

The underlying algorithm that we'll implement is given in a paper by Podelski and Rybalchenko called "A complete method for the synthesis of linear ranking functions". The crux of the paper is in Fig. 1:

In short, the paper encourages us to think of a relation R as a matrix of coefficients applied to the pre- and post-variables. Think of A as the coefficients that effect the pre-variables in R, and A' the coefficients that effect the post-variables (i.e. the variables with 's) . The paper says that if we can find a couple of vectors (lambda 1 and lambda 2) that meet some constraints, then we have proved the relation well-founded (i.e., the program that the relation represents will terminate under all conditions).

Consider a simple example of a loop like "while(x>0 and y>0) { x := x-1; y := y + posnum(); }". As a mathematical relation we can think of this as "x>0 and x'=x-1 and y>0 and y'>=y + 1". R is matrix representing this relation:

```let R  = [ [-1;0;0;0;1]      // x > 0
; [-1;0;1;0;1]      // x'>= x-1
; [1;0;-1;0;-1]     // x'<= x-1
; [0;-1;0;0;1]      // y > 0
; [0;1;0;-1;1]      // y' > y
]
```

Each line in R represents an inequality. The first line, in this example, represents "(x * -1) + (y * 0) + (x' * 0) + (y'* 0) + 1 <= 0", i.e. "x>0". We would like to build a tool that takes relations like R and proves them well-founded via the synthesis of a ranking function and bound. For example, imagine we have:

```let run (R:int list list) =
begin match synthesis R with
| None -> Printf.printf "Not well-founded\n"
| Some(rf,b) ->  Printf.printf "Rank function: %s\n" (rf.ToString());
Printf.printf "Bound: %s\n" (b.ToString())
end;
Printf.printf "\n"
```

If we execute this function run on our relation R from above, we'd like for our procedure to print something like:

```Rank function: [1;0]
Bound: 0
```

Where the rank function [1;0] represents the coefficients on the variables (x and y). In this case the ranking function is (x * 1) + (y * 0). That is, on every iteration of R to pairs of points on an infinite sequence we know that (x * 1) + (y * 0) would go down and that (x * 1) + (y * 0) would still be larger than or equal to the bound 0. For this reason we know that R is well-founded.

Implementation

We now define a working implementation of the algorithm. To get things started we need symbolic matrices. Thus we use generic matrices of type Matrix<z3e>, where z3e are the symbolic expressions used in Z3. This is where things get a little bit interesting, as the Matrix<T> library expects that T is an arithmetic type. Thus we need to instantiate a Z3 type such that operations like + and * are defined on it. The code to get Z3 loaded up, and to declare its type to be arithmetic is as follows:

```#light

#I "c:/Program Files/Microsoft Research/Z3-1.1/bin"
#r @"Microsoft.Z3.dll"

open Microsoft.Z3
open Microsoft.FSharp.Math
open Matrix.Generic

// Load up Z3.  Note: we need to dispose of z3 at the end.
let p = new Config()
p.SetParamValue("MODEL","true")
let z3 = new Context(p)
p.Dispose()

// We make a type z3e
type z3e = Z of nativeint
let z3ePtr (Z p) = p

let integer (i:int) = Z (z3.MkNumeral(i,z3.MkRealType()))
let zero  = integer 0
let one   = integer 1
let sub (Z a) (Z b) = Z(z3.MkSub(a,b))
let mul (Z a) (Z b) = Z(z3.MkMul(a,b))
let neg z = mul (integer (-1)) z

// Register arithmetic on z3e, allowing me to use overloaded +, *, etc
let FormNumerics =
{ new INumeric<_> with
member ops.Zero          = zero
member ops.One           = one
member ops.Subtract(a,b) = sub a b
member ops.Multiply(a,b) = mul a b
member ops.Negate(a)     = neg a

// The remainder of these ops are left unfinished for now.
member ops.Abs(a)  = failwith "no abs"
member ops.Sign(a) = failwith "no sign"
member ops.ToString((x:z3e),fmt,fmtprovider) = "formula"
member ops.Parse(s,numstyle,fmtprovider) = failwith "no parsing"
}

Math.GlobalAssociations.RegisterNumericAssociation (FormNumerics)
```
To call Z3 on z3e queries we simply call the following function
```let solve (Z q) =
z3.Push();
z3.AssertCnstr(q) ;
let model = ref null in
let ans = z3.CheckAndGetModel(model) in
z3.Pop();
match ans with
| LBool.True      -> Some (!model)
| LBool.False     -> None
| _ -> failwith "Error!\n"
```

Before we implement the algorithm we also need a few helper functions on z3e expressons:

```let make_var (name:string) =  Z (z3.MkConst(name,z3.MkRealType()))
let make_const (x:int) = integer x
let conjunction xs = Z(z3.MkAnd (Array.of_list (List.map (fun (Z x) -> x) xs))) ;;

let make_constraint f (v:RowVector<z3e>) =
let vl = Vector.Generic.fold (fun x (Z y) -> Z y::x) [] (v.Transpose) in
let cs = List.map f vl in
conjunction cs

// Some infix comparison operators that take vectors * integers
let (=*) v (k:int) = make_constraint (fun (Z x) ->
Z(z3.MkEq(x,z3.MkNumeral(k,z3.MkRealType())))) v
let (>=*) v (k:int) = make_constraint (fun (Z x) ->
Z(z3.MkGe(x,z3.MkNumeral(k,z3.MkRealType())))) v
let (<*) v (k:int) = make_constraint (fun (Z x) ->
Z(z3.MkLt(x,z3.MkNumeral(k,z3.MkRealType())))) v
```

We also need some code for generating fresh variables:

```// naturals = 0,1,2.......
let naturals = Seq.unfold (fun i -> Some (i,i+1I)) 0I

// vars = v0,v1,v2,v3,...
let vars = { for i in naturals -> make_var (sprintf "v%O" i) }

let varEnumerator = vars.GetEnumerator()
let fresh_var () = assert(varEnumerator.MoveNext()); varEnumerator.Current
```

From there we can almost literally implement the algorithm as it's found in the figure in the paper. Note especially that the definition of query comes directly from Fig 1.

```let synthesis R =

// Cut the input relation into its components
// (as described in Podelski & Rybalchenko)
// R defined as "coefs (X/X') + b <= 0" where
// coefs is a r*c matrix
let coefs = Matrix.Generic.of_list R in
let r,c = coefs.Dimensions in
let b  = - coefs.[ 0 .. r-1 , c-1        .. c-1     ] in
let A  =   coefs.[ 0 .. r-1 , 0          .. (c-2)/2 ] in
let A' =   coefs.[ 0 .. r-1 , (c-2)/2 +1 .. c-2     ] in

// Make z3e versions of matrices A, A', and b
// NOTE: the current version of F#'s generic map
// doesn't work.  See below for my version
let b_symb  = Matrix.Generic.map make_const b  in
let A_symb  = Matrix.Generic.map make_const A  in
let A_symb' = Matrix.Generic.map make_const A' in

// Create z3e-vectors lambda1 and lambda2
let lambda1 = RowVector.Generic.init r (fun i -> fresh_var ())
let lambda2 = RowVector.Generic.init r (fun i -> fresh_var ())

// Construct the query as described in Podelski & Rybalchenko, Fig. 1
let query = conjunction [ lambda1 >=* 0
; lambda2 >=* 0
; lambda1 * A_symb' =* 0
; (lambda1 - lambda2) * A_symb =* 0
; lambda2 * (A_symb + A_symb') =* 0
; lambda2 * b_symb <* 0
]
in

// Search for a model via Z3
match solve query with

// Case: relation not proved well-founded, thus no ranking
// function is returned
| None ->  None

// Case: relation is well-founded, construct a ranking function and bound
| Some (m as model) ->
// Get concrete instance for lambda1 and lambda2
let get_val (Z x) = m.GetNumeralValueInt(m.Eval(x)) in
let lambda1_inst  = RowVector.Generic.map get_val lambda1
let lambda2_inst  = RowVector.Generic.map get_val lambda2
m.Dispose();

// Compute ranking function "r" and bound "delta_0" as in
// Podelski & Rybalchenko, Fig. 1
let r = lambda2_inst * A' in
let delta_0 = - lambda1_inst * b in
Some(r,delta_0)
```

Examples

Now for some examples:

```// Simple example.  R0 represents "x'=x-1 && x>0"
let R0 = [ [-1;0;0]          // x >= 0
; [-1;1;1]          // x'>= x-1
; [1;-1;-1]         // x'<= x-1
]

// Example, slide 17 from in
// http://research.microsoft.com/TERMINATOR/l2.pps
// note: relation is A'A in the slides, here it is AA'
let R1 = [ [-1;0;0;0;1]      // x > 0
; [-1;0;1;0;1]      // x'>= x-1
; [1;0;-1;0;-1]     // x'<= x-1
; [0;-1;0;0;1]      // y > 0
; [0;1;0;-1;1]      // y' > y
]

// Example 1 from Podelski and Rybalchenko.
// R2 represents loop "while (i-j>=1) do (i,j) := (i-Nat(),j+Pos()); od"
// where "Nat()" could return any natural number and
// "Pos()" could return any positive number
let R2 = [ [-1;1;0;0;1]      // -i + j + 1 <= 0
;  [-1;0;1;0;0]     // -i + i'+ 0 <= 0
;  [0;1;0;-1;1]     // j - j' + 1 <= 0
]

// Example of termination bug.  Shouldn't be well-founded.
// Represents "x'=x+1 && x>0"
let R3 = [ [-1;0;1]         // x > 0
;  [-1;1;-1]       // x'>= x+1
;  [1;-1;1]        // x'<= x+1
]
```
When we call the sequence run R0;; run R1;; run R2;; run R3;; we get the following output:
```Rank function: rowvec [1]
Bound: rowvec [0]

Rank function: rowvec [1;0]
Bound: rowvec [1]

Rank function: rowvec [1;-1]
Bound: rowvec [1]

Not well-founded
```

One final thing

If you're wanting to make this example work today using the current version of F# you'll need to use this code:

```//  #'s Generic Matrix map is not general enough in type
// Additions to library -- lifted out of example
module Matrix =
module Generic =
let map f (m:Matrix<'a>) = Matrix.Generic.init m.NumRows m.NumCols
(fun r c -> f m.[r,c])

module RowVector =
module Generic =
let map f (rv:RowVector<'a>) = RowVector.Generic.init rv.Length
(fun i -> f rv.[i])
```

Acknowledgements

Nikolaj Bjorner, James Margetson, Leonardo de Moura, and Don Syme all looked at previous versions of the code above and made helpful suggestions.