计算几何总结

目录

• 凸包
• 极角排序
• 旋转卡壳
  • 最远点对
  • 最小外接矩形
• 最小外接圆
  • 代码
• 半平面交
  • 代码
• Simpson积分
  • 代码
• 其他常用操作
  • 向量运算
  • 求交点
  • 求垂足（投影）
  • 旋转
  • 三点确定一个圆
  • 求公切线（一条）
  • 求出所有公切线
• 注意事项

凸包

1、找到最左面的点，如有多个，取最下面的的点。
注：这一步是为了极角排序时方便直接叉积。
2、以这个点为原点，将所有点按与x轴正半轴的夹角为第一关键字（从大到小，顺时针），到原点距离为第二关键字（从小到大）排序。

比较的代码：

int left(SVe a,SVe b)
{
	return sgn(a*b);
}
int cmp(const void*a,const void *b)
{
	int t=left(*(SPt*)a-o,*(SPt*)b-o);
	if(t!=0)
		return t;
	return sgn(getlen(*(SPt*)a-o)-getlen(*(SPt*)b-o));
}

3、用栈维护，保证后一个向量在前一个的右面。

o=sz[0];
for(int i=1;i<n;i++)
{
	if(sgn(sz[i].x-o.x)==-1||(sgn(sz[i].x-o.x)==0&&sgn(sz[i].y-o.y)==-1))
		o=sz[i];
}
qsort(sz,n,sizeof(SPt),cmp);
st[tp++]=o;
for(int i=1;i<n;i++)
{
	while(tp>=2&&left(st[tp-1]-st[tp-2],sz[i]-st[tp-1])==1)
		tp-=1;
	st[tp++]=sz[i];
}

带曲线的凸包：[SHOI2012]信用卡凸包

通过平移，可以发现：所有的曲线正好形成一个圆，所有直线分别向内平移可以形成一个完整凸包。
因此，求出所有圆心组成的凸包即可。
再加上一个圆就是答案。

极角排序

极角：

atan2(y,x);

从x轴负半轴开始，逆时针为增加，从(-π)到(π)。其中上是(π)。
处理时可以类似环的处理方法，将范围变为2倍，从(-π到3π)。

使用(atan2)求极角通常会有精度损失。
另一种方法就是使用叉积。但是注意叉积对于转角超过(180)度时会错。
解决方案：比较时，若分别在(x)轴两侧，直接判断。
否则，使用叉积。因为这是可以保证转角严格小于(180)度。
代码：

int xix(SVe a)
{
    if(a.y<0||(a.y==0&&a.x<0))
        return 0;
    return 1;
}
int cmpj(SVe a,SVe b)
{
    int xa=xix(a),xb=xix(b);
    if(xa!=xb)return xa-xb;
    return cross(a,b);
}

旋转卡壳

可以求凸包上的最值等问题。
例如：最远点对，最小外界矩形，最大四边形面积等。

对于最远点对，就是先枚举一条边，然后用与其平行的直线切割凸包的另一面。
如果按照顺时针枚举这条边，那么另一个点也是顺时针转动的，可以使用双指针。

平行的切线可以使用叉积判断（面积最大）。

最远点对

double getmax(int n,SPt sz[50010])
{
	double ma=0;
	for(int i=0,j=1/*注意此处*/;i<n;i++)
	{
		while(mianj(sz[(j+1)%n],sz[i],sz[(i+1)%n])>mianj(sz[j],sz[i],sz[(i+1)%n]))
			j=(j+1)%n;
		double t=getlen(sz[i]-sz[j]);
		if(t>ma)ma=t;
		t=getlen(sz[(i+1)%n]-sz[(j+1)%n]);//要判断两种情况
		if(t>ma)ma=t;
	}
	return ma;
}

最小外接矩形

平行的边和之前一样。
与枚举的边垂直的边可以使用点积判断（最大化点的投影到端点的距离）。
注意只能垂直于枚举边，不能用对面的边（对面的与凸包仅有一个交点）。

