My Weblog

Blog about programming and math

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

September 9, 2011 Posted by | Programming | , , , | Leave a Comment

Converting Wikipedia html files in pdf

I want to convert html files from Wikipedia to pdf for off line reading purpose . After bit of searching , Wikipedia itself provides a link on left side [ Print/export ] of every article to convert it into pdf . After couple of clicks , we can download the pdf but I want to write Haskell script. This script generates the rendering url. Rendering url return empty tags while copy and pasting the rendering url to web browser generates the pdf file. After asking on Haskell-cafe revealed that the link is generated by javascript and i have to script an actual browser to generated pdf from this code. Technically this is still unfinished project :( but first time I played with some sort of web programming.

import Network.HTTP
import Text.HTML.TagSoup
import Data.Maybe
 
parseHelp :: Tag String -> Maybe String 
parseHelp ( TagOpen _ y ) = if any ( \( a , b ) -> b == "Download a PDF version of this wiki page" ) y 
                             then Just $  "http://en.wikipedia.org" ++   snd (   y !!  0 )
                              else Nothing
 
 
parse :: [ Tag String ] -> Maybe String
parse [] = Nothing 
parse ( x : xs ) 
   | isTagOpen x = case parseHelp x of 
                         Just s -> Just s 
                         Nothing -> parse xs
   | otherwise = parse xs
 
 
main = do 
        x <- getLine 
        tags_1 <-  fmap parseTags $ getResponseBody =<< simpleHTTP ( getRequest x ) --open url
        let lst =  head . sections ( ~== "<div class=portal id=p-coll-print_export>" ) $ tags_1
            url =  fromJust . parse $ lst  --rendering url
        putStrLn url
        tags_2 <-  fmap parseTags $ getResponseBody =<< simpleHTTP ( getRequest url )
        print tags_2
 

My second choice was obviously python and it finished the job perfectly . Python script for this purpose and in fact it can convert any html file to pdf. Its like opening a html file in web browser and printing it to pdf file.
import sys
from PyQt4.QtCore import *
from PyQt4.QtGui import *
from PyQt4.QtWebKit import *

#http://www.rkblog.rk.edu.pl/w/p/webkit-pyqt-rendering-web-pages/
#http://pastebin.com/xunfQ959
#http://bharatikunal.wordpress.com/2010/01/31/converting-html-to-pdf-with-python-and-qt/
#http://www.riverbankcomputing.com/pipermail/pyqt/2009-January/021592.html

def convertFile( ):
                web.print_( printer )
                print "done"
                QApplication.exit()


if __name__=="__main__":
        url = raw_input("enter url:")
        filename = raw_input("enter file name:")
        app = QApplication( sys.argv )
        web = QWebView()
        web.load(QUrl( url ))
        #web.show()
        printer = QPrinter( QPrinter.HighResolution )
        printer.setPageSize( QPrinter.A4 )
        printer.setOutputFormat( QPrinter.PdfFormat )
        printer.setOutputFileName(  filename + ".pdf" )
        QObject.connect( web ,  SIGNAL("loadFinished(bool)"), convertFile  )
        sys.exit(app.exec_())
~                              

September 9, 2011 Posted by | Programming | , , | Leave a Comment

SRM 385 – Division II, Level Three

Gettc is excellent software for practicing topcoder problems offline and specially in Haskell. I just wrote solution for SRM 385 Problem [ requires registration but its worth and excellent place for programmers ]. You can see the editorial for SRM 385. Haskell source code.

module SummingArithmeticProgressions where  

canRep :: Int -> Int -> Int -> Bool 
canRep a b n = canRephelp 1 a b n where 
	canRephelp x a b n 
	  | x > b = False
	  | n <= a * x = False
	  | mod ( n - a * x ) b == 0 = True
	  | otherwise = canRephelp  ( succ x ) a b n 


howMany :: Int -> Int -> Int -> Int
howMany left right k =  ans   where 
	a = k 
	b = div ( k * ( k - 1 ) ) 2
	d = gcd a b
	a' = div a d 
	b' = div b d
	left' = div ( left + d - 1 ) d 
	right' = div right d 
	ans = helphowMany a' b' left' right'  0
	

helphowMany :: Int -> Int -> Int -> Int -> Int -> Int
helphowMany a b left right cnt
	| and [ left <= right , left < 2 * a * b ] = if canRep a b left 
						      then helphowMany a b ( succ left ) right ( succ cnt )  
						       else helphowMany a b ( succ left ) right cnt 
	| otherwise = cnt + right - left + 1 

September 7, 2011 Posted by | Programming | , , | 1 Comment

   

Follow

Get every new post delivered to your Inbox.

Join 138 other followers