My Weblog

Blog about programming and math

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 ( 2n ). 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.
Advertisements

September 21, 2013 Posted by | Haskell, Programming | , , , , , , , | Leave a comment

Counting Inversions

First week of algorithm course programming assignment is to compute the number of inversion count of given array. We have a given array A[0…n – 1] of n distinct positive integers and we have to count all the pairs ( i , j ) such that i < j and A[i] > A[j].
1. Brute Force
Brute force idea is very simple. Run two loops and keep counting the respective pairs for which i < j and A[i] > A[j]. The complexity of brute force is O( n2 ).

bruteForce :: [ Int ] -> Int
bruteForce [] = 0
bruteForce [_] = 0
bruteForce ( x : xs ) = ( length . filter ( x > ) $ xs ) + bruteForce xs

*Main> bruteForce [ 3 , 1 , 2]
Loading package array-0.4.0.1 ... linking ... done.
Loading package deepseq-1.3.0.1 ... linking ... done.
Loading package bytestring-0.10.0.0 ... linking ... done.
2
*Main> bruteForce [ 2,3,8,6,1]
5

2. Divide and Conquer
We can modify the merge sort to compute the inversion count. The idea is to compute the inversion count of during the merge phase. Inversion count of array C resulting from merging two sorted array A and B is sum of inversion count of A, inversion count of B and their cross inversion. The inversion count of single element is 0 ( base case ). The only trick is compute the cross inversion. What will be the inversion count of following array if we merge them ?
A = [ 1 , 2 , 3 ]
B = [ 4 , 5 , 6 ]
Clearly 0 because all element of A is less than the first element of B. What is the inversion count of C = [ 4 , 5 , 6 ] and D = [ 1 , 2 , 3 ] ? First element of D ( 1 ) is less than first element of C ( 4 ) so we have total three pairs ( 4 , 1 ) ( 5 , 1 ) ( 6 , 1). Similarly for 2 and 3 ( total 9 inversion count ). It reveals that if we have any element y in D is less than certain element x in C then all the element of A after x will also be greater than y and we will have those many cross inversion counts. More clear explanation. Wrote a quick Haskell code and got the answer.

import Data.List
import Data.Maybe ( fromJust )
import qualified Data.ByteString.Lazy.Char8 as BS

inversionCnt :: [ Int ] -> [ Int ] -> [ Int ] -> Int -> ( Int , [ Int ] )
inversionCnt [] ys  ret  cnt = ( cnt , reverse ret ++ ys ) 
inversionCnt xs []  ret  cnt = ( cnt , reverse ret ++ xs )
inversionCnt u@( x : xs ) v@( y : ys )  ret  cnt
      | x <= y = inversionCnt xs v  (  x : ret   )   cnt
      | otherwise = inversionCnt u ys  ( y : ret )  ( cnt + length u )  

merge ::  ( Int , [ Int ] ) -> ( Int , [ Int ] ) -> ( Int , [ Int ] )
merge ( cnt_1 , xs ) ( cnt_2 , ys ) = ( cnt_1 + cnt_2 + cnt , ret ) where
        ( cnt , ret ) = inversionCnt xs ys []   0

mergeSort ::  [ Int ] -> ( Int , [ Int ] )
mergeSort [ x ] = ( 0 , [ x ] )
mergeSort  xs   = merge ( mergeSort xs' ) ( mergeSort ys' ) where
           ( xs' , ys' ) = splitAt ( div ( length xs ) 2 ) xs
          
main =  BS.interact $ ( \x -> BS.pack $ show x ++ "\n" ) . fst 
        .  mergeSort . map ( fst . fromJust .  BS.readInt :: BS.ByteString -> Int  ) 
        . BS.lines 


Mukeshs-MacBook-Pro:Puzzles mukeshtiwari$ ghc  -fforce-recomp -O2  Inversion.hs 
[1 of 1] Compiling Main             ( Inversion.hs, Inversion.o )
Linking Inversion ...
Mukeshs-MacBook-Pro:Puzzles mukeshtiwari$ time ./Inversion <  IntegerArray.txt 
2407905288

