旋转卡壳的经典应用，实现时直接比较角度。
1 /************************************************************** 2 Problem: 1185 3 User: idy002 4 Language: C++ 5 Result: Accepted 6 Time:172 ms 7 Memory:2580 kb 8 ****************************************************************/ 9 10 #include <cstdio> 11 #include <cmath> 12 #include <algorithm> 13 #define oo 1e20 14 #define N 50010 15 #define eps 1e-10 16 using namespace std; 17 18 int sg( double x ) { return (x>-eps)-(x<eps); } 19 struct Vector { 20 double x, y; 21 Vector(){} 22 Vector( double x, double y ):x(x),y(y){} 23 Vector operator+( const Vector & b ) const { return Vector(x+b.x,y+b.y); } 24 Vector operator-( const Vector & b ) const { return Vector(x-b.x,y-b.y); } 25 Vector operator*( double b ) const { return Vector(x*b,y*b); } 26 Vector operator/( double b ) const { return Vector(x/b,y/b); } 27 double operator^( const Vector & b ) const { return x*b.y-y*b.x; } 28 double operator&( const Vector & b ) const { return x*b.x+y*b.y; } 29 void rotate( double ang ) { 30 double cosa=x, sina=y; 31 double cosx=cos(ang), sinx=sin(ang); 32 x = cosa*cosx-sina*sinx; 33 y = sina*cosx+cosa*sinx; 34 } 35 double ang() const { return atan2(y,x); } 36 double len() const { return sqrt(x*x+y*y); } 37 bool operator<( const Vector & b ) const { 38 return x<b.x || (x-b.x<eps && y<b.y); 39 } 40 }; 41 typedef Vector Point; 42 43 int n, cn; 44 int nxt[N]; 45 Point pts[N], cvx[N], rct[4]; 46 47 bool onleft( Point &A, Point &B, Point &C ) { 48 return sg((B-A)^(C-A)) > 0; 49 } 50 int convex( int n, Point *p, Point *c ) { 51 if( n==1 ) { 52 c[0] = p[0]; 53 return 1; 54 } 55 int m; 56 sort( p, p+n ); 57 c[m=0] = p[0]; 58 for( int i=1; i<n; i++ ) { 59 while( m>0 && !onleft(c[m-1],c[m],p[i]) ) m--; 60 c[++m] = p[i]; 61 } 62 int k=m; 63 for( int i=n-2; i>=0; i-- ) { 64 while( m>k && !onleft(c[m-1],c[m],p[i]) ) m--; 65 c[++m] = p[i]; 66 } 67 return m; 68 } 69 inline double min( double a, double b ) { 70 return a<b ? a : b; 71 } 72 inline double max( double a, double b ) { 73 return a>b ? a : b; 74 } 75 inline double min( double a, double b, double c, double d ) { 76 double ab = min(a,b); 77 double cd = min(c,d); 78 return min(ab,cd); 79 } 80 inline double pdis( const Point &A, const Point &B, Point P ) { 81 return fabs( (B-A)^(P-A) ) / (A-B).len(); 82 } 83 inline Point inter( const Point &P, const Vector &u, const Point &Q, const Vector &v ) { 84 return P+u*(((Q-P)^v)/(u^v)); 85 } 86 inline void round_ang( double &ang ) { 87 while( ang>=M_PI ) ang-=M_PI; 88 while( ang<0 ) ang+=M_PI; 89 } 90 double minrect( int n, Point *p, Point *r ) { 91 if( n==1 ) { 92 r[0] = r[1] = r[2] = r[3] = p[0]; 93 return 0.0; 94 } else if( n==2 ) { 95 r[0] = r[1] = p[0]; 96 r[2] = r[3] = p[1]; 97 return 0.0; 98 } 99 for( int i=0; i<n; i++ ) nxt[i]=i+1; 100 nxt[n-1]=0; 101 102 double rt, suma; 103 int pa=0, pb=0, pc=0, pd=0; 104 Vector va(1,0), vb(0,1), vc(-1,0), vd(0,-1); 105 for( int i=1; i<n; i++ ) { 106 if( p[i].y<p[pa].y ) pa=i; 107 if( p[i].x>p[pb].x ) pb=i; 108 if( p[i].y>p[pc].y ) pc=i; 109 if( p[i].x<p[pd].x ) pd=i; 110 } 111 rt = (p[pc].y-p[pa].y)*(p[pb].x-p[pd].x); 112 r[0] = inter(p[pa],va,p[pb],vb); 113 r[1] = inter(p[pb],vb,p[pc],vc); 114 r[2] = inter(p[pc],vc,p[pd],vd); 115 r[3] = inter(p[pd],vd,p[pa],va); 116 117 /* 118 fprintf( stderr, "init bound: " ); 119 fprintf( stderr, "(%lf,%lf) (%lf,%lf) " 120 "(%lf,%lf) (%lf,%lf) ", 121 p[pa].x, p[pa].y, p[pb].x, p[pb].y, 122 p[pc].x, p[pc].y, p[pd].x, p[pd].y ); 123 */ 124 suma = 0.0; 125 while( suma <= M_PI_2*1.1 ) { 126 double da = (p[nxt[pa]]-p[pa]).ang()-va.ang(); 127 double db = (p[nxt[pb]]-p[pb]).ang()-vb.ang(); 128 double dc = (p[nxt[pc]]-p[pc]).ang()-vc.ang(); 129 double dd = (p[nxt[pd]]-p[pd]).ang()-vd.ang(); 130 round_ang(da); 131 round_ang(db); 132 round_ang(dc); 133 round_ang(dd); 134 double dm = min(da,db,dc,dd); 135 /* 136 va.rotate( dm ); 137 vb.rotate( dm ); 138 vc.rotate( dm ); 139 vd.rotate( dm ); 140 */ 141 if( dm==da ) { 142 va=p[nxt[pa]]-p[pa]; 143 vb=vc=vd=va; 144 vb.rotate(M_PI_2); 145 vc.rotate(M_PI); 146 vd.rotate(M_PI_2+M_PI); 147 pa=nxt[pa]; 148 } else if( dm==db ) { 149 vb=p[nxt[pb]]-p[pb]; 150 vc=vd=va=vb; 151 vc.rotate(M_PI_2); 152 vd.rotate(M_PI); 153 va.rotate(M_PI_2+M_PI); 154 155 pb=nxt[pb]; 156 } else if( dm==dc ) { 157 vc=p[nxt[pc]]-p[pc]; 158 vd=va=vb=vc; 159 vd.rotate(M_PI_2); 160 va.rotate(M_PI); 161 vb.rotate(M_PI_2+M_PI); 162 pc=nxt[pc]; 163 } else { 164 vd=p[nxt[pd]]-p[pd]; 165 va=vb=vc=vd; 166 va.rotate(M_PI_2); 167 vb.rotate(M_PI); 168 vc.rotate(M_PI_2+M_PI); 169 pd=nxt[pd]; 170 } 171 double area = pdis(p[pa],p[pa]+va,p[pc]) * pdis(p[pb],p[pb]+vb,p[pd]); 172 if( area<rt ) { 173 rt=area; 174 r[0] = inter(p[pa],va,p[pb],vb); 175 r[1] = inter(p[pb],vb,p[pc],vc); 176 r[2] = inter(p[pc],vc,p[pd],vd); 177 r[3] = inter(p[pd],vd,p[pa],va); 178 } 179 suma += dm; 180 } 181 return rt; 182 } 183 184 int main() { 185 scanf( "%d", &n ); 186 for( int i=0; i<n; i++ ) 187 scanf( "%lf%lf", &pts[i].x, &pts[i].y ); 188 cn = convex( n, pts, cvx ); 189 190 /* 191 fprintf( stderr, "There are %d points on convex ", cn ); 192 for( int i=0; i<cn; i++ ) 193 fprintf( stderr, "(%.5lf,%.5lf) ", cvx[i].x, cvx[i].y ); 194 */ 195 196 double marea = minrect( cn, cvx, rct ); 197 int c=0; 198 for( int i=1; i<4; i++ ) 199 if( rct[i].y<rct[c].y || (rct[i].y-rct[c].y<eps && rct[i].x<rct[c].x) ) 200 c=i; 201 rotate( rct, rct+c, rct+4 ); 202 printf( "%.5lf ", marea ); 203 for( int i=0; i<4; i++ ) 204 printf( "%.5lf %.5lf ", rct[i].x, rct[i].y ); 205 }
这是求最远点对的旋转卡壳：
1 #include <cstdio> 2 #include <cmath> 3 #include <iostream> 4 #include <algorithm> 5 #define N 500010 6 #define eps 1e-10 7 using namespace std; 8 9 10 inline int sg( double x ) { return (x>-eps)-(x<eps); } 11 12 struct Vector { 13 double x, y; 14 Vector(){} 15 Vector( double x, double y ):x(x),y(y){} 16 Vector operator+( Vector &b ) { return Vector(x+b.x,y+b.y); } 17 Vector operator-( Vector &b ) { return Vector(x-b.x,y-b.y); } 18 Vector operator*( double b ) { return Vector(x*b,y*b); } 19 Vector operator/( double b ) { return Vector(x/b,y/b); } 20 double operator^( const Vector &b ) const { return x*b.y-y*b.x; } 21 double operator&( const Vector &b ) const { return x*b.x+y*b.y; } 22 double len() { return sqrt(x*x+y*y); } 23 bool operator<( const Vector & b ) const { 24 return x<b.x || (x==b.x&&y<b.y); 25 } 26 }; 27 typedef Vector Point; 28 29 int n, cn; 30 Point pts[N], cvx[N]; 31 int nxt[N]; 32 33 inline bool onleft( Point &A, Point &B, Point &P ) { 34 return sg((B-A)^(P-A)) > 0; 35 } 36 int convex( int n, Point *p, Point *c ) { 37 int m; 38 sort( p, p+n ); 39 c[m=0] = p[0]; 40 for( int i=1; i<n; i++ ) { 41 while( m>0 && !onleft(c[m-1],c[m],p[i]) ) m--; 42 c[++m] = p[i]; 43 } 44 int k=m; 45 for( int i=n-2; i>=0; i-- ) { 46 while( m>k && !onleft(c[m-1],c[m],p[i]) ) m--; 47 c[++m] = p[i]; 48 } 49 return m+(n==1); 50 } 51 double area( Point &A, Point &B, Point &C ) { 52 return (B-A)^(C-A); 53 } 54 double maxdis( int n, Point *p ) { 55 for( int i=0; i<n; i++ ) nxt[i]=i+1; 56 nxt[n-1] = 0; 57 58 if( n==1 ) return 0.0; 59 if( n==2 ) return (p[0]-p[1]).len(); 60 61 double rt=-1e20; 62 int a, b, c=2, d; 63 for( int i=0; i<n; i++ ) { 64 a = i; 65 b = nxt[i]; 66 while( area(p[a],p[b],p[c])<area(p[a],p[b],p[nxt[c]]) ) c=nxt[c]; 67 d = nxt[c]; 68 rt = max( rt, (p[a]-p[c]).len() ); 69 rt = max( rt, (p[a]-p[d]).len() ); 70 rt = max( rt, (p[b]-p[c]).len() ); 71 rt = max( rt, (p[b]-p[d]).len() ); 72 } 73 return rt; 74 } 75 int main() { 76 scanf( "%d", &n ); 77 for( int i=0; i<n; i++ ) 78 scanf( "%lf%lf", &pts[i].x, &pts[i].y ); 79 cn = convex( n, pts, cvx ); 80 printf( "%.5lf ", maxdis(cn,cvx) ); 81 }