# 向量的板子

#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,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);
}


### 直线交点&线段交点

$\frac{\vec b×\vec c}{\vec b×\vec a}=\frac{S_{B_1B_2EA_1}}{S_{B_1B_2DC}}=\frac{|A_1F|}{|CG|}=\frac{|A_1P|}{|\vec a|}$

$$P=A_1+\frac{|A_1P|}{|\vec a|}\vec{a}$$，然后就做出来了。

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;
}


## 凸包

### 求凸包

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;
}


### 判断点是否在凸包内

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;
}


### 多边形面积

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;
}


# 进阶算法

## 凸包

### 旋转卡壳

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)$$

1. 加入直线时，先弹队尾，再弹队首。
2. 最后还要检查队尾交点是否在队首直线的右侧，如果是也要弹掉。
3. 特判平行直线，在右侧的要弹掉。
4. 如果题目给出的半平面不一定有限制边界，则应该手动加入一个INF边界。
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))<=0;}
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;
}


### 闵可夫斯基和

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;
}



## 其它

### 最小圆覆盖

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;
}


### 自适应辛普森积分

$\int_l^r f(x)dx\approx\int_l^r(ax^2+bx+c)dx\\ =\frac a3(r-l)^3+\frac b2(r-l)^2+c(r-l)\\ =\frac{r-l}6[2a(r^2-rl+l^2)+3b(r-l)+6c]\\ =\frac{r-l}6[f(l)+f(r)+4f(\frac{l+r}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;
}

posted @ 2019-01-12 09:23  Flash_Hu  阅读(1611)  评论(1编辑  收藏  举报