real	0m8.002s
user	0m7.881s
sys	0m0.038s

I tried to compare this solution with GeeksforGeeks‘s solution written in C and it was a complete surprise to me.

Mukeshs-MacBook-Pro:Puzzles mukeshtiwari$ time ./a.out <  IntegerArray.txt 
 Number of inversions are 2407905288 

real	0m0.039s
user	0m0.036s
sys	0m0.002s

Almost 200 times faster! My Haskell intuition told me that I am doing certainly some wrong. I started analyzing the code and got the point. In inversionCnt function, I am doing extra work for computing the length of list u ( O ( n ) ) so rather than doing this extra work, I can pass the length.

import Data.List
import Data.Maybe ( fromJust )
import qualified Data.ByteString.Lazy.Char8 as BS


inversionCnt :: [ Int ] -> [ Int ] -> [ Int ] -> Int -> Int -> ( Int , [ Int ] )
inversionCnt [] ys  ret _ cnt = ( cnt , reverse ret ++ ys ) 
inversionCnt xs []  ret _ cnt = ( cnt , reverse ret ++ xs )
inversionCnt u@( x : xs ) v@( y : ys )  ret n cnt
      | x <= y = inversionCnt xs v  (  x : ret   )  ( pred n ) cnt
      | otherwise = inversionCnt u ys  ( y : ret ) n ( cnt + n )  

merge ::  ( Int , [ Int ] ) -> ( Int , [ Int ] ) -> ( Int , [ Int ] )
merge ( cnt_1 , xs ) ( cnt_2 , ys ) = ( cnt_1 + cnt_2 + cnt , ret ) where
        ( cnt , ret ) = inversionCnt xs ys [] ( length xs )  0

mergeSort ::  [ Int ] -> ( Int , [ Int ] )
mergeSort [ x ] = ( 0 , [ x ] )
mergeSort  xs   = merge ( mergeSort xs' ) ( mergeSort ys' ) where
           ( xs' , ys' ) = splitAt ( div ( length xs ) 2 ) xs
          
main =  BS.interact $ ( \x -> BS.pack $ show x ++ "\n" ) . fst 
        .  mergeSort . map ( fst . fromJust .  BS.readInt :: BS.ByteString -> Int  ) 
        . BS.lines 

Mukeshs-MacBook-Pro:Puzzles mukeshtiwari$ time ./Inversion <  IntegerArray.txt 
2407905288

real	0m0.539s
user	0m0.428s
sys	0m0.031s

Now 10 times slower and still lot of room for improvement but for now I am happy with this :). You can also try to solve INVCNT. Run this code on Ideone.

July 4, 2013 Posted by | Haskell, Programming | , , , , | 2 Comments

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

While following the Agda tutorial, I wrote this code to compute an. 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.

June 12, 2013 Posted by | Agda, Programming | , , , | 2 Comments

Missionaries and cannibals problem

In the missionaries and cannibals problem, three missionaries and three cannibals must cross a river using a boat which can carry at most two people, under the constraint that, for both banks, if there are missionaries present on the bank, they cannot be outnumbered by cannibals (if they were, the cannibals would eat the missionaries.) The boat cannot cross the river by itself with no people on board. We can solve this problem using depth first search. Representing the states as number of missionaries, cannibals and boat location in following way.


data Side = Side { missionaries :: Int , cannibals :: Int } deriving ( Show , Eq  )
data Loc = L | R deriving ( Show , Eq  )
data State = State { left :: Side , right :: Side , bloc :: Loc } deriving ( Show , Eq  )

Target is to move every missionaries and cannibals from left to right so our initial state is initialState and final state is finalState.


initialState = State { left =  Side 3 3  , right =  Side 0 0  , bloc = L }
finalState = State  { left =  Side 0 0  , right = Side 3 3  , bloc = R }

