My Weblog

Propositional Logic course ( Coursera )

Recently I took a logic course on Coursera. Although I am bit familiar with logic but this course was great. I got more clear grasp on Proposition logic, Proof tree, bit of Prolog and came to know about Quine–McCluskey algorithm ( in my TODO list of implementation ) while learning about circuit minimisation. In one of the assignments, we have to tell if the given formula is tautology, contradiction or contingent so I wrote some quick Haskell code to solve the assignment

import Text.Parsec.Prim
import Text.Parsec.Expr
import Text.Parsec.Char
import Text.Parsec.String ( Parser )
import Control.Applicative hiding ( ( <|> ) , many )
import qualified Data.Map as M
import Data.Maybe
import Data.List
import Control.Monad

data LExpr = Lit Char
| Not LExpr
| And LExpr LExpr
| Or LExpr LExpr
| Imp LExpr LExpr  -- ( P =>  Q )
| Red LExpr LExpr  -- ( P <=  Q )
| Eqi LExpr LExpr  -- ( P <=> Q )
deriving Show

exprCal :: Parser LExpr
exprCal = buildExpressionParser table atom

table = [  [  Prefix ( Not <$string "~" ) ] , [ Infix ( And <$ string  "&"  ) AssocLeft ]
, [  Infix  ( Or  <$string "|" ) AssocLeft ] , [ Infix ( Eqi <$ try ( string "<=>" ) ) AssocLeft
, Infix  ( Imp <$string "=>" ) AssocLeft , Infix ( Red <$  string "<="         ) AssocLeft
]
]

atom :: Parser LExpr
atom =  char '(' *>  exprCal   <* char ')'
<|> ( Lit <$> letter ) assignment :: LExpr -> [ M.Map Char Bool ] assignment expr = map ( M.fromList . zip vs ) ps where vs = variables expr ps = replicateM ( length vs ) [ True, False] variables :: LExpr -> [ Char ] variables expr = map head . group . sort . vars expr$ []  where
vars ( Lit c )    xs = c : xs
vars ( Not expr ) xs = vars expr xs
vars ( And exprf exprs ) xs = vars exprf xs ++ vars exprs xs
vars ( Or  exprf exprs ) xs = vars exprf xs ++ vars exprs xs
vars ( Imp exprf exprs ) xs = vars exprf xs ++ vars exprs xs
vars ( Red exprf exprs ) xs = vars exprf xs ++ vars exprs xs
vars ( Eqi exprf exprs ) xs = vars exprf xs ++ vars exprs xs

expEval :: LExpr -> M.Map Char  Bool -> Bool
expEval ( Lit v   )         mp =  fromMaybe False ( M.lookup v mp )
expEval ( Not expr  )       mp =  not . expEval  expr $mp expEval ( And exprf exprs ) mp = expEval exprf mp && expEval exprs mp expEval ( Or exprf exprs ) mp = expEval exprf mp || expEval exprs mp expEval ( Imp exprf exprs ) mp = ( not . expEval exprf$ mp  ) || expEval exprs mp
expEval ( Red exprf exprs ) mp =  expEval ( Imp exprs exprf  ) mp
expEval ( Eqi exprf exprs ) mp =  expEval exprf mp == expEval exprs mp

values :: LExpr -> [ Bool ]
values expr = map ( expEval expr ) ( assignment expr )

isTautology :: LExpr -> Bool
isTautology = and . values

isContradiction :: LExpr -> Bool
isContradiction  = not . or . values

isContingent :: LExpr -> Bool
isContingent expr = not ( isTautology expr || isContradiction expr )

calculator :: String -> LExpr
calculator expr = case parse  exprCal ""  expr of
Left msg -> error "failed to parse"
Right val -> val

