0. 前置芝士
0.1. 误差
由于浮点数绝对值越大精度越不精确,需要尽量避免 大数与小数的加减。
尽量少用三角函数,除法,开方,求幂,取对数运算。
对于判断浮点数的大小关系,有以下两种方案:
-
定义
sgn(const double x)
计算浮点数的符号:const double eps=1e-9; inline int sgn(const double x) { return x<-eps?-1:x>eps; }
利用这个函数就有:
inline int dcmp(const double x,const double y) { return sgn(x-y); }
-
乘上 (10) 的幂转化成整数。
0.2. 反三角函数
需要注意 定义域!
double x=1.000001;
if(fabs(x-1.0)<eps || fabs(x+1.0)<eps)
x=round(x);
double acx=acos(x);
对于反正切函数,最好使用 cmath::atan2(y,x)
。
1. 向量与线
1.1. 转角公式
将向量 ((x,y)=(rcos alpha,rsin alpha)) 逆时针旋转 (eta)。
1.2. 向量外积
(|a imes b|) 等于由向量 (a) 和向量 (b) 构成的平行四边形的面积。可以理解成一个行列式 (egin{vmatrix} a.x & b.x \ a.y & b.y end{vmatrix})。
根据 (a imes b) 的正负还可以判断 (a,b) 的位置关系:
- (a imes b=0)。矩阵行列式为 (0),那么可以相互线性表示。(a,b) 共线。
- (a imes b>0)。将 (b) 转到 (a) 的角度顺时针方向较小。
- (a imes b<0)。将 (b) 转到 (a) 的角度逆时针方向较小。
1.3. 点到直线与线段距离
![](https://img2020.cnblogs.com/blog/1889894/202108/1889894-20210806165002850-115713978.png)
用外积求出四边形面积,再除以 (p_1p_2)。
线段需要特判一下 (p_0) 对 (p_1p_2) 的垂足是否落在线段 (p_1p_2) 上。
1.4. 求直线交点
先利用跨立实验判断是否相交,再利用面积比计算。图是我嫖的。
![](https://img2020.cnblogs.com/blog/1889894/202108/1889894-20210806173228015-1031417436.png)
需要注意的是,在跨立实验中是否传入两点重合的直线或线段。
2. 多边形
2.1. 多边形面积
任意选点进行三角剖分,即加上多边形相邻点和选点 (O) 的叉积。感性理解一下当 (p_i) 在 (p_{i+1}) 的右手边时加入面积,反之减去面积。
![](https://img2020.cnblogs.com/blog/1889894/202108/1889894-20210807094926706-1246259930.png)
2.2. 射线法
用于判断点 (p) 是否在多边形内部。从点 (p) 引一条射线,若与多边形有奇数个交点则在内部,反之在外部。它还可以处理复杂多边形,就像这样:
![](https://img2020.cnblogs.com/blog/1889894/202108/1889894-20210807100801788-1727870135.jpg)
不过有引出射线与多边形端点有交点的情况。事实上像这样:
bool inAngle(vec a,vec b,const vec p) {
if(sgn(cross(a,b))<0) swap(a,b);
return sgn(cross(a,p))>=0 and sgn(cross(b,p))<=0;
}
bool InPolygon_1(const vec p,const vector <vec> poy) {
int n=poy.size();
double r=(rand()/double(RAND_MAX)-0.5)*2*pi;
vec v(cos(r),sin(r));
bool ret=0;
for(int i=0;i<n;++i) {
if(onSeg(poy[i],poy[(i+1)%n],p))
return 1;
ret=inAngle(poy[i]-p,poy[(i+1)%n]-p,v)?(!ret):ret;
}
return ret;
}
我们将交点两端的边都计算进去,所以不会有问题。
但是如果射线与多边形的边重合就会出错… 不过这种情况很难发生就是了,因为出题人并不知道我是怎么随机射线的…
2.3. 回转数算法
用于判断点 (d) 是否在多边形内部。沿多边形走一圈,累计绕点 (d) 转了多少角度。
- (0)。多边形外。
- (±pi)。多边形上。
- (±2pi)。多边形内。
在实现中并不需要计算转角,只用转象限。最后保证多边形外的点答案为 (0) 即可。考虑同样从 (p_i) 转到 (p_{i+1}),内部和外部的点有什么区别?(p_i-d) 与 (p_{i+1}-d) 的外积是相反的。利用这个性质就可以规定一个方向来加减象限了。
bool InPolygon_2(const vec p,const vector <vec> poy) {
int n=poy.size(),ret=0;
for(int i=0;i<n;++i) {
if(onSeg(poy[i],poy[(i+1)%n],p))
return 1;
double c=cross(poy[i]-p,poy[(i+1)%n]-p);
double d1=poy[i].y-p.y;
double d2=poy[(i+1)%n].y-p.y;
if(sgn(c)<0 and sgn(d1)<=0 and sgn(d2)>0) ++ret;
if(sgn(c)>0 and sgn(d2)<=0 and sgn(d1)>0) --ret;
}
return ret^0;
}
当然你还可以朴素地直接加象限,就像这样:
int getD(const vec v) {
if(v.x>=0 && v.y>=0) return 0;
if(v.x>=0 && v.y<0) return 3;
if(v.x<0 && v.y<0) return 2;
return 1;
}
bool InPolygon_3(const vec p,const vector <vec> poy) {
int n=poy.size(),ret=0;
for(int i=0;i<n;++i) {
double c=cross(poy[i]-p,poy[(i+1)%n]-p);
if(sgn(c)==0 and sgn(dot(poy[i]-p,poy[(i+1)%n]-p))<=0)
return 1;
int d1=getD(poy[i]-p),d2=getD(poy[(i+1)%n]-p);
if(d2==(d1+1)%4) ++ret;
else if(d2==(d1+3)%4) --ret;
else if(d2==(d1+2)%4) {
if(sgn(c)>=0) ret+=2;
else ret-=2;
}
}
return ret^0;
}
2.4. 动态凸包
#include <cstdio>
#define print(x,y) write(x),putchar(y)
template <class T>
inline T read(const T sample) {
T x=0; char s; bool f=0;
while((s=getchar())>'9' or s<'0')
f|=(s=='-');
while(s>='0' and s<='9')
x=(x<<1)+(x<<3)+(s^48),
s=getchar();
return f?-x:x;
}
template <class T>
inline void write(const T x) {
if(x<0) {
putchar('-'),write(-x);
return;
}
if(x>9) write(x/10);
putchar(x%10^48);
}
#include <map>
#define X first
#define Y second
using namespace std;
typedef map <int,int> Map;
typedef Map :: iterator It;
Map up,down;
long long cross(It o,It a,It b) {
return 1ll*(a->X-o->X)*(b->Y-o->Y)-
1ll*(a->Y-o->Y)*(b->X-o->X);
}
bool In(Map &t,int x,int y) {
if(t.empty()) return 0;
if(x<t.begin()->X or x>t.rbegin()->X)
return 0;
if(t.count(x)) return y>=t[x];
t[x]=y;
It it,i,j;
it=t.lower_bound(x);
i=j=it,--i,++j;
bool ret=cross(i,it,j)<=0;
t.erase(it);
return ret;
}
void add(Map &t,int x,int y) {
if(In(t,x,y)) return;
It it,i,j;
t[x]=y;
it=t.lower_bound(x);
for(i=it,--i,j=i,--j;it!=t.begin() and i!=t.begin();i=j--)
if(cross(it,i,j)>=0) t.erase(i);
else break;
for(i=it,++i,j=i,++j;j!=t.end() and i!=t.end();i=j++)
if(cross(it,i,j)<=0) t.erase(i);
else break;
}
int main() {
int op,x,y;
for(int T=read(9);T;--T) {
op=read(9),x=read(9);
y=read(9);
if(op==1)
add(up,x,-y),add(down,x,y);
else if(In(up,x,-y) and In(down,x,y))
puts("YES");
else puts("NO");
}
return 0;
}
2.5. 最小圆覆盖
戳这。
#include <bits/stdc++.h>
using namespace std;
template <class T>
inline T read(const T sample) {
T x=0; char s; bool f=0;
while((s=getchar())>'9' or s<'0')
f|=(s=='-');
while(s>='0' and s<='9')
x=(x<<1)+(x<<3)+(s^48),
s=getchar();
return f?-x:x;
}
const double eps=1e-2,pi=acos(-1.0);
inline int sgn(const double x) {
return x<-eps?-1:x>eps;
}
inline int dcmp(const double x,const double y) {
return sgn(x-y);
} // if x>y, then return 1
struct vec {
double x,y;
vec(const double X=0,const double Y=0):x(X),y(Y){}
vec operator + (const vec t) const {
return vec(x+t.x,y+t.y);
}
vec operator - (const vec t) const {
return vec(x-t.x,y-t.y);
}
vec operator * (const double t) const {
return vec(x*t,y*t);
}
vec operator / (const double t) const {
return vec(x/t,y/t);
}
bool operator < (const vec p) const {
int c=dcmp(x,p.x);
if(c) return c==-1;
return dcmp(y,p.y)==-1;
}
bool operator == (const vec p) const {
return !dcmp(x,p.x) && !dcmp(y,p.y);
}
friend double dot(const vec a,const vec b) {
return a.x*b.x+a.y*b.y;
}
friend double cross(const vec a,const vec b) {
return a.x*b.y-a.y*b.x;
}
double len() const {
return sqrt(dot(*this,*this));
}
vec rot(const double alpha) const {
return vec(x*cos(alpha)-y*sin(alpha),x*sin(alpha)+y*cos(alpha));
}
vec normal() const {
double l=len();
return vec(-y/l,x/l);
}
friend double angle(const vec a,const vec b) {
return acos(dot(a,b)/a.len()/b.len());
}
friend double dis(const vec a,const vec b) {
return sqrt(dot(a-b,a-b));
}
void Read() {scanf("%lf %lf",&x,&y);}
void Print() {printf("%.2f %.2f ",x,y);}
} O;
struct line {
vec p,v;
line(const vec P,const vec V):p(P),v(V){}
friend double dis(const line &l,const vec &p) {
vec v=p-l.p;
return fabs(cross(v,l.v))/(l.v).len();
}
friend bool onLine(const line &l,const vec &p) {
vec v=p-l.p;
return sgn(cross(v, l.v))==0;
}
};
int n;
double r;
vec p[(int)1e5+5],o;
vec GetPoint(const line &l1,const line &l2) {
vec c=l2.p-l1.p;
double t=cross(l2.v,c)/cross(l2.v,l1.v);
return l1.p+l1.v*t;
}
vec Circle(vec &a,vec &b,vec &c) {
return GetPoint(line((a+b)/2,(b-a).normal()),line((a+c)/2,(c-a).normal()));
}
int main() {
while(233) {
n=read(9);
if(!n) break;
for(int i=0;i<n;++i)
p[i].Read();
srand(time(NULL));
random_shuffle(p,p+n);
for(int i=0;i<n;++i)
if(dcmp(dot(p[i]-o,p[i]-o),r)>0) {
o=p[i],r=0;
for(int j=0;j<i;++j)
if(dcmp(dot(p[j]-o,p[j]-o),r)>0) {
o=(p[i]+p[j])/2,r=dot(p[j]-o,p[j]-o);
for(int k=0;k<j;++k)
if(dcmp(dot(p[k]-o,p[k]-o),r)>0)
o=Circle(p[i],p[j],p[k]),
r=dot(p[i]-o,p[i]-o);
}
}
o.Print(),printf("%.2f
",sqrt(r));
}
return 0;
}
3. 板子
#include <cmath>
#include <cstdio>
#include <vector>
#include <cstdlib>
#include <iostream>
#include <algorithm>
using namespace std;
template <class T>
inline T read(const T sample) {
T x=0; char s; bool f=0;
while((s=getchar())>'9' or s<'0')
f|=(s=='-');
while(s>='0' and s<='9')
x=(x<<1)+(x<<3)+(s^48),
s=getchar();
return f?-x:x;
}
const double eps=1e-9,pi=acos(-1.0);
inline int sgn(const double x) {
return x<-eps?-1:x>eps;
}
inline int dcmp(const double x,const double y) {
return sgn(x-y);
} // if x>y, then return 1
struct vec {
double x,y;
vec(const double X=0,const double Y=0):x(X),y(Y){}
vec operator + (const vec t) const {return vec(x+t.x,y+t.y);}
vec operator - (const vec t) const {return vec(x-t.x,y-t.y);}
vec operator * (const double t) const {return vec(x*t,y*t);}
vec operator / (const double t) const {return vec(x/t,y/t);}
friend double dot(const vec a,const vec b) {return a.x*b.x+a.y*b.y;}
friend double cross(const vec a,const vec b) {return a.x*b.y-a.y*b.x;}
double len() const {return sqrt(dot(*this,*this));}
bool operator < (const vec p) const {
int c=dcmp(x,p.x);
if(c) return c==-1;
return dcmp(y,p.y)==-1;
}
bool operator == (const vec p) const {
return !dcmp(x,p.x) && !dcmp(y,p.y);
}
vec rot(const double alpha) const {
return vec(x*cos(alpha)-y*sin(alpha),x*sin(alpha)+y*cos(alpha));
}
vec normal() const {
double l=len();
return vec(-y/l,x/l);
}
friend double angle(const vec a,const vec b) {
return acos(dot(a,b)/a.len()/b.len());
}
void Read() {scanf("%lf %lf",&x,&y);}
void Print() {printf("(%f,%f)
",x,y);}
} O;
struct line {
vec p,v;
line(const vec P,const vec V):p(P),v(V){}
friend double dis(const line l,const vec p) {
vec v=p-l.p;
return fabs(cross(v,l.v))/(l.v).len();
}
friend bool onLine(const line l,const vec p) {
vec v=p-l.p;
return sgn(cross(v, l.v))==0;
}
};
double dis(const vec l1,const vec l2,const vec p) {
double d1=dot(p-l1,l2-l1),d2=dot(p-l2,l1-l2);
if(sgn(d1)>=0 && sgn(d2)>=0)
return dis(line(l1,l2-l1),p);
if(sgn(d1)<0) return (p-l1).len();
return (p-l2).len();
}
bool onSeg(const vec l1,const vec l2,const vec p) {
bool flag=onLine(line(l1,l2-l1),p);
if(!flag) return 0;
double d=dot(l1-p,l2-p);
if(sgn(d)<=0) return 1;
return 0;
}
bool SegIntersection(const vec a1,const vec a2,const vec b1,const vec b2) {
double c1=cross(a2-a1,b1-a1),c2=cross(a2-a1,b2-a1);
double c3=cross(b2-b1,a1-b1),c4=cross(b2-b1,a2-b1);
int r1=sgn(c1)*sgn(c2),r2=sgn(c3)*sgn(c4);
if(r1>0 || r2>0) return 0;
if(r1==0 && r2==0) {
if(onSeg(a1,a2,b1) || onSeg(a1,a2,b2))
return 1;
return 0;
}
return 1;
}
vec GetPoint(const line l1,const line l2) {
vec c=l2.p-l1.p;
if(fabs(cross(l1.v,l2.v))<eps)
return vec(-1e9,-1e9);
double t=cross(l2.v,c)/cross(l1.v,c);
return l1.p+l1.v*t;
}
int getD(const vec v) {
if(v.x>=0 && v.y>=0) return 0;
if(v.x>=0 && v.y<0) return 3;
if(v.x<0 && v.y<0) return 2;
return 1;
}
double getArea(const vector <vec> poy) {
double ans=0;
int n=poy.size();
for(int i=0;i<n;++i)
ans+=cross(poy[i]-O,poy[(i+1)%n]-O);
return fabs(ans)/2;
}
bool inAngle(vec a,vec b,const vec p) {
if(sgn(cross(a,b))<0) swap(a,b);
return sgn(cross(a,p))>=0 and sgn(cross(b,p))<=0;
}
bool InPolygon_1(const vec p,const vector <vec> poy) {
int n=poy.size();
double r=(rand()/double(RAND_MAX)-0.5)*2*pi;
vec v(cos(r),sin(r));
bool ret=0;
for(int i=0;i<n;++i) {
if(onSeg(poy[i],poy[(i+1)%n],p))
return 1;
ret=inAngle(poy[i]-p,poy[(i+1)%n]-p,v)?(!ret):ret;
}
return ret;
}
bool InPolygon_2(const vec p,const vector <vec> poy) {
int n=poy.size(),ret=0;
for(int i=0;i<n;++i) {
if(onSeg(poy[i],poy[(i+1)%n],p))
return 1;
double c=cross(poy[i]-p,poy[(i+1)%n]-p);
double d1=poy[i].y-p.y;
double d2=poy[(i+1)%n].y-p.y;
if(sgn(c)<0 and sgn(d1)<=0 and sgn(d2)>0) ++ret;
if(sgn(c)>0 and sgn(d2)<=0 and sgn(d1)>0) --ret;
}
return ret^0;
}
bool InPolygon_3(const vec p,const vector <vec> poy) {
int n=poy.size(),ret=0;
for(int i=0;i<n;++i) {
double c=cross(poy[i]-p,poy[(i+1)%n]-p);
if(sgn(c)==0 and sgn(dot(poy[i]-p,poy[(i+1)%n]-p))<=0)
return 1;
int d1=getD(poy[i]-p),d2=getD(poy[(i+1)%n]-p);
if(d2==(d1+1)%4) ++ret;
else if(d2==(d1+3)%4) --ret;
else if(d2==(d1+2)%4) {
if(sgn(c)>=0) ret+=2;
else ret-=2;
}
}
return ret^0;
}
vec getCentre(const vector <vec> poy) {
double area=0; int n=poy.size();
vec ret(0,0);
for(int i=0;i<n;++i) {
double t=cross(poy[i]-O,poy[(i+1)%n]-O);
area+=t;
ret.x+=(poy[i].x+poy[(i+1)%n].x)*t;
ret.y+=(poy[i].y+poy[(i+1)%n].y)*t;
}
ret.x/=3,ret.x/=area;
ret.y/=3,ret.y/=area;
return ret;
}
vec cox[1000];
int tp;
void Andrew(const vector <vec> poy) {
int n=poy.size();
vec t[1000];
for(int i=0;i<n;++i) t[i]=poy[i];
sort(t,t+n);
tp=0;
for(int i=0;i<n;++i) {
while(tp>1 and sgn(cross(cox[tp-1]-cox[tp-2],t[i]-cox[tp-1]))<=0)
--tp;
cox[tp++]=t[i];
}
int lim=tp;
for(int i=n-2;i>=0;--i) {
while(tp>lim and sgn(cross(cox[tp-1]-cox[tp-2],t[i]-cox[tp-1]))<=0)
--tp;
cox[tp++]=t[i];
}
if(n>1) --tp;
}
// 计算圆的交的面积
double calc(const vec &o,const double r1,const int i) {
double d=sqrt(dot(o-a[i],o-a[i]));
if(dcmp(d,r1+r[i])>0) return 0;
if(dcmp(d,fabs(r1-r[i]))<0)
return min(r1,r[i])*min(r1,r[i])*pi;
double a1=acos((r1*r1+d*d-r[i]*r[i])/(r1*d*2));
double a2=acos((r[i]*r[i]+d*d-r1*r1)/(r[i]*d*2));
return a1*r1*r1+a2*r[i]*r[i]-sin(a1)*r1*d;
}
int main() {
return 0;
}
/*
point line_intersection(point a,point a0,point b,point b0)
{
double a1,b1,c1,a2,b2,c2;
a1 = a.y - a0.y;
b1 = a0.x - a.x;
c1 = cross(a,a0);
a2 = b.y - b0.y;
b2 = b0.x - b.x;
c2 = cross(b,b0);
double d = a1 * b2 - a2 * b1;
return point((b1 * c2 - b2 * c1) / d,(c1 * a2 - c2 * a1) / d);
}
*/
4. 一些题
例 1.
( ext{Segments}):给定 (nle 100) 条线段,判断是否存在一条直线使得所有线段在直线上的投影有交点。
问题可以转化为是否存在一条直线经过所有线段,这条直线与我们要求的直线垂直。一定存在这样的直线经过线段的端点,所以可以直接枚举。
例 2.
( ext{That Nice Euler Circuit})
题目有一些比较重要的性质:
- (p_i eq p_{i-1})。
- 欧拉机器绝对不会画出任何与其他已经画出的线重叠的线。然而,这两条线也可能相交。
欧拉定理:设 (G) 为任意的连通的平面图,则 (v-e+f=2),(v) 是 (G) 的顶点数,(e) 是 (G) 的边数,(f) 是 (G) 的面数。
首先枚举给出线段求出所有交点,不过交点会有重合,用 std::unique()
即可。枚举所有交点 (i),给出线段 (j),如果 (i) 在 (j) 上(不包含端点)就增加一条边。
为什么?题目保证没有重叠的线,所以交点不会划分已划分的线段。
例 3.
( ext{CodeForces - 1C Ancient Berland Circus})
首先多边形一定在三角形的外接圆上,否则三个端点不可能和多边形端点重合。其次应使多边形端点尽量少,从而最小化面积。
利用正弦定理可知 (r=frac{abc}{4S})。三角形每条边对应的圆心角度数是 (arccos (frac{2r^2-l^2}{2r^2}))。
将圆心角度数取 (gcd) 可以最大化中心角,从而最小化多边形顶点数。第三个圆心角可以直接取 (2pi-alpha-eta) 因为它要么是 (alpha+eta) 要么是 (2pi-alpha-eta),
例 4.
( ext{UVA - 10256 The Great Divide})
判断两个凸包是否相交。
- 是否有相交线段。特判凸包退化成线段的情况。
- 点是否在另一凸包内。这是为了判断凸包相互包含的情况。