double getans(int n,SPt sz[50010],SPt dd[5])
{
	double ma=1e100;sz[n]=sz[0];
	for(int i=0,j=1,a=1,b=1;i<n;i++)
	{
		while(mianj(sz[j+1],sz[i],sz[i+1])>=mianj(sz[j],sz[i],sz[i+1]))
			j=(j+1)%n;
		while(dianj(sz[a+1]-sz[i],sz[i+1]-sz[i])>=dianj(sz[a]-sz[i],sz[i+1]-sz[i]))
			a=(a+1)%n;
		if(i==0)b=j;
		while(dianj(sz[b+1]-sz[i],sz[i]-sz[i+1])>=dianj(sz[b]-sz[i],sz[i]-sz[i+1]))
			b=(b+1)%n;
		SPt A=getH(sz[i],sz[i+1]-sz[i],sz[a]);
		SPt B=getH(A,getcz(sz[i+1]-sz[i]),sz[j]);
		SPt C=getH(B,getcz(A-B),sz[b]);
		SPt D=getH(sz[i],sz[i+1]-sz[i],sz[b]);
		double S=fabs(mianj(A,B,D));
		if(S<ma)
		{
			dd[0]=D;dd[1]=C;
			dd[2]=B;dd[3]=A;
			ma=S;
		}
	}
	return ma;
}

最小外接圆

有如下步骤：
首先，随机打乱，防止被卡。
将点按照1～n的顺序依次添加。
若当前点i在圆内，则跳过这个点。
若i不在圆内，那么这个点一定在外接圆上，以这个点为圆心，半径为0作为新圆。
枚举1～i-1，若j在圆内则跳过，否则以i，j为直径作为新圆。
枚举1～j-1，若k在圆内则跳过，否则用i，j，k确定一个新圆。
其实就是暴力，但期望(O(n))。

代码

SCr c(sz[0],0);
for(int i=1;i<n;i++)
{
	if(inc(c,sz[i]))
		continue;
	c=SCr(sz[i],0);
	for(int j=0;j<i;j++)
	{
		if(inc(c,sz[j]))
			continue;
		c=SCr(mid(sz[i],sz[j]),getle2(sz[i]-sz[j])/4);//为了方便，用的是平方
		for(int k=0;k<j;k++)
		{
			if(!inc(c,sz[k]))
				c=getcr(sz[i],sz[j],sz[k]);
		}
	}
}

半平面交

就是给出若干半平面，求交集。（一定是凸多边形）
用途很多，比如：

求一些凸多边形的交集 [CQOI2006]凸多边形
求二元一次不等式组的解集 [SCOI2015]小凸想跑步
求多边形的核：即多边形内的一个点集，点集内任意一点与多边形边上任意一点的连线都在多边形内。

求法（假设在直线左面）：

首先，按照极角从小到大排序，平行的线，靠左的排前面。

struct SLi
{
	SPt s,t;double ji;
	SLi(){}
	SLi(SPt S,SPt T)
	{
		s=S;t=T;
		ji=atan2(t.y-s.y,t.x-s.x);
	}
};
int cmp(const void*a2,const void*b2)
{
	SLi a=*(SLi*)a2,b=*(SLi*)b2;
	int t=sgn(a.ji-b.ji);
	if(t!=0)return t;
	return right(b,a.s);
}

然后，依次添加，和之前平行的舍去。
维护两个双端队列，一个保存直线，一个保存交点。
添加分三步：
1、从队尾，把交点在当前直线右面的交点去掉，同时也去掉队尾直线。
2、从队首，把交点在当前直线右面的交点去掉，同时也去掉队首直线。
3、添加当前直线，添加该直线和队列中上一条直线的交点。

最后，用队首直线排除队尾直线，方法同之前。

代码

double bpmj(SLi sz[200010],int n)
{
	qsort(sz,n,sizeof(SLi),cmp);
	int he=0,ta=0;
	for(int i=0;i<n;i++)
	{
		if(i>0&&pingx(sz[i],sz[i-1]))
			continue;
		while(he+1<ta&&right(sz[i],jd[ta-1])==1)
			ta-=1;
		while(he+1<ta&&right(sz[i],jd[he+1])==1)
			he+=1;
		if(he<ta)
			jd[ta]=jiaod(sz[i],dl[ta-1]);
		dl[ta++]=sz[i];
	}
	while(he+1<ta&&right(dl[he],jd[ta-1])==1)
		ta-=1;
	jd[he]=jiaod(dl[he],dl[ta-1]);
	jd[ta]=jd[he];
	double ans=0;
	for(int i=he;i<ta;i++)
		ans+=jd[i]*jd[i+1];
	return ans/2;
}

