计算几何的特点:
- 要么全场会要么全场不会。
- 一道题100行,其中90行是模板。所以模板必须可靠,写的时候必须注意结构。
- 注意精度,能用整数尽量用整数。
基本定义:
1.向量:参考初中数学教材。
意义:计算几何的单位不是点而是向量。
2.点积:$acdot b=|a||b|cos<a,b>=a.x*b.x+a.y*b.y$。
意义:若$acdot b=0$则两向量垂直,若$acdot b<0$则两向量夹角$>frac{pi}{2}$,若$acdot b>0$则两向量夹角$<frac{pi}{2}$。
3.叉积:$a imes b=|a||b|sin<a,b>=a.x*b.y-a.y*b.x$。
意义:若$a imes b=0$则两向量共线,若$a imes b>0$则b在a的逆时针方向,若$a imes b<0$则b在a的顺时针方向。
基础操作:
1.向量旋转:$(rcosA,rsinA) ightarrow (rcos(A+B),rsin(A+B))$
有$(rcos(A+B),rsin(A+B))=(rcosAcosB-rsinAsinB,rsinAcosB+rcosAsinB)$
将坐标带入得到$(x,y) ightarrow (xcosB-ysinB,ycosB+xsinB)$。
2.极角排序:将所有点按顺时针/逆时针排序。(有些时候要分象限讨论)
方法一:对于点$(x,y)$使用$atan2(y,x)$即可求出极角(弧度制),但精度堪忧。
方法二:使用叉积判断,较为常用。
3.两直线夹角:求直线$a,b$相交形成的角度$ heta$。
由于$frac{a imes b}{acdot b}=frac{sin heta }{cos heta }=tan heta$,直接$atan2(a imes b,acdot b)$即可。
4.两直线交点:求直线$a,b$的交点$p$。
设点$a1,a2$确定直线$a$,点$b1,b2$确定直线$b$。
画个图发现如果能确定$frac{|a1p|}{|a1a2|}$的值就可以了。
利用叉积代表面积的性质可以得到做法:
vec ins(line a,line b){ vec u=b1-a1,v=a2-a1,w=b2-b1; double k=(w*u)/(v*w); return a1+k*v; }
5.两线段交点:两次跨立实验判一下是不是相交。
跨立实验:$(b1-a1) imes (b1-a2)$与$(b2-a1) imes (b2-a2)$不同号,则线段$b$跨越直线$a$。
多边形相关:
1.多边形面积:$S=frac{1}{2}sum limits_{i=1}^{n}P_i imes P_{i+1}$。
画一下可以发现多边形外到原点的面积恰好被算了一次正的的和一次负的,于是抵消。
2.判断点是否在多边形内:从该点出发引一条射线,如果该射线与多边形边界相交奇数次则该点在多边形内。
3.求凸包:$Graham$扫描法。
步骤:
- 找到左下角的点,以该点为原点建系,将其他点极角排序。用队列维护构成当前凸包的点。
- 依次加入点,如果新加入的点与队列右端点构成的向量在上一条凸包边的顺时针方向,那么需要弹出上一个点。
- 最终将队列中所有点顺次连接构成凸包。
代码:
bool cmp1(node a,node b){return a.x==b.x?a.y<b.y:a.x<b.x;} bool cmp2(node a,node b){return a*b>0||(a*b==0&&dis(a)<dis(b));} void Graham(){ sort(A+1,A+1+n,cmp1); for(int i=n;i>=1;i--) A[i]=A[i]-A[1]; sort(A+2,A+1+n,cmp2); node sta[maxn]; int top=0; for(int i=1;i<=n;sta[++top]=A[i++]) while(top>=2 && (A[i]-sta[top])*(sta[top]-sta[top-1])>=0) top--; n=top; for(int i=1;i<=top;i++) A[i]=sta[i]; }
4.求若干条有向直线左侧平面的交:半平面交。
步骤:
- 以斜率为第一关键字,截距为第二关键字将直线排序。用双端队列维护构成当前半平面交的直线。
- 依次加入直线,如果队列右侧/左侧两条直线的交点不在新加入直线的左侧,那么需要弹出队列右端点/左端点。
- 最终将队列中所有直线顺次拼接构成半平面交。
几何意义:多边形内核(多边形内部的一块区域,站在该区域内任何一个点都可以看到多边形的任何一个角落)。
代码([CQOI2006]凸多边形):
#include<bits/stdc++.h> #define maxn 200005 #define maxm 500005 #define inf 0x7fffffff #define eps 1e-8 #define ll long long #define rint register int #define debug(x) cerr<<#x<<": "<<x<<endl #define fgx cerr<<"--------------"<<endl #define dgx cerr<<"=============="<<endl using namespace std; struct node{ double x,y; node operator+(const node t)const{return (node){x+t.x,y+t.y};} node operator-(const node t)const{return (node){x-t.x,y-t.y};} node operator^(const double t)const{return (node){x*t,y*t};} double operator*(const node t)const{return x*t.y-y*t.x;} }A[maxn],P[maxn]; struct line{ node a,b; double k,d; bool operator<(const line t)const{return (abs(k-t.k)<eps)?((b-a)*(t.b-a)>eps):(k-t.k<eps);} }que[maxn],L[maxn]; inline int read(){ int x=0,f=1; char c=getchar(); for(;!isdigit(c);c=getchar()) if(c=='-') f=-1; for(;isdigit(c);c=getchar()) x=x*10+c-'0'; return x*f; } inline node ins(line l1,line l2){ node u=l2.a-l1.a,v=l1.b-l1.a,w=l2.b-l2.a; double k=(u*v)/(v*w); return l2.a+(w^k); } inline bool onr(line l,node c){return (l.b-c)*(l.a-c)>eps;} int main(){ int T=read(),n=0; while(T--){ int m=read(); for(int i=1;i<=m;i++) A[i].x=read(),A[i].y=read(); for(int i=1;i<=m;i++){ L[++n].a=A[i],L[n].b=(i==m)?A[1]:A[i+1]; L[n].k=atan2((L[n].b-L[n].a).y,(L[n].b-L[n].a).x); L[n].d=L[n].a.y-L[n].a.x*L[n].k; } } sort(L+1,L+1+n); int l=1,r=0,tot=0,cnt=0; for(int i=1;i<=n;i++) if(abs(L[i].k-L[i+1].k)>eps) L[++cnt]=L[i]; for(int i=1;i<=cnt;i++){ while(l<r && onr(L[i],ins(que[r-1],que[r]))) r--; while(l<r && onr(L[i],ins(que[l],que[l+1]))) l++; que[++r]=L[i]; } while(l<r && onr(que[l],ins(que[r-1],que[r]))) r--; while(l<r && onr(que[r],ins(que[l],que[l+1]))) l++; for(int i=l;i<=r;i++){ line l1=que[i],l2=(i==r)?que[l]:que[i+1]; P[++tot]=ins(l1,l2); } double S=0; for(int i=1;i<=tot;i++) S+=P[i]*((i==tot)?P[1]:P[i+1]); printf("%.3lf ",(tot<=2)?0:(S*0.5)); return 0; }
5.求凸包直径(凸包上最远点对的距离):旋转卡壳。
步骤:
- 求出凸包。
- 按极角序枚举每条边,找到某个顶点使得该边与该点构成的三角形面积最大(注意特判面积相等时的情况),用该顶点与该边两端点的距离更新答案。
- 容易证明这样一定能求出上界,且复杂度为$O(nlog{n})$(单调的边对应的顶点也是单调的)。
代码([USACO03FALL]Beauty Contest G):
#include<bits/stdc++.h> #define maxn 200005 #define maxm 500005 #define inf 0x7fffffff #define ll long long #define rint register ll #define debug(x) cerr<<#x<<": "<<x<<endl #define fgx cerr<<"--------------"<<endl #define dgx cerr<<"=============="<<endl using namespace std; struct node{ ll x,y; inline node operator+(const node t)const{return (node){x+t.x,y+t.y};} inline node operator-(const node t)const{return (node){x-t.x,y-t.y};} inline node operator^(const ll t)const{return (node){x*t,y*t};} inline ll operator*(const node t)const{return x*t.y-y*t.x;} }A[maxn],S[maxn]; ll n; inline ll read(){ ll x=0,f=1; char c=getchar(); for(;!isdigit(c);c=getchar()) if(c=='-') f=-1; for(;isdigit(c);c=getchar()) x=x*10+c-'0'; return x*f; } inline bool cmp1(node a,node b){return (a.x==b.x)?(a.y<b.y):(a.x<b.x);} inline bool cmp2(node a,node b){return (a*b==0)?(a.x<b.x):(a*b>0);} inline ll dis(node a,node b){return (a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y);} inline void Graham(){ sort(A+1,A+1+n,cmp1); for(ll i=n;i>=1;i--) A[i]=A[i]-A[1]; sort(A+2,A+1+n,cmp2); ll top=0; for(ll i=1;i<=n;i++){ while(top>1 && (A[i]-S[top])*(S[top]-S[top-1])>=0) top--; S[++top]=A[i]; } n=top; for(ll i=1;i<=n;i++) A[i]=S[i]; } int main(){ n=read(); for(ll i=1;i<=n;i++) A[i].x=read(),A[i].y=read(); Graham(); if(n==2){printf("%lld ",dis(A[1],A[2]));return 0;} A[++n]=A[1]; ll res=0; for(ll i=1,j=1;i<n;i++){ while((A[i+1]-A[i])*(A[j]-A[i])<=(A[i+1]-A[i])*(A[j+1]-A[i])) j=(j%(n-1))+1; res=max(res,max(dis(A[j],A[i]),dis(A[j],A[i+1]))); } printf("%lld ",res); return 0; }
面积相关:
1.求$f(x)$积分:自适应$Simpson$法(使用二次函数拟合$f(x)$)。
公式:$int_{a}^{b}f(x)dx approx frac{(b-a)}{6}(f(a)+f(b)+4f(frac{a+b}{2}))$。(可通过简单展开整理得到该式,在此不再进行推导)
步骤:假设当前处理的状态是$(l,r,eps)$(即计算$[l,r]$区间,所需精度为$eps$)。
- 将区间$[l,r]$分成左右两半$[l,mid]$与$[mid,r]$。
- 套用上述公式计算这三个区间的积分,设计算结果为$S,S1,S2$。
- 若$|S1+S2-S|<eps$则返回$S1+S2$,否则继续递归处理$(l,mid,frac{eps}{2})$与$(mid,r,frac{eps}{2})$。
代码:
#include<bits/stdc++.h> #define maxn 200005 #define maxm 500005 #define inf 0x7fffffff #define ll long long #define rint register int #define debug(x) cerr<<#x<<": "<<x<<endl #define fgx cerr<<"--------------"<<endl #define dgx cerr<<"=============="<<endl using namespace std; double a,b,c,d,L,R; inline int read(){ int x=0,f=1; char c=getchar(); for(;!isdigit(c);c=getchar()) if(c=='-') f=-1; for(;isdigit(c);c=getchar()) x=x*10+c-'0'; return x*f; } inline double f(double x){return (c*x+d)/(a*x+b);} inline double simpson(double l,double r){return (r-l)*(f(r)+f(l)+4.0*f((l+r)/2.0))/6.0;} inline double solve(double l,double r,double eps){ double mid=(l+r)/2.0; double S=simpson(l,r),SL=simpson(l,mid),SR=simpson(mid,r); if(abs(SL+SR-S)<eps) return SL+SR; else return solve(l,mid,eps/2.0)+solve(mid,r,eps/2.0); } int main(){ scanf("%lf%lf%lf%lf%lf%lf",&a,&b,&c,&d,&L,&R); printf("%.6lf ",solve(L,R,1e-7)); return 0; }
2.求面积并:扫描线。
大概就是选定一个方向扫,把所有边转换成平行于该方向的入边和出边,然后用个数据结构维护。
比如求矩形的面积并就是线段树区间加减和区间求非0个数。
注意区间求非0个数可以写成标记永久化的形式,以及不要随便用map离散化,慢得离谱。
代码:
#include<bits/stdc++.h> #define maxn 400005 #define maxm 500005 #define inf 0x7fffffff #define ll long long #define rint register int #define debug(x) cerr<<#x<<": "<<x<<endl #define fgx cerr<<"--------------"<<endl #define dgx cerr<<"=============="<<endl using namespace std; struct node{ll l,r,k;}; struct square{ll x1,y1,x2,y2;}S[maxn]; ll T[maxn],mp[maxn],trcv[maxn<<2],trsz[maxn<<2]; vector<node> L[maxn]; inline ll read(){ ll x=0,f=1; char c=getchar(); for(;!isdigit(c);c=getchar()) if(c=='-') f=-1; for(;isdigit(c);c=getchar()) x=x*10+c-'0'; return x*f; } inline void pushup(ll l,ll r,ll k){ if(trcv[k]) trsz[k]=(mp[r+1]-mp[l]); else{ if(l!=r) trsz[k]=trsz[k<<1]+trsz[k<<1|1]; else trsz[k]=0; } } inline void upd(ll x,ll y,ll z,ll l,ll r,ll k){ if(x<=l && r<=y){trcv[k]+=z;pushup(l,r,k);return;} ll mid=l+r>>1; if(x<=mid) upd(x,y,z,l,mid,k<<1); if(y>mid) upd(x,y,z,mid+1,r,k<<1|1); pushup(l,r,k); } int main(){ freopen("1.in","r",stdin); ll n=read(),m=0,ans=0; for(ll i=1;i<=n;i++){ rint x1=read(),y1=read(),x2=read(),y2=read(); T[++m]=x1,T[++m]=y1,T[++m]=x2,T[++m]=y2; S[i].x1=x1,S[i].y1=y1,S[i].x2=x2,S[i].y2=y2; } sort(T+1,T+1+m); for(rint i=1;i<=m;i++) mp[i]=T[i]; mp[m+1]=mp[m]; for(rint i=1;i<=n;i++){ S[i].x1=lower_bound(T+1,T+1+m,S[i].x1)-T; S[i].y1=lower_bound(T+1,T+1+m,S[i].y1)-T; S[i].x2=lower_bound(T+1,T+1+m,S[i].x2)-T; S[i].y2=lower_bound(T+1,T+1+m,S[i].y2)-T; L[S[i].x1].push_back((node){S[i].y1,S[i].y2,1}); L[S[i].x2].push_back((node){S[i].y1,S[i].y2,-1}); } for(rint x=1;x<=m;x++){ ans+=(mp[x]-mp[x-1])*trsz[1]; for(rint i=0;i<L[x].size();i++) {node tp=L[x][i];upd(tp.l,tp.r-1,tp.k,1,m,1);} } printf("%lld ",ans); return 0; }