# My Weblog

## 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.

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 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 )
close s
destroy c
return ()



Client

import System.ZMQ3
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 )
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 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
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 sender
destroy c



Sink for collecting the results.

import System.ZMQ3
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

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
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 ## 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 )
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 ## Improved Elliptic Curve prime factorisation Currently i am trying to understand monad. I wrote elliptic curve prime factorisation some times ago and now using monads, i wrote a bit improved elliptic curve prime factorisation . This code factorise the product of prime numbers in range [1000..2000] withing less than second on my slow computer.Algorithm explained is previous post. Modified Haskell code. --y^2= x^3+ax+b mod n import Random import Control.Monad import Data.List import Data.Bits data Elliptic = Conelliptic Integer Integer Integer deriving ( Eq , Show ) data Point = Conpoint Integer Integer | Identity deriving ( Eq , Show ) powM :: Integer -> Integer -> Integer -> IO Integer powM a d n = return$ 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 -> IO ( Integer , Integer )
calSd n = return ( s , d ) where
s = until ( \x -> testBit ( n - 1) ( fromIntegral x ) )  ( +1 ) 0
d = div ( n - 1 ) (  shiftL 1 ( fromIntegral s )  )

isProbable::Integer->IO Bool
isProbable n
| n <= 1 = return False
| n == 2 = return True
| even n = return False
| otherwise	= calSd n >>= (\( s , d ) -> rabinMiller 0 n s d )