对于不等式(ax+by+cgeq0)，也可以转化为直线。

SLi::SLi(ll a,ll b,ll c)
{
	SPt s0,t0;
	if(b==0)
	{
		s0=SPt(double(-c)/a,0);
		t0=SPt(double(-c)/a,1);
	}
	else
	{
		s0=SPt(0,double(-c)/b);
		t0=SPt(1,double(-c-a)/b);
	}
	SPt jc=s0;
	if(a>0)jc.x+=1;
	else jc.x-=1;
	if(b>0)jc.y+=1;
	else jc.y-=1;
        //似乎是糟糕的做法？ --> 我也不知道是啥了，降智了
	if(right(SLi(s0,t0),jc)==-1)
		*this=SLi(s0,t0);
	else
		*this=SLi(t0,s0);
}

Simpson积分

用途：求面积时，若组合图形不易分割，但可以算出每个横坐标上的纵坐标值，就可用此方法。
例题 [NOI2005]月下柠檬树

基本思想就是把复杂的函数近似成二次函数。
先考虑二次函数的积分：

这样子直接计算我们发现误差非常大，毕竟原图像可能不能很好的用二次函数逼近
所以，有自适应 Simpson 积分：

发现当前区间分成两份和不分相差不大时，直接返回，否则分两半递归。

代码

double sim(double a,double b)
{
	return (b-a)*(F(a)+F(b)+4*F((a+b)/2))/6;
}
double asr(double a,double b,double z)
{
	double m=(a+b)/2,cl=sim(a,m),cr=sim(m,b);
	if(fabs(cl+cr-z)<eps)
		return z;
	return asr(a,m,cl)+asr(m,b,cr);
}

其他常用操作

向量运算

struct SPt
{
	double x,y;
	SPt(){}
	SPt(double X,double Y)
	{
		x=X;y=Y;
	}
};
typedef SPt SVe;
int sgn(double t)//符号
{	
	if(t>eps)return 1;
	else if(t<-eps)
		return -1;
	else return 0;
}
SVe operator-(SPt a,SPt b)
{
	return SVe(a.x-b.x,a.y-b.y);
}
SVe operator*(SVe a,double b)//数乘
{
	return SVe(a.x*b,a.y*b);
}
SPt operator+(SPt a,SVe b)
{
	return SPt(a.x+b.x,a.y+b.y);
}
double operator*(SVe a,SVe b)//叉积
{
	return a.x*b.y-a.y*b.x;
}
double dianj(SPt a,SPt b)//点积
{
	return a.x*b.x+a.y*b.y;
}
SVe getcz(SVe a)//垂直向量
{
	return SVe(a.y,-a.x);
}
double getlen(SVe a)
{
	return sqrt(a.x*a.x+a.y*a.y);
}

求交点

SPt jiaod(SPt a,SVe va,SPt b,SVe vb)
{
	double t=((b-a)*vb)/(va*vb);
	return a+va*t;
}

求垂足（投影）

SPt getH(SPt O,SVe a,SPt N)
{
	double t=dianj(a,N-O)/dianj(a,a);
	SPt M;
	M.x=O.x+a.x*t;M.y=O.y+a.y*t;
	return M;
}

旋转

SVe rotate(SVe a,double r)
{
	double tx=cos(r),ty=sin(r);
	SVe rt;
	rt.x=a.x*tx-a.y*ty;
	rt.y=a.x*ty+a.y*tx;
	return rt;
}

三点确定一个圆

struct SCr
{
	SPt o;
	double r;
	SCr(){}
	SCr(SPt O,double R)
	{
		o=O;r=R;
	}
};
SCr getcr(SPt a,SPt b,SPt c)
{
	SPt p1=mid(a,b),p2=mid(b,c);
	SPt o=jiaod(p1,getcz(p1-a),p2,getcz(p2-b));
	return SCr(o,getle2(o-a));
}

求公切线（一条）

利用相似

for(int i=0;i<n;i++)
{
	double d=x[i+1]-x[i],a=r[i+1]-r[i],c=sqrt(d*d-a*a);
	double x0=x[i]-r[i]/d*a,y0=r[i]/d*c;
	double x1=x[i+1]-r[i+1]/d*a,y1=r[i+1]/d*c;
	k[i]=(y1-y0)/(x1-x0);
	b[i]=y0-k[i]*x0;
}

