Proving the correctness of code using SMT sovler.

SMT solvers are great tool in computer science. One use of these solvers to prove the correctness of code. I wrote Haskell code to compute the inverse of polynomial using extended euclidean algorithm in Galois Field GF ( 2ⁿ ). If we try to analyze the code of inversePolynomial then every thing is normal except passing the depth parameter. It is unnecessary because we know that recursion is going to terminate. Initially I wrote the code without depth and the program was not terminating so I wrote mail to Levent Erkok and according to him

Your code suffers from what’s known as the “symbolic termination”. See Section 7.4 of: http://goo.gl/0E7wkm for a discussion of this issue in the context of Cryptol; but the same comments apply here as well. Briefly, when you recurse in inversePoly, SBV does not know when to stop the recursion of both branches in the ite: Thus it keeps expanding both branches ad infinitum.

Luckily, symbolic termination is usually easy to deal with once you identify what the problem is. Typically, one needs to introduce a recursion depth counter, which would stop recursion when it reaches a pre-determined depth that is determined to be safe. In fact, SBV comes with the regular GCD example that talks about this issue in detail:

http://hackage.haskell.org/packages/archive/sbv/2.10/doc/html/Data-SBV-Examples-CodeGeneration-GCD.html

I’d recommend reading the source code of the above; it should clarify the problem. You’ll need to find out a recursion-depth that’ll work for your example (try running it on concrete values and coming up with a safe number); and then prove that the bound is actually correct as explained in the above. The recursion depth bound just needs to be an overapproximation: If the algorithm needs 20 iterations but you use 30 as your recursion depth, nothing will go wrong; it’ll just produce bigger code and thus be less efficient and maybe more complicated for the backend prover. Of course, if you put a depth of 18 when 20 is needed, then it will be wrong and you won’t be able to establish correctness.

My verification condition is inverse of inverse of polynomial is same as polynomial. Right now I have no idea about proving the upper bound on depth so I just took 20. If you have any suggestion/comments then please let me know.

import Data.SBV.Bridge.Z3
import Data.Maybe


inversePolynomial ::  [ Int ] -> [ Int ] ->  SWord16
inversePolynomial poly reducer =  inversePoly  reducer' rem' first' second' depth'  where
                        poly' =  polynomial poly :: SWord16
                        reducer' =  polynomial reducer :: SWord16
                        first' =  0 :: SWord16
                        second' = 1 :: SWord16
                        depth' = 20 :: SWord16
                        ( quot', rem' ) = pDivMod poly' reducer'

inversePoly :: SWord16 -> SWord16 -> SWord16 -> SWord16  ->  SWord16 -> SWord16
inversePoly  reducer'' poly'' first'' second'' depth'' =
               ite ( depth'' .== 0 ||| rem'' .== ( 0 :: SWord16 ) ) (  second'' )  ( inversePoly  poly'' rem'' second'' ret ( depth'' - 1 ) )  where
                                   ( quot'', rem'' ) = pDivMod reducer'' poly''
                                   ret =  pAdd ( pMult ( quot'', second'', [] ) ) first''

isInversePolyCorrect :: SWord16 -> SWord16 -> SBool
isInversePolyCorrect poly reducer = inversePoly reducer ( inversePoly reducer rem first second depth ) first second depth .== rem where
                        ( quot, rem ) = pDivMod  poly reducer
                        first = 0 :: SWord16
                        second = 1 :: SWord16
                        depth = 20 :: SWord16

*Main> showPoly ( inversePolynomial  [ 6,4,1,0] [ 8,4,3,1,0] )
"x^7 + x^6 + x^3 + x"
*Main> proveWith z3 { verbose=False }  $ forAll ["x"] $ \x ->  pMod x ( 283 :: SWord16 ) ./=0 ==> isInversePolyCorrect x ( 283 :: SWord16 ) 
Q.E.D.
*Main> proveWith cvc4 { verbose=False }  $ forAll ["x"] $ \x ->  pMod x ( 283 :: SWord16 ) ./=0 ==> isInversePolyCorrect x ( 283 :: SWord16 ) 
*** Exception: fd:10: hFlush: resource vanished (Broken pipe)
*Main> proveWith yices { verbose=False }  $ forAll ["x"] $ \x ->  pMod x ( 283 :: SWord16 ) ./=0 ==> isInversePolyCorrect x ( 283 :: SWord16 ) 
Q.E.D.

Computing a^n using binary representation of natural number in Agda

While following the Agda tutorial, I wrote this code to compute aⁿ. Not very elegant example probably because I am not expert and still learning about dependent types. Maybe after getting some more experience and knowledge, I will try to prove the correctness and complexity of this algorithm. See the post of Twan van Laarhoven about proving correctness and complexity of merge sort in Agda.

module PowerFunction where

data ℕ⁺ : Set where
  one : ℕ⁺ 
  double : ℕ⁺ → ℕ⁺
  double⁺¹ : ℕ⁺ → ℕ⁺


add : ℕ⁺ → ℕ⁺ → ℕ⁺ 
add one one = double one
add one ( double x ) = double⁺¹ x
add one ( double⁺¹ x ) = double ( add x one )
add ( double x ) one =  double⁺¹ x
add ( double x ) ( double y ) = double ( add x y )
add ( double x ) (  double⁺¹ y ) =  double⁺¹ (  add x y )
add (  double⁺¹ x ) one =   double (  add x one )
add (  double⁺¹ x ) (  double y ) =  double⁺¹ (  add x y  )
add (  double⁺¹ x ) (  double⁺¹ y ) = double ( add (  add x y ) one )


mult : ℕ⁺ → ℕ⁺ → ℕ⁺
mult x one = x
mult x ( double y ) = double ( mult x y )
mult x ( double⁺¹ y ) = add x ( double ( mult x y ) )


pow : ℕ⁺ → ℕ⁺ → ℕ⁺
pow one _ = one
pow x  one =  x
pow x ( double one ) = mult x x
pow x ( double y ) =  pow ( pow x y ) ( double one )
pow x ( double⁺¹ y ) =  mult x ( pow ( pow x y ) ( double one ) ) 

Computing 2 ^ 6
pow ( double one ) ( double ( double⁺¹ one ))
double (double (double (double (double (double one)))))
If you have suggestion then please let me know.

SPOJ 8456. PRIMITIVEROOTS

SPOJ 8456. PRIMITIVEROOTS is related to number theory. Product of primitive roots of p ( prime number ) mod p is equal to 1.The product of all primitive roots of prime p 3 is ( mod p ). History of the theory of numbers By Leonard Eugene Dickson on page 183.This problem is not allowed in Haskell. Accepted C++ code.

#include<cstdio>
#include<iostream>
#include<cmath>
using namespace std;

bool p[10010];

void prime()
	{
		for(int i=2 ; i*i <=10000 ; i++ ) 
		 if(!p[i])
		 for(int j = i*i ; j <=10000 ; j += i ) p[j]=1;
		//for(int i=2; i <= 20 ; i++) if(!p[i]) cout<<i<<" ";cout<<endl;
	}

int main()
	{
		prime();
		int n , t , x=1 ;
		for( scanf("%d",&t) ; t>0 ; t--) 
		{
			scanf("%d",&n);
			if( n == 3 ) printf("%d:2\n", x++);
			else if ( !p[n] ) printf("%d:1\n",x++);
			else printf("%d:NOTPRIME\n",x++);
			/*if( t != 1 ) printf("\n");*/
		}

	}

SPOJ 9161. Legendre symbol

SPOJ 9161. Legendre symbol is related to Jacobian symbol.The Jacobi symbol is a generalization of the Legendre symbol. I just submitted my code from Solovay–Strassen primality test and accepted.Accepted Haskell code

import Data.List
import qualified Data.ByteString.Char8 as BS
import Data.Maybe

readInt :: BS.ByteString -> Integer
readInt = fst . fromJust . BS.readInteger


jacobi::Integral a => a -> a -> a
jacobi a n
 	| a == 0 = if n == 1 then 1 else 0
	| a == 2 && ( mod n 8 == 1 || mod n 8 == 7 ) = 1
	| a == 2 && ( mod n 8 == 3 || mod n 8 == 5 ) = -1
	| a >= n = jacobi ( mod a n ) n
	| even a = jacobi 2 n * jacobi ( div a 2 ) n
	| mod a 4 == 3 && mod n 4 == 3 = -jacobi n a
	| otherwise = jacobi n a


legend :: [Integer]  -> BS.ByteString
legend (a:p:_) = BS.pack.show $  b where 
		b = jacobi a p 

main = BS.interact $ BS.unlines . map ( legend . map readInt . BS.words ) . tail . BS.lines