My Weblog

Blog about programming and math

find the farthest pair of points

We have n points in the plane, and we have to find the distance between farthest pair of points.Cleary solution space points will belong to convex hull. So try to find out convex hull in O(nlogn) and using rotating calipers methods you can find distance between two farthest point in O(n).This problem is very well known diameter of convex polygon. See this problem on SPOJ
http://www.spoj.pl/problems/TFOSS/


#include<iostream>
#include<vector>
#include<cstdio>
typedef long long i64;
using namespace std;

	i64 cross(pair<i64,i64>P,pair<i64,i64>Q,pair<i64,i64> R)//if return >0 then point R is upper of lineseg PQ else down
		{
			return 	(Q.first-P.first)*(R.second-P.second)-(R.first-P.first)*(Q.second-P.second);
		}
	void hull(vector<pair<i64,i64> > &P,vector<pair<i64,i64> > &L,vector<pair<i64,i64> > &U)
		{
			int j=0,k=0,n=P.size();
			sort(P.begin(),P.end());
			U.resize(2*n);
			L.resize(2*n);
			for(int i=0;i<n;i++)
			 {
				while(j>=2 && cross(L[j-2],L[j-1],P[i])<=0)//p[i] is making right turn we need left turn
					j--;
				while(k>=2 && cross(U[k-2],U[k-1],P[i])>=0)//p[i] is making left turn we need right
					k--;
				L[j++]=P[i];
			 	U[k++]=P[i];
			}
			U.resize(k);
			L.resize(j);
		}
	i64 fun(pair<i64,i64> P ,pair<i64,i64> Q)
		{
			return (P.first-Q.first)*(P.first-Q.first)+(P.second-Q.second)*(P.second-Q.second);
		}

	int main()
		{
			i64 t,n,k1,k2;
			vector<pair<i64,i64> >v,U,L;
			cin>>t;
			while(t--)
			 {
				v.clear();
				U.clear();
				L.clear();
				cin>>n;
				for(int i=0;i<n;i++)
				 cin>>k1>>k2,v.push_back(make_pair(k1,k2));
				if(n==1)
				{
					cout<<"0\n";
					continue;
				}
				hull(v,L,U);
			 
				/*cout<<"lower hull\n";
				for(int i=0;i<L.size();i++)
					cout<<L[i].first<<" "<<L[i].second<<endl;
				cout<<"upper hull"<<endl;
				for(int i=0;i<U.size();i++)
					cout<<U[i].first<<" "<<U[i].second<<endl;*/
				int i=0,j,m;
				j=L.size()-1;
				m=U.size()-1;
				i64 dist=-1;
				while(i<m || j>0)
			 	{
					dist=max(dist,fun(U[i],L[j]));
					if(i==m)
						j--;
					else if(j==0)
						i++;
					else
				   	{
						if ( (U[i+1].second-U[i].second) * (L[j].first-L[j-1].first) > (L[j].second-L[j-1].second) * (U[i+1].first-U[i].first) )
						i++;
					else
						j--;
				   	}
				}
				cout<<dist<<endl;
			}
		}	



Tried some Haskell code but getting time limit exceed. If you are accepted in Haskell then please let me know.

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


About these ads

May 27, 2008 - Posted by | Programming

No comments yet.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 249 other followers

%d bloggers like this: