1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
|
-----------------------------------------------------------------------------
-- |
-- Module : Documentation.SBV.Examples.CodeGeneration.GCD
-- Copyright : (c) Levent Erkok
-- License : BSD3
-- Maintainer: erkokl@gmail.com
-- Stability : experimental
--
-- Computing GCD symbolically, and generating C code for it. This example
-- illustrates symbolic termination related issues when programming with
-- SBV, when the termination of a recursive algorithm crucially depends
-- on the value of a symbolic variable. The technique we use is to statically
-- enforce termination by using a recursion depth counter.
-----------------------------------------------------------------------------
{-# OPTIONS_GHC -Wall -Werror #-}
module Documentation.SBV.Examples.CodeGeneration.GCD where
import Data.SBV
import Data.SBV.Tools.CodeGen
-- $setup
-- >>> -- For doctest purposes only:
-- >>> import Data.SBV
-- >>> import Data.SBV.Tools.CodeGen
-----------------------------------------------------------------------------
-- * Computing GCD
-----------------------------------------------------------------------------
-- | The symbolic GCD algorithm, over two 8-bit numbers. We define @sgcd a 0@ to
-- be @a@ for all @a@, which implies @sgcd 0 0 = 0@. Note that this is essentially
-- Euclid's algorithm, except with a recursion depth counter. We need the depth
-- counter since the algorithm is not /symbolically terminating/, as we don't have
-- a means of determining that the second argument (@b@) will eventually reach 0 in a symbolic
-- context. Hence we stop after 12 iterations. Why 12? We've empirically determined that this
-- algorithm will recurse at most 12 times for arbitrary 8-bit numbers. Of course, this is
-- a claim that we shall prove below.
sgcd :: SWord8 -> SWord8 -> SWord8
sgcd a b = go a b 12
where go :: SWord8 -> SWord8 -> SWord8 -> SWord8
go x y c = ite (c .== 0 .|| y .== 0) -- stop if y is 0, or if we reach the recursion depth
x
(go y y' (c-1))
where (_, y') = x `sQuotRem` y
-----------------------------------------------------------------------------
-- * Verification
-----------------------------------------------------------------------------
{- $VerificationIntro
We prove that 'sgcd' does indeed compute the common divisor of the given numbers.
Our predicate takes @x@, @y@, and @k@. We show that what 'sgcd' returns is indeed a common divisor,
and it is at least as large as any given @k@, provided @k@ is a common divisor as well.
-}
-- | We have:
--
-- >>> prove sgcdIsCorrect
-- Q.E.D.
sgcdIsCorrect :: SWord8 -> SWord8 -> SWord8 -> SBool
sgcdIsCorrect x y k = ite (y .== 0) -- if y is 0
(k' .== x) -- then k' must be x, nothing else to prove by definition
(isCommonDivisor k' .&& -- otherwise, k' is a common divisor and
(isCommonDivisor k .=> k' .>= k)) -- if k is a common divisor as well, then k' is at least as large as k
where k' = sgcd x y
isCommonDivisor a = z1 .== 0 .&& z2 .== 0
where (_, z1) = x `sQuotRem` a
(_, z2) = y `sQuotRem` a
-----------------------------------------------------------------------------
-- * Code generation
-----------------------------------------------------------------------------
{- $VerificationIntro
Now that we have proof our 'sgcd' implementation is correct, we can go ahead
and generate C code for it.
-}
-- | This call will generate the required C files. The following is the function
-- body generated for 'sgcd'. (We are not showing the generated header, @Makefile@,
-- and the driver programs for brevity.) Note that the generated function is
-- a constant time algorithm for GCD. It is not necessarily fastest, but it will take
-- precisely the same amount of time for all values of @x@ and @y@.
--
-- > /* File: "sgcd.c". Automatically generated by SBV. Do not edit! */
-- >
-- > #include <stdio.h>
-- > #include <stdlib.h>
-- > #include <inttypes.h>
-- > #include <stdint.h>
-- > #include <stdbool.h>
-- > #include "sgcd.h"
-- >
-- > SWord8 sgcd(const SWord8 x, const SWord8 y)
-- > {
-- > const SWord8 s0 = x;
-- > const SWord8 s1 = y;
-- > const SBool s3 = s1 == 0;
-- > const SWord8 s4 = (s1 == 0) ? s0 : (s0 % s1);
-- > const SWord8 s5 = s3 ? s0 : s4;
-- > const SBool s6 = 0 == s5;
-- > const SWord8 s7 = (s5 == 0) ? s1 : (s1 % s5);
-- > const SWord8 s8 = s6 ? s1 : s7;
-- > const SBool s9 = 0 == s8;
-- > const SWord8 s10 = (s8 == 0) ? s5 : (s5 % s8);
-- > const SWord8 s11 = s9 ? s5 : s10;
-- > const SBool s12 = 0 == s11;
-- > const SWord8 s13 = (s11 == 0) ? s8 : (s8 % s11);
-- > const SWord8 s14 = s12 ? s8 : s13;
-- > const SBool s15 = 0 == s14;
-- > const SWord8 s16 = (s14 == 0) ? s11 : (s11 % s14);
-- > const SWord8 s17 = s15 ? s11 : s16;
-- > const SBool s18 = 0 == s17;
-- > const SWord8 s19 = (s17 == 0) ? s14 : (s14 % s17);
-- > const SWord8 s20 = s18 ? s14 : s19;
-- > const SBool s21 = 0 == s20;
-- > const SWord8 s22 = (s20 == 0) ? s17 : (s17 % s20);
-- > const SWord8 s23 = s21 ? s17 : s22;
-- > const SBool s24 = 0 == s23;
-- > const SWord8 s25 = (s23 == 0) ? s20 : (s20 % s23);
-- > const SWord8 s26 = s24 ? s20 : s25;
-- > const SBool s27 = 0 == s26;
-- > const SWord8 s28 = (s26 == 0) ? s23 : (s23 % s26);
-- > const SWord8 s29 = s27 ? s23 : s28;
-- > const SBool s30 = 0 == s29;
-- > const SWord8 s31 = (s29 == 0) ? s26 : (s26 % s29);
-- > const SWord8 s32 = s30 ? s26 : s31;
-- > const SBool s33 = 0 == s32;
-- > const SWord8 s34 = (s32 == 0) ? s29 : (s29 % s32);
-- > const SWord8 s35 = s33 ? s29 : s34;
-- > const SBool s36 = 0 == s35;
-- > const SWord8 s37 = s36 ? s32 : s35;
-- > const SWord8 s38 = s33 ? s29 : s37;
-- > const SWord8 s39 = s30 ? s26 : s38;
-- > const SWord8 s40 = s27 ? s23 : s39;
-- > const SWord8 s41 = s24 ? s20 : s40;
-- > const SWord8 s42 = s21 ? s17 : s41;
-- > const SWord8 s43 = s18 ? s14 : s42;
-- > const SWord8 s44 = s15 ? s11 : s43;
-- > const SWord8 s45 = s12 ? s8 : s44;
-- > const SWord8 s46 = s9 ? s5 : s45;
-- > const SWord8 s47 = s6 ? s1 : s46;
-- > const SWord8 s48 = s3 ? s0 : s47;
-- >
-- > return s48;
-- > }
genGCDInC :: IO ()
genGCDInC = compileToC Nothing "sgcd" $ do
x <- cgInput "x"
y <- cgInput "y"
cgReturn $ sgcd x y
|
