bzoj 1185

旋转卡壳的经典应用,实现时直接比较角度。

 

  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:\n" );
119     fprintf( stderr, "(%lf,%lf) (%lf,%lf)\n"
120                      "(%lf,%lf) (%lf,%lf)\n",
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\n", cn );
192     for( int i=0; i<cn; i++ )
193         fprintf( stderr, "(%.5lf,%.5lf)\n", 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\n", marea );
203     for( int i=0; i<4; i++ )
204         printf( "%.5lf %.5lf\n", rct[i].x, rct[i].y );
205 }
View Code

 

 

这是求最远点对的旋转卡壳:

 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\n", maxdis(cn,cvx) );
81 }
View Code

 

posted @ 2015-04-11 15:49  idy002  阅读(165)  评论(0编辑  收藏  举报