Mukeshs-MacBook-Pro:Compilers mukeshtiwari$ghci LogicPraser.hs GHCi, version 7.8.1: http://www.haskell.org/ghc/ 😕 for help Loading package ghc-prim ... linking ... done. Loading package integer-gmp ... linking ... done. Loading package base ... linking ... done. [1 of 1] Compiling Main ( LogicPraser.hs, interpreted ) Ok, modules loaded: Main. *Main> calculator "p=>q|q=>p" Loading package transformers-0.3.0.0 ... linking ... done. Loading package array-0.5.0.0 ... linking ... done. Loading package deepseq-1.3.0.2 ... linking ... done. Loading package containers-0.5.5.1 ... linking ... done. Loading package bytestring-0.10.4.0 ... linking ... done. Loading package mtl-2.1.3.1 ... linking ... done. Loading package text-1.1.1.1 ... linking ... done. Loading package parsec-3.1.5 ... linking ... done. Imp (Imp (Lit 'p') (Or (Lit 'q') (Lit 'q'))) (Lit 'p') *Main> calculator "(p=>q)|(q=>p)" Or (Imp (Lit 'p') (Lit 'q')) (Imp (Lit 'q') (Lit 'p')) *Main> isTautology ( Or (Imp (Lit 'p') (Lit 'q')) (Imp (Lit 'q') (Lit 'p')) ) True *Main> calculator "(p=>q)&(q=>p)" And (Imp (Lit 'p') (Lit 'q')) (Imp (Lit 'q') (Lit 'p')) *Main> isTautology ( And (Imp (Lit 'p') (Lit 'q')) (Imp (Lit 'q') (Lit 'p')) ) False *Main> isCont ( And (Imp (Lit 'p') (Lit 'q')) (Imp (Lit 'q') (Lit 'p')) ) isContingent isContradiction *Main> isContingent ( And (Imp (Lit 'p') (Lit 'q')) (Imp (Lit 'q') (Lit 'p')) ) True *Main> isCont ( And (Imp (Lit 'p') (Lit 'q')) (Imp (Lit 'q') (Lit 'p')) ) isContingent isContradiction *Main> isContradiction ( And (Imp (Lit 'p') (Lit 'q')) (Imp (Lit 'q') (Lit 'p')) ) False Proof tree for propositional logic was something which I learned from this course. Why proof tree method is great ? We will take a example to understand it more clearly. $p \Rightarrow ( q \vee r ) \newline s \Rightarrow ( \neg r) \newline \line(1,0){70}\newline ( p\wedge s ) \Rightarrow q$ A set of premises Δ logically entails a conclusion ϕ (written as Δ |= ϕ) if and only if every interpretation that satisfies the premises also satisfies the conclusion. Logical formulas above the line are premises and below is conclusion. $p \hspace{5 mm} q \hspace{5 mm} r \hspace{5 mm} s \hspace{5 mm} p \Rightarrow ( q \vee r ) \hspace{5 mm} s \Rightarrow ( \neg r) \hspace{5 mm}( p\wedge s ) \Rightarrow q$ $\begin{tabular}{ | l | l | l | l | p{2cm} | p{2cm} | p{2cm} | } \hline 0 & 0 & 0 & 0 & 1 & 1 & 1 \\ 0 & 0 & 0 & 1 & 1 & 1 & 1 \\ 0 & 0 & 1 & 0 & 1 & 1 & 1 \\ 0 & 0 & 1 & 1 & 1 & 0 & 1 \\ 0 & 1 & 0 & 0 & 1 & 1 & 1 \\ 0 & 1 & 0 & 1 & 1 & 1 & 1 \\ 0 & 1 & 1 & 0 & 1 & 1 & 1 \\ 0 & 1 & 1 & 1 & 1 & 0 & 1 \\ 1 & 0 & 0 & 0 & 0 & 1 & 1 \\ 1 & 0 & 0 & 1 & 0 & 1 & 0 \\ 1 & 0 & 1 & 0 & 1 & 1 & 1 \\ 1 & 0 & 1 & 1 & 1 & 0 & 0 \\ 1 & 1 & 0 & 0 & 1 & 1 & 1 \\ 1 & 1 & 0 & 1 & 1 & 1 & 1 \\ 1 & 1 & 1 & 0 & 1 & 1 & 1 \\ 1 & 1 & 1 & 1 & 1 & 0 & 1 \\ \hline \end{tabular}$ We can see from the truth table that all the interpretation which satisfy premises also satisfy conclusion. Looks great but lengthy. Proof tree rules $A \wedge B \hspace*{10mm} A \vee B \hspace*{10mm} A \Rightarrow B\hspace*{10mm} A \equiv B \hspace*{10mm} A \newline \hspace*{5mm}| \hspace*{15mm} / \hspace*{8mm} \backslash \hspace*{10mm} / \hspace*{8mm}\backslash \hspace*{10mm} / \hspace*{8mm} \backslash \hspace*{7mm} \neg A\newline \hspace*{4mm} A \hspace*{13mm} A \hspace*{8mm} B \hspace*{5mm} \neg A \hspace*{8mm}B \hspace*{6mm} A \hspace*{10mm} \neg A \hspace*{5mm} \times \newline \hspace*{4mm} B \hspace*{55mm} B\hspace*{10mm} \neg B\newline$ $\neg ( A \wedge B ) \hspace*{5mm} \neg ( A \vee B) \hspace*{5mm} \neg ( A \Rightarrow B) \hspace*{5mm} \neg ( A \equiv B) \hspace*{5mm} \neg \neg A \newline \hspace*{4mm} / \hspace*{8mm} \backslash \hspace*{15mm} | \hspace*{23mm} | \hspace*{15mm} /\hspace*{8mm} \backslash \hspace*{15mm} | \newline \hspace*{2mm} \neg A \hspace*{6mm} \neg B \hspace*{8mm} \neg A \hspace*{20mm} A \hspace*{13mm} A\hspace*{8mm} \neg A \hspace*{10mm} A \newline \hspace*{28mm} \neg B \hspace*{18mm} \neg B \hspace*{11mm} \neg B\hspace*{8mm} B \newline$ Theorem: Δ |= ϕ if and only if Δ ∪ {¬ϕ} is unsatisfiable. Following the theorem, our problem translates into proving that $p \Rightarrow ( q \vee r ) \newline s \Rightarrow ( \neg r) \newline \neg ( ( p\wedge s ) \Rightarrow q )$ is unsatisfiable. We can see that tree is closed and there is no interpretation which makes it satisfiable. Chose any branch and start moving up the tree and if it contains both A and ~A then we close the branch because it’s not possible to make A and ~A true simultaneously. You can see the course notes. . See Propositional logic library from hackage. Source code on github Advertisements May 2, 2014 While Interpreter Now a days I am trying to learn about static program analysis and started reading the Principal of Program Analysis. In order to understand the concepts, the book introduces a small programming language While. I thought about writing interpreter for While to improve my Haskell skills. The very first task is in compilation/interpretation is breaking the source code into keywords, identifiers, operators, numbers and other tokens, known as lexical analysis. These tokens are passed to syntax analysis phase to combine the tokens into well formed expressions, statements and programs according to the grammar. The output of syntax analysis is abstract syntax tree which gives structural representation of of the input. Once we have abstract syntax tree, we can interpret or manipulate it for code generation. If you are interested in code generation then see Stephen Dielh blog. The grammar for while program is Grammar for expression a ::= x | n | - a | a opa a b ::= true | false | not b | b opb b | a opr a opa ::= + | - | * | / opb ::= and | or opr ::= > | < Grammar for statements S ::= x := a | skip | S1; S2 | ( S ) | if b then S1 else S2 | while b do S The very first task is to define abstract syntax tree and Haskell algebraic data type makes this task almost trivial. module AST ( Opa (..), Opb (..), Opr (..), AExpr (.. ) , BExpr ( .. ), Stmt ( .. ) ) where data Opa = Add | Sub | Mul | Div deriving Show data Opb = And | Or deriving Show data Opr = Greater | Less deriving Show data AExpr = Var String | Num Integer | Neg AExpr | ABin Opa AExpr AExpr deriving Show data BExpr = Con Bool | Not BExpr | BBin Opb BExpr BExpr | AL Opr AExpr AExpr deriving Show data Stmt = List [ Stmt ] | Assing AExpr AExpr | If BExpr Stmt Stmt | While BExpr Stmt | Skip deriving Show The next task is design the lexical analyzer. We will use Parsec for this purpose. module Lexer where import Text.Parsec import qualified Text.Parsec.Token as T import Text.Parsec.Language ( emptyDef ) import Text.Parsec.String ( Parser ) lexer :: T.TokenParser () lexer = T.makeTokenParser emptyDef { T.commentStart = "{-" , T.commentEnd = "-}" , T.reservedOpNames = [ "+", "-", "*", "/", ":=", ">", "<", "not", "and", "or" ] , T.reservedNames = [ "if", "then", "else", "while", "do", "skip", "true", "false" ] } identifier :: Parser String identifier = T.identifier lexer whiteSpace :: Parser () whiteSpace = T.whiteSpace lexer reserved :: String -> Parser () reserved = T.reserved lexer reservedOp :: String -> Parser () reservedOp = T.reservedOp lexer parens :: Parser a -> Parser a parens = T.parens lexer integer :: Parser Integer integer = T.integer lexer semi :: Parser String semi = T.semi lexer semiSep :: Parser a -> Parser [ a ] semiSep = T.semiSep lexer Now lexical analyzer will assist the parser for creating the abstract syntax tree. A while program is list of statements and statement consists of expressions. We can build parse the expression by Parsec expression builder by giving the table of operators and associativity. module Parser ( whileParser ) where import Text.Parsec import Text.Parsec.Expr import Text.Parsec.String ( Parser ) import Control.Applicative hiding ( (<|>) ) import Lexer import AST aTable = [ [ Prefix ( Neg <$ reservedOp "-" )                ]
, [   Infix  ( ABin Mul <$reservedOp "*" ) AssocLeft , Infix ( ABin Div <$ reservedOp "/" ) AssocLeft ]
, [   Infix  ( ABin Add <$reservedOp "+" ) AssocLeft , Infix ( ABin Sub <$ reservedOp "-" ) AssocLeft ]
]