Possible movements
1. move ( 2 , 0 )
2. move ( 1 , 0 )
3. move ( 1 , 1 )
4. move ( 0 , 1 )
5. move ( 0 , 2 )
where move ( M , C ) is movement of M missionaries and C cannibals from one side to other side depending on boat location. If boat location is left then move ( M , C ) is moving M missionaries and C cannibals from left to right and vice versa. Using depth first search, we can explore every possible movement

path :: [ State ] -> [ State ] -> [ State ]
path []   vis = reverse vis
path  ( x : xs )  vis 
     | x  == finalState  = reverse ( x : vis )
     | elem x vis = path xs vis 
     | otherwise = path xs'  vis'  where 
              vis'  = x : vis
              xs' =   ( filter  isLegal  ( move x ) ) ++ xs

Encoding the moves and testing the legality

             
move :: State -> [ State ]
move ( State ( Side ml cl ) ( Side mr cr ) L ) = 
        [ State ( Side ( ml - 2 ) cl ) ( Side ( mr + 2 ) cr ) R 
        , State ( Side ( ml - 1 ) cl ) ( Side ( mr + 1 ) cr ) R
        , State ( Side ( ml - 1 ) ( cl - 1 ) ) ( Side ( mr + 1 ) ( cr + 1 ) ) R 
        , State ( Side ml ( cl - 2 ) ) ( Side mr ( cr + 2 ) ) R
        , State ( Side ml ( cl - 1 ) ) ( Side mr ( cr + 1 ) ) R        
        ] 
move  ( State ( Side ml cl ) ( Side mr cr ) R ) = 
        [ State ( Side ( ml + 2 ) cl ) ( Side ( mr - 2 ) cr ) L 
        , State ( Side ( ml + 1 ) cl ) ( Side ( mr - 1 ) cr ) L
        , State ( Side ( ml + 1 ) ( cl + 1 ) ) ( Side ( mr - 1 ) ( cr - 1 ) ) L 
        , State ( Side ml ( cl + 2 ) ) ( Side mr ( cr - 2 ) ) L
        , State ( Side ml ( cl + 1 ) ) ( Side mr ( cr - 1 ) ) L        
        ]           

isLegal :: State -> Bool 
isLegal ( State ( Side ml cl ) ( Side mr cr ) y ) 
        | ml == 0 = mr >= cr  && cr >= 0 
        | mr == 0 = ml >= cl  && cl >= 0
        | otherwise = ml >= cl && mr >= cr && cl >= 0 && cr >= 0  
                        

Complete source code ( In case if you are lazy 🙂 )

import Data.List


data Side = Side { missionaries :: Int , cannibals :: Int } deriving ( Show , Eq  )
data Loc = L | R deriving ( Show , Eq  )
data State = State { left :: Side , right :: Side , bloc :: Loc } deriving ( Show , Eq  )

initialState = State { left =  Side 3 3  , right =  Side 0 0  , bloc = L }
finalState = State  { left =  Side 0 0  , right = Side 3 3  , bloc = R }


path :: [ State ] -> [ State ] -> [ State ]
path []   vis = reverse vis
path  ( x : xs )  vis 
     | x  == finalState  = reverse ( x : vis )
     | elem x vis = path xs vis 
     | otherwise = path xs'  vis'  where 
              vis'  = x : vis
              xs' =   ( filter  isLegal  ( move x ) ) ++ xs
             
move :: State -> [ State ]
move ( State ( Side ml cl ) ( Side mr cr ) L ) = 
        [ State ( Side ( ml - 2 ) cl ) ( Side ( mr + 2 ) cr ) R 
        , State ( Side ( ml - 1 ) cl ) ( Side ( mr + 1 ) cr ) R
        , State ( Side ( ml - 1 ) ( cl - 1 ) ) ( Side ( mr + 1 ) ( cr + 1 ) ) R 
        , State ( Side ml ( cl - 2 ) ) ( Side mr ( cr + 2 ) ) R
        , State ( Side ml ( cl - 1 ) ) ( Side mr ( cr + 1 ) ) R        
        ] 