求出所有公切线

以一个圆心为原点，将坐标进行变换，另一个圆心放到x轴上。
就可以类似上一个的方法了，还是相似。
最后别忘变换回去。

#include <stdio.h> 
#include <stdlib.h> 
#include <math.h> 
#define eps 1e-6 
struct SPt {
	double x,y;
	SPt() {}
	SPt(double X, double Y) {
		x = X;y = Y;
	}
};
SPt rotate(SPt a, double co, double si) {
	SPt rt;
	rt.x = a.x * co - a.y * si;
	rt.y = a.x * si + a.y * co;
	return rt;
}
struct SLi {
	SPt s,t;
};
void geta2(double r1, double d, double r2, SLi & l1, SLi & l2) {
	double t = r2 - r1,l = sqrt(d * d - t * t);
	double si = l / d,co = t / d;
	l1.s = SPt( - r1 * co, r1 * si);
	l1.t = SPt(d - r2 * co, r2 * si);
	l2.s = SPt( - r1 * co, -r1 * si);
	l2.t = SPt(d - r2 * co, -r2 * si);
}
void getb2(double r1, double d, double r2, SLi & l1, SLi & l2) {
	double t = r2 + r1,l = sqrt(d * d - t * t);
	double si = l / d,co = t / d;
	l1.s = SPt(r1 * co, -r1 * si);
	l1.t = SPt(d - r2 * co, r2 * si);
	l2.s = SPt(r1 * co, r1 * si);
	l2.t = SPt(d - r2 * co, -r2 * si);
}
int getjd(double r1, double d, double r2, SLi li[5]) {
	if (fabs(d) < eps && fabs(r1 - r2) < eps) return - 1;
	double a1 = -r1,b1 = r1,a2 = d - r2,b2 = d + r2;
	if (a1 - a2 > eps || b1 - b2 > eps) return 0;
	if (fabs(a1 - a2) < eps) {
		li[0].s = li[0].t = SPt( - r1, 0);
		return 1;
	}
	if (fabs(b1 - b2) < eps) {
		li[0].s = li[0].t = SPt(r1, 0);
		return 1;
	}
	geta2(r1, d, r2, li[0], li[1]);
	if (b1 - a2 > eps) return 2;
	if (fabs(b1 - a2) < eps) {
		li[2].s = li[2].t = SPt(r1, 0);
		return 3;
	}
	getb2(r1, d, r2, li[2], li[3]);
	return 4;
}
SLi li[5];
int sgn(double x) {
	if (x > eps) return 1;
	if (x < -eps) return - 1;
	return 0;
}
int cmp(const void * a, const void * b) {
	int x = sgn(((SLi * ) a) ->s.x - ((SLi * ) b) ->s.x);
	int y = sgn(((SLi * ) a) ->s.y - ((SLi * ) b) ->s.y);
	if (x != 0) return x;
	return y;
}
int main() {
	while (1) {
		double x1,y1,r1,x2,y2,r2;
		scanf("%lf%lf%lf%lf%lf%lf", &x1, &y1, &r1, &x2, &y2, &r2);
		if (fabs(r1) < eps) break;
		double d = sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2));
		double a = x2 - x1,b = y2 - y1,co = a / d,si = b / d;
		int n = getjd(r1, d, r2, li);
		printf("%d
", n);
		if (n == -1) n = 0;
		for (int i = 0; i < n; i++) {
			li[i].s = rotate(li[i].s, co, si);
			li[i].t = rotate(li[i].t, co, si);
		}
		qsort(li, n, sizeof(SLi), cmp);
		for (int i = 0; i < n; i++) {
			double a = li[i].s.x - li[i].t.x;
			double b = li[i].s.y - li[i].t.y;
			printf("%.5lf %.5lf %.5lf %.5lf %.5lf
", li[i].s.x + x1, li[i].s.y + y1, li[i].t.x + x1, li[i].t.y + y1, sqrt(a * a + b * b));
		}
	}
	return 0;
}

注意事项

1、若比较的符号带等于，则需要使用eps判断。
2、注意输出-0的情况。

void output(double x)
{
	if(fabs(x)<eps)
		x=0;
	printf("%.5lf ",x);
}

3、若出现奇怪错误，可能是数组开小，eps过大/过小。

总之，注意精度。
万不得已可以手写分数类提高精度。