rabinMiller::Integer->Integer->Integer->Integer->IO Bool
rabinMiller cnt n s d
| cnt>=5 = return True
| otherwise =
randomRIO ( 2 , n - 2 ) >>=
(\a ->  powM a d n >>=
(\x -> case x == 1 of
True  -> rabinMiller ( cnt + 1 ) n s d
_   -> if any ( == pred n ) . take ( fromIntegral s ) . iterate (\e -> mod ( e^2 ) n ) $x then rabinMiller ( cnt + 1) n s d else return False ) ) --add points of elliptic curve addPoints::Elliptic->Point->Point-> Either Point Integer addPoints _ Identity p_2 = Left p_2 addPoints _ p_1 Identity = Left p_1 addPoints ( Conelliptic a b n ) ( Conpoint x_p y_p ) ( Conpoint x_q y_q ) | x_p /= x_q = let [ u , v , d ] = extended_gcd (x_p-x_q) n s = mod ( ( y_p - y_q ) * u ) n x_r = mod ( s * s - x_p - x_q ) n y_r= mod ( -y_p - s * ( x_r - x_p ) ) n in case ( ( Conpoint x_r y_r ) , d ) of ( _ , 1 ) -> Left ( Conpoint x_r y_r ) ( _ , d' ) -> Right d' | otherwise = if mod ( y_p + y_q ) n == 0 then Left Identity else let [ u , v , d ] = extended_gcd ( 2*y_p ) n s = mod ( ( 3 * x_p * x_p + a ) * u ) n x_r = mod ( s * s - 2 * x_p ) n y_r = mod ( -y_p - s * ( x_r - x_p ) ) n in case ( ( Conpoint x_r y_r ) , d ) of ( _ , 1 )-> Left (Conpoint x_r y_r) ( _ , d' ) -> Right d' extended_gcd::Integer->Integer->[Integer] extended_gcd u v= helpGcd u v 0 [ [ 1 , 0 ] , [ 0 , 1 ] ] where helpGcd u v n m @ ( [ [ a , b ] , [ c , d ] ] ) | v == 0 = if u<0 then [ - ( ( -1 ) ^ n ) * ( m !! 1 !! 1 ) , - ( ( -1 ) ^ ( n + 1 ) ) * ( m !! 0 !! 1 ) , -u ] else [ ( ( -1 ) ^ n ) * ( m !! 1 !! 1 ) , ( ( -1 ) ^ ( n + 1 ) ) * ( m !! 0 !! 1 ) , u ] | otherwise = helpGcd u' v' ( n + 1 ) m' where ( q , v' ) = quotRem u v t = [ [q , 1 ] , [ 1 , 0 ] ] m' = [ [ q * a + b , a ] , [ q * c + d , c ] ] --mult m t u' = v multiEC :: Elliptic -> Point -> Integer -> IO ( Either Point Integer ) multiEC _ _ 0 = return$ Left Identity
multiEC ecurve point k  = return $helpEC Identity point k where helpEC p _ 0 = Left p helpEC p q n = case (p',q') of ( Left p'' , Left q'' ) -> helpEC p'' q'' (div n 2) ( Right p'' , _ ) -> Right p'' ( _ , Right q'' ) -> Right q'' where p' = if odd n then addPoints ecurve p q else Left p q' = addPoints ecurve q q dscrmntEC a b = 4 * a * a * a + 27 * b * b randomCurve :: Integer -> IO [Integer] randomCurve n = randomRIO ( 1 , n ) >>= ( \a -> randomRIO ( 1 , n ) >>= ( \u -> randomRIO (1 , n) >>= (\v -> return [a , mod ( v*v - u*u*u - a*u ) n , n , u , v ] ) ) ) factor :: Integer -> Integer -> IO [Integer] factor 1 _ = return [] factor n m = isProbable n >>= (\x -> if x then return [n] else randomCurve n >>= (\[ a , b , n , u , v ] -> multiEC ( Conelliptic a b n ) ( Conpoint u v ) m >>= ( \p -> case p of Left p' -> factor n m Right p'-> factor ( div n p' ) m >>= ( \x -> factor p' m >>= (\y -> return$ x ++ y   ) ) ) ) )

solve :: Integer -> IO [ Integer ]
solve n = factor n ( foldl lcm 1 [1..10000] )

main = liftM read getLine >>=( \n -> solve n  ) >>= print

[user@localhost Haskell]$time ./Elliptic_Len 14771291995722167329218042370113255785101138530361517229308539087571387624569184860194674818433642463913959873342852255543111135151457209398370030102228858128758106619882182645353368949658870563576166323307466834284670054906895828663228915811739807554988946409285249055824035034910383070200873326731464685909980630202533115332294811841112677175465781631555282945350555352712075708686158031835188673039526232167090966030801645637 [1759,1997,1291,1637,1669,1889,1907,1823,1231,1579,1597,1327,1427,1499,1831,1951,1213,1753,1993,1487,1429,1601,1319,1447,1721,1847,1697,1873,1699,1663,1061,1879,1087,1093,1117,1423,1933,1493,1303,1021,1109,1973,1361,1811,1123,1543,1249,1783,1867,1301,1789,1901,1657,1693,1801,1091,1609,1913,1277,1787,1483,1283,1511,1381,1877,1297,1181,1571,1103,1229,1039,1471,1289,1373,1861,1069,1987,1481,1279,1559,1193,1129,1747,1949,1733,1489,1063,1201,1307,1187,1871,1031,1741,1433,1321,1013,1567,1627,1549,1399,1553,1523,1051,1439,1723,1259,1667,1621,1777,1151,1367,1019,1097,1451,1163,1223,1999,1979,1931,1583,1619,1217,1153,1009,1409,1613,1459,1171,1531,1453,1237,1607,1049,1709,1033] real 0m16.424s user 0m1.905s sys 0m0.009s  July 18, 2011 ## SPOJ 6560. N DIV PHI_N (Hard) SPOJ 6560. N DIV PHI_N (Hard) also accepted using the same idea accept using a bit memoization. Precalculation of product of primes in prime table $<=$ 10^25000 and linear search is enough to get accepted.Its first 6000 primes whose products are sufficient . Accepted Haskell code. I also tried a slightly improved Rabin Miller primality testing with idea give topcoder forum for faster modular exponentiation for large numbers and code with isProbable function got accepted with 3.95 and time limit exceed other time. import Data.List import Data.Maybe import Data.Bits import qualified Data.ByteString.Char8 as BS powM::Integer->Integer->Integer->Integer 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 powL :: Integer -> Integer -> Integer -> Integer powL a d n | d <= 10^10 = powM a d n | otherwise = mod ( powM ( powL a ( div d 1000000 ) n ) 1000000 n * powM a ( mod d 1000000 ) n ) n isProbable::Integer-> Bool isProbable n |n<=1 = False |n==2 = True |even n = False |otherwise = rabinMiller [2,3] n s d where s = until ( \x -> testBit ( n - 1) ( fromIntegral x ) ) ( +1 ) 0 d = div ( n - 1 ) ( shiftL 1 ( fromIntegral s ) ) --d=until (\x->mod x 2/=0) (div 2) (n-1) --(n-1=2^s*d) --s=until (\x->d*2^x==pred n) (+1) 0 --counts the power rabinMiller::[Integer]->Integer->Integer->Integer-> Bool rabinMiller [] _ _ _ = True rabinMiller (a:xs) n s d = if a == n then True else case x==1 of True -> rabinMiller xs n s d _ -> if any ( == pred n) . take ( fromIntegral s ) . iterate ( \e -> mod (e^2) n )$  x then rabinMiller xs n s d
else  False
where x=powL a d n

prime :: [Integer]
prime = 2:3:filter isPrime [2*k+1 | k<-[2..]]

isPrime :: Integer -> Bool
isPrime n = all ((/=0). mod n  ).takeWhile ( <= ( truncate.sqrt.fromIntegral $n) )$ prime

anstable ::[Integer]
anstable =  scanl1 (*)   prime

helpSolve :: [Integer] -> Integer -> Integer -> Integer
helpSolve (x:xs) n ret
| x > n = ret
| otherwise = helpSolve xs n x

solve :: Integer -> BS.ByteString
solve 1 = BS.pack "1"
solve n = BS.pack.show $a where a = helpSolve ( tail anstable ) n ( head anstable ) readInt :: BS.ByteString -> Integer readInt = fst.fromJust.BS.readInteger main = BS.interact$ BS.unlines.map ( solve .readInt ) .  tail . BS.lines


July 17, 2011

## Solovay–Strassen primality test

The Solovay–Strassen primality test, developed by Robert M. Solovay and Volker Strassen, is a probabilistic test to determine if a number is composite or probably prime. Wikipedia suggests that for any odd prime p and integer a $a^{\frac {p-1} 2 } \equiv(\frac a p )\mod p$ where $\left(\tfrac{a}{p}\right)$ is the Legendre symbol. The Jacobi symbol is a generalisation of the Legendre symbol to $\left(\tfrac{a}{n}\right)$, where n can be any odd integer.Given an odd number n

$a^{\frac {n-1} 2} \equiv \left(\frac{a}{n}\right) \pmod n$

holds for various values of a. If n is prime then this congruence is true for all a.Jacobi code is take from here.

import Data.List
import System.Random

powM :: Integer -> Integer -> Integer -> Integer
powM a b m
| b == 0 = 1
| b == 1 = mod a m
| odd b  =  mod ( a * powM ( mod ( a^2 ) m ) ( div b 2 ) m ) m
| otherwise =  mod ( powM ( mod ( a^2) m  ) ( div b 2 ) m ) m

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

isProPrime :: Integer -> Integer -> Bool
isProPrime n a
| gcd a n /= 1 = False
| otherwise = powM a ( div ( n - 1 ) 2 ) n == mod ( jacobi a n ) n

solovayStrassen :: Integer -> IO Bool
solovayStrassen n
| n < 2 = return False
| n == 2 = return True
| even n = return False
| otherwise = return.isProPrime n =<< randomRIO ( 2 , n-1 )


June 23, 2011