向量的板子
#include<bits/stdc++.h>
#define I inline
using namespace std;
typedef double DB;
struct Vec{
DB x,y;
I Vec(){x=y=0;}
I Vec(DB a){x=a;y=0;}
I Vec(DB a,DB b){x=a;y=b;}
I friend istream&operator>>(istream&cin,Vec&a){return cin>>a.x>>a.y;}
I friend ostream&operator<<(ostream&cou,Vec a){return cou<<a.x<<' '<<a.y;}
I Vec operator-(){return Vec(-x,-y);}
I Vec operator+(Vec a){return Vec(x+a.x,y+a.y);}
I Vec operator-(Vec a){return Vec(x-a.x,y-a.y);}
I Vec operator*(DB a){return Vec(x*a,y*a);}
I Vec operator/(DB a){return Vec(x/a,y/a);}
I friend Vec operator*(DB a,Vec b){b.x*=a,b.y*=a;return b;}
I friend Vec operator/(DB a,Vec b){b.x/=a,b.y/=a;return b;}
I Vec&operator+=(Vec a){x+=a.x,y+=a.y;return*this;}
I Vec&operator-=(Vec a){x-=a.x,y-=a.y;return*this;}
I Vec&operator*=(DB a){x*=a,y*=a;return*this;}
I Vec&operator/=(DB a){x/=a,y/=a;return*this;}
I DB operator*(Vec a){return x*a.x+y*a.y;}
I DB operator^(Vec a){return x*a.y-y*a.x;}
I friend bool operator<(Vec a,Vec b){
return (a^b)>0||((a^b)==0&&Len(a)<Len(b));
}
I friend DB Len(Vec a){
return sqrt(a.x*a.x+a.y*a.y);
}
I friend Vec Turn(Vec a,DB r){
const DB c=cos(r),s=sin(r);
return Vec(a.x*c-a.y*s,a.y*c+a.x*s);
}
};
初階
向量運算
點積：\(x_1x_2+y_1y_2\)，是一個向量在另一個向量上的投影
叉積：\(x_1y_2-x_2y_1\)，是兩個向量形成的平行四邊形的有向面積
用途很廣，搬一張圖
旋轉公式
\((x,y)\)轉\(r\)弧度\(\rightarrow(x\cos r-y\sin r,y\cos r+x\sin r)\)
點到直線距離
算一個叉積，除以底邊長就得到了高，也就是點到直線距離。
DB Dis(Vec&a,Vec&b1,Vec&b2){
return fabs((b1-a)^(b2-a))/Len(b1-b2);
}
直線交點&線段交點
點斜式要各種特判，不是理想的實現方法
考慮這樣的轉化方式：叉積->面積比->高的比->長度比->坐標
如下圖，兩直線相交於\(P\)，\(F\)和\(G\)是垂足，有\(\triangle A_1PF\sim\triangle CB_1G\)，可以推出
而\(P=A_1+\frac{|A_1P|}{|\vec a|}\vec{a}\)，然后就做出來了。
蒟蒻到現在好像還是背不了代碼，只能背下面這張圖了TAT
如何判斷線段相交呢？先把直線交點求出來，再點積判一下位置關系不就好啦
Vec LineCross(Vec a1,Vec a2,Vec b1,Vec b2){
a2-=a1;b2-=b1;
return b2^a2?a1+(b2^(b1-a1))/(b2^a2)*a2:Vec(NAN,NAN);
}
Vec SegCross(Vec&a1,Vec&a2,Vec&b1,Vec&b2){
Vec c=LineCross(a1,a2,b1,b2);
return (a1-c)*(a2-c)>0||(b1-c)*(b2-c)>0?Vec(NAN,NAN):c;
}
判斷點是否在多邊形內
過點任作一條射線（當然平行坐標軸的最方便啦），與多邊形交奇數次則在多邊形內。
凸包
極角排序
可以用atan2，也可以直接用叉積判斷。
把最下面的點的坐標找出來，把所有點整體平移，這時候所有點都不在\(x\)軸下面了，在\(180°\)內叉積可以起到比較方向的作用。
代碼見凸包部分。
求凸包
單調棧維護。如圖，\((P_i-P_{t-1})×(P_t-P_{t-1})\ge0\)，則\(P_t\)不在凸包上，需要彈掉。
洛谷P2742 【模板】二維凸包 / [USACO5.1]圈奶牛Fencing the Cows
I bool Polar(Vec&a,Vec&b){return (a^b)>0||((a^b)==0&&Len(a)<Len(b));}
//即Vec模板里重載的小於號。共線向量把短的放前面，求凸包的時候方便彈掉
I int Convex(Vec*a,Vec*e){
R n=e-a,k=0;
for(R i=1;i<n;++i)
if(a[i].y<a[k].y||(a[i].y==a[k].y&&a[i].x<a[k].x))k=i;
swap(a[0],a[k]);
const Vec tmp=a[0];
for(R i=0;i<n;++i)a[i]-=tmp;
sort(a+1,a+n,Polar);
R*st=new int[n],p=0;
for(R i=1;i<n;st[++p]=i++)
while(p&&((a[i]-a[st[p-1]])^(a[st[p]]-a[st[p-1]]))>=0)--p;
for(R i=0;i<=p;++i)a[i]=a[st[i]]+tmp;
return p+1;
}
int main(){
R n;cin>>n;
for(R i=1;i<=n;++i)cin>>a[i];
n=Convex(a+1,a+n+1);
DB ans=Len(a[1]-a[n]);
for(R i=1;i<n;++i)ans+=Len(a[i+1]-a[i]);
printf("%.2lf\n",ans);
return 0;
}
判斷點是否在凸包內
需要令凸包最下面的點為\((0,0)\)，不要還原。
把點的坐標也做相應的變換之后，丟進數組里lower_bound找到與它極角最接近的兩個點形成的線段，叉積判即可。
bool Inside(Vec*a,Vec*e,Vec v){
R i=lower_bound(a,e,v)-a-1;
return ((a[i+1]-a[i])^(v-a[i]))>=0;
}
多邊形面積
對於凸多邊形的理解：將多邊形划分成若干個三角形，每個三角形的面積是兩個向量叉積的\(\frac12\)。
update:其實任意多邊形面積都可以用三角剖分法求，詳見洛谷日報 計算幾何初步 by wjyyy。
DB PolygonArea(Vec*a,Vec*e){
R n=e-a;DB s=0;
for(R i=2;i<n;++i)s+=(a[i-1]-a[0])^(a[i]-a[0]);
return s/2;
}
進階算法
凸包
旋轉卡殼
利用凸包的單調性，在枚舉點的時候更新決策點（對踵點）並更新答案。
洛谷P1452 Beauty Contest
int main(){
R n;cin>>n;
for(R i=1;i<=n;++i)cin>>a[i];
n=Convex(a+1,a+n+1);
DB ans=0;
a[++n]=a[1];//最后加一個點
for(R i=1,p=1;i<n;++i){
while(Dis(a[p],a[i],a[i+1])<Dis(a[p+1],a[i],a[i+1]))p=p%n+1;
ans=max(ans,max(Len2(a[p]-a[i]),Len2(a[p]-a[i+1])));
}
printf("%.0lf\n",ans);
return 0;
}
半平面交
XZY巨佬提到了一種\(O(n^2)\)的算法。其實，如果不是在線插入的話，可以做到\(O(n\log n)\)。
我們用有向直線（一個點和一個方向向量）表示半平面，以下默認半平面在有向直線的左側。
對有向直線按方向向量的極角排序，維護一個雙端隊列，存儲當前構成半平面的直線以及相鄰兩直線的交點。
每次加入一條有向直線，如果隊首/隊尾的交點在直線右側（用叉積判）則彈掉隊首/隊尾的直線。
為什么這樣是對的呢？因為加入直線的單調性，所以要被彈出的直線一定在隊首或隊尾。感興趣的話可以自己手畫一些例子來理解。
需要注意的細節：
- 加入直線時，先彈隊尾，再彈隊首。
- 最后還要檢查隊尾交點是否在隊首直線的右側，如果是也要彈掉。
- 特判平行直線，在右側的要彈掉。
- 如果題目給出的半平面不一定有限制邊界，則應該手動加入一個INF邊界。
另有洛谷P3222 [HNOI2012]射箭 的題解
模板題：洛谷P4196 [CQOI2006]凸多邊形
struct Line{
Vec p,v;DB ang;
I Line(){}
I Line(Vec a,Vec b){p=a,v=b-a,ang=atan2(v.y,v.x);}
I bool operator<(Line&a){return ang<a.ang;}
I bool Right(Vec&a){return (v^(a-p))<-EPS;}
I friend Vec Cross(Line&a,Line&b){return a.p+(b.v^(b.p-a.p))/(b.v^a.v)*a.v;}
}a[N],q[N];
DB HalfPlane(Line*a,Line*e){
R n=e-a,h=0,t=0;
sort(a,e);q[0]=a[0];
for(R i=1;i<n;++i){
while(h<t&&a[i].Right(k[t-1]))--t;
while(h<t&&a[i].Right(k[h]))++h;
if(q[t].ang!=a[i].ang)q[++t]=a[i];
else if(a[i].Right(q[t].p))q[t]=a[i];
if(h<t)k[t-1]=Cross(q[t-1],q[t]);
}
while(h<t&&q[h].Right(k[t-1]))--t;
k[t]=Cross(q[t],q[h]);
return PolygonArea(k+h,k+t+1);
}
int main(){
R n,m,t=0;
cin>>n;
for(R i=1;i<=n;++i){
Vec fst,lst,now;
cin>>m>>fst;lst=fst;
for(R j=2;j<=m;++j)
cin>>now,a[++t]=Line(lst,now),lst=now;
a[++t]=Line(lst,fst);
}
printf("%.3lf\n",HalfPlane(a+1,a+t+1));
return 0;
}
閔可夫斯基和
對於歐氏空間的兩個點集\(A,B\)，其閔可夫斯基和為點集\(C=\{a+b|a\in A,b\in B\}\)
其中加法就是向量加法。
接下來我們只討論凸包的閔可夫斯基和。
一個四邊形，一個三角形，兩個凸包的閔可夫斯基和會長什么樣子呢？
因為凸包的特性，我們只需要取其中一個凸集所有最外層的點，將另一個凸多邊形沿這個點向量移動，就可以得到閔可夫斯基和。
觀察一下它的特點：外面正好有\(7\)個點，\(7\)條邊，有沒有注意到\(7=4+3\)呢？
實際上，對於任意兩個凸包來說，它們的閔可夫斯基和與它們的每一條邊按極角序順次相連所得的圖形都是全等的，也是一個凸包。
那豈不是很好求？二路歸並即可。需要處理三點共線的情況，了撇一點的話直接扔進Convex函數里再求一遍凸包就好了。
洛谷P4557 [JSOI2018]戰爭
typedef long long DB;//不涉及小數運算，直接開longlong
int Convex(Vec*a,Vec*e,Vec&bs){
R n=e-a,k=0,p=0;
for(R i=1;i<n;++i)
if(a[i].y<a[k].y||(a[i].y==a[k].y&&a[i].x<a[k].x))k=i;
swap(a[0],a[k]);bs+=a[0];
for(R i=n-1;~i;--i)a[i]-=a[0];
sort(a,e);
for(R i=1;i<n;st[++p]=i++)
while(p&&((a[i]-a[st[p-1]])^(a[st[p]]-a[st[p-1]]))>=0)--p;
for(R i=0;i<=p;++i)a[i]=a[st[i]];//做出來沒還原
return a[p+1]=0,p+1;
}
void Minkowski(R n,R m){
for(R i=n;i;--i)a[i]-=a[i-1];
for(R j=m;j;--j)b[j]-=b[j-1];
for(R i=1,j=1,k=0;i<=n||j<=m;++k)
c[k]=c[k-1]+(i<=n&&(j>m||a[i]<b[j])?a[i++]:b[j++]);
}
int main(){
ios::sync_with_stdio(0);
R n=in(),m=in(),q=in();
Vec bs=0,v;
for(R i=0;i<n;++i)cin>>a[i];
for(R i=0;i<m;++i)cin>>b[i],b[i]=-b[i];
n=Convex(a,a+n,bs);
m=Convex(b,b+m,bs);
Minkowski(n,m);
n=Convex(c,c+n+m+1,v);
while(q--)
cin>>v,cout<<Inside(c,c+n,v-bs)<<'\n';
return 0;
}
其它
最小圓覆蓋
把點集random_shuffle
第一層枚舉一個點，如果不在當前圓中，令其為圓心，半徑設為\(0\)。
第二層枚舉另一個點，如果不在當前圓中，將圓更新為以當前兩個點的連線為直徑的圓。
第三層枚舉另一個點，如果不在當前圓中，將圓更新為當前三個點的外接圓。
期望復雜度\(O(n)\)。
洛谷P1742 最小圓覆蓋
typedef long double DB;//此題略卡精度
I void PerpLine(Vec a1,Vec a2,Vec&b1,Vec&b2){//中垂線
b1=(a1+a2)/2;b2=b1+Turn(a1-a2);
}
I void CircleCoverage(Vec*a,Vec*e,Vec&O,DB&r){
R n=e-a;
for(R i=0;i<n;++i){
if(Len(a[i]-O)<=r)continue;
O=a[i],r=0;
for(R j=0;j<i;++j){
if(Len(a[j]-O)<=r)continue;
O=(a[i]+a[j])/2,r=Len(a[i]-a[j])/2;
Vec a2=O+Turn(a[i]-a[j]),b1,b2;
for(R k=0;k<j;++k){
if(Len(a[k]-O)<=r)continue;
PerpLine(a[i],a[k],b1,b2);
O=LineCross(O,a2,b1,b2),r=Len(a[i]-O);
}
}
}
}
int main(){
R n;cin>>n;
for(R i=1;i<=n;++i)cin>>a[i];
srand(20020307);
random_shuffle(a+1,a+n+1);
Vec O;DB r=0;
CircleCoverage(a+1,a+n+1,O,r);
cout<<fixed<<setprecision(10)<<r<<endl<<O<<endl;
return 0;
}
自適應辛普森積分
求函數的定積分，應用於算不規則圖形的面積。采用二次函數來擬合。
不斷分割區間並迭代，直到誤差控制在EPS以內。
因為誤差和區間長度有關，所以可以將EPS也進行迭代，每次除以\(2\)。
洛谷P4526 【模板】自適應辛普森法2
#include<bits/stdc++.h>
#define DB double
using namespace std;
DB a,EPS=1e-6;
inline DB f(DB x){
return pow(x,a/x-x);
}
inline DB Calc(DB l,DB r){
return (r-l)*(f(l)+f(r)+4*f((l+r)/2))/6;
}
DB Simpson(DB l,DB r,DB ans,DB EPS){
DB m=(l+r)/2,al=Calc(l,m),ar=Calc(m,r);
return fabs(al+ar-ans)<EPS?ans:Simpson(l,m,al,EPS/2)+Simpson(m,r,ar,EPS/2);
}
int main(){
cin>>a;
if(a<0)puts("orz");
else printf("%.5lf\n",Simpson(EPS,20,Calc(0,20),EPS));
return 0;
}