Verification of the SRT division table

A linear arithmetic friendly SRT division table

Verify it
SPEC

% A (slight) variant of the SRT division circuit presented by
% Clarke, German and Zhao in CAV 96.

% It is my (Nikolaj's) understanding that this model of the circuit
% does not reflect a realistic implementation very closely.
% In particular the lookup table has been modified by collapsing
% consecutive row-elements into one entry. This enables us to
% represent the lookup table in a concise manner, but quite 
% remotely from a realistic circuit design.
% The lemmas checked represent the computationally most intensive parts.
% Correctness of the division algorithm follows from additional inductive
% arguments that have not been included here. These require an application 
% of mathematical induction. 

% This benchmark is an excellent test for linear arithmetic decision procedures
% and in the way we present it here also tuples.

local rin,rin1: rat
local md,md1 : rat
in d,d1  : rat
local qdigit : int
local rout, rout1 : rat
local qsign : bool
local qpos, qneg : rat
local quotient : int

value power : rat * int --> rat
variable x : rat
SIMPLIFY : power(x,0) ---> 1
SIMPLIFY : power(x,1) ---> x
SIMPLIFY : power(x,2) ---> x * power(x,1)
SIMPLIFY : power(x,3) ---> x * power(x,2)
SIMPLIFY : power(x,4) ---> x * power(x,3)
SIMPLIFY : power(x,5) ---> x * power(x,4)
SIMPLIFY : power(x,6) ---> x * power(x,5)


macro table = 
	(( -7/2, -13/8,  -1/2, 1/2, 3/2,   3),
	 ( -7/2, -7/4,   -1/2, 1/2, 7/4,   7/2),
	 ( -4,   -2,     -3/4, 1/2, 2,     4),
	 ( -9/2, -9/4,   -3/4, 1/2, 2,     4),
	 ( -9/2, -5/2,     -1, 3/4, 9/4,   9/2),
	 ( -5, 	  -5/2,    -1, 1,   5/2,   5),
	 ( -11/2, -11/4,   -1, 1,   5/2,   5),
	 ( -11/2, -3,      -1, 1,   3,     11/2)
	)

macro lookup(d:rat) = 
      let
	 t = table;
      in
         (if d =  8/8 then #1(t) else
          if d =  9/8 then #2(t) else
          if d = 10/8 then #3(t) else
	  if d = 11/8 then #4(t) else
	  if d = 12/8 then #5(t) else
	  if d = 13/8 then #6(t) else
	  if d = 14/8 then #7(t) else #8(t))


AXIOM : rin1 <= rin /\ rin < rin1 + 1/(power(2,6))
AXIOM : md1 <= md /\ md < md1 +  1/(power(2,6))
AXIOM : (rin1 <= rin /\ rin < rin1 + 1/(power(2,6)))'
AXIOM : (md1 <= md /\ md < md1 +  1/(power(2,6)))'
AXIOM : d1 <= d /\ d < d1 + 1/power(2,3)

AXIOM : d1 = 8/8 \/ d1 = 9/8 \/ d1 = 10/8 \/ d1 = 11/8 \/ 
	d1 = 12/8 \/ d1 = 13/8 \/ d1 = 14/8 \/ d1 = 15/8 

AXIOM mux : md = d * qdigit

AXIOM mux' : (md = d * qdigit)'


AXIOM galu : 
	rout1 = if qsign 
		  then rin1 + md1 
		  else rin1 - md1 - 1/(power(2,6))

AXIOM dalu:   rout = if qsign then rin + md else rin - md 

AXIOM dalu': (rout = if qsign then rin + md else rin - md)' 


AXIOM qpos:  qpos' = if qsign then 4 * qpos else 4 * qpos + qdigit

AXIOM qneg:  qneg' = if qsign then 4 * qneg + qdigit else 4 * qneg

AXIOM quotient:   quotient = qpos - qneg

AXIOM quotient': (quotient = qpos - qneg)'


AXIOM remainder:
	rin' = 4 * rout
	
AXIOM quo logic: 
	(qdigit',qsign') = 
		let
		   row = lookup(d1);
		   b1  = #1 row;
		   b2  = #2 row;
		   b3  = #3 row;
		   b4  = #4 row;
		   b5  = #5 row;
		   b6  = #6 row;
		   r   = 4 * rout1
		in
		   (if b1 <= r /\ r <= b2-1/power(2,4) \/
		       b5 <= r /\ r <= b6-1/power(2,4)
		       then 2
		       else
		    if b2 <= r /\ r <= b3-1/power(2,4) \/
		       b4 <= r /\ r <= b5-1/power(2,4) 
		       then 1
		       else
		    if b3 <= r /\ r <= b4-1/power(2,4)
		       then 0
		       else 
		    if r < b1 \/ b6 <= r
		       then 3
		       else 4, r < b3)
 

% 0.5600 seconds
% The format of axioms mux and mux' is essential for the verification of
% this verification condition using only linear arithmetic.
% modeling the mux as a conditional statement standard linear arithmetic
% cannot establish this. The range assumption in PROPERTY 3 and 4 will
% also be relevant in this case.
PROPERTY 1 : (4 * quotient * d + rin)' = 4 * (4 * quotient * d + rin)

% This is originally used as a lemma for PROPERTY 4.
% In our verification we will not make use of it, but
% verify it anyway.
% 0.8200 seconds.
PROPERTY 2 : rout1 <= rout /\ rout < rout1 + 1/(power(2,5))


% 4660 leaves in 44.8700 seconds
% A different version of this lemma verified in 30 seconds
% and was reported in the CADE-97 paper.
PROPERTY 3 : -2/3 * d <= rout /\ rout < 2/3 *d /\ qdigit' != 4 --> 
		(qdigit = 1 \/ qdigit = 2 \/ (qdigit = 0 /\ !qsign))'


% 9320 leaves in 90.1900 seconds
PROPERTY 4 : -2/3 * d <= rout /\ rout < 2/3 *d /\ qdigit' != 4 --> 
		(-2/3 * d <= rout /\ rout < 2/3 * d)'




A discrete SRT division table

Verify it
SPEC

% This lookup table was also presented by Clarke, German and Zhao.
% It was not verified in that paper for various computational reasons.
% It was verified in the paper by Ruess, Shankar and Srivas
% where it reportedly took the better half of 3 hours. 
% PVS has been optimized later to improve on these figures.

% A very rudimentary way to model bit-vectors: as tuples of booleans
type d = bool * bool * bool * bool
type g = bool * bool * bool * bool * bool * bool * bool

% Soft-coded support for exponentiation.
value power : rat * int --> rat
variable x : rat
SIMPLIFY : power(x,0) ---> 1
SIMPLIFY : power(x,1) ---> x
SIMPLIFY : power(x,2) ---> x * power(x,1)
SIMPLIFY : power(x,3) ---> x * power(x,2)
SIMPLIFY : power(x,4) ---> x * power(x,3)
SIMPLIFY : power(x,5) ---> x * power(x,4)
SIMPLIFY : power(x,6) ---> x * power(x,5)



macro lookup(g:g,d:d) = 
      let
	  g1 = #1 g;
	  g2 = #2 g;
	  g3 = #3 g;
	  g4 = #4 g;
	  g5 = #5 g;
	  g6 = #6 g;
	  g7 = #7 g;

	  d2 = #2 d;
	  d3 = #3 d;
	  d4 = #4 d;

	  A = -(2 - g2 * g1);
	  B = -(2 - g2);
	  C = 1 + g2;
	  D = -(1-g2);
	  E = g2;

	  table = 
	   (( 3, 3, 3, 3,-2,-2,-2, A,-1,-1,0,0,1,1,2,2,2,3,3,3,3,3),
	    ( 3, 3, 3, 3,-2,-2,-2, B,-1,-1,0,0,1,1,C,2,2,2,3,3,3,3),
	    ( 3, 3, 3,-2,-2,-2,-2,-1,-1, D,0,0,1,1,1,2,2,2,2,3,3,3),
	    ( 3, 3,-2,-2,-2,-2, B,-1,-1, D,0,0,1,1,1,2,2,2,2,3,3,3),
	    ( 3, 3,-2,-2,-2,-2,-1,-1,-1, 0,0,0,E,1,1,C,2,2,2,2,3,3),
	    ( 3,-2,-2,-2,-2,-2,-1,-1,-1, 0,0,0,0,1,1,1,2,2,2,2,2,3),
	    (-2,-2,-2,-2,-2, B,-1,-1,-1, 0,0,0,0,1,1,1,2,2,2,2,2,3),
	    (-2,-2,-2,-2,-2,-1,-1,-1,-1, 0,0,0,0,1,1,1,1,2,2,2,2,2)
	   );

	 row = if d2 
		  then
		  if d3
		     then 
		     if d4 then #8 table else #7 table
		     else
		     if d4 then #6 table else #5 table
		  else
		  if d3
		     then 
		     if d4 then #4 table else #3 table
		     else
		     if d4 then #2 table else #1 table;

      in
	 (if g7
	    then
	    if g6
	       then
	       if g5 
		  then
		  if g4
		     then
		     if g3 then #11 row else #10 row
		     else
		     if g3 then #9 row else #8 row
		  else
		  if g4
		     then
		     if g3 then #7 row else #6 row
		     else
		     if g3 then #5 row else #4 row
	       else
	       if g5
		  then 
		  if g4
		     then 
		     if g3 then #3 row else #2 row
		     else
		     if g3 then #1 row else 3
		  else 3
	    else
	    if g6
	       then
	       if g5 
		  then 3
		  else
		  if g4
		     then
		     if g3 then 3 else #22 row
		     else
		     if g3 then #21 row else #20 row
	       else
	       if g5
		  then 
		  if g4
		     then 
		     if g3 then #19 row else #18 row
		     else
		     if g3 then #17 row else #16 row
		  else 
		  if g4
		     then 
		     if g3 then #15 row else #14 row
		     else
		     if g3 then #13 row else #12 row)

(* most significant bit is first *)
macro valOfD(d:d) = 1 + (#2 d)/2 + (#3 d)/4 + (#4 d)/8 
(* two's complement. most significant bit is last *)
macro valOfG(g:g) = (#1 g)/32 + (#2 g)/16 + (#3 g)/8 + (#4 g)/4 
		      + (#5 g)/2 + (#6 g) - 2*(#7 g) 

macro isApprox(x:rat,y:rat,z:int) = 
      y <= x /\ x < y + 1/power(2,z)


in    d1     : d
in    d      : rat
local rout1  : g
local qdigit : int
local rin    : rat

% Wires
macro md   = qdigit * d
macro rout = rin - md

% Latches
AXIOM : rin' = 4 * rout
AXIOM : qdigit' = lookup(rout1,d1)

% This is actually a derived property in a more elaborate
% modeling of the circuit. Here we concentrate on the 
% lookup-table where this is needed.
AXIOM : isApprox(rout,valOfG(rout1),5) 

% d1 is the truncation of the first 4 bits of d.
AXIOM : isApprox(d,valOfD(d1),3)

% 1024 leaves 156.9900 seconds
PROPERTY : -2/3 * d <= rout /\ rout < 2/3 * d --> qdigit' != 3

% 84.4200 seconds
PROPERTY : 
	-2/3 * d <= rout /\ rout < 2/3 * d --> 
		(-2/3 * d <= rout /\ rout < 2/3 * d)'

% This is the combined statement of correct table-lookup.
% qdigit' != 3 means that we never get out of bounds (the pentium bug)
% |rout'| <= 2/3 * d means that the remainder always stays within 
% the allowed approximation of d.
% 1955 leaves in 219.9100 seconds.
PROPERTY Correct estimate and in bounds : 
	-2/3 * d <= rout /\ rout < 2/3 * d --> 
		(-2/3 * d <= rout /\ rout < 2/3 * d /\ qdigit != 3)'