bTable = [  [  Prefix ( Not  <$reservedOp "not" ) ] , [ Infix ( BBin And <$ reservedOp "and" ) AssocLeft ]
, [  Infix  ( BBin Or  <$reservedOp "or" ) AssocLeft ] ] aExpression :: Parser AExpr aExpression = buildExpressionParser aTable aTerm where aTerm = parens aExpression <|> Var <$> identifier
<|> Num <$> integer bExpression :: Parser BExpr bExpression = buildExpressionParser bTable bTerm where bTerm = parens bExpression <|> ( Con True <$ reserved "true"  )
<|> (  Con False  <$reserved "false" ) <|> try ( AL Greater <$>  ( aExpression  <* reserved ">" )
<*>    aExpression )
<|> (  AL Less    <$> ( aExpression <* reserved "<" ) <*> aExpression ) whileParser :: Parser Stmt whileParser = whiteSpace *> stmtParser <* eof where stmtParser :: Parser Stmt stmtParser = parens stmtParser <|> List <$> sepBy stmtOne semi
stmtOne :: Parser Stmt
stmtOne =  ( Assing <$> ( Var <$> identifier )
<*> ( reserved ":=" *> aExpression ) )
<|> ( If <$> ( reserved "if" *> bExpression <* reserved "then" ) <*> stmtParser <*> ( reserved "else" *> stmtParser ) ) <|> ( While <$> ( reserved "while" *> bExpression <*  reserved "do" )
<*>   stmtParser )
<|> ( Skip <$reserved "nop" ) We have abstract syntax tree so we can interpret our program. You can think of a program as collection of commands which manipulates some memory location. module Interpreter ( evalProgram ) where import qualified Data.Map as M import AST type Store = M.Map String Integer evalA :: AExpr -> Store -> Integer evalA ( Var v ) s = M.findWithDefault 0 v s evalA ( Num n ) _ = n evalA ( Neg e ) s = - ( evalA e s ) evalA ( ABin Add e1 e2 ) s = evalA e1 s + evalA e2 s evalA ( ABin Sub e1 e2 ) s = evalA e1 s - evalA e2 s evalA ( ABin Mul e1 e2 ) s = evalA e1 s * evalA e2 s evalA ( ABin Div e1 e2 ) s = div ( evalA e1 s ) ( evalA e2 s ) evalB :: BExpr -> Store -> Bool evalB ( Con b ) _ = b evalB ( Not e ) s = not ( evalB e s ) evalB ( BBin And e1 e2 ) s = ( && ) ( evalB e1 s ) ( evalB e2 s ) evalB ( BBin Or e1 e2 ) s = ( || ) ( evalB e1 s ) ( evalB e2 s ) evalB ( AL Greater e1 e2 ) s = ( evalA e1 s ) > ( evalA e2 s ) evalB ( AL Less e1 e2 ) s = ( evalA e1 s ) < ( evalA e2 s ) interpret :: Stmt -> Store -> Store interpret ( Assing ( Var v ) expr ) s = M.insert v ( evalA expr s ) s interpret ( List [] ) s = s interpret ( List ( x : xs ) ) s = interpret ( List xs ) ( interpret x s ) interpret ( If e st1 st2 ) s | evalB e s = interpret st1 s | otherwise = interpret st2 s interpret ( While e st ) s | not t = s | otherwise = interpret ( While e st ) w where t = evalB e s w = interpret st s evalProgram :: Stmt -> Store evalProgram st = interpret st M.empty Now every thing is complete so let’s add main and write some while program. module Main where import System.Environment import Text.Parsec import Parser import Interpreter main = do ( file : _ ) <- getArgs program <- readFile file case parse whileParser "" program of Left e -> print e >> fail "Parse Error" Right ast -> print ( evalProgram ast ) While program to compute the greatest common divisor Mukeshs-MacBook-Pro:whileinterpreter mukeshtiwari$ cat Gcd.while
a := 10 ;
b := 100  ;
while ( b > 0 ) do
(
t := b ;
b := a - ( a / b ) * b ;
a := t
)

