计算几何随笔
Lifting the Stone
对于多边形的重心，将多边形三角剖分，求出每一个三角形的重心，和三角形的面积加权以后求平均值，三角形的重心即为三个点的平均值。
循环除6容易爆精度，应该把除数放在循环外，
而且判断double是否为0，应该用==0，<eps是不行的，有负数情况。
复杂度O（n）
#include<bits/stdc++.h> using namespace std; const int N=1e6+5000; typedef double db; db eps=1e-6; int n; struct Point{ db x,y; Point(double X=0,double Y=0){x=X;y=Y;} Point operator +(Point B){return Point(x+B.x,y+B.y);} Point operator -(Point B){return Point(x-B.x,y-B.y);} Point operator *(double k){return Point(x*k,y*k);} Point operator /(double k){return Point(x/k,y/k);} }; typedef Point Vector; db Cross(Vector A,Vector B){return A.x*B.y-A.y*B.x;} // Point P[N]; vector<Point>P(N); db Polygon_area(){ double ans=0; for(int i=0;i<n;i++){ ans+=Cross(P[i],P[(i+1)%n]); } return ans/2; } Point Polygon_center(){ Point ans(0,0); if(Polygon_area()==0)return ans; for(int i=0;i<n;i++){ ans=ans+(P[i]+P[(i+1)%n])*Cross(P[i],P[(i+1)%n]); } return ans/Polygon_area()/6; } int main(){ int t;cin>>t; while(t--){ scanf("%d",&n); for(int i=0;i<n;i++)scanf("%lf %lf",&P[i].x,&P[i].y); Point ans=Polygon_center(); printf("%.2f %.2f\n",ans.x,ans.y); } // system("pause"); return 0; }
Surround the Trees
求凸包周长
将多边形的点按x坐标排序，扫描上凸包和下凸包，把偏向右的点剔除，求凸包一个周长。
复杂度O（nlogn）
#include<bits/stdc++.h> using namespace std; const int N=1e3+5000; typedef double db; db eps=1e-6; int sgn(db x){ if(fabs(x)<eps)return 0; else return x<0?-1:1; } struct Point{ db x,y; Point(double X=0,double Y=0){x=X;y=Y;} Point operator +(Point B){return Point(x+B.x,y+B.y);} Point operator -(Point B){return Point(x-B.x,y-B.y);} Point operator *(double k){return Point(x*k,y*k);} Point operator /(double k){return Point(x/k,y/k);} // Point operator <(Point B){ // return sgn(x-B.x)<0||(sgn(x-B.x)==0&&sgn(y-B.y)<0); // } // Point operator ==(Point B){ // return sgn(x-B.x)==0&&sgn(y-B.y)==0; // } friend bool operator<(Point A,Point B){ return sgn(A.x-B.x)<0||(sgn(A.x-B.x)==0&&sgn(A.y-B.y)<0); } friend bool operator==(Point A,Point B){ return sgn(A.x-B.x)==0&&sgn(A.y-B.y)==0; } }; int n; bool cmp(Point A,Point B){ return sgn(A.x-B.x)<0||(sgn(A.x-B.x)==0&&sgn(A.y-B.y)<0); } Point P[N],ch[N]; typedef Point Vector; db Cross(Vector A,Vector B){return A.x*B.y-A.y*B.x;} db Distance(Point A,Point B){return hypot(A.x-B.x,A.y-B.y);} int convexHull(){ sort(P,P+n); n=unique(P,P+n)-P; int cnt=0; for(int i=0;i<n;i++){ while(cnt>1&&sgn(Cross(ch[cnt-1]-ch[cnt-2],P[i]-ch[cnt-2]))<=0)cnt--; ch[cnt++]=P[i]; } int j=cnt; for(int i=n-2;i>=0;i--){ while(cnt>j&&sgn(Cross(ch[cnt-1]-ch[cnt-2],P[i]-ch[cnt-2]))<=0)cnt--; ch[cnt++]=P[i]; } if(n>1)cnt--; return cnt; } int main(){ while(~scanf("%d",&n),n){ for(int i=0;i<n;i++)scanf("%lf %lf",&P[i].x,&P[i].y); int cnt=convexHull(); double ans=0; if(cnt==1)ans=0; else if(cnt==2)ans=Distance(ch[0],ch[1]); else { for(int i=0;i<cnt;i++)ans+=Distance(ch[i],ch[(i+1)%cnt]); } printf("%.2lf\n",ans); } // system("pause"); return 0; }
Run
一些人跑步，给定若干个人的s=s0+vt，求有多少个人在某一个时刻在第一。
建立坐标系，添加一个沿y轴向下，沿x轴向左的向量，对若干个向量求一个半平面交即可。
复杂度O（log（n））
#include<bits/stdc++.h> using namespace std; const int N=1e4+5000; typedef double db; const db INF=1e16; const db eps=1e-6; int sgn(db x){ if(fabs(x)<eps)return 0; else return x<0?-1:1; } struct Point{ db x,y; Point(double X=0,double Y=0){x=X;y=Y;} Point operator +(Point B){return Point(x+B.x,y+B.y);} Point operator -(Point B){return Point(x-B.x,y-B.y);} Point operator *(double k){return Point(x*k,y*k);} Point operator /(double k){return Point(x/k,y/k);} }; typedef Point Vector; struct Line{ Point p; Vector v; db ang; Line(){} Line(Point p,Vector v):p(p),v(v){ang=atan2(v.y,v.x);} bool operator <(Line &L){return ang<L.ang;} }; db Cross(Vector A,Vector B){return A.x*B.y-A.y*B.x;} db Distance(Point A,Point B){return hypot(A.x-B.x,A.y-B.y);} //判断点在向量左边，即半平面里面 bool OnLeft(Line L,Point p){return sgn(Cross(L.v,p-L.p))>0;} Point Cross_point(Line a,Line b){ Vector u=a.p-b.p; db t=Cross(b.v,u)/Cross(a.v,b.v); return a.p+a.v*t; } vector<Point> HPI(vector<Line>L){ int n=L.size(); sort(L.begin(),L.end());//极角排序 int first,last; vector<Point>p(n); vector<Line>q(n); q[first=last=0]=L[0]; for(int i=1;i<n;i++){ // cout<<"que :"<<endl; //覆盖队尾 while(first<last&&!OnLeft(L[i],p[last-1]))last--; //覆盖队头 while(first<last&&!OnLeft(L[i],p[first]))first++; q[++last]=L[i]; // 极角相同保留左边 if(fabs(Cross(q[last].v,q[last-1].v))<eps){ last--; if(OnLeft(q[last],L[i].p))q[last]=L[i]; } if(first<last)p[last-1]=Cross_point(q[last-1],q[last]); } // cout<<"out :"<<endl; // 剔除无用的队尾元素 while(first<last&&!OnLeft(q[first],p[last-1]))last--; vector<Point>ans; if(last-first<=1)return ans; p[last]=Cross_point(q[last],q[first]); for(int i=first;i<=last;i++)ans.push_back(p[i]); return ans; } int main(){ int t;cin>>t; while(t--){ int n;scanf("%d",&n); vector<Line>lines; lines.push_back(Line(Point(0,0),Vector(0,-1))); lines.push_back(Line(Point(0,INF),Vector(-1,0))); // cout<<"test :"<<endl; for(int i=1;i<=n;i++){ db s,v;scanf("%lf %lf",&s,&v); lines.push_back(Line(Point(0,s),Vector(1,v))); } vector<Point>ans=HPI(lines); // cout<<"test :"<<endl; // for(int i=0;i<ans.size();i++){ // printf("test :%.2lf %.2lf \n",ans[i].x,ans[i].y); // } printf("%d\n",ans.size()-2); } // system("pause"); return 0; }
Problem G. Interstellar Travel
给出若干个点，从1到n，每次走的费用为叉积，求字典序最小的，费用最小的方案。
沿着上凸包走必然最小，把点排序出来，求一个上凸包即可。
【排序开始写的有点问题一直wa】
#include<bits/stdc++.h> using namespace std; const int N=2e5+5000; typedef double db; db eps=1e-6; int sgn(db x){ if(fabs(x)<eps)return 0; else return x<0?-1:1; } struct Point{ db x,y; int id; Point(double X=0,double Y=0){x=X;y=Y;} Point operator +(Point B){return Point(x+B.x,y+B.y);} Point operator -(Point B){return Point(x-B.x,y-B.y);} Point operator *(double k){return Point(x*k,y*k);} Point operator /(double k){return Point(x/k,y/k);} friend bool operator<(Point A,Point B){ return sgn(A.x-B.x)<0||(sgn(A.x-B.x)==0&&sgn(A.y-B.y)<0); } friend bool operator==(Point A,Point B){ return sgn(A.x-B.x)==0&&sgn(A.y-B.y)==0; } }; int n; bool cmp(Point A,Point B){ // return sgn(A.x-B.x)<0||(sgn(A.x-B.x)==0&&sgn(A.y-B.y)<0); return A==B?A.id<B.id:A<B; } Point P[N],ch[N]; typedef Point Vector; db Cross(Vector A,Vector B){return A.x*B.y-A.y*B.x;} db Distance(Point A,Point B){return hypot(A.x-B.x,A.y-B.y);} // vector<int>ans; int convexHull(){ sort(P,P+n,cmp); n=unique(P,P+n)-P; int cnt=0; for(int i=0;i<n;i++){ while(cnt>=2){ int d=sgn(Cross(ch[cnt-2]-ch[cnt-1],ch[cnt-1]-P[i])); if(d>0|| (d==0&&P[i].id<ch[cnt-1].id)) cnt--; else break; } ch[cnt++]=P[i]; } for(int i=0;i<cnt;i++){ printf("%d%c",ch[i].id," \n"[i==cnt-1]); } // for(int i=0;i<n;i++){ // while(cnt>1&&sgn(Cross(ch[cnt-1]-ch[cnt-2],P[i]-ch[cnt-2]))<=0)cnt--; // ch[cnt++]=P[i]; // } // int j=cnt; // for(int i=n-2;i>=0;i--){ // while(cnt>j&&sgn(Cross(ch[cnt-1]-ch[cnt-2],P[i]-ch[cnt-2]))<0)cnt--; // ch[cnt++]=P[i]; // ans.push_back(P[i].id); // } // reverse(ans.begin(),ans.end()); // ans.push_back(n); // if(n>1)cnt--; // return cnt; } int main(){ int t;cin>>t; while(t--){ scanf("%d",&n); // ans.clear(); for(int i=0;i<n;i++)scanf("%lf %lf",&P[i].x,&P[i].y),P[i].id=i+1; // cout<<convexHull()<<endl; convexHull(); // for(int i=0;i<ans.size();i++)printf("%d%c",ans[i]," \n"[i==ans.size()-1]); // printf("%.2lf\n",ans); } // system("pause"); return 0; }
Strange fuction
给出一个函数，求最小值。
对函数求导，再求导，发现函数是个先减后增的函数，直接三分。
#include<bits/stdc++.h> using namespace std; const int N=2e5+5000; typedef double db; db eps=1e-6; // int sgn(db x){ // if(fabs(x)<eps)return 0; // else return x<0?-1:1; // } // struct Point{ // db x,y; // Point(double X=0,double Y=0){x=X;y=Y;} // Point operator +(Point B){return Point(x+B.x,y+B.y);} // Point operator -(Point B){return Point(x-B.x,y-B.y);} // Point operator *(double k){return Point(x*k,y*k);} // Point operator /(double k){return Point(x/k,y/k);} // friend bool operator<(Point A,Point B){ // return sgn(A.x-B.x)<0||(sgn(A.x-B.x)==0&&sgn(A.y-B.y)<0); // } // friend bool operator==(Point A,Point B){ // return sgn(A.x-B.x)==0&&sgn(A.y-B.y)==0; // } // }; // int n; // bool cmp(Point A,Point B){ // return sgn(A.x-B.x)<0||(sgn(A.x-B.x)==0&&sgn(A.y-B.y)<0); // } // Point P[N],ch[N]; // typedef Point Vector; // db Cross(Vector A,Vector B){return A.x*B.y-A.y*B.x;} // db Distance(Point A,Point B){return hypot(A.x-B.x,A.y-B.y);} db y; db func(db x){ return 6*pow(x,7)+8*pow(x,6)+7*pow(x,3)+5*pow(x,2)-y*x; } int main(){ int t;cin>>t; while(t--){ scanf("%lf",&y); db ans=1e18,l=0,r=100; while(r>l+eps){ db k=(r-l)/3; db mid1=l+k,mid2=r-k; if(func(mid1)>=func(mid2)){ l=mid1; ans=min(ans,func(mid2)); } else r=mid2; } printf("%.4lf\n",func(l)); } // system("pause"); return 0; }
模拟退火做法，
每次随机跳一个点，更优则加入答案，否则设概率为$p=e^{\frac{E}{T}}$，协调E，只要p<1,即可，随机跳进去。
神奇的是降温系数r=0.98，其他是奇奇怪怪的东西。
#include<bits/stdc++.h> using namespace std; const int N=2e5+5000; typedef double db; db eps=1e-8; db y; db func(db x){ return 6*pow(x,7)+8*pow(x,6)+7*pow(x,3)+5*pow(x,2)-y*x; } int dir[]={1,-1}; db solve(){ db T=100; db r=0.98; db now=50,next,ans=1e18,nowans=func(50); while(T>eps){ next=now+dir[rand()%2]*T/10; if(next>100||next<0)continue; db nextans=func(next); db E=nowans-nextans; if(E>eps){ nowans=nextans; now=next; ans=min(ans,nowans); } else if(exp(E/T)*RAND_MAX>rand()){ nowans=nextans; now=next; ans=min(ans,nowans); } T*=r; } return ans; } int main(){ int t; cin>>t; while(t--){ scanf("%lf",&y); printf("%.4lf\n",solve()); } system("pause"); return 0; }
Segment set
每次判断两个线段是否相交，加入并查集即可。注意端点相交的情况。
#include<bits/stdc++.h> using namespace std; #define pb push_back const int N=1e3+500; typedef double db; const db eps=1e-6; const db PI=acos(-1.0); int sgn(db x){ if(fabs(x)<eps)return 0; else return x<0?-1:1; } struct Point{ db x,y; Point(double X=0,double Y=0){x=X;y=Y;} Point operator +(Point B){return Point(x+B.x,y+B.y);} Point operator -(Point B){return Point(x-B.x,y-B.y);} Point operator *(double k){return Point(x*k,y*k);} Point operator /(double k){return Point(x/k,y/k);} friend bool operator<(Point A,Point B){ return sgn(A.x-B.x)<0||(sgn(A.x-B.x)==0&&sgn(A.y-B.y)<0); } friend bool operator==(Point A,Point B){ return sgn(A.x-B.x)==0&&sgn(A.y-B.y)==0; } }; typedef Point Vector; struct Line{ Point p1,p2;//线上的两个点 Line(){} Line(Point p1,Point p2):p1(p1),p2(p2){} // Line(Point x,Point y){p1 = x;p2 = y;} // Point(double x,double y):x(x),y(y){} //根据一个点和倾斜角 angle 确定直线,0<=angle<pi Line(Point p,double angle){ p1 = p; if(sgn(angle - PI/2) == 0){p2 = (p1 + Point(0,1));} else{p2 = (p1 + Point(1,tan(angle)));} } //ax+by+c=0 Line(double a,double b,double c){ if(sgn(a) == 0){ p1 = Point(0,-c/b); p2 = Point(1,-c/b); } else if(sgn(b) == 0){ p1 = Point(-c/a,0); p2 = Point(-c/a,1); } else{ p1 = Point(0,-c/b); p2 = Point(1,(-c-a)/b); } } }; typedef Line Segment; //定义线段，两端点是Point p1,p2 double Cross(Vector A,Vector B){return A.x*B.y-A.y*B.x;} db Distance(Point A,Point B){return hypot(A.x-B.x,A.y-B.y);} // bool Cross_segment(Point a,Point b,Point c,Point d){//Line1:ab, Line2:cd // double c1=Cross(b-a,c-a),c2=Cross(b-a,d-a); // double d1=Cross(d-c,a-c),d2=Cross(d-c,b-c); // return sgn(c1)*sgn(c2)<0 && sgn(d1)*sgn(d2)<0;//注意交点是端点的情况不算在内 // } bool Cross_segment(Point a,Point b,Point c,Point d) { Point p1=a,p2=b,p3=c,p4=d; if( min(p1.x,p2.x)>max(p3.x,p4.x) || min(p1.y,p2.y)>max(p3.y,p4.y) || min(p3.x,p4.x)>max(p1.x,p2.x) || min(p3.y,p4.y)>max(p1.y,p2.y) ) return 0; double k1,k2,k3,k4; k1=(p2.x-p1.x)*(p3.y-p1.y) - (p2.y-p1.y)*(p3.x-p1.x); k2=(p2.x-p1.x)*(p4.y-p1.y) - (p2.y-p1.y)*(p4.x-p1.x); k3=(p4.x-p3.x)*(p1.y-p3.y) - (p4.y-p3.y)*(p1.x-p3.x); k4=(p4.x-p3.x)*(p2.y-p3.y) - (p4.y-p3.y)*(p2.x-p3.x); return (k1*k2<=eps && k3*k4<=eps); } // Point Cross_point(Line a,Line b){ // Vector u=a.p-b.p; // db t=Cross(b.v,u)/Cross(a.v,b.v); // return a.p+a.v*t; // } Line lines[N]; int fa[N],size[N]; int find(int x){return fa[x]==x?x:fa[x]=find(fa[x]);} void build(int x,int y){ int dx=find(x),dy=find(y); if(dx!=dy){ fa[dx]=dy; size[dy]+=size[dx]; } } int main(){ int t; cin>>t; while(t--){ int n;scanf("%d",&n); int cnt=1; for(int i=1;i<=n;i++)fa[i]=i,size[i]=1; while(n--){ char op[10];scanf("%s",&op); if(op[0]=='P'){ scanf("%lf %lf %lf %lf",&lines[cnt].p1.x,&lines[cnt].p1.y,&lines[cnt].p2.x,&lines[cnt].p2.y); for(int i=0;i<cnt;i++){ if(Cross_segment(lines[i].p1,lines[i].p2,lines[cnt].p1,lines[cnt].p2)){ if(find(cnt)==find(i))continue; build(find(i),cnt); } } cnt++; } else { int x;scanf("%d",&x); printf("%d\n",size[find(x)]); } } if(t)printf("\n"); } // system("pause"); return 0; }
想的太多，做的太少；
