# My Weblog

## 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

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 ## Monotone Chain Convex Hull Monotone Chain Convex Hull algorithm does not need angular sorting. Wikipedia suggest ” The same basic idea works also if the input is sorted on x-coordinate instead of angle, and the hull is computed in two steps producing the upper and the lower parts of the hull respectively. This modification was devised by A. M. Andrew and is known as Andrew’s Monotone Chain Algorithm. It has the same basic properties as Graham’s Scan but eschews costly comparisons between polar angles “. Also tested this code on SPOJ 3421. Garden Hull and accepted . Also have look at Algorithmist. import Data.List data Point a = P a a deriving ( Show , Ord , Eq ) data Turn = S | L | R deriving ( Show , Ord , Eq , Enum ) --start of monotone 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 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 monotone

--from here on testing part for SPOJ

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

helpSqrt :: (  Floating  a ) => Point a -> Point a ->  a
helpSqrt ( P x0 y0 ) ( P x1 y1 ) =  sqrt  $( x0 - x1 ) ^ 2 + ( y0 - y1 ) ^ 2 solve :: ( Num a , RealFrac a , Floating a ) => [ Point a ] -> String solve xs = ( show . fromIntegral . truncate . snd . foldl ( \( P x0 y0 , s ) ( P x1 y1 ) -> ( P x1 y1 , s + helpSqrt ( P x0 y0 ) ( P x1 y1 ) ) ) ( head xs , 0 )$  ( tail xs  ++ [ head xs ] ) ) ++ "\n"

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

main = interact $solve . convexHull . format . map ( map readInt . words ) . tail . lines  July 20, 2011 ## Convex Hull in Haskell Implementation Convex hull using Graham scan algorithm and tested it on SPOJ 3421. Garden Hull. Accepted Haskell code.  import Data.List import qualified Data.Sequence as DS data Point a = P a a deriving ( Show , Eq , Ord ) data Turn = S | L | R deriving ( Show , Eq , Ord , Enum ) -- straight left right 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 --from here on testing part for SPOJ format::(Num a , Ord a ) => [[a]] -> [Point a] format xs = map (\[x0 , y0] -> P x0 y0 ) xs helpSqrt :: ( Floating a ) => Point a -> Point a -> a helpSqrt ( P x0 y0 ) ( P x1 y1 ) = sqrt$  ( x0 - x1 ) ^ 2 + ( y0 - y1 ) ^ 2

solve :: ( Num a , RealFrac a , Floating a  ) => [ Point a ] -> String
solve xs = ( show . fromIntegral . truncate  .  snd . foldl ( \(  P x0 y0  , s )  ( P x1 y1 ) -> ( P x1 y1   , s + helpSqrt  ( P  x0 y0  ) ( P x1 y1 ) ) )  ( head xs , 0 ) $( tail xs ++ [ head xs ] ) ) ++ "\n" readInt ::( Num a , Read a ) => String -> a readInt = read main = interact$ solve . convexHull . format . map  ( map readInt . words ) . tail . lines



July 20, 2011

## Improved Elliptic Curve prime factorisation

Currently i am trying to understand monad. I wrote elliptic curve prime factorisation some times ago and now using monads, i wrote a bit improved elliptic curve prime factorisation . This code factorise the product of prime numbers in range [1000..2000] withing less than second on my slow computer.Algorithm explained is previous post. Modified Haskell code.

--y^2= x^3+ax+b mod n
import Random
import Data.List
import Data.Bits

data Elliptic = Conelliptic Integer Integer Integer deriving ( Eq , Show )
data Point = Conpoint Integer Integer | Identity deriving ( Eq , Show )