Factorial program

Mukeshs-MacBook-Pro:whileinterpreter mukeshtiwari$cat Fact.while x := 10 ; y := x ; z := 1 ; while ( y > 1 ) do ( z := z * y ; y := y - 1 ); y := 0 Since our language doesn’t have IO instruction so we will have to see which variable store the result. In gcd program, the variable t stores the result and variable z stores the factorial of number. Mukeshs-MacBook-Pro:src mukeshtiwari$ ghc -fforce-recomp Main.hs
[1 of 5] Compiling AST              ( AST.hs, AST.o )
[2 of 5] Compiling Lexer            ( Lexer.hs, Lexer.o )
[3 of 5] Compiling Interpreter      ( Interpreter.hs, Interpreter.o )
[4 of 5] Compiling Parser           ( Parser.hs, Parser.o )
[5 of 5] Compiling Main             ( Main.hs, Main.o )
Linking Main ...
Mukeshs-MacBook-Pro:src mukeshtiwari$./Main ../Gcd.while fromList [("a",10),("b",0),("t",10)] Mukeshs-MacBook-Pro:src mukeshtiwari$ ./Main ../Fact.while
fromList [("x",10),("y",0),("z",3628800)]

I am not expert in either Haskell or compiler so if you have any comments then please let me know. The complete source code on github.

January 30, 2014

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:

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.

