这是一个小透明辣鸡 的【计算几何】屑模板
杂七杂八（头文件啥的）
1 #include <iostream> 2 #include <cmath> 3 #include <cstring> 4 #include <cstdio> 5 #include <cstdlib> 6 #include <algorithm> 7 #include <vector> 8 #include <deque> 9 #include <queue> 10 #include <cassert> 11 #include <set> 12 #include <map> 13 #define mp make_pair 14 #define fi first 15 #define se second 16 #define pb push_back 17 #define f(c,a,b) for(int c=a; c<=b; c++) 18 19 using namespace std; 20 typedef double db; 21 const db eps=1e-8; 22 const db pi=acos(-1); 23 struct P{ 24 db x, y; 25 P operator - (P a){ return (P){x-a.x, y-a.y}; } 26 db dot(P a){ return a.x*x + a.y*y; } 27 P operator * (db a) { return (P){a*x, a*y}; } 28 P operator + (P a){ return (P){x+a.x, y+a.y}; } 29 P operator / (db a) { return (P) {x/a, y/a}; } 30 db times(P a){ return x*a.y-y*a.x; } 31 db le() { return sqrt( x*x+y*y ); } 32 db le2() { return x*x+y*y; } 33 P rot90(){ return (P){-y,x}; } 34 db rad(P a) { return atan2(times(a),dot(a)); } 35 void pri(){ printf("%.1lf %.1lf ", x, y); } 36 void swap(P& a){ 37 db tx = x, ty = y; 38 x = a.x; y = a.y; 39 a.x = tx; a.y = ty; 40 } 41 bool operator < (const P a) const { return (x==a.x)?(y<a.y):(x<a.x); } 42 bool operator == (const P a) const { return (x==a.x)&&(y==a.y); } 43 }po[nmax]; 44 struct L{ 45 P p, q; //p->q取左面 46 db arc; 47 bool inc(P a){ return (q-p).times(a-p)>0; } //必须严格大于，不然会RE（待填） 48 void carc(){ arc = atan2((q-p).y, (q-p).x); } 49 }l[nmax]; 50 int sign(db x) { return (x<-eps) ? -1 : x>eps; } 51 bool operator < (L a, L b){ return sign(a.arc-b.arc)? a.arc<b.arc: b.inc(a.p) ; }
基础操作
求点到直线的投影点
1 P projection(P a, L l){ //求得点a在l上的投影点坐标 2 P tl = l.q - l.p; //表示l的两个点形成的线段 3 db ll = tl.len(); //线段长度 4 db k = (a-l.p).dot(tl)/ll/ll; //投影长度与线段长度的比值 5 return P( l.p.x + k*tl.x, l.p.y + k*tl.y ); //最终结果 6 }
求点关于直线的对称点
1 P reflection(P a, L l){ //点a关于点l的对称点 2 P tmp = projection(a, l); 3 return a+(tmp-a)*2; 4 }
求距离
点和直线：
1 db disLP(P p1, P p2, P p){ return abs( (p1-p).times(p2-p) ) / (p1-p2).le(); } 2 db disLP(L l, P p){ return disLP(l.p, l.q, p); }
点和线段：
1 db disSP(P p1, P p2, P p){ 2 if((p1-p2).dot(p-p2) < eps) return (p-p2).le(); 3 if((p2-p1).dot(p-p1) < eps) return (p-p1).le(); 4 return disLP(p1, p2, p); 5 }
线段和线段
1 db disSS(P p1, P q1, P p2, P q2){ 2 if( intersection(p1, q1, p2, q2) ) return 0; 3 db a = disSP(p1,q1,p2); 4 db b = disSP(p1,q1,q2); 5 db c = disSP(p2,q2,p1); 6 db d = disSP(p2,q2,q1); 7 return min(min(a,b),min(c,d) ); 8 }
求交点
判断两线段是否相交
1 bool intersection(P p1, P q1, P p2, P q2){ 2 if(max(p1.x,q1.x)+eps < min(p2.x,q2.x)) return false; 3 if(max(p1.y,q1.y)+eps < min(p2.y,q2.y)) return false; 4 if(max(p2.x,q2.x)+eps < min(p1.x,q1.x)) return false; 5 if(max(p2.y,q2.y)+eps < min(p1.y,q1.y)) return false; 6 int t1=sign( (p2-p1).times(q1-p1) )*sign( (q1-p1).times(q2-p1) ), 7 t2=sign( (p1-p2).times(q2-p2) )*sign( (q2-p2).times(q1-p2) ); 8 return (t1>=0) && (t2>=0); 9 }
求两线段交点
1 P cpSS(P p1, P q1, P p2, P q2){ 2 P p=p2-p1, l2=q2-p2, l1=q1-p1; 3 if( !sign( p.times(l2) ) ) return p1; 4 P tmp = l1* abs( p.times(l2)/l2.times(l1) ) ; 5 return p1+tmp; 6 }
两直线交点
1 P isLL(P p1, P q1, P p2, P q2){ return p1+(q1-p1)*( (q2-p1).times(p2-q2)/(q1-p1).times(p2-q2) ); } 2 P isLL(L l1, L l2){ return isLL(l1.p,l1.q,l2.p,l2.q); }
直线和圆交点
1 void isCL(P c1, db r1, P p1, P q1){ 2 if(q1<p1) p1.swap(q1); //p1的x值小于q1 3 db dis = (p1-c1).times(p1-q1)/(p1-q1).clen(); 4 dis *= (p1-c1).times(p1-q1); 5 dis = sqrt( r1*r1-dis ); 6 P pj = projection(c1, p1, q1); 7 P a1, a2; 8 a1 = pj+(p1-q1)*( dis/sqrt((p1-q1).clen()) ); 9 a2 = pj+(q1-p1)*( dis/sqrt((q1-p1).clen()) ); 10 printf("%.8lf %.8lf %.8lf %.8lf\n", a1.x, a1.y, a2.x, a2.y); 11 }
两圆交点
1 void isCC(P c1, db r1, P c2, db r2){ 2 if( (c2-c1).y>0 ) { c1.swap(c2); swap(r1,r2); } 3 db tl = (c1-c2).l2(); 4 db xk = ( (r1*r1-r2*r2)/tl + 1)/2; 5 db yk = r1*r1/tl-xk*xk; 6 if(yk<0) yk=0; 7 P tp = c1+(c2-c1)*xk; 8 P lp = (c2-c1).rot90()*sqrt(yk); 9 printf("%.8lf %.8lf %.8lf %.8lf\n", (tp-lp).x, (tp-lp).y, (tp+lp).x, (tp+lp).y ); 10 }
平面最近点对
1 P p[nmax], t[nmax]; 2 int n; 3 bool c_x(P& a, P& b) { return a.x < b.x; } 4 bool c_y(P& a, P& b) { return a.y < b.y; } 5 6 db closestP(int l, int r){ 7 db ta = inf; 8 if(r-l<=4){ 9 f(i,l,r) f(j,i+1,r) ta = min(ta, (p[i]-p[j]).clen() ); 10 sort(p+l, p+r+1, c_y); 11 }else{ 12 int mid = (l+r)>>1; 13 db mx = p[mid].x, d; 14 ta = d = min( closestP(l,mid), closestP(mid+1,r) ); 15 merge(p+l, p+mid+1, p+mid+1, p+r+1, t, c_y); 16 copy(t, t+r-l+1, p+l); 17 int jsz=0; 18 f(i,l,r) if( abs(p[i].x-mx) < d ) { 19 for (int j=jsz; j>=1 && p[i].y-t[j].y<d; j--) ta = min(ta, (t[j]-p[i]).clen() ); 20 t[++jsz] = p[i]; 21 } 22 } 23 return ta; 24 }
半平面交
1 vector <L> getHL(vector<L>& l){ 2 sort(l.begin(), l.end()); 3 deque <L> q; 4 f(i,0,l.size()-1){ 5 if(i && sign(l[i].arc-l[i-1].arc)==0) continue; 6 while(q.size()>1 && !l[i].inc( isLL(q[q.size()-1],q[q.size()-2]) ) ) q.pop_back(); 7 while(q.size()>1 && !l[i].inc( isLL(q[0],q[1]) ) ) q.pop_front(); 8 if(q.size()==1 && sign( (q[0].q-q[0].p).times(l[i].q-l[i].p) ) <=0 ) break; 9 q.push_back(l[i]); 10 } 11 while(q.size()>2 && !q[0].inc( isLL(q[q.size()-1],q[q.size()-2]) ) ) q.pop_back(); 12 vector<L> ans; 13 if(q.size()>1) f(i,0,q.size()-1) ans.push_back(q[i]); 14 return ans; 15 }
多边形相关
点在多边形内外判断
1 int contain(P a){ 2 bool flag=false; 3 f(i,1,n){ 4 P b = po[i]-a, c = po[i-1]-a; 5 if( !sign(b.times(c)) && b.dot(c)<eps ) return 1; 6 if(b.y < c.y) { P t=b; b=c; c=t; } 7 if(c.times(b)>eps && b.y>eps && c.y<eps ) flag=!flag; 8 } 9 return flag?2:0; 10 }
求面积
1 db area(){ 2 db ans = 0; 3 f(i,3,n) ans += (po[i-1]-po[1]).times(po[i]-po[1]); 4 return ans/2; 5 }
判断是不是凸多边形
1 int inc(P c1, db r1, P c2, db r2){ 2 db dis = (c1-c2).clen(); 3 int x1 = sign( (r1+r2)*(r1+r2)-dis ), x2=sign( (r1-r2)*(r1-r2)-dis ); 4 if(x1 == 0) return 3; 5 else if(x1 == -1) return 4; 6 else { 7 if(x2 == 0) return 1; 8 else if(x2 == 1) return 0; 9 else return 2; 10 } 11 }
凸包相关
求凸包
1 void andrew(){ 2 ps[++is]=po[1]; 3 f(i,2,n){ 4 while( is>=2 && (ps[is]-ps[is-1]).times(po[i]-ps[is])<0 ) is--; 5 ps[++is]=po[i]; 6 } 7 for(int i=n-1; i>=1; i--){ 8 while( is>=2 && (ps[is]-ps[is-1]).times(po[i]-ps[is])<0 ) is--; 9 ps[++is]=po[i]; 10 } 11 }
直线切凸包求剩下部分面积
1 db convexcut(P p, P q){//p->q 右边舍弃 2 int al = -1; 3 f(i,0,n-1){ 4 int t1 = sign( (q-p).times(po[i]-p) ); 5 int t2 = sign( (q-p).times(po[(i+1)%n]-p) ); 6 if(t1 >= 0) pa[++al] = po[i]; 7 if(t1*t2 < 0) pa[++al] = isLL(po[i],po[(i+1)%n],p,q); 8 } 9 return area(al,pa); 10 }
凸包直径（旋转卡壳）
1 db Diameter(){ 2 int lm=0, rm=1; 3 f(i,0,n-1){ 4 if(po[lm].x > po[i].x) lm=i; 5 if(po[rm].x < po[i].x) rm=i; 6 } 7 int i=rm, j=lm; 8 db ans = 0; 9 while(i!=lm || j!=rm){ 10 ans = max(ans, (po[j]-po[i]).clen()); 11 if( (po[(i+1)%n]-po[i]).times(po[(j+1)%n]-po[j])<=0 ) i=(i+1)%n; else j=(j+1)%n; 12 } 13 return ans; 14 }
圆相关
过点求圆切线
1 void tanCP(P c1, db r1, P p){ 2 db x = (c1-p).l2(); 3 db d = x - r1*r1; 4 P tp = p + (c1-p)*(d/x); 5 P lp = (c1-p).rot90()*( r1*sqrt(d)/x ); 6 P a1 = tp-lp, a2 = tp+lp; 7 if(a2 < a1) a1.swap(a2); 8 printf("%.10lf %.10lf\n%.10lf %.10lf\n", a1.x, a1.y, a2.x, a2.y ); 9 }
两圆公切线数量
1 int inC(P c1, db r1, P c2, db r2){ 2 db dis = (c1-c2).le2(); 3 int x1 = sign( (r1+r2)*(r1+r2)-dis ), x2=sign( (r1-r2)*(r1-r2)-dis ); 4 if(x1 == 0) return 3; 5 else if(x1 == -1) return 4; 6 else { 7 if(x2 == 0) return 1; 8 else if(x2 == 1) return 0; 9 else return 2; 10 } 11 }
多边形和圆相交面积
三角形和圆相交面积之和两圆位置关系（外切，相交等）判断
输出切线数量表示位置关系求两圆公切线（返回点坐标）
1 vector<P> tanCC(P c1, db r1, P c2, db r2, int qx){ 2 vector<P> ta; 3 if(!qx) return ta; 4 if( sign(r1-r2)==0 ) { 5 P tp = ( (c2-c1)/(c2-c1).le()*r1 ).rot90(); 6 ta = { c1+tp,c1-tp }; 7 }else{ 8 P tp = c1 + (c2-c1)*(r1)/(r1-r2); 9 ta = tanCP(c1, r1, tp); 10 } 11 //cout<<qx<<endl; 12 if(qx<=2) return ta; 13 P tp = c1 + (c2-c1)*r1/(r1+r2); 14 vector<P> tmp = tanCP(c1, r1, tp); 15 ta.insert(ta.end(), tmp.begin(), tmp.end()); 16 return ta; 17 }
三角形
三角形外接圆半径和圆心
1 db triarea(db a, db b, db c) { db t=(a+b+c)/2; return sqrt(t*(t-a)*(t-b)*(t-c)); } //三角形面积 2 db triccr(db a, db b, db c) { return a*b*c/4/triarea(a,b,c); } //外接圆半径 3 P triccp(P a, P b, P c) { //外接圆圆心 4 db a1 = 2*(b.x-a.x), b1 = 2*(b.y-a.y), c1 = (b.x*b.x+b.y*b.y-a.x*a.x-a.y*a.y); 5 db a2 = 2*(c.x-b.x), b2 = 2*(c.y-b.y), c2 = (c.x*c.x+c.y*c.y-b.x*b.x-b.y*b.y); 6 return (P){ (c1*b2-c2*b1)/(a1*b2-a2*b1) , (a1*c2-a2*c1)/(a1*b2-a2*b1) }; 7 }
反演
有个圆，圆心为$P$，半径为$R$，平面上某点对于这个圆的反演点为：设该点到圆心的距离反演前为$lbefore$，反演后为$lafter$。则有$lbefore*lafter=R^2$
直线反演为过P点的圆：
1 C fyltoc(L l, C fy) { 2 P tp = projection(fy.c, l); 3 db dis = (fy.c-tp).le(); 4 db fr = fy.r*fy.r/(2*dis); 5 return (C){ fy.c+(tp-fy.c)/dis*fr, fr }; 6 }
过P点的圆反演为直线：
1 L fyctol(C c1, C fy){ return isCC(c1, fy); }
不过P点的圆反演为不过P点的圆：
1 C fyctoc(C c, C fy){ 2 db l=(c.c-fy.c).le2(); 3 db yx = sqrt(l)*fy.r*fy.r/(l-c.r*c.r); 4 return (C){ fy.c+(c.c-fy.c)/(c.c-fy.c).le()*yx , c.r*fy.r*fy.r/(l-c.r*c.r) }; 5 }