powM :: Integer -> Integer -> Integer -> IO Integer
powM a d n = return $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 -> IO ( Integer , Integer ) calSd n = return ( s , d ) where s = until ( \x -> testBit ( n - 1) ( fromIntegral x ) ) ( +1 ) 0 d = div ( n - 1 ) ( shiftL 1 ( fromIntegral s ) ) isProbable::Integer->IO Bool isProbable n | n <= 1 = return False | n == 2 = return True | even n = return False | otherwise = calSd n >>= (\( s , d ) -> rabinMiller 0 n s d ) rabinMiller::Integer->Integer->Integer->Integer->IO Bool rabinMiller cnt n s d | cnt>=5 = return True | otherwise = randomRIO ( 2 , n - 2 ) >>= (\a -> powM a d n >>= (\x -> case x == 1 of True -> rabinMiller ( cnt + 1 ) n s d _ -> if any ( == pred n ) . take ( fromIntegral s ) . iterate (\e -> mod ( e^2 ) n )$ x
then rabinMiller ( cnt + 1) n s d
else return False  ) )

addPoints _ Identity p_2 = Left p_2
addPoints _ p_1 Identity = Left p_1
addPoints ( Conelliptic a b n ) ( Conpoint x_p y_p ) ( Conpoint x_q y_q )
| x_p /= x_q = let
[ u , v , d ] = extended_gcd (x_p-x_q) n
s = mod  ( ( y_p - y_q ) * u ) n
x_r = mod ( s * s - x_p - x_q ) n
y_r= mod ( -y_p - s * ( x_r - x_p ) ) n
in case ( ( Conpoint x_r y_r ) , d ) of
( _ , 1 ) -> Left ( Conpoint x_r y_r )
( _ , d' ) -> Right d'
| otherwise = if mod ( y_p + y_q ) n == 0 then Left Identity
else  let
[ u , v , d ] = extended_gcd ( 2*y_p ) n
s = mod  ( ( 3 * x_p * x_p + a ) * u ) n
x_r = mod ( s * s - 2 * x_p ) n
y_r = mod ( -y_p - s * ( x_r - x_p ) ) n
in case ( ( Conpoint x_r y_r ) , d ) of
( _ , 1 )-> Left (Conpoint x_r y_r)
( _ , d' ) -> Right d'

extended_gcd::Integer->Integer->[Integer]
extended_gcd u v= helpGcd u v 0 [ [ 1 , 0 ] , [ 0 , 1 ] ] where
helpGcd u v n m @ ( [ [ a , b ] , [ c , d ] ] )
| v == 0 = if u<0 then [ - ( ( -1 ) ^ n ) * ( m !! 1 !! 1 ) , - ( ( -1 ) ^ ( n + 1 ) ) * ( m !! 0 !! 1 ) , -u ]
else [ ( ( -1 ) ^ n ) * ( m !! 1 !! 1 ) , ( ( -1 ) ^ ( n + 1 ) ) * ( m !! 0 !! 1 ) , u ]
| otherwise = helpGcd u' v' ( n + 1 ) m' where
( q , v' ) = quotRem u v
t = [ [q , 1 ] , [ 1 , 0 ] ]
m' = [ [ q * a + b , a ] , [ q * c + d , c ] ] --mult m t
u' = v

multiEC :: Elliptic -> Point -> Integer -> IO ( Either Point Integer )
multiEC _ _ 0 = return $Left Identity multiEC ecurve point k = return$ helpEC Identity point k
where
helpEC p _ 0 = Left p
helpEC p q n =
case (p',q') of
( Left p'' , Left q'' ) -> helpEC p'' q'' (div n 2)
( Right p'' , _ ) -> Right p''
( _ , Right q'' ) -> Right q''
where
p' = if odd n then addPoints ecurve p q else Left p
q' =  addPoints ecurve q q

dscrmntEC a b = 4 * a * a * a + 27 * b * b

randomCurve :: Integer -> IO [Integer]
randomCurve n =
randomRIO ( 1 , n ) >>=
( \a -> randomRIO ( 1 , n ) >>=
( \u -> randomRIO (1 , n) >>=
(\v -> return [a , mod ( v*v - u*u*u - a*u ) n , n , u , v ] ) ) )

factor :: Integer -> Integer -> IO [Integer]
factor 1 _   = return []
factor n m   =
isProbable n >>= (\x -> if x then  return [n]
else
randomCurve n >>= (\[ a , b , n , u , v ] ->
multiEC ( Conelliptic a b n ) ( Conpoint u v ) m >>=
( \p -> case p of
Left p' -> factor n m
Right p'->  factor ( div n p' ) m  >>= ( \x -> factor p' m  >>= (\y  -> return $x ++ y ) ) ) ) ) solve :: Integer -> IO [ Integer ] solve n = factor n ( foldl lcm 1 [1..10000] ) main = liftM read getLine >>=( \n -> solve n ) >>= print  [user@localhost Haskell]$ time ./Elliptic_Len
14771291995722167329218042370113255785101138530361517229308539087571387624569184860194674818433642463913959873342852255543111135151457209398370030102228858128758106619882182645353368949658870563576166323307466834284670054906895828663228915811739807554988946409285249055824035034910383070200873326731464685909980630202533115332294811841112677175465781631555282945350555352712075708686158031835188673039526232167090966030801645637
[1759,1997,1291,1637,1669,1889,1907,1823,1231,1579,1597,1327,1427,1499,1831,1951,1213,1753,1993,1487,1429,1601,1319,1447,1721,1847,1697,1873,1699,1663,1061,1879,1087,1093,1117,1423,1933,1493,1303,1021,1109,1973,1361,1811,1123,1543,1249,1783,1867,1301,1789,1901,1657,1693,1801,1091,1609,1913,1277,1787,1483,1283,1511,1381,1877,1297,1181,1571,1103,1229,1039,1471,1289,1373,1861,1069,1987,1481,1279,1559,1193,1129,1747,1949,1733,1489,1063,1201,1307,1187,1871,1031,1741,1433,1321,1013,1567,1627,1549,1399,1553,1523,1051,1439,1723,1259,1667,1621,1777,1151,1367,1019,1097,1451,1163,1223,1999,1979,1931,1583,1619,1217,1153,1009,1409,1613,1459,1171,1531,1453,1237,1607,1049,1709,1033]

real	0m16.424s
user	0m1.905s
sys	0m0.009s


July 18, 2011

## 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 )

main = BS.interact $BS.unlines.map ( solve .readInt ) . tail . BS.lines  July 17, 2011 ## SPOJ 8456. PRIMITIVEROOTS SPOJ 8456. PRIMITIVEROOTS is related to number theory. Product of primitive roots of p ( prime number ) mod p is equal to 1.The product of all primitive roots of prime p $\not=$ 3 is $\equiv 1$ ( mod p ). History of the theory of numbers By Leonard Eugene Dickson on page 183.This problem is not allowed in Haskell. Accepted C++ code. #include<cstdio> #include<iostream> #include<cmath> using namespace std; bool p[10010]; void prime() { for(int i=2 ; i*i <=10000 ; i++ ) if(!p[i]) for(int j = i*i ; j <=10000 ; j += i ) p[j]=1; //for(int i=2; i <= 20 ; i++) if(!p[i]) cout<<i<<" ";cout<<endl; } int main() { prime(); int n , t , x=1 ; for( scanf("%d",&t) ; t>0 ; t--) { scanf("%d",&n); if( n == 3 ) printf("%d:2\n", x++); else if ( !p[n] ) printf("%d:1\n",x++); else printf("%d:NOTPRIME\n",x++); /*if( t != 1 ) printf("\n");*/ } }  July 17, 2011 ## SPOJ Power Sums SPOJ Power Sums accepted after couple of attempts. This problem is similar to SRM 397 SumOfPowers . This discussion on topcoder forum shows couple of ideas to solve this problem using Gaussian elimination , recursive sum and Bernoulli numbers . First i implemented Gaussian elimination mod prime . The idea is $\sum_{i=1}^x i^k = a_{k+1}x^{k+1}+a_kx^k+....+a_1x+a_0$ . For this equation , we have k+2 unknown coefficients so we need k+2 equations to solve this. Putting x from 0 to ( k + 1 ) and we can have k + 2 equation . Solve this system of equation using Gaussian elimination mod prime . $\begin{pmatrix} x_0^{k+1} & x_0^k & . & . & x_0 & 1 \\ x_1^{k+1} & x_1^k & . & . & x_1 & 1 \\ x_2^{k+1} & x_2^k & . & . & x_2 & 1 \\ . & . & . & . & . & . \\ . & . & . & . & . & . \\ x_{k+1}^{k+1} & x_{k+1}^k & . & . & x_{k+1} & 1\end{pmatrix} * \begin{pmatrix} a_{k+1} \\ a_k \\ . \\ . \\ a_1 \\ a_0 \end{pmatrix} = \begin{pmatrix} \sum_{i=1}^{x_0} i^k \\ \sum_{i=1}^{x_1} i^k \\ \sum_{i=1}^{x_2} i^k \\ . \\ .\\ \sum_{i=1}^{x_{k+1}} i^k \end{pmatrix}{\mod p }$ Haskell code for this implementation but it got time limit exceed . This code can also be useb for computation of determinant mod primes and combine the result using chinese remainder theorem. import Data.List extendedGcd :: Integer -> Integer -> ( Integer , Integer ) extendedGcd a b | b == 0 = ( 1 , 0 ) | otherwise = ( t , s - q * t ) where ( q , r ) = quotRem a b ( s , t ) = extendedGcd b r modInv :: Integer -> Integer -> Integer modInv a b | gcd a b /= 1 = error " gcd is not 1 " | otherwise = d where d = until ( > 0 ) ( + b ) . fst.extendedGcd a$ b

adjustwitH::[ Integer ] -> [ Integer ] -> [ Integer ]
adjustwitH ( m:ms ) ( n:ns ) = map ( flip mod 10007 )  $zipWith ( - ) ms$ map ( ( m * modInv n 10007 )* )  ns --m/n mod p modular inverse

echeloN::[[ Integer ]]->[[ Integer ]]
echeloN rs
| null rs || null ( head rs) = rs
| null rs_2 = map (0:) ( echeloN $map tail rs ) --haskell road to logic | otherwise = piv : map ( 0:) ( echeloN rs' ) where rs'= map ( adjustwitH piv ) ( rs_1 ++ rs_3 ) ( rs_1 , rs_2 ) = span ( \( x:xs )-> x == 0 ) rs ( piv:rs_3 ) = rs_2 rows , cols :: [[Integer ]] -> Int rows m = length m cols m = length.head$ m

backS :: [[Integer]] -> [Integer]
backS rs = backS' rs [] where
backS' [] ps = ps
backS' rs ps = backS' rs' ( p:ps ) where
a = ( last rs ) !! ( ( cols rs ) - 2 )
c = ( last rs ) !! ( ( cols rs ) - 1 )
p = mod ( c * modInv a 10007 ) 10007
rs' = eliminate p ( init rs )

eliminate :: Integer -> [[Integer]] -> [[Integer]]
eliminate p rs = map ( simplify p ) rs where
simplify p row = init ( init row ) ++ [ mod ( d - b * p) 10007 ] where
d = last row
b = last ( init row )

helpMod::Integer->Integer->Integer->Integer
helpMod a d n= fun a d where
fun a 0 = 1
fun a 1 = mod a n
fun a 2 = mod (a^2) n
fun a d = let
p=fun (mod (a^2) n) (div d 2)
q=if odd d then mod (a*p) n else p
in mod q n

modS :: Integer -> Integer -> Integer
modS p k = mod ( sum [ helpMod i k 10007 | i <- [1..p]] )  10007

solve :: Integer -> [Integer]
solve k =  backS . echeloN $[ [ helpMod i n 10007 | n <- [ ( k + 1 ) , k .. 0 ] ] ++ [ modS i k ] | i <- [1 .. (k+2)] ] format ::Integer -> [Integer] -> String format k [x] = show x ++ "x^" ++ show k format k t@(x:xs) = ( show x ++ "x^" ++ show k ++ " + ") ++ format ( k - 1 ) xs final :: Integer -> String final k = format ( k + 1) ( solve k ) main = interact$ unlines. map ( final.read ) . lines


My second attempt was recursive formula $( m + 1 )\sum_{k=0}^m = ( n + 1 ) ^{m+1} - \sum_{r=2}^{m+1} ( {{n}\choose{r}} *\sum_{k=0}^n k^{m+1-r} )$ but unfortunately this one also time limit exceed .

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

polySum :: [Integer] -> [Integer] -> [Integer]
polySum ( x : xs ) ys =map ( flip mod 10007 ) $x : zipWith ( + ) xs ys polySub :: [Integer] -> [Integer] -> [Integer] polySub (x:xs) ys = map ( flip mod 10007 )$ x : zipWith ( - ) xs ys

binomM :: Integer -> Integer -> Integer -> Integer
binomM n k p =  pas !! ( fromIntegral n) !! ( fromIntegral k )

pas :: [[Integer]]
pas = iterate ( \x -> map ( flip mod 10007 ) $1 : zipWith (+) x ( tail x ) ++ [1] ) [1] extendedGcd :: Integer -> Integer -> ( Integer , Integer ) extendedGcd a b | b == 0 = ( 1 , 0 ) | otherwise = ( t , s - q * t ) where ( q , r ) = quotRem a b ( s , t ) = extendedGcd b r modInv :: Integer -> Integer -> Integer modInv a b | gcd a b /= 1 = error " gcd is not 1 " | otherwise = d where d = until ( > 0 ) ( + b ) . fst.extendedGcd a$ b

sum' ::[[Integer]] -> [Integer]
sum' xs = foldr ( \x y -> polySum x y ) [0] xs

solve :: Integer -> [[Integer]]
solve m = foldl ( \xs y -> ( polySub ( [ modInv ( y + 1) 10007 * binomM ( y + 1 ) i 10007 | i <- [ 0 .. ( y + 1 ) ] ] ) ( sum' $zipWith map [ ( ( modInv ( y + 1 ) 10007 * binomM ( y + 1 ) r 10007 )* ) | r <- [ 2 .. ( y + 1 ) ] ] xs ) ) : xs ) [[1,1]] [ 1 .. m ]  Finally attempt was using Bernoulli numbers and bingo. Sum is given by $S_m(n) = {1\over{m+1}}\sum_{k=0}^m {m+1\choose{k}} B_k n^{m+1-k}$ and Bernoulli numbers can be computed using $B_m= 1 - \sum_{k=0}^{m-1}\binom mk\frac{B_k}{m-k+1}$.Accepted Haskell code . Learn lot while solving this problem. import Data.List import Data.Maybe import qualified Data.ByteString.Char8 as BS readInt :: BS.ByteString -> Integer readInt = fst . fromJust . BS.readInteger binomM :: Integer -> Integer -> Integer -> Integer binomM n k p = pas !! ( fromIntegral n ) !! ( fromIntegral k ) pas :: [[Integer]] pas = iterate ( \x -> 1 : zipWith (\ a b -> mod ( a + b ) 10007 ) x ( tail x ) ++ [1] ) [1] extendedGcd :: Integer -> Integer -> ( Integer , Integer ) extendedGcd a b | b == 0 = ( 1 , 0 ) | otherwise = ( t , s - q * t ) where ( q , r ) = quotRem a b ( s , t ) = extendedGcd b r modInv :: Integer -> Integer -> Integer modInv a b | gcd a b /= 1 = error " gcd is not 1 " | otherwise = d where d = until ( > 0 ) ( + b ) . fst.extendedGcd a$ b

ber ::  [Integer]
ber  = foldl ( \ xs y ->  if and [ odd y , y /= 1 ] then xs ++ [0] else  xs ++ [ mod ( 1 -  ( sum $zipWith (*) [ binomM y k 10007 * modInv ( y - k + 1 ) 10007 | k <-[ 0 .. pred y ] ] xs ) ) 10007 ] ) [1] [1..300] {--ber is fast than berIt berIt :: [( [Integer] , Integer )] berIt = iterate ( \( xs , y) -> ( map ( flip mod 10007 )$ if and [ odd y , y /= 1 ] then xs ++ [0] else  xs ++ [ 1 -  ( sum $zipWith (*) [ binomM y k 10007 * modInv ( y - k + 1 ) 10007 | k <-[ 0 .. pred y ] ] xs ) ] , y + 1 ) ) ( [1] , 1 ) --} solve :: Integer -> [Integer] solve m = [ mod ( binomM ( m + 1 ) k 10007 * ber !! ( fromIntegral k ) * modInv ( m + 1 ) 10007 ) 10007 | k <- [0..m] ] final :: Integer -> BS.ByteString final k = format ( k + 1 ) ( solve k ) format::Integer -> [Integer] -> BS.ByteString format _ [] = BS.pack "0x^0" format k ( x : xs ) = BS.append ( BS.append ( BS.append ( BS.append ( BS.pack.show$ x ) ( BS.pack  "x^" ) ) ( BS.pack.show $k ) ) ( BS.pack " + " ) ) ( format ( k - 1 ) xs ) main = BS.interact$ BS.unlines . map ( final . readInt )  . BS.lines


Here is more faster way to generate Bernoulli numbers. Replace ber !! ( fromIntegral k ) in solve m with bernoulli k for improved run time.

bernoulli :: Integer -> Integer
bernoulli = genericIndex bernoulli_table
bernoulli_table = map bernoulli' [0..]

bernoulli' :: Integer ->  Integer
bernoulli' 0 = 1
bernoulli' 1 = 1 * modInv 2 10007
bernoulli' n | odd n          = 0
| n mod 6 == 0 = let val = (n+3) * modInv 3 10007 - bernoulliA n (n div 6) in
mod ( val * modInv ( binomM ( n + 3 ) n 10007 ) 10007 ) 10007 --fromIntegral (choose (n+3) n)
| n mod 6 == 2 = let val = (n+3) * modInv 3 10007 - bernoulliA n ((n-2) div 6) in
mod ( val * modInv  ( binomM ( n + 3 ) n 10007 ) 10007 ) 10007 --fromIntegral (choose (n+3) n)
| n mod 6 == 4 = let val = -(n+3) * modInv 6 10007 - bernoulliA n ((n-4) div 6) in
mod (  val * modInv ( binomM ( n + 3 ) n 10007 ) 10007 ) 10007 -- fromIntegral (choose (n+3) n)

bernoulliA :: Integer -> Integer ->  Integer
bernoulliA n m = mod ( sum (map foo [1..m]) ) 10007
where foo j = binomM  (n+3) (n-6*j) 10007  *  bernoulli (n-6*j)


July 15, 2011

## SPOJ 9161. Legendre symbol

SPOJ 9161. Legendre symbol is related to Jacobian symbol.The Jacobi symbol is a generalization of the Legendre symbol. I just submitted my code from Solovay–Strassen primality test and accepted.Accepted Haskell code

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

jacobi::Integral a => a -> a -> a
jacobi a n
| a == 0 = if n == 1 then 1 else 0
| a == 2 && ( mod n 8 == 1 || mod n 8 == 7 ) = 1
| a == 2 && ( mod n 8 == 3 || mod n 8 == 5 ) = -1
| a >= n = jacobi ( mod a n ) n
| even a = jacobi 2 n * jacobi ( div a 2 ) n
| mod a 4 == 3 && mod n 4 == 3 = -jacobi n a
| otherwise = jacobi n a

legend :: [Integer]  -> BS.ByteString
legend (a:p:_) = BS.pack.show $b where b = jacobi a p main = BS.interact$ BS.unlines . map ( legend . map readInt . BS.words ) . tail . BS.lines



July 12, 2011

## SPOJ 7897. Skyline

SPOJ 7897. Skyline is related to Catalan numbers.From wiki ” $C_n$ is the number of permutations of {1, …, n} that avoid the pattern 123 (or any of the other patterns of length 3); that is, the number of permutations with no three-term increasing subsequence”. Accepted Haskell code .

import Data.List
import Data.Ratio
import Data.Maybe
import qualified   Data.ByteString.Char8 as BS
import Data.Array

catlan :: [Integer]
catlan = 1 : zipWith (\c n ->  div ( c * 2 * ( 2 * n + 1 ) )   ( n + 2 )  ) catlan [1..]

catArray_1 ::  Array Int Integer
catArray_1  = listArray ( 0 , 1000 )  catlan

cat :: Array Integer Integer
cat = array (0,1000) ( [( 0 , 1 ) ] ++ [ ( i , div ( cat ! ( i - 1 ) *  ( 4 * i + 2 )  )  ( i + 2 )  ) | i <- [1..1000] ] )

solve :: Integer -> BS.ByteString
solve n = BS.pack.show $mod ( cat ! ( n - 1 ) ) 1000000 main = BS.interact$ BS.unlines. map (  solve.readInt ) . init . BS.lines



July 11, 2011

## SPOJ 1754. Divisor Summation (Hard)

SPOJ 1754. Divisor Summation (Hard) finally accepted in Haskell. I tried this problem in 2007 and almost my 20 solutions in C/C++ and python got time limit exceed. Nothing special , just implemented the pollard rho prime factorisation in Haskell and accepted. Accepted Haskell code.

import Data.List
import Random

helpMod::Integer->Integer->Integer->Integer
helpMod a d n= fun a d where
fun a 1 = mod a n
fun a 2 = mod (a^2) n
fun a d = let
p=fun (mod (a^2) n) (div d 2)
q=if odd d then mod (a*p) n else p
in mod q n

isProbable::Integer-> Bool
isProbable n |n<=1 = False
|n==2 = True
|even n = False
|otherwise	= rabinMiller [2,3,5,7,11,13,17,19] n s d where
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=helpMod a d n

fn :: Integer -> Integer -> Integer -> Integer
fn x c n = mod ( x^2 + c ) n

pollard :: Integer -> IO  Integer
pollard n
| even n = return  2
| otherwise =
randomRIO ( 2 , n - 1 ) >>=
(\x' -> randomRIO ( 2 , n - 1 ) >>=
(\c ->  return ( until ( \( _ , _ , d ) -> d /= 1 || d == n )  ( \( x , y , d ) -> ( fn x c n , fn ( fn y c n ) c n , gcd ( fn x c n - fn ( fn y c n ) c n  ) n ) )  ( x' , x' , 1 ) ) >>=
( \( _ , _ , k ) -> return k  ) ) )

factor :: Integer -> IO [ Integer ]
factor 1 = return []
factor n =
(return . isProbable $n ) >>= (\x -> case x of True -> return [n] _ -> pollard n >>= (\d -> case d of 1 -> factor n _ -> factor d >>= (\l -> ( factor.div n$ d ) >>= ( \m -> return $l ++ m ) ) ) ) divSum :: [[Integer]] -> Integer divSum [] = 1 divSum (x:xs) = div ( ( head x )^( genericLength x + 1 ) -1 ) ( head x - 1 ) * divSum xs main = do n <- liftM read getLine let helpfun n cnt | cnt == n = return [] | otherwise = do m <- liftM read getLine lst <- factor m putStrLn.show$ ( divSum.groupBy (\x y -> x == y ).sort \$ lst ) - m
helpfun n (cnt+1)

helpfun n 0



July 8, 2011