My Weblog

Blog about programming and math

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

State Monads

Finally i got little bit about Haskell monads. There are lot of articles on net but Mike’s World-O-Programming is really excellent . Thanks to Divyanshu for pointing me to the link. Any way i am not going to write another monad tutorial because i am not fit for this topic. I wrote couple of codes using state monads.

import Control.Monad
import Control.Monad.State
import System.Random

--from  wiki books about Haskell
rollDice :: Int -> State StdGen Int 
rollDice n = get >>= ( \s ->  ( return $ randomR ( 1 , n ) s ) >>= ( \ ( v , n ) -> put n >> return v ) )

gcD :: State ( Integer , Integer ) Integer 
gcD = get >>= ( \( a , b ) -> case b == 0 of 
				True -> put ( a , a ) >> return a  --change the state before returning the result  just for fun 
				_    -> put ( b , mod a b ) >> gcD )


fibo :: Int -> State ( Integer , Integer , Int )  Integer 
fibo n = get >>=( \( a , b , t ) -> case t == n of 
					True -> return a 
					_    -> put ( b , a + b , t + 1 ) >> fibo n )
fib :: State ( Integer , Integer ) () 
fib = get >>= ( \( a , b ) -> put ( b , a + b ) ) 

binaryGcd :: State ( Integer , Integer , Int ) Integer 
binaryGcd = get >>=( \( u , v , k ) ->  
			case () of 
			  _ | u == 0 -> return ( v *  2 ^ k  )
			    | v == 0 -> return ( u *  2 ^ k  ) 
			    | and [ even u , even v ] -> put ( div u 2 , div v 2 , succ k ) >> binaryGcd 
			    | and [ even u , odd v ] -> put ( div u 2 , v , k ) >> binaryGcd 
			    | and [ odd u , even v ] -> put ( u , div v 2 , k ) >> binaryGcd 
			    | otherwise -> if u >= v then put ( div ( u - v ) 2 , v , k ) >> binaryGcd 
					    else put ( div ( v - u ) 2 , u , k ) >> binaryGcd  )


findMax :: Int -> State Int () 
findMax n = get >>= ( \x -> case () of 
				_ | x > n -> put ( x ) 
				  | otherwise ->  put ( n )  )


sumList :: [ Integer ] -> State Integer Integer
sumList [] = get >>=(\t -> return t )
sumList ( x : xs ) = 
     get >>= ( \ t -> put ( t + x ) >> sumList xs ) 

sum' :: State ( Integer , Integer ) Integer
sum' = get >>=( \( a , b ) -> return $ a + b ) --here we don't need to change the state 


main = do 
	print $ evalState ( rollDice 1000 ) ( mkStdGen 1000) 
	print $ evalState gcD ( 12345 , 3123123 )
	print $ evalState ( fibo 10 ) ( 0 , 1 , 1 ) -- first two term 0 and 1 and count 0 
	print $ evalState  binaryGcd ( 12345 , 3123123  , 0 )
	print $ execState ( mapM findMax [ 1.. 10 ] ) 0 --remember our maximum value will be kept in state variable
	print $ fst . execState ( replicateM ( 10 - 1 )  fib ) $ ( 0 , 1 )  --to get nth decrase it one [ replicateM ( n - 1 ) fib ]
	print $ evalState ( sumList [ 1..10 ] ) 0
	print $ evalState sum' ( 10 , 100 )  

August 24, 2011 Posted by | Programming | , | 1 Comment

SPOJ 3931. Maximum Triangle Area

SPOJ 3931. Maximum Triangle Area is related to geometry and maximum area triangle will lie on convex hull of given points . My first thought was rotating calipers and I assumed that maximum area triangle’s one side will overlap with side of convex polygon which is totally incorrect. A bit of searching led me to this stackoverflow post. In short , the algorithm is [ borrowed from spoj user vipul ]
Required triangle’s edges may not coincide with that of convex hull but the 3 points will coincide with its vertices.
1. Choose a,b,c as first three points of convex hull( initial area = area of triangle abc )
2. Keep selecting next point to c as new c till area increases.
3. Now, select next point to b as new b, if area increases go to step 2.
4. Repeat 2,3 for all a and keep a track of maximum area found so far.
Accepted Haskell code

import Data.List
import Data.Array
import Data.Maybe
import qualified Data.ByteString.Lazy.Char8 as BS
import Text.Printf


--monotone

data Point a = P a a deriving ( Show , Ord , Eq ) 
data Turn = RS | L deriving ( Show , Ord  , Eq , Enum ) 

calAngle :: ( Num a , Ord a , Eq a ) => Point a -> Point a -> Point a -> Ordering
calAngle ( P x0 y0 ) ( P x1 y1 ) ( P x2 y2 ) = compare ( ( y1 - y0 ) * ( x2 - x0 ) ) ( ( y2 - y0 ) * ( x1 - x0 ) )


