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 add (Z a) (Z b) = Z(z3.MkAdd(a,b))
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.Add(a,b)      = add a b
      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)

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
         ]

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.