September 21, 2013

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 ZeroMQ for distributed computing. This post is influenced by ØMQ – The Guide By Pieter Hintjens and translation of codes in Haskell. I suggest you to read The Guide By Pieter Hintjens and if you are interested in Haskell code then you see these codes. Request-Reply The client sends “Accept hello from Client. ” to the server, which replies with “I got you.”. This Haskell server opens a ØMQ socket on port 5555, reads requests on it, and replies with “I got you” to each request.The REQ-REP socket pair is in lockstep. The client issues send and then receive, in a loop (or once if that’s all it needs). Doing any other sequence (e.g., sending two messages in a row) will result in a return code of -1 from the send or receive call. import System.ZMQ3 import Control.Monad import qualified Data.ByteString.Char8 as BS import Data.ByteString.Lazy.Internal as BL import Data.IORef import Control.Concurrent ( threadDelay ) main = do c <- context s <- socket c Rep bind s "tcp://127.0.0.1:5555" counter <- newIORef 0 forever$ do
t <- readIORef counter
res <- receive s
print res
send' s [] ( BL.packChars $"I got you. " ++ show t ) modifyIORef counter ( +1 ) threadDelay 10000 close s destroy c return () Client import System.ZMQ3 import Control.Monad import qualified Data.ByteString.Char8 as BS import qualified Data.ByteString.Lazy.Internal as BL import Data.IORef import Control.Concurrent ( threadDelay ) main = do c <- context s <- socket c Req connect s "tcp://127.0.0.1:5555" counter <- newIORef 0 forever$ do
t <- readIORef counter
send' s [] ( BL.packChars $"Accept hello from Client. " ++ show t ) msg <- receive s print msg modifyIORef counter ( +1 ) threadDelay 10000 return () Running a server with two clients. Mukeshs-MacBook-Pro:ZMQ mukeshtiwari$ ./ZeroMqServer
"Accept hello from Client. 0"
"Accept hello from Client. 1"
"Accept hello from Client. 2"
"Accept hello from Client. 3"
"Accept hello from Client. 4"
"Accept hello from Client. 5"
"Accept hello from Client. 6"
"Accept hello from Client. 7"
"Accept hello from Client. 8"
"Accept hello from Client. 9"
"Accept hello from Client. 10"
"Accept hello from Client. 11"
"Accept hello from Client. 12"
"Accept hello from Client. 13"
"Accept hello from Client. 14"
"Accept hello from Client. 15"

"Accept hello from Client. 0"
"Accept hello from Client. 228"
"Accept hello from Client. 1"
"Accept hello from Client. 229"
"Accept hello from Client. 2"
"Accept hello from Client. 230"
"Accept hello from Client. 3"
"Accept hello from Client. 231"
"Accept hello from Client. 4"
"Accept hello from Client. 232"

First client.
Mukeshs-MacBook-Pro:ZMQ mukeshtiwari$./ZeroMqClient "I got you. 0" "I got you. 1" "I got you. 2" "I got you. 3" "I got you. 4" "I got you. 5" "I got you. 6" "I got you. 7" "I got you. 8" "I got you. 9" "I got you. 10" "I got you. 11" "I got you. 12" Second client Mukeshs-MacBook-Pro:ZMQ mukeshtiwari$ ./ZeroMqClient
"I got you. 228"
"I got you. 230"
"I got you. 232"
"I got you. 234"
"I got you. 236"
"I got you. 238"
"I got you. 240"
"I got you. 242"
"I got you. 244"
"I got you. 246"
"I got you. 248"
"I got you. 250"
"I got you. 252"
"I got you. 254"
"I got you. 256"
"I got you. 258"
"I got you. 260"
"I got you. 262"
"I got you. 264"
"I got you. 266"

Publish-Subscribe

Data publishing server which publishes weather data for zip codes in range 500 and 2000.

import System.ZMQ3
import Control.Monad
import qualified Data.ByteString.Char8 as BS hiding ( putStrLn )
import qualified Data.ByteString.Lazy.Internal as BL
import Control.Concurrent ( threadDelay )
import System.Random

main = do
c <- context
publisher <- socket c Pub
bind publisher "tcp://127.0.0.1:5556"
bind publisher "ipc://weather.ipc"
forever $do zipcode <- randomRIO ( ( 500 , 2000 ) :: ( Int , Int ) ) temp <- randomRIO ( 10 , 45 ) :: IO Int relhum <- randomRIO ( 0 , 100 ) :: IO Int putStrLn$ show zipcode ++ " " ++ show temp ++ " " ++ show relhum
send' publisher [] ( BL.packChars $show zipcode ++ " " ++ show temp ++ " " ++ show relhum ) close publisher destroy c return () Client who is only interested in two zip codes. import System.ZMQ3 import Control.Monad import qualified Data.ByteString.Char8 as BS import qualified Data.ByteString.Lazy.Internal as BL import Control.Concurrent ( threadDelay ) import System.Random main = do c <- context subscriber <- socket c Sub connect subscriber "tcp://127.0.0.1:5556" subscribe subscriber ( BS.pack "1000" ) subscribe subscriber ( BS.pack "1010" ) forever$ do
update <- receive subscriber
print update

close subscriber
destroy c

Our server is publishing lot of data but client is only interested in two zip codes.

Mukeshs-MacBook-Pro:ZMQ mukeshtiwari$./ZeroMqWeatherPubServer 1568 27 85 924 41 46 1461 15 46 1867 44 28 1013 23 100 1052 13 6 1720 45 6 1480 12 94 1852 33 4 1295 20 6 925 18 77 935 37 94 1670 11 6 1285 39 38 1613 44 99 1888 26 62 1011 21 45 993 45 45 1402 26 86 1639 13 65 1285 18 40 1960 38 18 1160 27 39 1374 16 59 665 25 22 Mukeshs-MacBook-Pro:ZMQ mukeshtiwari$ ./ZeroMqWeatherClient
"1000 34 89"
"1010 19 70"
"1010 15 86"
"1000 38 57"
"1010 12 42"
"1000 25 1"
"1000 28 78"
"1000 25 16"
"1000 28 82"
"1010 28 98"
"1000 12 77"
"1010 11 16"
"1010 44 14"
"1010 12 89"
"1000 32 8"
"1010 37 55"
"1010 17 21"
"1000 13 96"
"1000 18 51"
"1010 16 38"
"1000 18 21"
"1010 37 60"
"1010 17 25"
"1000 33 43"
"1000 34 44"
"1010 33 78"
"1010 35 63"
"1000 39 50"
"1000 45 70"
"1000 26 3"
"1010 34 12"
"1010 26 3"
"1000 32 75"
"1010 14 68"
"1000 44 75"
"1010 27 54"
"1000 21 39"
"1010 12 65"
"1010 43 29"
"1010 25 60"

Divide and Conquer

For this problem we will calculate the number of primes less than $10^{7}$. Our ventilator will push the task to workers and they will perform the task. After finishing the job, they will send the result back to sink.

Ventilator

import System.ZMQ3
import Control.Monad
import qualified Data.ByteString.Char8 as BS
import Data.ByteString.Lazy.Internal as BL
import Data.IORef
import Control.Concurrent ( threadDelay )

main = do
c <- context
sender <- socket c Push
bind sender "tcp://127.0.0.1:5557"

sink <- socket c Push
connect sink "tcp://127.0.0.1:5558"

putStrLn "Press enter when workers are ready."
_ <- getChar
putStrLn "Sending the task to workers.\n"

send' sink [] ( BL.packChars . show $0 ) forM_ [ 0..9999 ]$ \ i -> do
putStrLn $"Sending the range [ " ++ show ( 1000 * i + 1 ) ++ " .. " ++ show ( 1000 * i + 1 + 999 ) ++ "] to worker for primality testing." send' sender [] ( BL.packChars . show$ 1000 * i + 1 )

close sink
close sender
destroy c

Worker for computing the prime number.

import System.ZMQ3
import Control.Monad
import qualified Data.ByteString.Char8 as BS
import qualified Data.ByteString.Lazy.Internal as BL
import Data.IORef
import Control.Concurrent ( threadDelay )
import Data.Bits

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 ]

main = do
c <- context
receiver <- socket c Pull
connect receiver "tcp://127.0.0.1:5557"

sender <- socket c Push
connect sender "tcp://127.0.0.1:5558"

forever $do num <- receive receiver let n = read . BS.unpack$ num  :: Integer
let len = length . filter ( == True ) . primeRange  n  $n + 999 putStrLn$ "Received range [ " ++  show n ++ " .. " ++ show ( n + 999 ) ++ "  and number of primes in this range is " ++ show len
send' sender [] ( BL.packChars . show $len ) close receiver close sender destroy c Sink for collecting the results. import System.ZMQ3 import Control.Monad import qualified Data.ByteString.Char8 as BS import qualified Data.ByteString.Lazy.Internal as BL import Data.IORef import Control.Concurrent ( threadDelay ) main = do c <- context receiver <- socket c Pull bind receiver "tcp://127.0.0.1:5558" m <- receive receiver print m sum <- newIORef 0 forM_ [ 1.. 10000 ]$ \i -> do
num <- receive receiver
let n = read . BS.unpack $num :: Integer putStrLn$ "Received a number " ++ show n ++ " from one of worker"
modifyIORef sum ( + n )

t <- readIORef sum
putStrLn $"Total number of primes in the range is " ++ show t close receiver destroy c Running Ventilator, 3 workers and sink Mukeshs-MacBook-Pro:ZMQ mukeshtiwari$ ./TaskVent
Sending the range [ 9992001 .. 9993000] to worker for primality testing.
Sending the range [ 9993001 .. 9994000] to worker for primality testing.
Sending the range [ 9994001 .. 9995000] to worker for primality testing.
Sending the range [ 9995001 .. 9996000] to worker for primality testing.
Sending the range [ 9996001 .. 9997000] to worker for primality testing.
Sending the range [ 9997001 .. 9998000] to worker for primality testing.
Sending the range [ 9998001 .. 9999000] to worker for primality testing.
Sending the range [ 9999001 .. 10000000] to worker for primality testing.

