凸包
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过大/过小。
总之，注意精度。
万不得已可以手写分数类提高精度。