Factorial Example

The factorial of a natural number is recursively defined as follows:

     fact(0) = 1
     fact(n) = n*fact(n-1),   for n>0

A procedure that calculates the factorial of x is then specified as:

     PROCEDURE Fact(x:Nat)
          x:=x' | x'=fact(x)        <--

We want to implement this with a program that accumulates fact(x) by iterating x times. We start by introducing a local variable which will be used to store the multiplier:

     VAR y:Nat .  y:=x ;
     { y=x }  x,y:=x',y' | x'=fact(x)        <--

In order to implement this using iteration, we will use the Introduce DO-statement law. But this law does not allow us to refer to the initial value of x in the assignment part of the specification. To remove the reference to x, we introduce a constant giving it the value of $x$:

CON N:Nat | N=x .
{ y=x /\ N=x } x:=x' | x'=fact(x) <--

and use this to rewrite:

     { y=x /\ N=x }  x,y:=x',y' | x'=fact(N)        

The computation will commece by initialising the accumulated value to 1. So we get, using the leading assignment law:

{y=x /\ N=x } x:=1; x,y:=x',y' | x'=fact(N) <--

We want to push the precondition through the leading assignment, so we first transform it:

     { y=N }  x:=1

then push it:

         x:=1  {y=N}

then introduce an extra assertion:

         x:=1  {x=1} {y=N}

and merge these assertions:

         x:=1  {x=1 /\ y=N}

to get:

     x:=1;
     { x=1 /\ y=N }  x,y:=x',y' | x'=fact(N)        <--

From now on we will introduce, push, and merge assertions automatically as required.

In order to introduce a DO-statement, we require a guard, an invariant such that

       x=1 /\ y=N  ==>  inv              (1)
       inv /\ ~G   ==>  x=fact(N)        (2)

As we accumulate fact(N) in the loop, we shall decrement the multiplier y. We should do this while y is greater than zero, so we choose the guard to be

       G  ==  y>0

Choosing an invariant requires care and insight. Chapter 16 of Gries' book, The Science of Programming, describes several ways of coming up with loop invariants. Basically we want to generalise the required postcondition, to get something that is true not only on termination, but also before entering the loop. By choosing

     inv  ==  x*fact(y) = fact(N)

both conditions (1) and (2) are clearly satisfied.

The obvious choice for a variant is the value of variable y. Now the law tells us that the refined program will have the form

     DO G ->
       { inv /\ G }
       x,y:=x',y' | inv'  /\  0 <= E < E[x,y\x',y']
     OD       

By instantiating G, inv, and E, we get:

     DO y>0 ->
       { x*fact(y) = fact(N)  /\  y>0 }                     <--
       x,y:=x',y' | x'*fact(y') = fact(N)  /\   0<=y'<y      <--
     OD       

In the loop body, we want to decrement y by 1, so we strengthen the postcondition:

{ x*fact(y) = fact(N) /\ y>0 }
x,y:=x',y' | x'*fact(y') = fact(N) /\ 0<= y'<y /\ y'=y-1

We rewrite using the fact that y'=y-1:

{ x*fact(y) = fact(N) /\ y>0 }
x,y:=x',y' | x'*fact(y-1) = fact(N) /\ 0<= y-1 <y /\ y'=y-1

Since

               y>0  ==>  0<= y-1<y

we can remove the highlighted conjunct (using Strengthen Relation):

{ x*fact(y) = fact(N) /\ y>0 }
x,y:=x',y' | x'*fact(y-1) = fact(N) /\ y'=y-1

We further transform the relation while assuming the precondition: (using Strengthen Relation):

          x'*fact(y-1) = fact(N)
     ==   x'*fact(y-1) = x*fact(y)           since fact(N) = x*fact(y)
     ==   x' = x*(fact(y) / fact(y-1))       since y>0 ==> fact(y-1)>0
     ==   x' = x*(y*fact(y-1) / fact(y-1))   since y>0 
     ==   x' = x*y                           

This yields:

     x,y:=x',y' | x' = x*y  /\  y'=y-1      

We could transform this directly into a simultaneous assignment:

     x,y:=x*y,y-1      

or into a sequential composition by proceeding as follows: First, introduce a trailing assignment:

     (x:=x' | x' = x*y )  ;            <--
     y:=y-1     

then simply introduce an assignment:

     x:=x*y        

Gathering all the refined pieces together, we get:

     PROCEDURE Factorial(x:Nat)
     VAR y:Nat .
     CON N:Nat |  N=x .
     y:=x ; x:=1;
     DO y>0 -> x:=x*y ; y:=y-1 OD       

Finally, since constant N has disappeared from the code, we can remove it:

     PROCEDURE Factorial(x:Nat)
     VAR y:Nat .
     y:=x ; x:=1;
     DO y>0 -> x:=x*y ; y:=y-1 OD