sortByangle :: ( Num a , Ord a , Eq a ) => [ Point a ] -> [ Point a ]
sortByangle ( x : xs ) = x : sortBy (\ p1 p2 -> calAngle x p1 p2 ) xs 

findTurn :: ( Num a , Ord a , Eq a ) => Point a -> Point a -> Point a -> Turn
findTurn ( P x0 y0 ) ( P x1 y1 ) ( P x2 y2 ) = 
   case compare ( ( y1 - y0 ) * ( x2 - x0 ) ) ( ( y2 - y0 ) * ( x1 - x0 ) ) of 
	    LT -> L 
	    _  -> RS 


computeHull :: ( Num a , Ord a , Eq a ) => [ Point a ] -> [ Point a ] -> [ Point a ] 
computeHull [ x ] ( z : ys ) = computeHull [ z , x ] ys 
computeHull ys [] =  ys 
computeHull u@( y : x : xs ) t@( z : ys )  
  | findTurn x y z == RS = computeHull ( x : xs ) t 
  | otherwise  = computeHull ( z : u ) ys 


convexHull :: ( Num a , Ord a , Eq a ) => [ Point a ] -> [ Point a ]
convexHull xs = final   where 
	t1@( x1 : y1 : xs1)  = sort xs 
	lhull = computeHull [ y1 , x1 ] xs1 
	t2@( x2 : y2 : xs2 ) = reverse t1 
	uhull = computeHull [ y2 , x2 ] xs2 
	final =  nub $  lhull ++  uhull
	


--end of monotone


caltmp :: ( Num a , Ord a , Floating a ) => Int -> Int -> Int -> Array Int ( Point a ) -> a 
caltmp a b c arr = area where 
	P x1 y1 = arr ! a 
	P x2 y2 = arr ! b 
	P x3 y3 = arr ! c 
	area = 0.5 * ( abs $ ( x1 * y2 + x2 * y3 + x3 * y1 )  -  ( y1 * x2 + y2 * x3 + y3 * x1 ) )
	 

calArea :: ( Num a , Ord a , Floating a ) => Int -> Int -> Int  -> Int  -> a -> Array Int ( Point a ) -> ( Int , Int , Int , a ) 
calArea a b c  n area arr
 | area' >= area = calArea a b ( mod ( c + 1 ) n )   n area' arr      --area a b  ( c + 1 )  is greater than area a b c
 | area'' >=  area  = calArea a ( mod ( b + 1 ) n )  c  n area'' arr   --check if area a ( b + 1 ) c is greater area a b c 
 | otherwise =  ( a , b , c , area ) 
	where 
	 area' = caltmp a b ( mod ( c + 1 ) n ) arr
	 area'' = caltmp a ( mod ( b + 1 ) n ) c arr 
		    
looPing :: ( Num a , Ord a , Eq a , Floating a ) => Int -> Int -> Int -> Int -> a -> a -> Array Int ( Point a ) -> a 
looPing a b c n area best arr  
 | a == n = max area best
 | otherwise = looPing a'' b'' c'' n area'' ( max area' best )  arr where 
	( a' , b' , c' , area' ) = calArea a b c n area arr 
	a'' = a' + 1 
	b'' = if a'' == b' then mod ( b' + 1 ) n else b'
	c'' = if b'' == c' then mod ( c' + 1 ) n else c'
 	area'' = caltmp ( mod a'' n ) b'' c'' arr

solve :: ( Num a , Ord a , Floating a ) => [ Point a ] ->  a 
solve [] = 0
solve [ p ] = 0
solve [ p1 , p2 ] = 0
solve arr =  
	looPing 0 1 2  n area area arr' where 
	n = length arr 
	arr' = listArray ( 0 , pred n ) arr
        area = caltmp 0 1 2 arr'	

 
final :: ( Num a , Ord a , RealFloat a ) => [ Point a ] -> a
final [] = 0
final [ p ] =  0 
final [ p1 , p2 ] =  0 
final arr = solve . convexHull $ arr 


