My Weblog

Blog about programming and math

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 
Advertisements

July 17, 2011 - Posted by | Programming | , , , ,

2 Comments »

  1. Looks a little too heavy in my opinion. Why not just use sieve of Eratosthenes from haskellwiki? Although anyway nice.

    Comment by gorlum0 | July 26, 2011 | Reply

    • Yes its bit heavy because of Rabin Miller code. Thanks for the suggestion for sieve of Eratosthenes. It did not occur to me to solve this problem using array.

      Comment by tiwari_mukesh | July 26, 2011 | Reply


Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: