Sieve of eratosthenes using Haskell IOArray

More about Haskell array . This one is simple Haskell code for Sieve of eratosthenes and Euler’s totient function . This is first time i played with mutable array. Also posted this code on ideone.

import Data.Array
import Data.Array.IO
--start of prime

loopB :: Int -> Int -> Int ->  IOArray Int Int -> IO ( IOArray Int Int )
loopB st a b arr
        | st > b = return arr
	| otherwise = 
	    do
		writeArray arr st 0 
		loopB ( st + a ) a b arr 


loopA :: Int -> Int ->  IOArray Int Int -> IO ( IOArray Int Int ) 
loopA a b arr 
	| a ^ 2 > b = return arr
	| otherwise = 
	   do
		t <- readArray arr a 
		arr' <- if t /= 0 then loopB ( a ^ 2 ) a  b arr else return arr 
		loopA ( a + 1 ) b arr' 

printerA :: Int -> Int -> IOArray Int Int -> IO ( )
printerA a b arr 
	| a > b = return ()
	| otherwise = do 
		t <- readArray arr a 
		print ( a , if t == 0 then "composite" else "prime" ) 
		printerA ( a + 1 ) b arr 

--end of prime generation
--start of euler phi function

helpEuler :: Int -> Int -> IOArray Int Int -> IO ( IOArray Int Int ) 
helpEuler a b phi 
	| a > b = return phi
	| otherwise = 
	   do 
		writeArray phi a a 
		helpEuler ( a + 1 ) b phi

loopInner :: Int -> Int -> Int -> IOArray Int Int -> IO ( IOArray Int Int ) 
loopInner st a b phi 
	| st > b = return phi
	| otherwise = 
	   do 
		t <- readArray phi st 
		let k = div ( t * ( a - 1 ) ) a 
		writeArray phi st k 
		loopInner ( st + a ) a b phi

loopOuter :: Int -> Int -> IOArray Int Int -> IO ( IOArray Int Int ) 
loopOuter a b phi
	| a > b = return phi
	| otherwise = 
	   do 
	     t <- readArray phi a 
	     phi' <- if t == a then loopInner (  2 * a ) a b phi else return phi
	     loopOuter ( a + 1 ) b phi' 
 
modifyValue :: Int -> Int -> IOArray Int Int -> IO ( IOArray Int Int ) 
modifyValue a b phi 
	| a > b = return phi
	| otherwise = 
	   do 
	      t <- readArray phi a 
	      if t == a then writeArray phi a ( a - 1 ) else writeArray phi a t 
	      modifyValue ( succ a ) b phi 


printerB :: Int -> Int -> IOArray Int Int -> IO (  )
printerB a b phi
       | a > b = return ()
       | otherwise = 
	   do 
		t <- readArray phi a 
		print ( a , t ) 
		printerB ( a + 1 ) b phi

--end of euler phi function


main = do
	{-- 
	arr <- newArray ( 2 , 1000 ) 1 :: IO ( IOArray Int Int ) 
	arr' <- loopA  2 1000 arr
	printerA 2 1000 arr'
	--}
	
	phi <- newArray ( 1 , 100 ) 1 :: IO ( IOArray Int Int ) 
	tmp <- helpEuler  2 100 phi
	phi' <- loopOuter 2 100 tmp	
	phi'' <- modifyValue 2 100 phi'
	printerB 1 100 phi''
	
	{--
	a <- readArray arr 1
	writeArray arr 1 64 
	b <- readArray arr 1
	print ( a , b ) 
	--}

Couple of more compact code to compute prime , totient , divisor sum and mobius of number.

import Control.Monad.ST
import Data.Array.ST
import Data.Array.Unboxed
import Control.Monad

prime :: Int -> UArray Int Int 
prime n = runSTUArray $ do 
	arr <- newArray ( 2 , n ) 1 :: ST s ( STUArray s Int Int ) 
	forM_ ( takeWhile ( \x -> x*x <= n ) [ 2 .. n ] ) $ \i -> do 
		ai <- readArray arr i 
		when ( ai == 1 ) $ forM_ [ i^2 , i^2 + i .. n ] $ \j -> do 
			writeArray arr j 0 

	return arr 

totient :: Int -> UArray Int Int 
totient n = runSTUArray $ do 
	arr <- newListArray ( 1 , n ) [ 0 .. pred n ] :: ST s (STUArray s Int Int ) 
	forM_ [ 2 .. n ] $ \i -> do 
		ai <- readArray arr i
		when ( ai == pred i ) $ forM_ [ 2*i , 3*i .. n ] $ \j -> do 
			aj <- readArray arr j 
			writeArray arr j ( aj - div aj i )
	writeArray arr 1 1  
	return arr

divisorSum :: Int -> UArray Int Int
divisorSum n = runSTUArray $ do
	arr <- newArray ( 1 , n ) 0 :: ST s ( STUArray s Int Int ) 
	forM_ [ 1 .. n ] $ \i-> do 
	  forM_ [ i , 2 * i .. n ] $ \j -> do
		aj <- readArray arr j 
		writeArray arr j ( aj + i ) 
	return arr


mobius :: Int -> UArray Int Int 
mobius n = runSTUArray $ do
	arr <- newArray ( 1 , n ) 0 :: ST s ( STUArray s Int Int ) 
	writeArray arr 1 1
	forM_ [ 1 .. n ] $ \i -> do
	  forM_ [ 2*i , 3*i .. n ] $ \j -> do
	    ai <- readArray arr i
	    aj <- readArray arr j
	    writeArray arr j ( aj - ai )
	return arr 	 


main = do
	
	print $ ( prime 100 ) ! 37
	print $ ( totient 100  ) ! 50
	print $ ( divisorSum 100 ) ! 50
	print $ ( mobius 100 ) ! 50

About these ads