move  ( State ( Side ml cl ) ( Side mr cr ) R ) = 
        [ State ( Side ( ml + 2 ) cl ) ( Side ( mr - 2 ) cr ) L 
        , State ( Side ( ml + 1 ) cl ) ( Side ( mr - 1 ) cr ) L
        , State ( Side ( ml + 1 ) ( cl + 1 ) ) ( Side ( mr - 1 ) ( cr - 1 ) ) L 
        , State ( Side ml ( cl + 2 ) ) ( Side mr ( cr - 2 ) ) L
        , State ( Side ml ( cl + 1 ) ) ( Side mr ( cr - 1 ) ) L        
        ]           

isLegal :: State -> Bool 
isLegal ( State ( Side ml cl ) ( Side mr cr ) y ) 
        | ml == 0 = mr >= cr  && cr >= 0 
        | mr == 0 = ml >= cl  && cl >= 0
        | otherwise = ml >= cl && mr >= cr && cl >= 0 && cr >= 0  
                        
                         


Main> path [ initialState ] []
[State {left = Side {missionaries = 3, cannibals = 3}, right = Side {missionaries = 0, cannibals = 0}, bloc = L},State {left = Side {missionaries = 2, cannibals = 2}, right = Side {missionaries = 1, cannibals = 1}, bloc = R},State {left = Side {missionaries = 3, cannibals = 2}, right = Side {missionaries = 0, cannibals = 1}, bloc = L},State {left = Side {missionaries = 3, cannibals = 0}, right = Side {missionaries = 0, cannibals = 3}, bloc = R},State {left = Side {missionaries = 3, cannibals = 1}, right = Side {missionaries = 0, cannibals = 2}, bloc = L},State {left = Side {missionaries = 1, cannibals = 1}, right = Side {missionaries = 2, cannibals = 2}, bloc = R},State {left = Side {missionaries = 2, cannibals = 2}, right = Side {missionaries = 1, cannibals = 1}, bloc = L},State {left = Side {missionaries = 0, cannibals = 2}, right = Side {missionaries = 3, cannibals = 1}, bloc = R},State {left = Side {missionaries = 0, cannibals = 3}, right = Side {missionaries = 3, cannibals = 0}, bloc = L},State {left = Side {missionaries = 0, cannibals = 1}, right = Side {missionaries = 3, cannibals = 2}, bloc = R},State {left = Side {missionaries = 1, cannibals = 1}, right = Side {missionaries = 2, cannibals = 2}, bloc = L},State {left = Side {missionaries = 0, cannibals = 0}, right = Side {missionaries = 3, cannibals = 3}, bloc = R}]
*Main> map bloc . path [ initialState ] $ []
[L,R,L,R,L,R,L,R,L,R,L,R]

If you have any suggestion or comment then please let me know.

April 19, 2013 Posted by | Haskell, Programming | , , | Leave a comment

SPOJ Pell (Mid pelling)

Pell (Mid pelling) is related to Pell’s equation. It is similar to Project Euler 66 and SPOJ EQU2. I just wrote the quick solution from mathworld but I found Chakravala method very interesting. Accepted Haskell code.

import qualified Data.ByteString.Char8 as BS
import Data.List
import Data.Maybe ( fromJust )

continuedFraction :: Integer -> [ Integer ] 
continuedFraction n = map ( \ ( a , _ , _ ) -> a ) . iterate fun $ ( d , 0 , 1 ) 
                  where 
                  d = truncate . sqrt . fromIntegral $ n 
                  fun ( a0 , p0 , q0 ) = ( a1 , p1 , q1 ) where
                      p1 = a0 * q0 - p0
                      q1 = div ( n - p1 ^ 2 ) q0
                      a1 = div ( d + p1 ) q1 

pellSolver :: Integer -> BS.ByteString
pellSolver n 
  | perfectSqr n = BS.pack. show $ ( -1 ) 
  | otherwise =  ( BS.pack . show $ p ) `BS.append` ( BS.pack " " ) 
                 `BS.append` ( BS.pack.show $ q ) where 
           d = truncate . sqrt . fromIntegral $ n
           lst = takeWhile ( /= 2 * d ) . continuedFraction $ n 
           len = length lst
           r@( x : y : xs ) = take ( if even len then len else 2 * len ) . 
                              continuedFraction $ n
           ( p0 , p1 , q0 , q1 ) = ( x , x * y + 1 , 1 , y )
           ( p , _ ) = foldl' compute ( p1 , p0 ) $ xs
           ( q , _ ) = foldl' compute ( q1 , q0 ) $ xs
           compute :: ( Integer , Integer ) -> Integer -> ( Integer , Integer )
           compute ( p1 , p0 ) a = ( a * p1 + p0 , p1 )


perfectSqr :: Integer -> Bool
perfectSqr n = d * d == n where 
           d = truncate . sqrt . fromIntegral $ n 

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

main = BS.interact $  BS.unlines . map ( pellSolver . readI ) . tail . BS.lines

March 23, 2013 Posted by | Haskell, Programming | , , , | 1 Comment

SPOJ EMAIL ID

Finally solved EMAIL ID using python though it was very hard for me to switch from Haskell to python. While doing this problem, I learned quite a lot about regular expressions, found some cool site like pythontutor and problem solving with python . Accepted code in python.


import re

if __name__ == "__main__":
    n = int ( raw_input() )
    c = 1
    while c <= n :
        email =  re.findall ( "[a-zA-Z0-9][a-zA-Z0-9._]{4,}@[a-zA-Z0-9]+\.(?:com|edu|org|co\.in)", raw_input() )  
        t = len ( email )
        print 'Case #' + str ( c ) + ': ' + str ( t )
        for i in xrange ( t ) : print email[i]
        c += 1





March 11, 2013 Posted by | Programming, python | , , , | Leave a comment

Generating primes in parallel

We can use Miller-Rabin primality testing to test a number is prime or not with certain probability. We can use this method on chunk of data to get the primes in parallel. I wrote this code to see how it goes and it’s quite interesting. My other post about data parallelism and Don Stewart answer regarding different options to improve parallelism. The crux of this code is

primeRange :: Integer -> Integer -> [ Bool ]
primeRange m n = ( map isPrime [ m .. n ] ) `using` parListChunk 10000 rdeepseq 

For parallelism in Haskell , see Haskell-for-multicores and Real World Haskell. Code on github.

import Data.Bits
import Control.Parallel.Strategies 

powM :: Integer -> Integer -> Integer ->  Integer
powM a d n = powM' a d n where 
  powM' a d n 
        | d == 0 = 1
        | d == 1 = mod a n
        | otherwise = mod q n  where
	       p = powM'  ( mod ( a^2 ) n ) ( shiftR d 1 ) n
	       q = if (.&.) d 1 == 1 then mod ( a * p ) n else p

calSd :: Integer -> ( Integer , Integer ) 
calSd n = ( s , d ) where 
      s = until ( \x -> testBit ( n - 1) ( fromIntegral x ) )  ( +1 ) 0
      d = div ( n - 1 ) (  shiftL 1 ( fromIntegral s )  )


rabinMiller::Integer->Integer->Integer->Integer-> Bool
rabinMiller  n s d a 
   | n == a = True 
   | otherwise =  case x == 1 of 
          True -> True 
          _ ->   any ( == pred n ) . take ( fromIntegral s ) 
                      . iterate (\e -> mod ( e^2 ) n ) $ x  
        where 
              x = powM a d n 
   

isPrime::Integer-> Bool
isPrime n 
   | n <= 1 = False
   | n == 2 = True
   | even n = False
   | otherwise	= all ( == True ) . map ( rabinMiller n s d ) $  [ 2 , 3 , 5 , 7 , 11 , 13 , 17 ] where 
                ( s , d ) = calSd n 



