# Enclosing Circle Problem implementation in Haskell

I implemented the algorithm by Pr. Chrystal described here in Haskell so could someone please tell me: Do i have implemented this algorithm correctly?
Initial calling of findPoints takes the first and second point of convex hull and the list of points on convex hull.

My second question is: I am trying to solve the problem on sphere online judge [link of problem in code below, because I can not post more than 2 links] using this algorithm but I'm getting the wrong answer.
The code for the problems is on ideone. Why am I getting the wrong answer.

--http://www.spoj.pl/problems/QCJ4/
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 * sb * sc ) ) where
b = ( x0 - x2 ) ^ 2 + ( y0 - y2 ) ^ 2
c = ( x1 - x2 ) ^ 2 + ( y1 - y2 ) ^ 2
a = ( x0 - x1 ) ^ 2 + ( y0 - y1 ) ^ 2
sb = sqrt b
sc = sqrt c

findPoints :: ( Num a , Ord a , Floating a ) => Point a -> Point a -> [ Point  a ] -> ( Point a , Point a , Point a , a )
findPoints u v xs
| 2 * theta >= pi   =     ( a , b , t , theta )
| and [ 2 * alpha <= pi , 2 * beta <= pi ]   = ( 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 ( \(_,_,_, t1 ) ( _ , _ , _ ,t2 ) -> compare  t1 t2 ) . map ( findAngle u v )  $xs ( _ , _ , _ , 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  • I have solved this problem. Accepted code of problem QCJ4 spoj on ideone. I hope Mihai don't need to give away his reputation points as bounty to any one :). Jul 29, 2011 at 16:05 • Can you post solution below? (I'll give you the bounty in this case) Aug 2, 2011 at 7:30 • @keep_learning: Now I got the bounty, which is not okay as I just posted some cosmetic changes, but AFAIK there is no way to give it back. If there is no other solution for that, at least I could start a 50 points bounty for a question of your choice... Aug 3, 2011 at 10:43 • @Mihai posted solution here Aug 4, 2011 at 14:00 • @Landel I think your code lead me to write more concise Haskell code so you can think of these bounties as taken from me :) Aug 4, 2011 at 14:02 ## 2 Answers As per request by Mihai , here is the accepted code of spoj problem . The only error in previous code was " The minimal enclosing circle is determined by the diametric circle of S ". Instead of finding a circle passing through two points , i was calculating circle from three points. 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 = ( init lhull ) ++ ( init 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  Not an answer, just a small improvement (not tested): import Data.Function(on) findAngle :: (Num a , Ord a , Floating a ) => Point a -> Point a -> Point a -> ( Point a , Point a , Point a , a ) findAngle u v t | 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 * (sqrt $b * c)) 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 , Point a , Point a , a ) findPoints u v xs | 2 * theta >= pi = ( a , b , t , theta ) | and [ 2 * alpha <= pi , 2 * beta <= pi ] = (a , b , t , theta) | otherwise = if 2 * alpha > pi then findPoints v t xs else findPoints u t xs where t4 (_ , _ , _ , t) = t (a , b , t , theta) = minimumBy (compare on t4) . map (findAngle u v)$ xs
alpha = t4 $findAngle v t u --angle between v u t angle subtended at u by v t beta = t4$ findAngle u t v   -- angle between u v t angle subtended at v by  u t