这段时间刷了一群计算几何题，呼！
做计算几何时一定要保证自己做到以下两点：
1.下手之前保证自己已经考虑到每一个细节，每一步的作法都要非常清晰。
2.敲代码时必须非常冷静，必须保持头脑清醒，保证自己敲下来的代码必然是对的（计算几何的代码几乎无法调试）。
接下来就是这一道道题了：
POJ2653
题目意思是一堆木条从天而降，要找到哪些木条是在最上面的（像一个古老的游戏）。
做法不是很难（计算几何做法都不是很难），直接上O(n2)就行了，考察叉积的应用。
1 #include <iostream> 2 #include <cstdio> 3 using namespace std; 4 const double eps=1e-10; 5 struct point 6 { 7 double x,y; 8 point(){x=y=0;} 9 point(double X,double Y){x=X,y=Y;} 10 }; 11 point V(point start,point end){return point(end.x-start.x,end.y-start.y);} 12 double crossP(point a,point b){return a.x*b.y-a.y*b.x;} 13 struct line 14 { 15 bool istop; 16 int prev,next; 17 point s,e; 18 }g[110000];int head; 19 bool have(line a,line b) 20 { 21 return (crossP(V(a.s,a.e),V(a.s,b.s))*crossP(V(a.s,a.e),V(a.s,b.e))<-eps) 22 &&(crossP(V(b.s,b.e),V(b.s,a.s))*crossP(V(b.s,b.e),V(b.s,a.e))<-eps); 23 } 24 int main(int argc, char *argv[]) 25 { 26 int n; 27 while(scanf("%d",&n)&&n) 28 { 29 head=0; 30 for(int i=1;i<=n;i++) 31 { 32 scanf("%lf%lf%lf%lf",&g[i].s.x,&g[i].s.y,&g[i].e.x,&g[i].e.y); 33 g[i].istop=true;g[i].next=0;g[i].prev=head;g[head].next=i; 34 for(int j=head;j;j=g[j].prev) 35 if(have(g[i],g[j])) 36 { 37 g[j].istop=false; 38 g[g[j].prev].next=g[j].next; 39 g[g[j].next].prev=g[j].prev; 40 } 41 head=i; 42 } 43 printf("Top sticks:"); 44 for(int i=1;i<n;i++) 45 if(g[i].istop)printf(" %d,",i); 46 printf(" %d. ",n); 47 } 48 return 0; 49 }
POJ2318
题意是将一个箱子用木板隔成好几格子，然后在箱子里撒一些点，求每个格子里点的个数。
直接枚举每个格子，用叉积判断点的是否在里面即可。
其实可以更快，比方说给点排个序，或者二分位置，可以把效率压下来。
1 #include <iostream> 2 #include <cstdio> 3 #include <cstring> 4 using namespace std; 5 struct point 6 { 7 int x,y; 8 point(){x=y=0;} 9 point(int X,int Y){x=X,y=Y;} 10 }up[5100],down[5100]; 11 inline point V(point s,point e){return point(e.x-s.x,e.y-s.y);} 12 inline int crossP(point a,point b){return a.x*b.y-a.y*b.x;} 13 int ans[5100]; 14 int main(int argc, char *argv[]) 15 { 16 int n,m;point ul,lr;//upper-left lower-right 17 while(scanf("%d",&n)&&n) 18 { 19 scanf("%d%d%d%d%d",&m,&ul.x,&ul.y,&lr.x,&lr.y); 20 memset(ans,0,sizeof ans); 21 up[0]=point(ul.x,ul.y);down[0]=point(ul.x,lr.y); 22 up[n+1]=point(lr.x,ul.y);down[n+1]=point(lr.x,lr.y); 23 for(int i=1;i<=n;i++) 24 { 25 scanf("%d%d",&up[i].x,&down[i].x); 26 up[i].y=ul.y;down[i].y=lr.y; 27 } 28 while(m--) 29 { 30 point t;scanf("%d%d",&t.x,&t.y); 31 for(int i=0;i<=n;i++) 32 if(crossP(V(down[i],up[i]),V(down[i],t))<0&&crossP(V(down[i+1],up[i+1]),V(down[i+1],t))>0) 33 { 34 ans[i]++; 35 break; 36 } 37 } 38 for(int i=0;i<=n;i++)printf("%d: %d ",i,ans[i]); 39 printf(" "); 40 } 41 return 0; 42 }
poj1113
好了凸包君出现了。
答案是凸包的周长+半径为l的周长的长度。
1 #include <iostream> 2 #include <cstdio> 3 #include <algorithm> 4 #include <cmath> 5 using namespace std; 6 const double pi=3.14159265358979323846; 7 struct point{ 8 int x,y; 9 point(){x=y=0;} 10 point(int X,int Y){x=X,y=Y;} 11 }p[1100]; 12 point V(point s,point e){return point(e.x-s.x,e.y-s.y);} 13 int crossP(point a,point b){return a.x*b.y-a.y*b.x;} 14 double dis(point a,point b){return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));} 15 bool cmp1(point a,point b){return a.x==b.x?a.y<b.y:a.x<b.x;} 16 bool cmp2(point a,point b){return a.x==b.x?a.y>b.y:a.x>b.x;} 17 point s[1100];int r; 18 int main(int argc, char *argv[]) 19 { 20 int n,l;scanf("%d%d",&n,&l); 21 for(int i=1;i<=n;i++)scanf("%d%d",&p[i].x,&p[i].y); 22 sort(p+1,p+n+1,cmp1); 23 r=0;s[++r]=p[1];s[++r]=p[2]; 24 double ans=2*pi*l; 25 for(int i=3;i<=n;i++) 26 { 27 while(r>1&&crossP(V(s[r-1],s[r]),V(s[r],p[i]))>0)r--; 28 s[++r]=p[i]; 29 } 30 for(int i=1;i<r;i++)ans+=dis(s[i],s[i+1]); 31 sort(p+1,p+n+1,cmp2); 32 r=0;s[++r]=p[1];s[++r]=p[2]; 33 for(int i=3;i<=n;i++) 34 { 35 while(r>1&&crossP(V(s[r-1],s[r]),V(s[r],p[i]))>0)r--; 36 s[++r]=p[i]; 37 } 38 for(int i=1;i<r;i++)ans+=dis(s[i],s[i+1]); 39 printf("%d ",(int)(ans+0.5)); 40 return 0; 41 }
POJ1556
题意是求(0,5)到(10,5)的最短路，只要枚举所有的特殊点，互相连边，SPFA即可。
P.S.这是第一次被奇葩POJ卡… G++输出double不能用%lf要用%f，不然会WA（C++提交不会）。
1 #include <iostream> 2 #include <cstdio> 3 #include <cmath> 4 #include <queue> 5 #include <cstring> 6 using namespace std; 7 const double eps=1e-10; 8 int n,pnum; 9 //point 10 struct point 11 { 12 double x,y; 13 int ord; 14 point(){x=y=0;} 15 point(double X,double Y,int Ord){x=X,y=Y,ord=Ord;} 16 }; 17 double dis(point a,point b) 18 { 19 return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y)); 20 } 21 //line 22 struct line 23 { 24 point s,e; 25 line(point S,point E){s=S,e=E;} 26 }; 27 inline double gety(line l,double x) 28 { 29 return (x-l.s.x)*(l.e.y-l.s.y)/(l.e.x-l.s.x)+l.s.y; 30 } 31 //wall 32 struct wall 33 { 34 double x,y[5]; 35 point p[5]; 36 }w[200]; 37 int head[600]; 38 //graph 39 struct edge{int to,next;double c;}g[80000];int gnum=0; 40 void addedge(int from,int to,double c) 41 { 42 g[++gnum].to=to;g[gnum].c=c;g[gnum].next=head[from];head[from]=gnum; 43 } 44 double d[500];bool isin[500]; 45 double SPFA() 46 { 47 queue<int> q; 48 for(int i=1;i<=pnum;i++)d[i]=1e10,isin[i]=false; 49 d[1]=0;isin[1]=true;q.push(1); 50 while(!q.empty()) 51 { 52 int now=q.front(); 53 for(int i=head[now];i;i=g[i].next) 54 if(d[now]+g[i].c+eps<d[g[i].to]) 55 { 56 d[g[i].to]=d[now]+g[i].c; 57 if(!isin[g[i].to]) 58 { 59 q.push(g[i].to); 60 isin[g[i].to]=true; 61 } 62 } 63 isin[now]=false;q.pop(); 64 } 65 return d[2]; 66 } 67 //main 68 bool check(line l) 69 { 70 for(int i=1;i<=n;i++) 71 { 72 if(w[i].x+eps<l.s.x||w[i].x>l.e.x+eps)continue; 73 double y=gety(l,w[i].x); 74 if(y+eps<w[i].y[1])return false; 75 else if(y-eps<w[i].y[2])continue; 76 else if(y+eps<w[i].y[3])return false; 77 else if(y-eps<w[i].y[4])continue ; 78 else return false; 79 } 80 return true; 81 } 82 void addedge(point s,point e) 83 { 84 if(check(line(s,e)))addedge(s.ord,e.ord,dis(s,e)); 85 } 86 int main(int argc, char *argv[]) 87 { 88 //freopen("1.in","r",stdin); 89 //freopen("1.out","w",stdout); 90 while(scanf("%d",&n)!=EOF&&n!=-1) 91 { 92 memset(head,0,sizeof head); 93 pnum=2;gnum=0; 94 point start(0,5,1),end(10,5,2); 95 for(int i=1;i<=n;i++) 96 { 97 scanf("%lf",&w[i].x); 98 for(int j=1;j<=4;j++) 99 { 100 scanf("%lf",&w[i].y[j]); 101 w[i].p[j]=point(w[i].x,w[i].y[j],++pnum); 102 } 103 } 104 addedge(start,end); 105 for(int i=1;i<=n;i++) 106 for(int j=1;j<=4;j++) 107 { 108 addedge(start,w[i].p[j]); 109 addedge(w[i].p[j],end); 110 } 111 for(int i=1;i<=n;i++) 112 for(int j=1;j<=4;j++) 113 for(int i1=i+1;i1<=n;i1++) 114 for(int j1=1;j1<=4;j1++) 115 addedge(w[i].p[j],w[i1].p[j1]); 116 printf("%.2f ",SPFA());//神奇的POJ告诉我如果这里用"%0.2lf "用G++就过不了 C++能过 117 } 118 return 0; 119 }
POJ1228
题意是判断给定的点所构成的凸包是否是稳定的。
稳定的定义是说，给出点集S,不存在一个点P使得
1.S的凸包≠S∪{P}的凸包
2.S中的点全部出现在S∪{P}的凸包上。
做这题需要一个结论：一个凸包是稳定的当且仅当这个凸包的每条边上都有数量大于等于3的点。
证明也很简单：如果有一条边上只有两个点，我们可以在凸包外侧离这条线极近的地方取一点P，P就能够使得1、2均成立。
1 #include <iostream> 2 #include <cstdio> 3 #include <algorithm> 4 #include <cmath> 5 using namespace std; 6 struct point{ 7 int x,y; 8 friend bool operator==(point a,point b){return a.x==b.x&&a.y==b.y;} 9 point(){x=y=0;} 10 point(int X,int Y){x=X,y=Y;} 11 }p[1100]; 12 point V(point s,point e){return point(e.x-s.x,e.y-s.y);} 13 int crossP(point a,point b){return a.x*b.y-a.y*b.x;} 14 bool cmp1(point a,point b){return a.x==b.x?a.y<b.y:a.x<b.x;} 15 bool cmp2(point a,point b){return a.x==b.x?a.y>b.y:a.x>b.x;} 16 bool inaline(point a,point b,point c) 17 { 18 return !crossP(V(a,b),V(b,c)); 19 } 20 point s[1100];int r; 21 int main(int argc, char *argv[]) 22 { 23 int T;scanf("%d",&T); 24 while(T--) 25 { 26 int n;scanf("%d",&n); 27 for(int i=1;i<=n;i++)scanf("%d%d",&p[i].x,&p[i].y); 28 if(n<=5){printf("NO ");continue;} 29 bool ans=true; 30 sort(p+1,p+n+1,cmp1); 31 r=0;s[++r]=p[1];s[++r]=p[2]; 32 for(int i=3;i<=n;i++) 33 { 34 while(r>1&&crossP(V(s[r-1],s[r]),V(s[r],p[i]))>0)r--; 35 s[++r]=p[i]; 36 } 37 if(r<=2||r>=n){printf("NO ");continue;} 38 for(int i=3;i<=r-2;i++) 39 if(!(inaline(s[i-1],s[i],s[i+1])||(inaline(s[i-2],s[i-1],s[i])&&inaline(s[i],s[i+1],s[i+2]))))ans=false; 40 if(!ans){printf("NO ");continue;} 41 point em1=s[r-1],em2=s[r-2],sa1=s[2],sa2=s[3]; 42 sort(p+1,p+n+1,cmp2); 43 r=0;s[++r]=p[1];s[++r]=p[2]; 44 for(int i=3;i<=n;i++) 45 { 46 while(r>1&&crossP(V(s[r-1],s[r]),V(s[r],p[i]))>0)r--; 47 s[++r]=p[i]; 48 } 49 if(r<=2||r>=n){printf("NO ");continue;} 50 for(int i=3;i<=r-2;i++) 51 if(!(inaline(s[i-1],s[i],s[i+1])||(inaline(s[i-2],s[i-1],s[i])&&inaline(s[i],s[i+1],s[i+2]))))ans=false; 52 if(!ans){printf("NO ");continue;} 53 if(!(inaline(s[r-1],s[r],sa1)||(inaline(s[r],sa1,sa2)&&inaline(s[r-2],s[r-1],s[r]))))ans=false; 54 if(!(inaline(em1,s[1],s[2])||(inaline(em2,em1,s[1])&&inaline(s[1],s[2],s[3]))))ans=false; 55 printf(ans?"YES ":"NO "); 56 } 57 return 0; 58 }
POJ1151
矩形面积并，经典线段树。
1 #include <iostream> 2 #include <cstdio> 3 #include <algorithm> 4 #include <cmath> 5 using namespace std; 6 const double eps=1e-10; 7 double sa[6100]; 8 struct node{int l,r,havecover;double len;}t[6100]; 9 void construct(int l,int r,int k) 10 { 11 t[k].l=l; 12 t[k].r=r; 13 t[k].havecover=t[k].len=0; 14 if(l<r) 15 { 16 int mid=(l+r)/2; 17 construct(l,mid,2*k); 18 construct(mid+1,r,2*k+1); 19 } 20 } 21 void add(int l,int r,int k,int num) 22 { 23 if(t[k].l==l&&t[k].r==r) 24 { 25 t[k].havecover+=num; 26 if(t[k].havecover)t[k].len=sa[r+1]-sa[l]; 27 else t[k].len=t[k].l<t[k].r?t[2*k].len+t[2*k+1].len:0; 28 return ; 29 } 30 int mid=(t[k].l+t[k].r)/2; 31 if(l>mid)add(l,r,2*k+1,num); 32 else if(r<=mid)add(l,r,2*k,num); 33 else { 34 add(l,mid,2*k,num); 35 add(mid+1,r,2*k+1,num); 36 } 37 if(!t[k].havecover)t[k].len=t[2*k].len+t[2*k+1].len; 38 } 39 struct line{ 40 double x; 41 int ys,ye,isleft; 42 friend bool operator<(line a,line b){return a.x<b.x;} 43 }l[6100];int lsum; 44 struct isort 45 { 46 int pos; 47 double num; 48 bool isstart; 49 friend bool operator<(isort a,isort b){return a.num<b.num;} 50 }s[6100];int snum; 51 int main(int argc, char *argv[]) 52 { 53 int n,time=0; 54 while(scanf("%d",&n)&&n) 55 { 56 lsum=snum=0; 57 for(int i=1;i<=n;i++) 58 { 59 double x1,y1,x2,y2;scanf("%lf%lf%lf%lf",&x1,&y1,&x2,&y2); 60 l[++lsum].x=x1;l[lsum].isleft=1; 61 s[++snum].num=y1;s[snum].pos=lsum;s[snum].isstart=true; 62 s[++snum].num=y2;s[snum].pos=lsum;s[snum].isstart=false; 63 l[++lsum].x=x2;l[lsum].isleft=-1; 64 s[++snum].num=y1;s[snum].pos=lsum;s[snum].isstart=true; 65 s[++snum].num=y2;s[snum].pos=lsum;s[snum].isstart=false; 66 } 67 sort(s+1,s+snum+1); 68 double last=s[1].num;int ord=1; 69 sa[ord]=s[1].num; 70 for(int i=1;i<=snum;i++) 71 { 72 if(s[i].num!=last)sa[++ord]=last=s[i].num; 73 if(s[i].isstart)l[s[i].pos].ys=ord; 74 else l[s[i].pos].ye=ord; 75 } 76 construct(1,ord-1,1); 77 sort(l+1,l+lsum+1); 78 double lastx=l[1].x,ans=0; 79 cerr<<sa[4]-sa[1]<<endl; 80 for(int i=1;i<=lsum;i++) 81 { 82 //cerr<<" == "<<t[1].len<<" "<<lastx<<" "<<l[i].x<<endl; 83 //cerr<<l[i].ys<<" "<<l[i].ye<<endl; 84 ans+=t[1].len*(l[i].x-lastx); 85 add(l[i].ys,l[i].ye-1,1,l[i].isleft); 86 lastx=l[i].x; 87 //for(int j=1;j<=10;j++)cerr<<t[j].l<<" "<<t[j].r<<" "<<t[j].havecover<<" "<<t[j].len<<endl; 88 } 89 printf("Test case #%d Total explored area: %.2f ",++time,ans);//POJ????%lf 90 } 91 return 0; 92 }
POJ2187
求最远点对，旋转卡壳。
这里有我以前一直理解错的一个地方：旋转卡壳的单调性应该是由对踵点所构成三角形的面积决定的而不是由其距离决定的。
1 #include <iostream> 2 #include <cstdio> 3 #include <algorithm> 4 #include <cmath> 5 using namespace std; 6 inline int imax(int a,int b){return a>b?a:b;} 7 struct point 8 { 9 int x,y; 10 point(){x=y=0;} 11 point(int x,int y){this->x=x;this->y=y;} 12 }; 13 point V(point s,point e){return point(e.x-s.x,e.y-s.y);} 14 int crossP(point a,point b){return a.x*b.y-a.y*b.x;} 15 int dis(point a,point b){return (a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y);} 16 point p[51000]; 17 bool cmp1(point a,point b){return a.x==b.x?a.y<b.y:a.x<b.x;} 18 bool cmp2(point a,point b){return a.x==b.x?a.y>b.y:a.x>b.x;} 19 20 point c[51000];int cnum; 21 inline int prev(int num){return num-1?num-1:cnum;} 22 inline int next(int num){return (num+1<=cnum)?num+1:1;} 23 point s[51000];int r; 24 int main(int argc, char *argv[]) 25 { 26 freopen("1.in","r",stdin); 27 freopen("1.out","w",stdout); 28 int n;scanf("%d",&n); 29 for(int i=1;i<=n;i++)scanf("%d%d",&p[i].x,&p[i].y); 30 31 sort(p+1,p+n+1,cmp1); 32 r=0;s[++r]=p[1];s[++r]=p[2]; 33 for(int i=3;i<=n;i++) 34 { 35 while(r>1&&crossP(V(s[r-1],s[r]),V(s[r],p[i]))>=0)r--; 36 s[++r]=p[i]; 37 } 38 for(int i=1;i<r;i++)c[++cnum]=s[i]; 39 sort(p+1,p+n+1,cmp2); 40 r=0;s[++r]=p[1];s[++r]=p[2]; 41 for(int i=3;i<=n;i++) 42 { 43 while(r>1&&crossP(V(s[r-1],s[r]),V(s[r],p[i]))>=0)r--; 44 s[++r]=p[i]; 45 } 46 for(int i=1;i<r;i++)c[++cnum]=s[i]; 47 int ans=0,d=2; 48 /* 49 for(int i=1;i<=cnum;i++) 50 for(d=i+1;d<=cnum;d++) 51 ans=imax(ans,dis(c[i],c[d])); 52 */ 53 for(int i=1;i<=cnum;i++) 54 { 55 while(abs(crossP(V(c[i],c[next(i)]),V(c[i],c[d])))<abs(crossP(V(c[i],c[next(i)]),V(c[i],c[next(d)])))) 56 d=next(d); 57 ans=imax(ans,dis(c[i],c[d]));ans=imax(ans,dis(c[next(i)],c[d])); 58 } 59 printf("%d ",ans); 60 return 0; 61 }
BZOJ1043: [HAOI2008]下落的圆盘
与圆有关。其实不论是并还是交，不论是周长还是面积，写法都是几乎完全一样的。
1 /************************************************************** 2 Problem: 1043 3 User: zhuohan123 4 Language: C++ 5 Result: Accepted 6 Time:280 ms 7 Memory:1392 kb 8 ****************************************************************/ 9 10 #include <iostream> 11 #include <cstdio> 12 #include <cmath> 13 #include <algorithm> 14 using namespace std; 15 const double pi=3.14159265358979323846,eps=1e-10; 16 inline double angout(double t){return t/pi*180;} 17 inline double theta(double y,double x) 18 { 19 double ang=atan2(y,x); 20 return ang>-eps?ang:2*pi+ang; 21 } 22 inline double tabs(double ang) 23 { 24 if(ang<-eps)return ang+2*pi; 25 if(ang>2*pi)return ang-2*pi; 26 return ang; 27 } 28 struct point 29 { 30 double x,y; 31 point(){} 32 point(double X,double Y){x=X,y=Y;} 33 }; 34 inline double dis(point a,point b){return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));} 35 inline double sqrdis(point a,point b){return (a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y);} 36 struct circle{point o;double r;}c[1100]; 37 int n; 38 struct T 39 { 40 double ang;int isstart; 41 friend bool operator<(const T& a,const T& b) 42 { 43 if(abs(a.ang-b.ang)<eps) 44 { 45 if(a.isstart==b.isstart)return false; 46 else if(a.isstart==0)return true; 47 else if(a.isstart==1)return true; 48 return false; 49 } 50 else return a.ang+eps<b.ang; 51 } 52 53 }a[7000];int anum; 54 int main(int argc, char *argv[]) 55 { 56 //freopen("1.in","r",stdin); 57 //freopen("1.out","w",stdout); 58 scanf("%d",&n); 59 for(int i=1;i<=n;i++)scanf("%lf%lf%lf",&c[i].r,&c[i].o.x,&c[i].o.y); 60 double ans=0; 61 for(int i=1;i<=n;i++) 62 { 63 bool isallcover=false; 64 for(int j=i+1;j<=n;j++) 65 if(dis(c[i].o,c[j].o)+c[i].r<c[j].r+eps){isallcover=true;break;} 66 if(isallcover)continue; 67 anum=0; 68 for(int j=i+1;j<=n;j++) 69 { 70 double d=dis(c[i].o,c[j].o); 71 if(d+eps>c[i].r+c[j].r)continue; 72 if(d+c[j].r<c[i].r+eps)continue; 73 double ang1=acos((c[i].r*c[i].r+sqrdis(c[i].o,c[j].o)-c[j].r*c[j].r)/(2*c[i].r*d)), 74 ang2=theta(c[j].o.y-c[i].o.y,c[j].o.x-c[i].o.x); 75 double eang=tabs(ang1+ang2),sang=tabs(ang2-ang1); 76 if(sang>eang+eps) 77 { 78 a[++anum].ang=0;a[anum].isstart=1; 79 a[++anum].ang=eang;a[anum].isstart=-1; 80 a[++anum].ang=sang;a[anum].isstart=1; 81 a[++anum].ang=2*pi;a[anum].isstart=-1; 82 } 83 else 84 { 85 a[++anum].ang=sang;a[anum].isstart=1; 86 a[++anum].ang=eang;a[anum].isstart=-1; 87 } 88 } 89 a[++anum].ang=0;a[anum].isstart=0; 90 sort(a+1,a+anum+1); 91 a[++anum].ang=2*pi;a[anum].isstart=0; 92 int k=0; 93 for(int j=1;j<anum;j++) 94 { 95 k+=a[j].isstart; 96 if(k==0)ans+=c[i].r*(a[j+1].ang-a[j].ang); 97 } 98 } 99 printf("%.3lf ",ans); 100 return 0; 101 }
POJ1981
这道题很有意思，求一个单位圆能够最多能够覆盖多少点。
做法很巧妙。分别以每个点为圆心，1为半径画圆，对于任意两个点，在这两个圆相交的区域中的任意点为圆心所画的单位圆必然同时包含这两个点。
所以，对于每一个圆，让它与其他所有的圆求交，看哪一段被其他圆覆盖的次数最多，更新答案。
1 #include <iostream> 2 #include <cstdio> 3 #include <algorithm> 4 #include <cmath> 5 using namespace std; 6 inline int imax(int a,int b){return a>b?a:b;} 7 const double pi=3.14159265358979323846,eps=1e-10; 8 struct point{ 9 //circle: R=1 10 double x,y; 11 point(){} 12 point(double X,double Y){x=X,y=Y;} 13 }c[310]; 14 inline double tabs(double ang){ 15 if(ang>2*pi+eps)return ang-2*pi; 16 if(ang<-eps)return ang+2*pi; 17 return ang; 18 } 19 inline double dis(point a,point b){return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));} 20 struct angle{ 21 double ang;int bra; 22 angle(){ang=0,bra=0;} 23 angle(double Ang,int Bra){ang=Ang,bra=Bra;} 24 friend bool operator<(angle a,angle b) 25 { 26 if(abs(a.ang-b.ang)<eps)return a.bra>b.bra; 27 return a.ang+eps<b.ang; 28 } 29 }a[3100];int anum; 30 int main(int argc, char *argv[]) 31 { 32 int n; 33 while(scanf("%d",&n)&&n) 34 { 35 for(int i=1;i<=n;i++)scanf("%lf%lf",&c[i].x,&c[i].y); 36 int ans=0; 37 for(int i=1;i<=n;i++) 38 { 39 anum=0; 40 a[++anum]=angle(0,0); 41 for(int j=1;j<=n;j++) 42 { 43 double d=dis(c[i],c[j]); 44 if(d<eps||d>2+eps)continue ; 45 double ang1=tabs(atan2(c[j].y-c[i].y,c[j].x-c[i].x)), 46 ang2=tabs(acos(d/2)); 47 double angs=tabs(ang1-ang2),ange=tabs(ang1+ang2); 48 if(angs<ange+eps) 49 { 50 a[++anum]=angle(angs,1); 51 a[++anum]=angle(ange,-1); 52 } 53 else 54 { 55 a[++anum]=angle(0,1); 56 a[++anum]=angle(ange,-1); 57 a[++anum]=angle(angs,1); 58 a[++anum]=angle(2*pi,-1); 59 } 60 } 61 sort(a+2,a+anum+1); 62 a[++anum]=angle(2*pi,0); 63 int k=0; 64 for(int i=1;i<anum;i++) 65 { 66 k+=a[i].bra; 67 ans=imax(k+1,ans); 68 } 69 } 70 printf("%d ",ans); 71 } 72 return 0; 73 }
BZOJ2178: 圆的面积并
做法如题…
1 /************************************************************** 2 Problem: 2178 3 User: zhuohan123 4 Language: C++ 5 Result: Accepted 6 Time:780 ms 7 Memory:1408 kb 8 ****************************************************************/ 9 10 #include <iostream> 11 #include <cstdio> 12 #include <cmath> 13 #include <algorithm> 14 //#define ONLINE_JUDGE 15 using namespace std; 16 const double pi=3.14159265358979323846,eps=1e-10; 17 double tabs(double ang) 18 { 19 if(ang>2*pi+eps)return ang-2*pi; 20 if(ang<-eps)return ang+2*pi; 21 return ang; 22 } 23 struct point 24 { 25 double x,y; 26 point(){x=y=0;} 27 point(double X,double Y){x=X,y=Y;} 28 }; 29 inline double crossP(point a,point b){return a.x*b.y-a.y*b.x;} 30 inline double dis(point a,point b){return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));} 31 inline double sqrdis(point a,point b){return (a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y);} 32 struct circle 33 { 34 point o;double r; 35 double area(double ang){return (ang-sin(ang))/2*r*r;} 36 point getpo(double ang){return point(cos(ang)*r+o.x,sin(ang)*r+o.y);} 37 }; 38 circle c[1100]; 39 bool cantuse[1100]; 40 struct angle 41 { 42 double ang;int bra; 43 angle(){ang=0,bra=0;} 44 angle(double Ang,int Bra){ang=Ang,bra=Bra;} 45 friend bool operator<(angle a,angle b) 46 { 47 if(abs(a.ang-b.ang)<eps)return a.bra>b.bra; 48 return a.ang<b.ang+eps; 49 } 50 }a[8100];int anum; 51 int main(int argc, char *argv[]) 52 { 53 #ifndef ONLINE_JUDGE 54 freopen("1.in","r",stdin); 55 freopen("1.out","w",stdout); 56 #endif 57 int n;scanf("%d",&n); 58 for(int i=1;i<=n;i++)scanf("%lf%lf%lf",&c[i].o.x,&c[i].o.y,&c[i].r); 59 for(int i=1;i<=n;i++) 60 for(int j=1;j<=n;j++) 61 if(i!=j&&!cantuse[j]&&(dis(c[i].o,c[j].o)+c[i].r<c[j].r+eps)) 62 { 63 cantuse[i]=true; 64 break; 65 } 66 double ans=0; 67 for(int i=1;i<=n;i++) 68 if(!cantuse[i]) 69 { 70 anum=0; 71 a[++anum]=angle(0,0); 72 for(int j=1;j<=n;j++) 73 if(!cantuse[j]&&i!=j) 74 { 75 double d=dis(c[i].o,c[j].o); 76 if(d+eps>c[i].r+c[j].r)continue ; 77 double ang1=tabs(atan2(c[j].o.y-c[i].o.y,c[j].o.x-c[i].o.x)), 78 ang2=acos((c[i].r*c[i].r+sqrdis(c[i].o,c[j].o)-c[j].r*c[j].r)/(2*c[i].r*d)); 79 double angs=tabs(ang1-ang2),ange=tabs(ang1+ang2); 80 if(angs>ange+eps) 81 { 82 a[++anum]=angle(0,1); 83 a[++anum]=angle(ange,-1); 84 a[++anum]=angle(angs,1); 85 a[++anum]=angle(2*pi,-1); 86 } 87 else 88 { 89 a[++anum]=angle(angs,1); 90 a[++anum]=angle(ange,-1); 91 } 92 } 93 sort(a+2,a+anum+1); 94 a[++anum]=angle(2*pi,0); 95 int k=0; 96 for(int j=1;j<anum;j++) 97 { 98 k+=a[j].bra; 99 if(!k) 100 { 101 ans+=c[i].area(a[j+1].ang-a[j].ang); 102 ans+=crossP(c[i].getpo(a[j].ang),c[i].getpo(a[j+1].ang))/2; 103 } 104 } 105 } 106 printf("%.3lf ",ans); 107 return 0; 108 }
POJ1279
题意是求多边形核的面积，直接上半平面交。
这里用的是ZZY的排序增量法，效率O(nlgn)，这个算法像极了水平可见直线的做法，只不过是用双端队列而不是栈来维护所有的半平面。
一开始我曾经尝试用ax+by+c≤0的表示法来弄半平面交，一下子就放弃了…计算几何还是必须上向量啊。
1 #include <iostream> 2 #include <cstdio> 3 #include <cmath> 4 #include <algorithm> 5 using namespace std; 6 const double eps=1e-10; 7 struct point 8 { 9 double x,y; 10 point(){x=y=0;} 11 point(double X,double Y){x=X,y=Y;} 12 friend bool operator==(point a,point b){return abs(a.x-b.x)<eps&&abs(a.y-b.y)<eps;} 13 friend bool operator!=(point a,point b){return !(a==b);} 14 }; 15 inline point V(point s,point e){return point(e.x-s.x,e.y-s.y);} 16 inline double crossP(point a,point b){return a.x*b.y-a.y*b.x;} 17 double area(point p[],int size)//逆时针 18 { 19 double ans=0; 20 for(int i=1;i<size;i++)ans+=crossP(p[i],p[i+1]); 21 ans+=crossP(p[size],p[1]); 22 return ans/2; 23 } 24 struct line 25 { 26 //halfplane:(在向量s->e的左侧) 27 point s,e; 28 double theta; 29 line(){} 30 line(point S,point E){s=S;e=E;theta=atan2(e.y-s.y,e.x-s.x);} 31 friend bool operator<(line a,line b) 32 { 33 if(abs(a.theta-b.theta)<eps)return crossP(V(b.s,a.e),V(b.s,b.e))<-eps; 34 else return a.theta<b.theta; 35 } 36 bool isout(point p){return crossP(V(s,e),V(s,p))<-eps;} 37 }; 38 inline point intersection(line a,line b) 39 { 40 double s1=crossP(V(a.s,b.e),V(a.s,b.s)),s2=crossP(V(a.e,b.s),V(a.e,b.e)); 41 return point((a.s.x*s2+a.e.x*s1)/(s2+s1),(a.s.y*s2+a.e.y*s1)/(s2+s1)); 42 } 43 line h[5000];int hnum; 44 void addhalfplane(point s,point e){h[++hnum]=line(s,e);} 45 point ans[5000];int pnum; 46 line q[5000];int ql=1,qr=0; 47 void halfplane_intersection() 48 { 49 int usenum=1;pnum=0; 50 sort(h+1,h+hnum+1); 51 for(int i=2;i<=hnum;i++) 52 if(abs(h[i].theta-h[i-1].theta)>eps)h[++usenum]=h[i]; 53 //for(int i=1;i<=usenum;i++)cerr<<"("<<h[i].s.x<<","<<h[i].s.y<<") -> ("<<h[i].e.x<<","<<h[i].e.y<<") "<<(h[i].theta/3.1415926535)*180<<endl; 54 ql=1;qr=0; 55 q[++qr]=h[1];q[++qr]=h[2]; 56 for(int i=3;i<=usenum;i++) 57 { 58 while(ql<qr&&h[i].isout(intersection(q[qr],q[qr-1])))qr--; 59 while(ql<qr&&h[i].isout(intersection(q[ql],q[ql+1])))ql++; 60 q[++qr]=h[i]; 61 } 62 while(ql<qr&&q[ql].isout(intersection(q[qr],q[qr-1])))qr--; 63 while(ql<qr&&q[qr].isout(intersection(q[ql],q[ql+1])))ql++; 64 if(qr-ql<=1)return ; 65 int tpnum=0; 66 for(int i=ql;i<qr;i++)ans[++tpnum]=intersection(q[i],q[i+1]); 67 ans[++tpnum]=intersection(q[ql],q[qr]); 68 for(int i=1;i<=tpnum;i++) 69 //printf("(%lf,%lf) ",ans[i].x,ans[i].y); 70 pnum=1; 71 for(int i=2;i<=tpnum;i++) 72 if(ans[i]!=ans[i-1])ans[++pnum]=ans[i]; 73 return ; 74 } 75 point p[5000]; 76 77 int main(int argc, char *argv[]) 78 { 79 //freopen("D.DAT","r",stdin); 80 //freopen("1.out","w",stdout); 81 int T;scanf("%d",&T); 82 while(T--) 83 { 84 hnum=0; 85 int n;scanf("%d",&n); 86 for(int i=1;i<=n;i++)scanf("%lf%lf",&p[i].x,&p[i].y); 87 if(area(p,n)<-eps) 88 { 89 for(int i=1;i<n;i++)addhalfplane(p[i+1],p[i]); 90 addhalfplane(p[1],p[n]); 91 } 92 else 93 { 94 for(int i=1;i<n;i++)addhalfplane(p[i],p[i+1]); 95 addhalfplane(p[n],p[1]); 96 } 97 halfplane_intersection(); 98 printf("%.2lf ",area(ans,pnum)); 99 } 100 return 0; 101 }
BZOJ1038: [ZJOI2008]瞭望塔
半平面交的经典题。山的每一个面可以被看到当且仅当瞭望塔那个面所对应的半平面之内。
1 /************************************************************** 2 Problem: 1038 3 User: zhuohan123 4 Language: C++ 5 Result: Accepted 6 Time:40 ms 7 Memory:4360 kb 8 ****************************************************************/ 9 10 #include <iostream> 11 #include <cstdio> 12 #include <cmath> 13 #include <algorithm> 14 using namespace std; 15 const double eps=1e-10,limit=1e13; 16 inline double imin(double a,double b){return a<b?a:b;} 17 struct point 18 { 19 double x,y; 20 point(){x=y=0;} 21 point(double X,double Y){x=X,y=Y;} 22 }; 23 point V(point s,point e){return point(e.x-s.x,e.y-s.y);} 24 double crossP(point a,point b){return a.x*b.y-a.y*b.x;} 25 struct line 26 { 27 //halfplane: s->e left 28 point s,e; 29 double ang; 30 line(){} 31 line(point S,point E){s=S;e=E;ang=atan2(e.y-s.y,e.x-s.x);} 32 friend bool operator<(line a,line b) 33 { 34 if(abs(a.ang-b.ang)<eps)return crossP(V(b.s,b.e),V(b.s,a.e))>eps; 35 else return a.ang<b.ang+eps; 36 } 37 bool isin(point po){return crossP(V(s,e),V(s,po))>-eps;} 38 double gety(double x){return (x-s.x)*((e.y-s.y)/(e.x-s.x))+s.y;} 39 }; 40 point intersection(line a,line b) 41 { 42 double s1=crossP(V(b.s,b.e),V(b.s,a.s)),s2=crossP(V(b.s,a.e),V(b.s,b.e)); 43 return point((s2*a.s.x+s1*a.e.x)/(s1+s2),(s2*a.s.y+s1*a.e.y)/(s1+s2)); 44 } 45 line h[5100];int hnum; 46 void addhalfplane(point s,point e){h[++hnum]=line(s,e);} 47 line q[51000];int ql,qr; 48 point hans[51000];int ansnum=0; 49 void halfplane_intersection() 50 { 51 ql=1;qr=0;ansnum=0; 52 sort(h+1,h+hnum+1); 53 int unum=1; 54 for(int i=2;i<=hnum;i++) 55 if(abs(h[i].ang-h[i-1].ang)>eps)h[++unum]=h[i]; 56 q[++qr]=h[1];q[++qr]=h[2]; 57 for(int i=3;i<=unum;i++) 58 { 59 while(ql<qr&&!h[i].isin(intersection(q[qr],q[qr-1])))qr--; 60 while(ql<qr&&!h[i].isin(intersection(q[ql],q[ql+1])))ql++; 61 q[++qr]=h[i]; 62 } 63 while(ql<qr&&!h[ql].isin(intersection(q[qr],q[qr-1])))qr--; 64 while(ql<qr&&!h[qr].isin(intersection(q[ql],q[ql+1])))ql++; 65 if(qr-ql<2)return ; 66 int tansnum=0; 67 hans[++tansnum]=intersection(q[ql],q[qr]); 68 for(int i=ql;i<qr;i++)hans[++tansnum]=intersection(q[i],q[i+1]); 69 ansnum=1; 70 for(int i=2;i<=tansnum;i++) 71 if(abs(hans[i].x-hans[i-1].x)>eps||abs(hans[i].y-hans[i-1].y)>eps) 72 hans[++ansnum]=hans[i]; 73 return ; 74 } 75 point p[5100]; 76 int main(int argc, char *argv[]) 77 { 78 //freopen("1.in","r",stdin); 79 //freopen("1.out","w",stdout); 80 int n;scanf("%d",&n); 81 for(int i=1;i<=n;i++)scanf("%lf",&p[i].x); 82 for(int i=1;i<=n;i++)scanf("%lf",&p[i].y); 83 addhalfplane(point(-limit,limit),point(-limit,-limit)); 84 addhalfplane(point(-limit,-limit),point(limit,-limit)); 85 addhalfplane(point(limit,-limit),point(limit,limit)); 86 addhalfplane(point(limit,limit),point(-limit,limit)); 87 for(int i=1;i<n;i++)addhalfplane(p[i],p[i+1]); 88 halfplane_intersection(); 89 double ans=1e10; 90 for(int i=1,j=1;i<=ansnum;i++) 91 { 92 if(hans[i].x+eps<p[1].x||hans[i].x+eps>p[n].x+eps)continue ; 93 while(j<n&&p[j+1].x<hans[i].x+eps)j++; 94 ans=imin(ans,hans[i].y-line(p[j],p[j+1]).gety(hans[i].x)); 95 } 96 for(int i=1,j=1;i<=n;i++) 97 { 98 while(j<ansnum&&hans[j+1].x<p[i].x+eps)j++; 99 ans=imin(ans,line(hans[j],hans[j+1]).gety(p[i].x)-p[i].y); 100 } 101 printf("%.3lf ",ans); 102 return 0; 103 }
poj1070
能做到这道题完全是意外…
题意很简单，做法很直接，就是很难搞！
考点只有向量的旋转。本来每次旋转的角度应该二分的，像我一样用三角函数是很不靠谱的做法……
这题有几个特殊的点要考虑的：
1.可能一开始是朝右转的。
2.不能输出-0.000。
1 #include <iostream> 2 #include <cstdio> 3 #include <cmath> 4 #include <algorithm> 5 using namespace std; 6 const double pi=3.14159265358979,eps=1e-9; 7 int n; 8 struct point 9 { 10 double x,y; 11 point(){x=y=0;} 12 point(double X,double Y){x=X,y=Y;} 13 point rotate(point po,double theta) 14 { 15 double tx=po.x-x,ty=po.y-y; 16 double Sin=sin(theta),Cos=cos(theta); 17 double ax=tx*Cos-ty*Sin,ay=ty*Cos+tx*Sin; 18 return point(x+ax,y+ay); 19 } 20 friend point operator+(point a,point b){return point(a.x+b.x,a.y+b.y);} 21 inline point normal(){return point(-y,x);} 22 inline double len(){return sqrt(x*x+y*y);} 23 }p[20],zx,hig; 24 inline point V(point s,point e){return point(e.x-s.x,e.y-s.y);} 25 inline double dotP(point a,point b){return a.x*b.x+a.y*b.y;} 26 inline double crossP(point a,point b){return a.x*b.y-a.y*b.x;} 27 inline double dis(point a,point b){return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));} 28 struct line 29 { 30 point s,e; 31 line(){} 32 line(point S,point E){s=S,e=E;} 33 inline bool ison(point a) 34 { 35 return a.x+eps>s.x&&a.x<e.x+eps&&fabs(crossP(V(s,a),V(a,e)))<eps; 36 } 37 }l[110]; 38 int lnum; 39 double len[110],k[110]; 40 41 point getleftintersection(point o,double r,line l,double lpdis) 42 { 43 double ang1=atan2(l.e.y-l.s.y,l.e.x-l.s.x); 44 double ang2=asin(lpdis/r); 45 double ang=ang1+ang2+pi; 46 return point(o.x+r*cos(ang),o.y+r*sin(ang)); 47 } 48 int goleft(int o) 49 { 50 double ang=2*pi; 51 for(int i=1;i<=n;i++) 52 if(i!=o) 53 { 54 double r=dis(p[i],p[o]); 55 for(int j=lnum;j;j--) 56 { 57 double lpdis=dotP(V(l[j].s,p[o]),V(l[j].s,l[j].e).normal())/V(l[j].s,l[j].e).len(); 58 if(lpdis>r+eps)continue; 59 point inte=getleftintersection(p[o],r,l[j],lpdis); 60 if(l[j].ison(inte)) 61 { 62 double tang=2*asin(dis(p[i],inte)/(2*r)); 63 if(tang+eps<ang&&tang>eps)ang=tang; 64 } 65 } 66 } 67 for(int i=1;i<=n;i++) 68 if(i!=o) p[i]=p[o].rotate(p[i],ang); 69 zx=p[o].rotate(zx,ang); 70 for(int i=n;i>=1;i--) 71 if(i!=o) 72 for(int j=1;j<=lnum;j++) 73 if(l[j].ison(p[i])) 74 return i; 75 } 76 77 point getrightintersection(point o,double r,line l,double lpdis) 78 { 79 double ang1=atan2(l.e.y-l.s.y,l.e.x-l.s.x); 80 double ang2=asin(lpdis/r); 81 double ang=ang1-ang2; 82 return point(o.x+r*cos(ang),o.y+r*sin(ang)); 83 } 84 int goright(int o) 85 { 86 //cerr<<o<<endl; 87 double ang=2*pi; 88 for(int i=1;i<=n;i++) 89 if(i!=o) 90 { 91 double r=dis(p[i],p[o]); 92 for(int j=lnum;j;j--) 93 { 94 double lpdis=dotP(V(l[j].s,p[o]),V(l[j].s,l[j].e).normal())/V(l[j].s,l[j].e).len(); 95 if(lpdis>r+eps)continue; 96 point inte=getrightintersection(p[o],r,l[j],lpdis); 97 if(l[j].ison(inte)) 98 { 99 double tang=2*asin(dis(p[i],inte)/(2*r)); 100 if(tang+eps<ang&&tang>eps)ang=tang; 101 break ; 102 } 103 } 104 } 105 for(int i=1;i<=n;i++) 106 if(i!=o) p[i]=p[o].rotate(p[i],2*pi-ang); 107 zx=p[o].rotate(zx,2*pi-ang); 108 for(int i=n;i>=1;i--) 109 if(i!=o) 110 for(int j=1;j<=lnum;j++) 111 if(l[j].ison(p[i])) 112 return i; 113 } 114 115 int main(int argc, char *argv[]) 116 { 117 //freopen("wheel.in","r",stdin); 118 //freopen("wheel.out","w",stdout); 119 int T;scanf("%d",&T); 120 while(T--) 121 { 122 scanf("%d",&n); 123 for(int i=1;i<=n;i++)scanf("%lf%lf",&p[i].x,&p[i].y); 124 scanf("%lf%lf",&zx.x,&zx.y); 125 lnum=0; 126 do 127 { 128 lnum++; 129 scanf("%lf%lf",&len[lnum],&k[lnum]); 130 } 131 while(abs(k[lnum])>eps); 132 scanf("%lf%lf",&hig.x,&hig.y); 133 l[1].e=hig; 134 for(int i=1;i<=lnum;i++) 135 { 136 double tx=len[i]/sqrt(1+k[i]*k[i]); 137 l[i].s=l[i+1].e=point(l[i].e.x-tx,l[i].e.y-k[i]*tx); 138 } 139 //for(int i=1;i<=lnum;i++)printf("(%.4lf,%.4lf)->(%.4lf,%.4lf) ",l[i].s.x,l[i].s.y,l[i].e.x,l[i].e.y); 140 int start[3]={0},ssum=0; 141 for(int i=1;i<=n;i++) 142 for(int j=1;j<=lnum;j++) 143 if(l[j].ison(p[i])) 144 start[++ssum]=i; 145 if(ssum==2&&p[start[2]].x<p[start[1]].x)swap(start[1],start[2]); 146 if(ssum==1)start[2]=start[1]; 147 if(p[start[1]].x+eps>zx.x) 148 { 149 while(p[start[1]].x+eps>zx.x) 150 start[1]=goleft(start[1]); 151 } 152 else if(p[start[2]].x+eps<zx.x) 153 { 154 while(p[start[2]].x+eps<zx.x) 155 start[2]=goright(start[2]); 156 } 157 printf("%.3lf %.3lf ",zx.x,zx.y+eps); 158 } 159 return 0; 160 }