primeRange :: Integer -> Integer -> [ Bool ]
primeRange m n = ( map isPrime [ m .. n ] ) `using` parListChunk 10000 rdeepseq 


main = do 
     let t = filter ( == True ) . primeRange 2 $ 10 ^ 6 
     print . length $ t 

Mukeshs-MacBook-Pro:Haskell mukeshtiwari$ ghc -O2 -rtsopts --make RabinMiller.hs -threaded  -fforce-recomp
[1 of 1] Compiling Main             ( RabinMiller.hs, RabinMiller.o )
Linking RabinMiller ...
Mukeshs-MacBook-Pro:Haskell mukeshtiwari$ time ./RabinMiller +RTS -N1
78498

real	0m4.841s
user	0m4.695s
sys	0m0.075s
Mukeshs-MacBook-Pro:Haskell mukeshtiwari$ time ./RabinMiller +RTS -N2
78498

real	0m3.115s
user	0m5.969s
sys	0m0.215s
Mukeshs-MacBook-Pro:Haskell mukeshtiwari$ time ./RabinMiller +RTS -N3
78498

real	0m2.687s
user	0m7.556s
sys	0m0.415s

I am no expert in parallelism or Haskell so if you have any comment or suggestion then please let me know.

February 5, 2013 Posted by | Haskell, Programming | , , , , , , | 3 Comments

Newton–Raphson method to solve Ax+Bsin(x)=C

I was reading Newton–Raphson method and realized that you can encode this whole algorithm using until function in Haskell.

Type signature of until 
ghci>:t until
until :: (a -> Bool) -> (a -> a) -> a -> a

To encode Newton-Raphson algorithm in until

until ( you testing condition ) ( function )  ( initial value of x )

TRIGALGE is related to Newton-Raphson . \sin x will be in between [ -1.0 , 1.0 ] so Ax - B \leq C \leq Ax + B . My initial guess was \frac{C-B}{A} . See this stackoverflow discussion for initial approximation. Accepted Haskell code.

import Data.ByteString.Lazy.Char8 as BS hiding  ( map , tail , filter  , null )
import Data.Maybe ( fromJust )
import Text.Printf ( printf )

diffEQ :: Double -> Double -> Double -> Double
diffEQ a b x = a + b * cos x

valEQ :: Double -> Double -> Double -> Double -> Double
valEQ a b c x = a * x   +  b * sin x - c