format :: ( Num a , Ord a , Floating a ) => [ Int ] -> [ [ Point a ]]
format [] = []
format (x:xs ) =  t : format b  where
        ( a , b ) = splitAt ( 2 * x ) xs
        t = helpFormat a where
            helpFormat [] = []
            helpFormat ( x' : y' : xs' ) =  P ( fromIntegral x' ) ( fromIntegral  y' )  : helpFormat xs'
 
readD :: BS.ByteString -> Int
readD = fst . fromJust . BS.readInt
 
main = BS.interact $ BS.unlines . map ( BS.pack . ( printf "%.2f" :: Double -> String ) . final ) . format . concat . map  ( map  readD . BS.words ) . init  . BS.lines  


--main = interact $ unlines . map ( show . convexHull ) . format . concat . map ( map read . words ) . init . lines 

August 17, 2011 Posted by | Programming | , , , , | Leave a Comment

SPOJ 515. Collatz

SPOJ 515. Collatz is bit archaic in expression. Its 3 * n + 1 problem . f(1) = N and f(K) = ( 0.5 + 2.5 * ( f(K-1) \bmod 2)) * f(K-1) + (f(K-1) \bmod 2) if K>1. If f(k-1) is odd then given expression f(K) = 3 * f(K-1) + 1 and if f(k-1) is even then f(K) = \frac {f(k-1)}{2} . Accepted Haskell code.

import Data.List
import Data.Maybe
import qualified Data.ByteString.Lazy.Char8 as BS


{--
collatz :: Integer -> Integer
collatz n = head . filter (\t -> helpCollatz t == 1 ) $  [ 1.. ] where 
	helpCollatz 1 = n 
	helpCollatz k = l where 
		t = helpCollatz ( k - 1 ) 
		l = if odd t then 3 * t + 1 else div t 2 
--}

solve :: Integer -> Integer
solve n = succ . fromIntegral . fromJust . findIndex ( \t -> t == 1 ) . scanl ( \ x y -> if odd x then 3 * x + 1 else div x 2 ) n $ [ 2..]

readInt :: BS.ByteString -> Integer
readInt = fst . fromJust . BS.readInteger

main = BS.interact $ BS.unlines. map ( BS.pack.show . solve . readInt ) . BS.lines 

August 8, 2011 Posted by | Programming | , , , | 3 Comments

SPOJ 4421. Irreducible polynomials over GF2

SPOJ 4421. Irreducible polynomials over GF2 is too easy in Haskell . For algorithm , mathworld explains it very well. My only concern was calculation of mobius function but under given time limit , mobius can be calculated using prime factorisation. Today i learn nice use of sequence from here.Before this , i used sequence only with IO monads. Accepted Haskell code

import Data.List
import Control.Monad
import Data.Maybe
import qualified Data.ByteString.Char8 as BS

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

pfactor :: Integer -> [ Integer ] -> [Integer]
pfactor n p@( x : xs ) 
   | n <= 1 = [] 
   | x ^ 2  > n = [ n ]
   | mod n x == 0 = x : pfactor ( div n x ) p 
   | otherwise = pfactor n xs 

mobius :: Integer -> Integer 
mobius 1 = 1 
mobius n 
    | facs == ( nub facs ) = ( -1 ) ^ len
    | otherwise = 0
	where 
	 facs = pfactor n prime 
	 len = genericLength facs 

--nice use of monads . https://github.com/bjin/fancy-walks/blob/master/spoj.pl/GF2.hs
divisor :: Integer -> Integer -> [Integer]
divisor n  p = sort $ map product $ sequence mul
  where
    fs = group $ pfactor n prime
    mul = map (scanl (*) 1) fs 

{--
divisor :: Integer -> Integer -> [ Integer ]
divisor n p 
  | n <= 1 = [1]
  | p * p > n = []
  | mod n p == 0 = p : ( div n p ) : divisor n (  p  + 1 ) 
  | otherwise = divisor n  ( p + 1 ) 
--}

solve :: Integer  -> String 
solve  n = show $ div ans n  where 
	ans = sum .  map ( \t -> mobius ( div n t ) * 2 ^ t ) . sort . nub . divisor  n $ 1


main = interact $ unlines . map ( solve . read ) . lines 

August 6, 2011 Posted by | Programming | , , , | Leave a Comment

SPOJ 8612. Penney Game

SPOJ 8612. Penney Game is an easy simulation problem . Accepted Haskell code , not very Haskellish but accepted .

import Data.List
import qualified Data.ByteString.Char8 as BS



countfN :: String -> ( Int , Int , Int , Int , Int , Int , Int , Int ) -> String
countfN [ x , y ]( a , b , c , d , e , f , g , h )  = foldl ( \ x y ->  x ++ " " ++ show y ) ( show a ) [ b , c , d , e , f , g , h ]   
countfN ( x : y : z : xs ) ( a , b , c , d , e , f , g , h ) 
	| x : y : z :[] == "TTT" = countfN ( y : z : xs ) ( a + 1 , b , c , d , e , f , g , h )
	| x : y : z :[] == "TTH" = countfN ( y : z : xs ) ( a  , b + 1 , c , d , e , f , g , h )
	| x : y : z :[] == "THT" = countfN ( y : z : xs ) ( a  , b , c + 1 , d , e , f , g , h )
	| x : y : z :[] == "THH" = countfN ( y : z : xs ) ( a  , b , c , d + 1 , e , f , g , h )
	| x : y : z :[] == "HTT" = countfN ( y : z : xs ) ( a  , b , c , d , e + 1 , f , g , h )
	| x : y : z :[] == "HTH" = countfN ( y : z : xs ) ( a  , b , c , d , e , f + 1 , g , h )
	| x : y : z :[] == "HHT" = countfN ( y : z : xs ) ( a  , b , c , d , e , f , g + 1 , h )
	| otherwise =  countfN ( y : z : xs ) ( a  , b , c , d , e , f , g , h + 1 )


solve :: String -> String 
solve xs = countfN xs ( 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ) 
 
format :: [ String ] -> [ String ] 
format [] = []
format ( x : y : xs ) = y : format xs 

main = interact $ unlines . zipWith (\x y -> show x ++ " " ++  y ) [1..] . map solve . format . tail . lines 

August 5, 2011 Posted by | Programming | , , , | Leave a Comment

Rotating calipers [ The diameter of a convex polygon ]

Rotating calipers is great tool in computational geometry similar to sweep line . According to Wikipedia ” In computational geometry, rotating calipers is a method used to construct efficient algorithms for a number of problems”. Given code is Haskell transformation of pseudo code described at Wikipedia and computes the square of the diameter of a convex polygon . I tried to test this code on TFOSS but time limit is very tight so its hard for Haskell to get accepted.I generated some random test cases and it worked fine. A excellent paper ” Solving geometric problems with the rotating calipers ” investigates some problems of computational geometry which can be solved by this tool.

import Data.List
import Data.Array
import Data.Maybe
import Data.Function
import Text.Printf
import qualified Data.ByteString.Char8 as BS

data Point a = P a a deriving ( Show , Ord , Eq )
data Vector a = V a a deriving ( Show , Ord , Eq )
data Turn = S | L | R deriving ( Show , Eq , Ord , Enum  )

--start of convex hull

compPoint :: ( Num  a , Ord a ) => Point a -> Point a -> Ordering
compPoint ( P x1 y1 ) ( P x2 y2 )
  | compare x1 x2 == EQ = compare y1 y2
  | otherwise = compare x1 x2 

findMinx :: ( Num a , Ord a ) => [ Point a ] -> [ Point a ]
findMinx xs = sortBy ( \x  y  -> compPoint  x y  ) xs

compAngle ::(Num a , Ord a ) => Point a -> Point a -> Point a -> Ordering
compAngle ( P x1 y1 ) ( P x2 y2 ) ( P x0 y0 ) = compare ( (  y1 - y0 ) * ( x2 - x0 )  ) ( ( y2 - y0) * ( x1 - x0 ) )

sortByangle :: ( Num a , Ord a ) => [ Point a ] -> [ Point a ]
sortByangle (z:xs) = z : sortBy ( \x y -> compAngle x y z ) xs 

findTurn :: ( Num a , Ord a , Eq a ) => Point a -> Point a -> Point a -> Turn
findTurn ( P x0 y0 ) ( P x1 y1 ) ( P x2 y2 )
 | ( y1 - y0 ) * ( x2- x0 ) < ( y2 - y0 ) * ( x1 - x0 ) = L
 | ( y1 - y0 ) * ( x2- x0 ) == ( y2 - y0 ) * ( x1 - x0 ) = S
 | otherwise = R 

findHull :: ( Num a , Ord a  )  => [ Point a ] ->   [ Point a ] -> [ Point a ]
findHull [x]  ( z : ys )  = findHull [ z , x ]  ys  --incase of second point  on line from x to z
findHull xs  [] = xs
findHull ( y : x : xs )  ( z : ys )
  | findTurn x y z == R = findHull (  x : xs )   ( z:ys )
  | findTurn x y z == S = findHull (  x : xs )   ( z:ys )
  | otherwise = findHull ( z : y : x : xs  )   ys


convexHull ::( Num a , Ord a )  => [ Point a ] -> [ Point a ]
convexHull xs = reverse . findHull [ y , x ]  $ ys where
        ( x : y : ys ) = sortByangle . findMinx $ xs


--end of convex hull 

--start of rotating caliper part http://en.wikipedia.org/wiki/Rotating_calipers
--dot product for getting angle

angVectors :: ( Num a , Ord a , Floating a ) => Vector a -> Vector a -> a
angVectors ( V ax ay ) ( V bx by ) = theta where
    dot = ax * bx + ay * by
    a = sqrt $ ax ^ 2 + ay ^ 2
    b = sqrt $ bx ^ 2 + by ^ 2
    theta = acos $ dot / a / b  

--rotate the vector x y by angle t

rotVector :: ( Num a , Ord a , Floating a ) => Vector a -> a -> Vector a
rotVector ( V x y ) t = V ( x * cos t - y * sin t ) ( x * sin t + y * cos t )  

--square of dist between two points 

distPoints :: ( Num a , Ord a  ) => Point a -> Point a -> a
distPoints ( P x1 y1 ) ( P x2 y2 ) =  ( x1 - x2 ) ^ 2 + ( y1 - y2 ) ^ 2 

--rotating caliipers 

rotCal :: ( Num a , Ord a , Floating a ) => Array Int ( Point a )  -> a -> Int -> Int -> Vector a -> Vector a -> a -> Int -> a
rotCal arr ang  pa pb ca@( V ax ay ) cb@( V bx by ) dia n
   | ang > pi = dia
   | otherwise = rotCal arr ang' pa' pb' ca' cb' dia' n where
	P x1 y1 = arr ! pa
	P x2 y2 = arr ! ( mod ( pa + 1 ) n )
	P x3 y3 = arr ! pb
	P x4 y4 = arr ! ( mod ( pb + 1 ) n )
	t1 = angVectors ca ( V ( x2 - x1 ) ( y2 - y1 ) )
	t2 = angVectors cb ( V ( x4 - x3 ) ( y4 - y3 ) )
	ca' = rotVector ca  $ min t1 t2
	cb' = rotVector cb  $ min t1 t2
	ang' = ang + min t1 t2
	( pa' , pb' )  = if t1 < t2 then ( mod ( pa + 1 ) n  , pb ) else ( pa , mod ( pb + 1 ) n  )
	dia' = max dia $ distPoints ( arr ! pa' ) ( arr ! pb' )

solve :: ( Num a , Ord a , Floating a ) => [ Point a ] -> a
solve [] = 0
solve [ p ] = 0
solve [ p1 , p2 ] =  distPoints p1 $ p2
solve arr =  rotCal arr' 0 pa pb ( V 1 0 ) ( V (-1) 0 ) dia   n where
	   y1 = minimumBy ( on  compare fN  ) arr
	   y2 = maximumBy ( on  compare fN  ) arr
	   pa = fromJust . findIndex (  == y1 ) $ arr
	   pb = fromJust . findIndex (  == y2 ) $ arr
	   dia = distPoints ( arr !! pa ) ( arr !! pb )
	   n = length arr
	   arr' = listArray ( 0 , n ) arr
	   fN ( P x y ) = y 

 --end of rotating caliper 

--spoj code for testing but time limit is very tight so hard to get accepted in Haskell 

final :: ( Num a , Ord a , Floating a ) => [ Point a ] -> a
final [] = 0
final [ p ] = 0
final [ p1 , p2 ] =  distPoints p1 $ p2
final arr = solve . convexHull $ arr

format :: ( Num a , Ord a , Floating a ) => [ Int ] -> [ [ Point a ]]
format [] = []
format (x:xs ) =  t : format b  where
	( a , b ) = splitAt ( 2 * x ) xs
	t = helpFormat a where
	    helpFormat [] = []
	    helpFormat ( x' : y' : xs' ) = ( P ( fromIntegral x' ) ( fromIntegral  y' ) ) : helpFormat xs'

readD :: BS.ByteString -> Int
readD = fst . fromJust . BS.readInt

main = BS.interact $ BS.unlines . map  ( BS.pack . ( printf "%.0f" :: Double -> String ) .final ) . format . concat . map  ( map  readD . BS.words ) . tail  . BS.lines  

--end of spoj code

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

SPOJ 4186. Break a New RSA system

SPOJ 4186. Break a New RSA system is related to searching a root in given interval. We have
n = p * q * r \hspace{ 20 pt} (1)
\phi_{n} = ( p - 1 ) * ( q - 1 ) * ( r - 1 ) \hspace{ 20 pt} (2)
\sigma_{ n }= ( p + 1 ) * ( q + 1 ) * ( r + 1 )  \hspace{ 20 pt} (3)
Adding 2 and 3
\frac {\phi_{n} + \sigma_{n}} {2} = p + q + r + n  \hspace{ 20 pt} (4)
Subtracting 2 and 3
\frac {\sigma_{n} - \phi_{n}  } {2} = qr + pr + pq + 1  \hspace{ 20 pt} (5)
In 5 , we can replace qr by \frac {n} {p} and take p common from second and third so our final equation will be
\frac {\sigma_{n} - \phi_{n}  } {2} = \frac {n} {p} + p(r + q ) + 1  \hspace{ 20 pt} (6)
From equation 4 , replacing the valuve r + q in 6 we get a cubic equation in p . Our final equation will be
2*p^3 + (2 * n - \phi_{n} - \sigma_{n})*p^2 + (  \sigma_{n} - \phi_{n} - 2 ) * p + (-2n) = 0 \hspace{ 20 pt} (7) .If we differentiate equation 7 and equate to 0 ,two roots of quadratic equation are either local minima and maxima or vice-verse of equation 7 and one of the root of equation 7 will be in this interval. We can search this root by binary search [ why ? because either the function will be increasing or decreasing in this interval ]. After getting one root , divide the equation 7 by this root [ see factorisation ] and get the other 2 roots of quadratic equation by closed form . The time limit for this problem is very tight so Haskell solution got time limit exceed while c++ passed. Accepted code

#include<cstdio>
#include<iostream>
#include<cmath>
using namespace std;
typedef long long ll;

ll fuN( ll a , ll b , ll c , ll d , ll x ) 
	{

		return ( ( a * x + b ) * x + c ) * x + d ;  
	}
  
ll bsearch( ll a , ll b , ll c , ll d , ll lo , ll hi )
	{
		ll mid ;
		while( lo < hi ) 
		{
			mid = ( lo + hi ) >> 1 ; 
			ll f = fuN ( a , b, c , d , lo );
			ll s = fuN ( a , b, c, d , mid ) ; 
			if( ( f > 0 && s > 0 ) || ( f < 0 && s < 0) ) lo = mid + 1 ;
			else hi = mid ;
		}

		return lo;
	}

void solve ( ll n , ll phi , ll sigma ) 
	{
		
		ll a = 2 , b = 2 * n - phi - sigma , c = sigma -  phi - 2 , d = -2 * n ;
		ll dis = ( ll ) sqrtl( ( 4 * b * b - 12 * a * c ) * 1.0 ) ;
		ll low =( ll )  ceill( ( -2 * b - dis ) / ( 6.0 * a ) );
		ll high = ( ll ) floor ( ( -2 * b + dis ) / ( 6.0 * a ) );
		ll r = bsearch (  a , b , c , d ,  low , high ) ; 
		ll d1 = -( b + r * a );
        	ll d2 = ( ll ) sqrtl ( b * b  - 4 * a * c - 2 * a * b * r - 3 * a * a * r * r ) ;
		ll r1 = ( ll ) ( ( d1 - d2 ) / ( 2 * a ) ) ;
        	ll r2 =  ( ll ) ( ( d1 + d2 ) / ( 2 * a ) ) ;
		printf("%lld %lld %lld\n", r1 , r , r2 );

	}

int main()
	{

		int t;
		ll n , phi , sigma ; 
		for(scanf("%lld",&t ) ; t > 0 ; t--) 
		 {
			scanf("%lld %lld %lld", &n , &phi , &sigma );
			solve ( n , phi , sigma );
		 }
	}

Haskell code which got time limit exceed. If your Haskell solution is accepted then please let me know.
import Data.List
import Data.Bits
import Data.Maybe
import qualified Data.ByteString.Char8 as BS

bSearch :: ( Num a , Integral a , Bits a ) => [ a ] -> a -> a -> a
bSearch [ a , b , c , d ] low high = helpBsearch low high where 
	helpBsearch lo hi 
          | lo >= hi = lo  
	  | fuN lo * fuN mid > 0 = helpBsearch ( succ mid ) hi
	  | otherwise = helpBsearch lo mid
	  where mid = shiftR ( lo + hi )  1
		fuN x = ( ( a * x + b ) * x + c ) * x + d  



solve :: ( Num a , Integral  a , Bits a ) => [ a ] ->  BS.ByteString
solve [n , phi , sigma ] = roots    where 
	( a , b , c , d )  = ( 2 , 2 * n - phi - sigma ,  sigma - phi - 2 , -2 * n )
	dis =   sqrt . fromIntegral $ 4 * b ^ 2 - 12 * a * c 
	low =ceiling $ ( ( fromIntegral $  -2 * b  )  - dis  ) / ( fromIntegral $  6 * a  )
	high = ceiling $  ( ( fromIntegral $  -2 * b  )   +  dis  ) / ( fromIntegral $  6 * a  )
	r = bSearch [ a , b , c , d ]  low  high 
	d1 = -( b + r * a ) 
	d2 = truncate.sqrt.fromIntegral $ b ^ 2 - 4 * a * c - 2 * a * b * r - 3 * a ^ 2 * r ^ 2 
	r1 = div ( d1 - d2 ) ( 2 * a ) 
	r2 = div ( d1 + d2 ) ( 2 * a )
	roots =  BS.append  ( BS.append ( BS.append ( BS.append ( BS.pack.show $ r1 ) ( BS.pack " " ) ) ( BS.pack.show $ r ) ) ( BS.pack " " ) ) ( BS.pack.show $ r2 ) 


readInt :: BS.ByteString -> Integer
readInt = fst . fromJust . BS.readInteger  


main = BS.interact $ BS.unlines . map ( solve . map readInt . BS.words ) . tail . BS.lines  

August 1, 2011 Posted by | Programming | , , , | Leave a Comment

SPOJ 6044. Minimum Diameter Circle [ Enclosing Circle Problem ]

SPOJ 6044. Minimum Diameter Circle problem is also know as enclosing circle problem. Accoriding to Wikipedia ” The smallest circle problem or minimum covering circle problem is a mathematical problem of computing the smallest circle that contains all of a given set of points in the Euclidean plane. The corresponding problem in n-dimensional space, the smallest bounding sphere problem, is to compute the smallest n-sphere that contains all of a given set of points. The smallest circle problem was initially proposed by the English mathematician James Joseph Sylvester in 1857″. Couple of algorithm explained here and implementation Pr. Chrystal algorithm. Accepted Haskell code. You can run this code on ideone.

This algorithm is due to Pr. Chrystal (1885).

The Algorithm

Foremost, observe that the minimal enclosing circle, MEC, is entirely determined by the Convex Hull of a given set of point. The reason is the points of the set touched by the minimal enclosing circle are always on the convex hull of the set.

First step of the algorithm is to compute the convex hull, H, of the points in linear time. Clearly, we can do this since points are kept ordered by x-coordinate. Call the convex hull H and the number of convex hull vertices h. Next step is to pick any side, say S, of the convex hull. The main loop of the algorithm is as follows:

Main Loop

For each vertex of H other than those of S, we compute the angle subtended by S. The minimum such angle occurs at vertex v, and the value of the minimum angle is α (alpha).
1 . IF  α  ≥  90 degrees THEN We are done!
    /* The minimal enclosing circle is determined  by the  diametric circle of S */
    IF α < 90 degrees THEN
      Goto step 2 below.
2.  Since α < 90 degrees, check the remaining vertices of the triangle formed by S and v.
    IF no vertices are obtuse THEN We are done!
/* The MEC is determined by the vertices of S and the vertex v */
3. If one of the other angles of the triangle formed by S and v is obtuse, then set S to be the side opposite the obtuse angle and repeat the main loop of the algorithm (i.e., Goto step 1 above).


import Data.List
import qualified Data.Sequence as DS
import Text.Printf
import Data.Function ( on ) 

data Point a = P a a deriving ( Show , Ord , Eq )
data Turn = S | L | R deriving ( Show , Eq , Ord , Enum  ) -- straight left right  

--start of convex hull  http://en.wikipedia.org/wiki/Graham_scan
{--
compPoint :: ( Num  a , Ord a ) => Point a -> Point a -> Ordering
compPoint ( P x1 y1 ) ( P x2 y2 )
  | compare x1 x2 == EQ = compare y1 y2
  | otherwise = compare x1 x2 

findMinx :: ( Num a , Ord a ) => [ Point a ] -> [ Point a ]
findMinx xs = sortBy ( \x  y  -> compPoint  x y  ) xs

compAngle ::(Num a , Ord a ) => Point a -> Point a -> Point a -> Ordering
compAngle ( P x1 y1 ) ( P x2 y2 ) ( P x0 y0 ) = compare ( (  y1 - y0 ) * ( x2 - x0 )  ) ( ( y2 - y0) * ( x1 - x0 ) )

sortByangle :: ( Num a , Ord a ) => [ Point a ] -> [ Point a ]
sortByangle (z:xs) = z : sortBy ( \x y -> compAngle x y z ) xs 

convexHull ::( Num a , Ord a )	=> [ Point a ] -> [ Point a ]
convexHull xs = reverse . findHull [y,x]  $ ys where
	(x:y:ys) = sortByangle.findMinx $ xs 

findTurn :: ( Num a , Ord a , Eq a ) => Point a -> Point a -> Point a -> Turn
findTurn ( P x0 y0 ) ( P x1 y1 ) ( P x2 y2 )
 | ( y1 - y0 ) * ( x2- x0 ) < ( y2 - y0 ) * ( x1 - x0 ) = L
 | ( y1 - y0 ) * ( x2- x0 ) == ( y2 - y0 ) * ( x1 - x0 ) = S
 | otherwise = R 

findHull :: ( Num a , Ord a  )  => [ Point a ] ->   [ Point a ] -> [ Point a ]
findHull [x]  ( z : ys )  = findHull [ z , x ]  ys  --incase of second point  on line from x to z
findHull xs  [] = xs
findHull ( y : x : xs )  ( z:ys )
  | findTurn x y z == R = findHull (  x : xs )   ( z:ys )
  | findTurn x y z == S = findHull (  x : xs )   ( z:ys )
  | otherwise = findHull ( z : y : x : xs  )   ys
--}

--start of monotone hull give .04 sec lead over graham scan 

compPoint :: ( Num  a , Ord a ) => Point a -> Point a -> Ordering
compPoint ( P x1 y1 ) ( P x2 y2 )
  | compare x1 x2 == EQ = compare y1 y2
  | otherwise = compare x1 x2 

sortPoint :: ( Num a , Ord a ) => [ Point a ] -> [ Point a ]
sortPoint xs = sortBy ( \ x y -> compPoint x y ) xs

findTurn :: ( Num a , Ord a , Eq a ) => Point a -> Point a -> Point a -> Turn
findTurn ( P x0 y0 ) ( P x1 y1 ) ( P x2 y2 )
 | ( y1 - y0 ) * ( x2- x0 ) < ( y2 - y0 ) * ( x1 - x0 ) = L
 | ( y1 - y0 ) * ( x2- x0 ) == ( y2 - y0 ) * ( x1 - x0 ) = S
 | otherwise = R 


hullComputation :: ( Num a , Ord a ) => [ Point a ] -> [ Point a ] -> [ Point a ]
hullComputation [x] ( z:ys ) = hullComputation [z,x] ys
hullComputation xs [] = xs 
hullComputation  ( y : x : xs ) ( z : ys ) 
  |  findTurn x y z == R = hullComputation ( x:xs ) ( z : ys ) 
  |  findTurn x y z == S = hullComputation ( x:xs ) ( z : ys ) 
  |  otherwise = hullComputation ( z : y : x : xs ) ys 


convexHull :: ( Num a , Ord a ) => [ Point a ] -> [ Point a ]
convexHull xs = final  where 
	txs = sortPoint xs 
	( x : y : ys  ) = txs 
        lhull = hullComputation [y,x] ys  
	( x': y' : xs' ) = reverse txs 
	uhull = hullComputation [ y' , x' ] xs' 
	final = nub $  ( reverse lhull )  ++  uhull 

--end of convex hull
--start of finding point algorithm http://www.personal.kent.edu/~rmuhamma/Compgeometry/MyCG/CG-Applets/Center/centercli.htm  Applet’s Algorithm 

findAngle :: ( Num a , Ord a , Floating a ) => Point a -> Point a -> Point  a  -> ( Point a , Point a , Point  a , a )
findAngle u@(P x0 y0 ) v@(P x1 y1 ) t@(P x2 y2)
	| u == t || v == t = ( u , v , t , 10 * pi )  -- two points are same so set the angle more than pi
	| otherwise =  ( u , v, t , ang ) where
          	ang = acos $  ( b + c - a ) / ( 2.0 * ( sqrt $ b * c ) )  where
			sqrDist ( P x y ) ( P x' y' ) = ( x - x' ) ^ 2 + ( y - y' ) ^ 2
			[ b , c , a ] = map ( uncurry sqrDist ) [ ( u , t ) , ( v , t ) , ( u , v ) ]

findPoints :: ( Num a , Ord a , Floating a ) => Point a -> Point a -> [ Point  a ] -> ( Point a  , a )
findPoints u v xs
  | 2 * theta >= pi   =  circle2Points a b 	--( a , b , t , theta )
  | and [ 2 * alpha <= pi , 2 * beta <= pi ]   = circle3Points a b t --( a , b , t , theta )
  | otherwise = if 2 * alpha > pi then findPoints v t xs else findPoints u t xs
     where
	( a , b , t , theta ) = minimumBy ( on compare fn  ) . map ( findAngle u v )  $ xs
	fn ( _ , _ , _ , t ) = t
        ( _ , _ , _ , alpha ) = findAngle v t u  --angle between v u t angle subtended at u by v t
	( _ , _ , _ , beta ) = findAngle u t v   -- angle between u v t angle subtended at v by  u t

--end of finding three points algorithm
--find the circle through three points http://paulbourke.net/geometry/circlefrom3/ 

circle2Points :: ( Num a , Ord a , Floating a ) => Point a -> Point a -> ( Point a , a )
circle2Points ( P x0 y0 ) ( P x1 y1 ) = ( P x y , r ) where
	x = ( x0 + x1 ) / 2.0
	y = ( y0 + y1 ) / 2.0
	r = sqrt $ ( x0 - x1 ) ^ 2 + ( y0 - y1 ) ^ 2 

circle3Points :: ( Num a , Ord a , Floating a ) => Point a -> Point a -> Point a -> ( Point a , a ) --( center , radius )
circle3Points u@(P x1 y1 ) v@(P x2 y2 ) t@(P x3 y3 )
	| x2 == x1 = circle3Points u t v
	| x3 == x2 = circle3Points v u t
	| otherwise =  ( P x y , 2 *  r )
	  where
		m1 = ( y2 - y1 ) / ( x2 - x1 )
		m2 = ( y3 - y2 ) / ( x3 - x2 )
		x = (  m1 * m2 * ( y1 - y3 )  +   m2 * ( x1 + x2 )  -  m1 * ( x2 + x3 )  ) / ( 2 * ( m2 - m1 ) )
		y = if y2 /= y1
                     then ( ( x1 + x2 -  2 * x  ) /  ( 2 * m1 ) ) + ( ( y1 + y2 ) / 2.0 )
		      else  ( ( x2 + x3 -  2 * x  ) / ( 2 * m2 ) ) + ( ( y2 + y3 ) / 2.0 )
		r = sqrt $ ( x - x1 ) ^2 + ( y - y1 ) ^ 2 

--start of SPOJ code 

format::(Num a , Ord a ) => [[a]] -> [Point a]
format xs = map (\[x0 , y0] -> P x0 y0 ) xs 

readInt  ::( Num a , Read a ) =>   String -> a
readInt  = read

solve :: ( Num a , Ord a , Floating a ) => [ Point a ] -> (  Point a , a )
solve [ P x0 y0 ] = ( P x0 y0 , 0 ) --in case of one point
solve [ P x0 y0 , P x1 y1 ] = circle2Points ( P x0 y0 ) ( P x1 y1 )   -- in case of two points
solve  xs = findPoints x y  t where
	t@( x : y : ys )  = convexHull xs  

final :: ( Num a , Ord a , Floating a ) => ( Point a , a ) -> a
final ( t , w )
	| w == 0 = 0
	| otherwise =  w

main = interact $   ( printf "%.2f\n" :: Double -> String ) . final . solve . convexHull . format . map  ( map readInt . words ) . tail . lines

July 24, 2011 Posted by | Programming | , , , , | Leave a Comment

Follow

Get every new post delivered to your Inbox.

Join 138 other followers