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
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
readInt = read
main = interact $ solve . convexHull . format . map ( map readInt . words ) . tail . lines
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
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 Control.Monad
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 ) )
--add points of elliptic curve
addPoints::Elliptic->Point->Point-> Either Point Integer
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
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 
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
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 
#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");*/
}
}
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
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 
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 
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)
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 readInt :: BS.ByteString -> Integer readInt = fst . fromJust . BS.readInteger 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
SPOJ 7897. Skyline
SPOJ 7897. Skyline is related to Catalan numbers.From wiki ” 
import Data.List import Data.Ratio import Control.Monad import Data.Maybe import qualified Data.ByteString.Char8 as BS import Data.Array readInt :: BS.ByteString -> Integer readInt = fst.fromJust.BS.readInteger 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
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
import Control.Monad
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
About
Hello All, My name is Mukesh Tiwari and I graduated from IIITM Gwalior in 2009. I am passionate about programming and mathematics. I love solving problems on SPOJ , UVA , Topcoder and Project Euler. I am curious about functional languages specially Haskell and it’s one the excellent language I encountered after python. I am also looking for some challenging functional programming job specially related to Haskell. You can reach me mukeshtiwariDOTiiitmATgmailDOTcom. Replace DOT with . and AT with @
Recent Comments
