My Weblog

Blog about programming and math

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 | , , , | Leave a 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 | , , | 1 Comment

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 | , , , , | 2 Comments

Project Euler 197

Project Euler 197 was totally based on intuition ( its still working although now a days I am not very active ). After doing bit of math, I saw the pattern and rest of work was done by Haskell ( iterate function ). I just wrote the function f and ran

ghci>take 4 . drop 0 .  iterate f $ ( -1 )
[-1.0,0.7100000000000001,1.001242148,0.708777686]
ghci>take 4 . drop 10 .  iterate f $ ( -1 )
[1.005832407,0.7042658820000001,1.0068861020000002,0.7032313570000001]
ghci>take 4 . drop 500 .  iterate f $ ( -1 )
[1.029461825,0.6811758910000001,1.029461827,0.681175889]

Haskell source code

import Data.List

f :: Double -> Double
f x = 1E-9 * ret where
   a = 30.403243784
   y = a - x ^ 2
   ret = fromIntegral ( floor $ 2 ** y )

recur :: [ Double ] -> Double
recur ( a : b : c : xs )
   | abs ( a - c ) <= 1E-10   = a + b
   | otherwise = recur ( b : c : xs )

solve :: Double
solve = recur . iterate  f $ ( - 1 )

main :: IO ()
main = print $ solve  

August 10, 2012 Posted by | Programming | , , | Leave a Comment

Project Euler Problem 107 ( Minimum Spanning Tree )

Project Euler Problem 107 is related to minimum spanning tree. Trying to solve some easy problems to cross the 200 mark. This is the first time I have implemented the graph algorithm in Haskell so it was kind of fun. Wrote Prim’s algorithm to calculated the minimum spanning tree. Replaced comma ( , ) in input file by space. See modified input file and total sum of inputs from ideone.

import Data.List
import qualified Data.IntMap.Lazy as M
import qualified Data.PQueue.Min as PQ
import qualified Data.Bits as B


buildGraph :: [ ( Int , Int , Int ) ]  -> M.IntMap [ ( Int , Int ) ]
buildGraph xs = M.fromListWith ( ++ ) . map f_1  $ xs where
          f_1 ( i , j , val ) = ( i , [ (  val , j ) ] )



visited :: Int -> Integer -> ( Bool , Integer )
visited x vis = ( t == 0 , vis' ) where
            t = ( B..&. ) ( B.shiftL ( 1 :: Integer ) x ) vis
            vis' = ( B..|. ) ( B.shiftL ( 1 :: Integer ) x ) vis


primAlgorithm ::  PQ.MinQueue ( Int , Int )  -> Integer -> Int -> M.IntMap [ ( Int , Int ) ] -> String
primAlgorithm  queue  vis cost m
       | PQ.null queue =  show cost  ++ "\n"
       | t == False = primAlgorithm  queue' vis' cost m
       | otherwise =   primAlgorithm queue'' vis' ( cost + c ) m where 
              ( ( c , x ) , queue' ) = PQ.deleteFindMin queue 
              ( t , vis' ) = visited x vis
              queue'' = PQ.union queue' ( PQ.fromList . ( M.! ) m  $ x )


test :: [ ( Int , Int , Int ) ] -> String
test xs = primAlgorithm  ( PQ.fromList [ ( 0 , 0 ) ] ) 0 0  ( buildGraph xs ) 

parseInput :: [ ( Int , String ) ] -> [ ( Int , Int ) ]
parseInput [] = []
parseInput ( ( k , x )  : xs ) 
   | x == "-" = parseInput xs 
   | otherwise = ( k , readD x ) : parseInput xs where 
        readD :: String -> Int
        readD = read

fun :: Int -> [ [ ( Int , Int ) ] ] -> [ [ ( Int , Int , Int ) ] ]
fun _ [] = []
fun k ( x : xs ) = ( map ( \( i , val ) -> ( k , i , val ) ) x ) : fun ( succ k ) xs 


main =  interact $  primAlgorithm ( PQ.fromList [ ( 0 , 0 ) ] ) 0 0 .  buildGraph . concat . fun 0 . map ( parseInput  . zip [ 0 , 1 .. ] . words ) . lines   

Macintosh-0026bb610428:Programming mukesh$ ghc-7.4.1 –make -O2 Prim.hs
[1 of 1] Compiling Main ( Prim.hs, Prim.o )
Linking Prim …
Macintosh-0026bb610428:Programming mukesh$ ./Prim < network.txt
2153

August 4, 2012 Posted by | Programming | , , , , , | Leave a Comment

Project Euler Problem 381

Project Euler Problem 381 is easy problem. Using Wilson theorem
1\cdot 2\cdots ( p - 1 )\ \equiv\  -1\  \pmod{p}
1\cdot 2\cdots ( p - 2 )\ \equiv\ \frac{-1}{p-1}   \pmod{p}
1\cdot 2\cdots ( p - 3 )\ \equiv\ \frac{-1}{( p-1) (p-2)}   \pmod{p}
1\cdot 2\cdots ( p - 4 )\ \equiv\ \frac{-1}{( p-1) (p-2)(p-3)}   \pmod{p}
1\cdot 2\cdots ( p - 5 )\ \equiv\ \frac{-1}{( p-1) (p-2)(p-3) ( p-4 )}   \pmod{p}

Haskell source code using monad par

import Data.List
import Data.Numbers.Primes
import Control.Monad.Par
import GHC.Conc.Sync

--use wilson theorem ( p - 1 ) ! = -1 mod p 

prime = drop 2 . takeWhile ( < 10 ^ 8 ) $ primes

extendedGcd :: Int -> Int -> ( Int , Int )                 
extendedGcd a b
   | b == 0 = ( 1 , 0 )
   | otherwise = ( t , s - q * t ) where
        ( q , r ) = quotRem a b
        ( s , t ) = extendedGcd b r


modInv :: Int -> Int -> Int
modInv  a b
   | gcd a b /= 1 = error " gcd is not 1 "
   | otherwise = d where
        d = until ( > 0 ) ( + b  ) . fst.extendedGcd a $ b


computeFinal :: Int -> Int
computeFinal p = computeS ( p - 1 ) 0 1 where
     computeS inv ret cnt
        | cnt >= 5 = mod ( ret + inv ) p
        | otherwise = computeS ( mod ( inv * inv' ) p  )  ( mod ( ret + inv ) p )  ( cnt + 1 )   where
         inv' = modInv ( p - cnt ) p

chunk :: [ Int ] -> Int -> [ [ Int ] ]
chunk xs n = ret where
    tot = numCapabilities
    brk = div n tot
    ret = unfoldr fun ( xs , 1 ) where
    fun ( xs' , cnt )
     | null xs' = Nothing
     | cnt >= tot = Just ( xs' , ( [] , cnt ) )
     | otherwise = Just ( a , ( b , cnt + 1 ) ) where
        ( a , b ) = splitAt brk xs'

ans =   ret'  where
  ret = runPar $ parMap ( map computeFinal )  ( chunk prime len )
  len = length prime
  ret' = sum . map sum $ ret


main = do
   print ans

Compile it “ghc-7.4.1 -O2 -rtsopts -threaded -fforce-recomp Euler_381.hs” and run it with as many cores you have. My run time is
[mukesh@user Euler]$ time ./Euler_381 +RTS -A16M -H2G -N8 -RTS
139602943319822

real 0m28.473s
user 0m58.204s
sys 0m21.713s

April 28, 2012 Posted by | Programming | , , , | Leave a Comment

Follow

Get every new post delivered to your Inbox.

Join 187 other followers