Task worker - 1
Mukeshs-MacBook-Pro:ZMQ mukeshtiwari$./TaskWorkder Received range [ 9978001 .. 9979000 and number of primes in this range is 60 Received range [ 9981001 .. 9982000 and number of primes in this range is 49 Received range [ 9984001 .. 9985000 and number of primes in this range is 69 Received range [ 9987001 .. 9988000 and number of primes in this range is 57 Received range [ 9990001 .. 9991000 and number of primes in this range is 60 Received range [ 9993001 .. 9994000 and number of primes in this range is 71 Received range [ 9996001 .. 9997000 and number of primes in this range is 58 Received range [ 9999001 .. 10000000 and number of primes in this range is 53 Task worker - 2 Mukeshs-MacBook-Pro:ZMQ mukeshtiwari$ ./TaskWorkder
Received range [ 9982001 .. 9983000  and number of primes in this range is 69
Received range [ 9985001 .. 9986000  and number of primes in this range is 62
Received range [ 9988001 .. 9989000  and number of primes in this range is 58
Received range [ 9991001 .. 9992000  and number of primes in this range is 57
Received range [ 9994001 .. 9995000  and number of primes in this range is 64
Received range [ 9997001 .. 9998000  and number of primes in this range is 67
Task Worker - 3
Received range [ 9986001 .. 9987000  and number of primes in this range is 59
Received range [ 9989001 .. 9990000  and number of primes in this range is 58
Received range [ 9992001 .. 9993000  and number of primes in this range is 58
Received range [ 9995001 .. 9996000  and number of primes in this range is 62
Received range [ 9998001 .. 9999000  and number of primes in this range is 64

Sink
Received a number 53 from one of worker
Received a number 67 from one of worker
Received a number 64 from one of worker
Received a number 52 from one of worker
Received a number 59 from one of worker
Received a number 58 from one of worker
Received a number 58 from one of worker
Received a number 62 from one of worker
Received a number 64 from one of worker
Total number of primes in the range is 664579

First start all the workers and sink and then run ventilator. See more about distributed computing in Haskell Eden and distributed-process. If you have any suggestion then please let me know. Some of the contents and images are taken from Pieter Hintjens‘s tutorial by his permission.

May 13, 2013

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 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 Parsing Email ID Now a days I am exploring one of the most awesome feature of Haskell, parsing. There are lot of parsing libraries but Parsec is the popular one. This parsing code is written for SPOJ EMAIL ID problem ( Unfortunately it’s getting time limit exceed. I have seen couple of python solution accepted so I am hopeful that there must be another algorithm probably using regular expressions to solve the problem ) but you can build more sophisticated email-id parser by adding more functionality. Tony Morris excellent parsing tutorial is must read. import Data.List import qualified Text.Parsec.ByteString as PB import Text.Parsec.Prim import Text.Parsec.Char import Text.Parsec.Combinator import qualified Data.ByteString.Char8 as BS import Control.Applicative hiding ( ( <|> ) , many ) validChars :: PB.Parser Char validChars = alphaNum <|> oneOf "._" dontCare :: PB.Parser Char dontCare = oneOf "~!@#$%^&*()<>?,."

{--
emailAddress :: PB.Parser  String
emailAddress = do
_ <- many dontCare
fi <- alphaNum
se <- validChars
th <- validChars
fo <- validChars
ft <- validChars
restAddr <- many validChars
let addr = fi : se : th : fo : ft : restAddr
char '@'
dom <- many1 alphaNum
rest <- try ( string ".com" <|> string ".org"
<|>  string ".edu" ) <|> try ( string ".co.in" )
_ <- many dontCare
return $addr ++ ( '@': dom ++ rest ) --} emailAddress :: PB.Parser String emailAddress = conCatfun <$> ( many dontCare *> alphaNum ) <*> validChars <*>
validChars <*> validChars <*> validChars <*> many alphaNum  <*>
( char '@' *> many1 alphaNum ) <*> ( try ( string ".com" <|>
string ".org" <|>  string ".edu" ) <|> try ( string ".co.in" )
<* many dontCare ) where
conCatfun a b c d e f dom rest =
( a : b : c : d : e : f ) ++ ( '@' : dom ) ++ rest

collectEmail :: BS.ByteString -> String
collectEmail email = case  parse emailAddress "" email of
Right addr -> addr
Left err ->  ""

process :: ( Int , [ String ] ) -> BS.ByteString
process ( k , xs ) = ( BS.pack "Case " ) BS.append ( BS.pack . show $k ) BS.append ( BS.pack ": " ) BS.append ( BS.pack . show . length$ xs )
BS.append ( BS.pack "\n" ) BS.append ( BS.pack
(  unlines . filter ( not . null ) $xs ) ) main = BS.interact$ BS.concat .  map process . zip [ 1.. ] . map ( map collectEmail .
BS.words ) . tail . BS.lines

Mukeshs-MacBook-Pro:Haskell mukeshtiwari$cat t.txt 2 svm11@gmail.com svm11@gmail.com svm12@gmail.co.in ~!@#$%^&*()<>?svm12@gmail.co.in~!@#$%^&*() Mukeshs-MacBook-Pro:Haskell mukeshtiwari$ ./Spoj_11105 < t.txt
Case 1: 1
svm11@gmail.com
Case 2: 2
svm11@gmail.com
svm12@gmail.co.in
svm12@gmail.co.in

Update

I tried again to solve this problem using regular expression. Following the tutorial, I wrote this code but I got compiler error because of old version of ghc on SPOJ ( GHC-6.10.4 ). It’s working fine on my system but I still have to test if it is correct and fast enough to get accepted.

import Data.List
import Text.Regex.Posix
import qualified Data.ByteString.Char8 as BS

pat :: BS.ByteString
pat = BS.pack "[^~!@#$%^&*()<>?,.]*[a-zA-Z0-9][a-zA-Z0-9._][a-zA-Z0-9._][a-zA-Z0-9._][a-zA-Z0-9._][a-zA-Z0-9._]*@[a-zA-Z0-9]+.(com|edu|org|co.in)[^~!@#$%^&*()<>?,.a-zA-Z0-9]*"

collectEmail :: BS.ByteString -> BS.ByteString
collectEmail email = ( =~ ) email pat

process :: ( Int , [ BS.ByteString ] ) -> BS.ByteString
process ( k , xs ) = ( BS.pack "Case " ) BS.append ( BS.pack . show $k ) BS.append ( BS.pack ": " ) BS.append ( BS.pack . show . length$ xs )
BS.append ( BS.pack "\n" ) BS.append ( BS.unlines xs )

main = BS.interact $BS.concat . map process . zip [ 1 .. ] . map ( filter ( not . BS.null ) . map collectEmail . BS.words ) . tail . BS.lines Mukeshs-MacBook-Pro:SPOJ mukeshtiwari$ cat t.txt
2
svm11@gmail.com
svm11@gmail.com svm12@gmail.co.in %&^%&%&%&%&^%&^%&^%&^mukeshtiwari.iiitm@gmail.com%&%^&^%&%&%&%&^%&%&^%&^%&^% %$%$#%#%#%#%$#%&%&%&tiwa@gmail.com Mukeshs-MacBook-Pro:SPOJ mukeshtiwari$ ./Spoj_11105 < t.txt
Case 1: 1
svm11@gmail.com
Case 2: 3
svm11@gmail.com
svm12@gmail.co.in
mukeshtiwari.iiitm@gmail.com

March 8, 2013

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

Parsing with Applicative Functors

While reading Write Yourself a Scheme in 48 Hours, I thought of converting the monadic codes of parsing into applicative style. Learn you a Haskell for Great Good has excellent explanation about functors and applicative functors. Real World Haskell has excellent chapter on using applicative functors for parsing. Below is my attempt to convert the monadic codes of parsing into applicative style for Scheme parsing. On parsing note, see the excellent tutorial of Tony Morris and Albert Y. C. Lie.

import Control.Applicative hiding ( many , ( <|> )  )
import Text.ParserCombinators.Parsec hiding ( spaces )

symbol :: Parser Char
symbol =  oneOf "!$%&|+-*/:<=?>@^_~#" spaces :: Parser () spaces = skipMany space parseString :: Parser LispVal parseString = String <$> ( char '"' *>  x <*  char '"' ) where
x = many ( noneOf "\"" )
--This is not completely applicative style because I have to take the result for testing conditions.
parseAtom :: Parser LispVal
parseAtom = do
atom <-  ( : ) <$> ( letter <|> symbol ) <*> many ( letter <|> digit <|> symbol ) return$ case atom of
"#t" -> Bool True
"#f" -> Bool False
_ -> Atom atom

parseNumber :: Parser LispVal
parseNumber = Number . read <$> many1 digit parseList :: Parser LispVal parseList = List <$> sepBy parseExpr space

parseDottedList :: Parser LispVal
parseDottedList = DottedList <$> ( endBy parseExpr space ) <*> ( char '.' *> spaces *> parseExpr ) parseQuoted :: Parser LispVal parseQuoted = f <$>   ( char '\'' *> parseExpr ) where
f x =   List [ Atom "quote" , x ]

parseExpr :: Parser LispVal
parseExpr =  parseAtom
<|> parseString
<|> parseNumber
<|> parseQuoted
<|> ( char '(' *> spaces *> ( ( try parseList ) <|> parseDottedList ) <* spaces  <* char ')' )

Complete code on github ( Apart from applicative parsing every other code is from Scheme parsing tutorial ). If you have any comments or suggestion then please let me know.

January 18, 2013