evalFun :: [ Int ]  -> Double
evalFun [ a' , b' , c' ]  = ret where
          ( a , b , c  ) = ( fromIntegral a' , fromIntegral b' , fromIntegral c' )
          ret = fst . until ( \ ( _ , cnt ) -> cnt  >=  500 )  ( \( x , cnt ) -> ( x - (  valEQ a b c x  /  diffEQ a b x ) , succ cnt ) ) $ (  ( c - b ) / a  , 0   )

readD :: BS.ByteString -> Int
readD = fst . fromJust . BS.readInt

main = BS.interact $ BS.unlines . map (  BS.pack . ( printf "%.6f" :: Double -> String ) . evalFun .  map readD ) . filter ( not.null ) . map BS.words . tail . BS.lines

Changing the initial approximation of \frac{C}{A} and 50 iteration improves the time by 0.60 seconds ( from 0.80 to 0.22 ).

January 4, 2013 Posted by | Haskell, Programming | , , , | Leave a comment

SPOJ Aritho-geometric Series (AGS)

AGS is rather simple problem. Write down couple of terms
a , a + d ,  ar + dr , ar + dr + d , ar^2 + dr^2 + dr . You can see
1. n is even
a*r^{\frac{n-1}{2}} + d*( 1 + r + r^2 \cdots + r^{\frac{n-1}{2}})
2. n is odd
a*r^{\frac{n-1}{2}} + d*( 1 + r + r^2 \cdots + r^{\frac{n-1}{2}}) - d
Accepted Haskell code.

import qualified Data.ByteString.Lazy.Char8 as BS

powM :: Integer -> Integer -> Integer -> Integer
powM a n m   -- a^n mod m
  | n == 0 = 1
  | n == 1 = mod a m
  | even n = ret
  | otherwise = mod ( a * ret ) m
  where
     ret = mod ( powM ( mod ( a ^ 2 ) m ) ( div n 2 ) m ) m

geoSum :: Integer -> Integer -> Integer -> Integer
geoSum r n m
  | n == 0 = mod 1 m
  | n == 1 = mod ( 1 + r ) m
  | odd n = mod ( ( 1 + powM r ( div ( n + 1 ) 2 ) m ) * mod ( geoSum r ( div ( n - 1 ) 2 ) m ) m ) m
  | otherwise = mod (  ( 1 + powM r ( div n 2 ) m ) * mod ( geoSum r ( div ( n - 1 ) 2 ) m ) m  +  powM r n m ) m

solve :: [ Integer ] -> [ BS.ByteString ]
solve [] = []
solve ( a : d : r : n : m : xs )
   | even n = ( BS.pack . show  $ ret )  :  solve xs
   | otherwise = (  BS.pack . show $ ( mod ( ret - d + m ) m ) ) :  solve xs
   where ret  =  mod ( mod ( a * powM r ( div ( n - 1 ) 2 ) m ) m  + mod ( d * geoSum r ( div ( n - 1 ) 2 ) m ) m ) m

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


main =  BS.interact $   BS.unlines . solve . map readD . concat . map BS.words . tail . BS.lines

December 12, 2012 Posted by | Haskell, Programming | , , | 6 Comments

SPOJ DIE HARD

DIEHARD I tried to find out greedy solution but could not come up so trying every possibility and memoization of each state. See more about dynamic programming. The basic idea is that you have to always go to air to increase your chance for survival because its increase the health and armor. We start in air with health increased by 3 and armor increased by 2. After that we have to find out which one, going to water or fire gives maximum survival. A simple Haskell recursive solution which works for small values. We start from funsolve_Air after that we chose the maxiumum value of two functions respectively for water and fire.

funsolve_WF :: Int -> Int -> Int -> Int
funsolve_WF h a cnt
     | h <= 0 || a <= 0 = cnt
     | otherwise = funsolve_Air h a ( cnt + 1 )

funsolve_Air :: Int -> Int -> Int -> Int
funsolve_Air h a cnt = max ( funsolve_WF ( h + 3 - 5 ) ( a + 2 - 10 ) cnt' ) ( funsolve_WF ( h + 3  - 20 )  ( a + 2  + 5 )  cnt' ) where
                      cnt' = cnt + 1

*Main> funsolve_Air 20 8 0
5
*Main> funsolve_Air 2 10 0
1
*Main> funsolve_Air 4 4  0
1

Accepted source code in C++.

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

int memo[1100][1100] ;

int recurse( int h , int a , int cnt , bool flag )
    {
      if ( h <= 0 || a <= 0 ) return cnt ;
      if ( memo[h][a] ) return memo[h][a] ;
      if ( flag ) memo[h][a] = max ( memo[h][a] , recurse ( h + 3 , a + 2 , cnt + 1 , !flag ) ) ;
      else
         memo[h][a] = max ( memo[h][a] ,  max ( recurse ( h - 5 , a - 10 , cnt + 1 , !flag ) , recurse ( h - 20 , a + 5 , cnt + 1 , !flag ) ) ) ;

     return memo[h][a];
   }

int main()
  {
    int n , a , b ;
    scanf( "%d", &n );
    for(int i = 0 ; i < n ; i++)
    {
     memset ( memo , 0 , sizeof memo ) ;
     scanf("%d%d", &a , &b );
     printf("%d\n" , recurse( a , b , -1 ,  1 ));
     if( i != ( n - 1 ) ) printf("\n");
    }

  }

December 11, 2012 Posted by | Programming | , , , , | 3 Comments