目錄
凸包
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π\)。
旋轉卡殼
可以求凸包上的最值等問題。
例如：最遠點對，最小外界矩形，最大四邊形面積等。
對於最遠點對，就是先枚舉一條邊，然后用與其平行的直線切割凸包的另一面。
如果按照順時針枚舉這條邊，那么另一個點也是順時針轉動的，可以使用雙指針。
平行的切線可以使用叉積判斷（面積最大）。
最遠點對
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+c\geq0\)，也可以轉化為直線。
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", 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\n", 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過大/過小。