Gay算几何

七月 09, 2021

讲课

前言

抱着负责的心态把要讲的东西重写了一遍。

[TOC]

知识

基础

在计算几何题中，各种图形的表示往往是以坐标形式出现的，因此我们需要对这些坐标进行一个好的存储，从而方便我们进行操作。

判断精度

我们有时会遇见精度问题，于是整个好玩的判一下就行了。

const double eps = 1e-8;
int ceps(double x)//处理精度，判断正负 
{
	if(fabs(x) < eps) return 0;
	if(x < 0) return -1;
	return 1;
}

表示

点和向量

这个十分简单，一个点只有 $x$ 和 $y$ 坐标，甚至向量也是这样，我们整个结构体，再整个 typedef 就 ok。

另一种好的方式是 pair。

struct poi//点或向量
{
	double x,y;
	poi(double X = 0, double Y = 0){x = X, y = Y;}
};
typedef poi vec
pair<double,double> poi;

两点间的距离

任意两点之间 $P_1(x_1,y_1),P_2(x_2,y_2)$ 的距离记为 $|P_1P_2| $，由勾股定理可知：

$|P_1P_2| = \sqrt{(x_2-x_1)^2 + (y_2-y_1)^2}$ .

double dis(poi a,poi b)//两点距离
{
	return (double) sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y))
}

直线

两点确定一条直线，所以我们存两个点。

另一种表示，可以是一个点加上一个方向向量。

但是这里只写点表示法了。

struct line
{
	poi A,B;
	line(poi a,poi b){A = a, B = b;}
}

线段

同样道理，左右端点就ok。

struct Segment
{
	poi A,B;
	Segment(poi a,poi b) {A = a, B = b;} 
}

多边形

考虑一下，emm，多边形不就是一堆点嘛。。。

poi a[N];

圆

思考，坐标系上的一个圆…

很显然会有圆心和半径，那么存个圆心存个半径就行了。

struct circle
{
	poi p;//圆心 
	double r;
	circle(point A,double B){p = A, r = B;}
}

至此，表示方法就结束了。

等等，是不是有另一种表示方法？

有的！

極角坐标系

选一个 $O$ 称之为极点，然后引出一条射线叫极轴，再选一个单位长度（数学中常为 $1$ ），一个角度（通常为弧度），以及其正方向（通常为逆时针），这样一个极角坐标系就出来了。

~~IG要无了，我shy哥还没上场~~

我们先在这里放一张图：

在极坐标系中，设 $A$ 为平面上一点，极点 $O$ 和 $A$ 之间的距离 $|OA|$ 为极径，记为 $\rho$ ,以极轴为始边， $OA$ 为终边的角 $\angle xOA$ 为极角，记为 $\beta$ ,那么有序序对 $(\rho,\beta)$ 即为 $A$ 的极坐标。

有些时候我们需要把极坐标系转化到平面直角坐标系下去研究，两者转化关系如下：

点 $A(\rho,\theta)$ 的直角坐标 $(x,y)$ 可以表示为： $\left\{\begin{aligned} \rho \cos\theta&=x\\\rho\sin\theta&=y\end{aligned}\right.$

进而可知： $\left\{\begin{aligned} \rho^2 &=x^2 + y^2\\\tan\theta&={y\over x}\end{aligned}\right.$

于是极角 $\theta = \arctan{y\over x} $ ,然后这样就可以求出极角了。

如果题目中要求反正切函数，尽量使用 $atan2({y,x})$ ,这个函数用途比 $atan(x)$ 要广。

~~搬了搬了~~

三角和差公式

$\sin(a + \theta) = \sin a\cos \theta + \cos a\sin \theta$

$\sin(a-\theta) = \sin a\cos \theta - \cos a\sin \theta$

$\cos(a + \theta) = \cos a\cos \theta - \sin a\sin \theta$

$\cos(a - \theta) = \cos a \cos \theta + \sin a\sin \theta$

矢量及其运算

矢量就是向量…

模长

若向量 $\vec a = (x,y)$ , 则 $|a| = \sqrt{x^2+y^2} = \sqrt{|\vec a|^2}$

double len(poi a)//向量的模长 
{
	return (double) sqrt(a.x*a.x+a.y*a.y);
}

加减法

对于 $\vec a = (x_1,y_1)$ , $\vec b = (x_2,y_2)$ , 则 $\vec a + \vec b = (x_1+x_2,y_1+y_2)$

对于 $\vec a = (x_1,y_1)$ , $\vec b = (x_2,y_2)$ , 则 $\vec a -\vec b = (x_1-x_2,y_1-y_2)$

1 2	poi operator + (poi a,poi b) {return poi(a.x+b.x,a.y+b.y);} poi operator - (poi a,poi b) {return poi(a.x-b.x,a.y-b.y);}

数乘

对于 $\vec a = (x,y) $ , $\lambda\vec a = (\lambda x,\lambda y)$

1	poi operator * (poi a,double k) {return poi(a.xk,a.yk);}//数乘

数量积

数量积也就是点积，结果是一个数，大小等于这两个矢量的模长的乘积，再乘以他们夹角的余弦，即： $a \cdot b = |a||b|\cos<a,b>$ 。

对于 $\vec a = (x_1,y_1)$ , $\vec b = (x_2,y_2)$ , $\vec a \cdot \vec b = (x_1x_2+y_1y_2)$

夹角 $\theta$ 与点积大小的关系 (可以结合三角函数理解一下)：

若 $\theta = 0$ , 则 $\vec a\cdot \vec b = |a||b|$
若 $\theta = 180^\circ$ ,则 $\vec a \cdot \vec b = -|a||b|$
若 $\theta < 90^\circ$ , 则 $\vec a\cdot \vec b > 0$
若 $\theta = 90^\circ$ , 则 $\vec a \cdot \vec b = 0$
若 $\theta > 90^\circ$ , 则 $\vec a\cdot\vec b < 0$

1	double Dot(poi a,poi b){return a.xb.x+a.yb.y;}//点积

满足以下几个性质：

交换律： $\vec a \cdot \vec b = \vec b\cdot \vec a$

分配律： $\vec a\cdot (\vec b + \vec c) = \vec a\cdot \vec b + \vec a \cdot \vec c$

结合律： $(m\vec a) \cdot \vec b = m (\vec a\cdot \vec b) = \vec a \cdot (m\vec b)$

矢量积

也叫叉积，其结果是一个向量。

矢量 $a$ 和矢量 $b$ 叉乘的几何意义为：他的模长等于由 $a$ 和 $b$ 组成的平行四边形的面积，方向与四边形所在平面垂直。即： $|\vec a \times \vec b| = |\vec a||\vec b| \sin<a,b>$ 。

对于 $\vec a = (x_1,y_1 )$ ， $\vec b = (x_2,y_2)$ , $\vec a\times \vec b = (x_1y_2-x_2y_1)$ (交叉相乘)

如图所示，向量位置与叉积模长大小的关系：

若 $\vec a\times \vec b = 0$ ，则 $\vec a$ 和 $\vec b$ 共线，方向相同或相反。
若 $\vec a\times \vec b > 0$ , 则 $\vec a$ 在 $\vec b$ 的顺时针方向。
若 $\vec a\times \vec b < 0$ , 则 $\vec a$ 在 $\vec b$ 的逆时针方向.

1	double Cro(poi a,poi b){return a.xb.y-a.yb.x;}//叉积的模长

注意 $\sin$ 函数的符号，如果 $\vec a$ 在 $\vec b$ 的顺时针方向，那么 $\sin<\vec a,\vec b>$ 就是正值，否则为负值，记住一句话：顺正逆负即可。

又因为前面乘的是绝对值，所以叉积的符号就是 $\sin$ 函数的符号。

小技巧

三角形面积公式： $2S_{\Delta ABC} = ab\sin C$

因此，三角形的面积就是任意两条临边的叉积绝对值的一半。

同理，平行四边行的面积就是两条临边对应向量的叉积的绝对值。

类似的，任意四边形的面积等于两条对角线的叉积的绝对值的一半。

向量的旋转

$v^\prime = (x\cos a - y\sin a, x\sin a + y\cos a)$ .

其中 $a$ 表示你逆时针旋转的弧度。这个公式比较容易忘，建议熟记。

如果实在记不住，就把直角坐标转化为极坐标，在把角度加上 $a$ 之后再转回来（其实这个公式就是这么来的）。

poi rtt(poi a,double altha)
{
	return poi(a.x * cos(altha) - a.y * sin(altha),a.x * sin(altha) + a.y *cos(altha));
}

例题： uva11178

图形之间的关系

点与线段

判断点是否在线段上

做法1：点 $p$ 在线段 $AB$ 上，当且仅当 $|PA| + |PB| = |AB|$

做法2：先利用叉积判断点是否在线段所在的直线上（叉积为 0），在用点积判断点是否在线段上（点积小于等于 0）即可。如果点在线段的端点上点积等于零，否则点积小于零。

一般来说，做法 $1$ 损失的精度比较大，所以我们一般采用做法 2.

bool pd_PS(poi p,poi A,poi B)//判断 p 是否在 线段 AB 上 
{
//	return dis(p,A) + dis(p,B) == dis(A,B);//这种做法精度损失较大。
	return (ceps(Cro(p-A,p-B)) == 0) &&(ceps(Dot(p-A,p-B) <= 0));
}

口胡：如果说 $p$ 在线段 $A B$ 的外部，显然 $\vec {AP}$ 和 $\vec {BP}$ 是同向的，点积为正，反之，逆向为负。

点到线段的距离

从 $P$ 向 $AB$ 直线，做一条垂线，设垂足为 $Q$ 。

若 $Q$ 在线段 $AB$ 上，则点 $P$ 到线段 $AB$ 的距离为 $PQ$ 。

如果 $Q$ 在射线 $BA$ 上，则所求距离为 $PA$ 的长度，否则为 $PB$ 的长度。

如果 $\angle PAB$ 和 $\angle PBA$ 均为锐角或直角，则 $Q$ 在线段 $AB$ 上，若 $\angle PAB$ 为钝角，则 $Q$ 在射线 $BA$ ,若 $\angle PBA$ 为钝角，则 $Q$ 在射线 $AB$ 上，利用点积可以很好的判断出来。

可以自己画一下图理解一下。

double dis_PS(point p,point A,point B)//点到直线的距离
{
	if(Dot(p-A,B-A) < 0) return dis(p,A);//垂足离A更近
	if(Dot(p-B,A-B) < 0) return dis(p,B);//垂足离B更近
	return abs(Cro(A-p,B-p))/dis(A,B);//垂足的距离最小
}

例题： UVA10263 Railway

判断两线段是否相交(不包括交点)

如果两个线段相交，则必然有一条线段的两个端点在另一条线段的两侧。

利用叉积判断即可。

bool pd_SS(point A,point B,point C,point D)//线段 AB 和线段 CD 是否相交 
{
	double f1 = Cro(C-A,B-A), f2 = Cro(D-A,B-A);//第一条线段的两个端点是否在第二条线段的两侧 
	double g1 = Cro(A-C,D-C), g2 = Cro(B-C,D-C);//第二条线段的两个端点是否在第一条线段的两侧 
	return ((f1 < 0) ^ (f2 < 0)) && ((g1 < 0) ^ (g2 < 0));
}

点与直线

点与直线的距离

如图所示我们要求点 $A$ 到直线 $BC$ 的距离：

我们连接 $AB$ 和 $AC$ ,在做 $A$ 关于直线 $BC$ 的垂线，设垂足为 $D$

则有 $2S_{\Delta ABC} = AD * BC$ .

点 $A$ 到直线的距离 $AD$ 则为 $2S_{\Delta ABC} \over BC$ ,用一个柿子表示出来就是： $AD = {|\vec {AB} \times \vec {AC}| \over |\vec{BC}|}$

double dis_PL(point p,point A,point B)//点到直线 AB 的距离 
{
	return abs(Cross(p-A,p-B)) / len(A-B);//面积除以底 
}

点在直线上的垂足（投影）

直线 $AB$ 的参数式: $A + tv$ , 点 $P$ 的垂足 $Q$ 的参数为 $t_0$ , 根据 $\vec {QP} \perp \vec {AB}$ ，则有：

Q = A + t_0 v\\ (P-Q)\cdot v = 0\\ (P - A - t_0v) \cdot v = 0\\ (P-A) \cdot v - t_0v \cdot v = 0\\ t_0 = {(P-A) \cdot v\over v \cdot v}\\

point GetLineRoot(point p,point A,point B)//p在直线AB上的投影 
{
	point v = B - a;//方向向量v 
	return A + v * (Dot(P-A,v) / Dot(v,v));
}

点在直线的那一边

假设直线上一点 $P$ 方向向量为 $v$ ,则 $Q$ 在直线的那一边可以用 $\vec {PQ} \times v$ 来判断。叉积为负， $Q$ 在直线上方，

叉积为零，在直线上，叉积为正， $Q$ 在直线下方。

点关于直线的对称点

先找到点 $P$ 关于直线的垂足 $Q$ , 然后把 $PQ$ 延长两倍即可。

point Symmetry_PL(point p,point A,point B)//p关于直线AB的对称点
{
	return p + (GetLineRoot(p,A,B)-p) * 2;//把Q延长两倍	
}

（这里原来有一道巨ex的题）

线与线

两直线交点

设两直线方程分别为： $l_1 = P + t\vec v$ , $l_2 = Q + t\vec w$ 。

则：交点在第一条直线的参数为 $ t $ ,则 $t = {(P-Q) \times \vec w \over \vec w \times \vec v}$ 。

因为是交点，所以同时满足两个直线的方程，令交点在两条直线上的参数分别为 $t$ 和 $\mu$ ,

则有： $P + t \vec v = Q + \mu \vec w$ 。

两边同时叉上 $\vec w$ 可得： $(P + t\vec v)\times \vec w = (Q + \mu\vec w)\times \vec w$ 。

去括号可得： $P\times \vec w + t\vec v \times\vec w = Q\times \vec w + \mu \vec w \times \vec w$ 。

又因为 $\vec w \times \vec w = 0$ , 所以可得： $P \times \vec w + t\vec v \times \vec w = Q \times \vec w$

移项可得： $t\vec v \times \vec w = (Q-P)\times \vec w$

所以： $t = {(P-Q)\times \vec w\over \vec w \times \vec v}$ 。

point jiao_LL(point A,point B,point C,point D)//线段AB和线段CD的交点 
{
	point v = B - A, w = D - C, PQ = A - C;//直接套用公式就ok了 
	double t = Cro(PQ,w) / Cro(w,v);
	return A + t * v;
}//用的时候记得要特判一下两直线是否平行

点与多边形

判断点是否在多边形内

设凸 $n$ 边形的顶点为 $P_0P_1....P_n$ ，以及点 $A$ 。记 $\theta_i = <\vec{AP_i},\vec{AP_{i+1}}>$ ，则点 $A$ 在多边形内部等价于 $\displaystyle\sum_{i=0}^{n-1} \theta_i = 2π$

需要特判一下在多边形顶点上的情况。

point p[1005];
int n;
bool InPloygon(point A)//精度损失会比较大一些 
{
	double alpha = 0;
	for(int i = 1; i <= n; i++)
	{
		if(i == n) alpha += fabs(Angle(p[i]-A,p[1]-A));
		else alpha += fabs(Angle(p[i]-A,p[i+1]-A));
	}
	return dcmp(alpha-2 * pai) == 0;
}

例题： UVA109 SCUD Busters

判断点是否在凸多边形内部

给你 $m$ 个点和由 $n$ 个点组成的凸多边形，问你这 $m$ 个点是否在凸多边形内部。

~~转角法和射线法？我也不知道，全家桶里的东西被偷吃了~~

如果用上面的转角法和射线法来判断的话，复杂度为 $O(nm)$ 的。

有没有一个更高效的算法呢？我们可以用二分来做到 $O(mlogn)$ 的。

首先任选一个点为起点，向凸多边形其他顶点做一条射线。如图：

显然，会把整个凸多边形分成 $n-2$ 个三角形。

第二步：判断这个点是否在所有射线的包含范围之内，即判断定点是否在最左侧射线的左边或者最右侧射线的右边。

第三部：二分出定点在那个三角形内部，然后判断定点是否在这个三角形内部即可。

int InPolygn(point p,point *a)
{
	if(Cro(p-a[1],a[2]-a[1]) < 0 || Cro(p-a[1],a[n]-a[1]) > 0) return 0;//在边界外 
	if(pd_PS(p,a[1],a[2]) || pd_PS(p,a[n],a[1])) return 2;//在边界上 
	int l = 2, r = n;
	while(r - l > 1)
	{
		int mid = (l + r)>>1;
		if(Cro(p-a[1],a[mid]-a[1]) < 0) l = mid;
		else r = mid;
	}
	if(Cro(p-a[l],a[r]-a[l]) > 0) return 0;//在边界外 
	if(pd_PS(p,a[l],a[r])) return 2;//在边界上 
	return 1;
}

圆

判断点是否在圆内

很简单，判断点到圆心的距离和圆的半径大小即可。

过定点做圆的切线

我们先求出 $AB$ 的长度和 $AC$ 的长度，然后可以根据这个求出来 $\angle CBA$ .

那么两条切线就是直线 $AB$ 顺时针（逆时针）旋转 $\angle CBA$ 形成的直线。

int get_PC(point p,Circle c,line *v)//返回交点个数，v存储切线 
{
	Vector u = c.p - p;//BA向量 
	double dis = dis(p,c.p);
	if(dis < r) return 0;
	else if(dcmp(dis,r) == 0)//点在圆上的情况 
	{
		v[0] = retate(u,PI/2);
		return 1;
	}
	else//点不在圆上的情况 
	{
		double ang = asin(c.r/dis);
		v[0] = retate(u,ang);
		v[1] = retate(u,-ang);
		return 2;
	}
}

两个圆的交点

point get(Circle c,double a)//返回c圆上圆心角为a的点的坐标 
{
	return point(c.x+cos(a)*c.r,c.y+sin(a)*r);
}
double Angle(point a){return atan2(a.x,a.y);}
int getCC(Cirle c1,Cirle c2,vector<point>& sta)//返回交点个数，sta存交点的坐标 
{
	double d = len(c1.c-c2.c);
	if(cdmp(d) == 0)
	{
		if(dcmp(c1.r-c2.r) == 0) return -1//两个圆重合 
		return 0;//没有交点 
	}
	if(dcmp(c1.r+c2.r-d) < 0) return 0;//没有交点的情况 
	if(dcmp(fabs(c1.r-c2.r)-d) > 0) return 0;
	double alpha = Angle(c2.c-c1.c);
	double deta = acos((c1.r*c1.r+d*d-c2.r*c2.r)/(2*c1.r*d));//余弦定理算夹角 
	point p1 = get(c1,alpha-deta);
	point p2 = get(c1,alpha+deta);//旋转 
	sta.push_back(p1);
	if(p1 == p2) return 1;
	sta.push_back(p2);
	return 2;
}

三角形重心

三点各坐标的平均值： $({x_1+x_2+x_3\over 3}, {y_1+y_2+y_3\over 3})$

图形面积

任意多边形面积

因为从第一个顶点出发可以把多边形分为 $n-2$ 个三角形，多边形的面积就是这 $n-2$ 个三角形的面积之和，直接叉积计算即可。

因为你叉积运算的时候是带符号的，所以你计算结果可能是面积也可能是面积的相反数，取个绝对值就行了。

在你计算复杂多边形的面积的时候，多算的面积就会被叉积的符号抵消掉。

double PloyArea(point *a,int n)
{
	double res = 0;
	for(int i = 2; i < n; i++) res += Cro(a[i]-a[1],a[i+1]-a[1]);
	return abs(res / 2);//如果你提前按极角排一下序，就可以省去取绝对值这一步
}

P1183 多边形的面积

自适应辛普森法

用来求定积分。

当函数极为鬼畜的时候，我们采用二次拟合的形式，当精度相等的时候，便视为相等。

假如说我们要求 $\int_{a}^{b} f(x) d(x)$ ,你现在要求出这个函数的近似值。

我们可以把这个函数近似成矩形和梯形和等等。我们为了保证准确性，这里选择二次函数来拟合曲线。

设 $f(x) \approx Ax^2 + Bx+ C$

$\int_{a}^{b} f(x)d(x) \approx \int_{a}^{b} Ax^2+Bx+C$

$= \int_{a}^{b} Ax^2 + \int_{a}^{b} Bx + \int_{a}^{b} C$

$={A\over 3}x^3 + {B\over 2}x^2 + Cx|_{a}^{b}$

$={(a-b)\over 6} (2A(a^2+ab+b^2) + 3B(a+b) + 6C)$

$={a-b\over 6} ((Aa^2+Ba+c) + (Ab^2+Bb+c) + 4(({a+b\over 2})^2 + {a+b\over 2} + c))$

$ = {a-b\over 6} (f(a) + f(b) + 4f({a+b\over 2}))$

然后就有个结论： $$\int_a^{bf(x)dx\approx\int_a}{b}(Ax^2+Bx+C)dx=\frac{(a-b)(f(a)+f(b)+4f(\frac{a+b}{2}))}{6}$$

这样直接算的话，我们的精度误差会很大，我们考虑吧 $[a,b]$ 从终点分为两部分 $[a,mid],[mid,b]$ 来计算。

当 $\int_{a}^{mid} f(x)d(x) + f_{mid}^{b}f(x)d(x) = f_{a}^{b}f(x)d(x) + eps$ 时，便停止递归。

double f(double x){return 1.0*(c*x+d)/(a*x+b);}
double calc(double l,double r)
{
	double mid = (l +r)/2.0;
	return (r-l)*(f(l)+f(r)+4*f(mid))/6;
}
double simpon(double l,double r,double last)
{
	double mid = (l + r)/2; 
	double lx = calc(l,mid), rx = calc(mid,r);//递归算左右侧的积分
	if(fabs(lx + rx - last) < eps) return lx+rx;//lx+rx-last<eps 说明我们l-r这一段的积分的值误差小到可以忽略
	else return simpon(l,mid,lx) + simpon(mid,r,rx);//这个你可以理解为simpon把l-r这一段划分成了若干个小段，l-r这一段的积分等于这些小段积分的和。
}
int main()
{
	scanf("%lf%lf%lf%lf%lf%lf",&a,&b,&c,&d,&l,&r);
	printf("%.6lf\n",simpon(l,r,calc(l,r)));
	return 0;
}

例题：板子1 板子2

其实这个算法挺冷门的，学了这个关键是做这道题（~~柠檬树上懵逼果，柠檬树下你和我~~）

面积并

全家桶里也没有，我稍微看了看，大概就是一条扫描线进行扫描，然后

离散化: 这些技巧都是老生常谈的了, 不然浮点数怎么建树, 离散化坐标就可以了
扫描线: 首先把矩形按轴分成两条边, 上边和下边, 对轴建树, 扫描线可以看成一根平行于轴的直线.
从开始往上扫, 下边表示要计算面积, 上边表示已经扫过了, 直到扫到最后一条平行于轴的边
但是真正在做的时候, 不需要完全模拟这个过程, 一条一条边地插入线段树就好了
线段树: 用于动态维护扫描线在往上走时, 轴哪些区域是有合法面积的
这种线段树是不用的, 因为不用, 为啥不用, 因为没有查询操作

不太懂www

不如我们再看看图？

没咋看的代码（不如先放放吧）：


using namespace std;
const int N = 205, INF = 0x3f3f3f3f, MOD = 1e9 + 7;
 
int n;
struct Seg {
    double l, r, h; int d;
    Seg() {}
    Seg(double l, double r, double h, int d): l(l), r(r), h(h), d(d) {}
    bool operator< (const Seg& rhs) const {return h < rhs.h;}
} a[N];
 
int cnt[N << 2]; //根节点维护的是[l, r+1]的区间
double sum[N << 2], all[N];
 
#define lson l, m, rt << 1
#define rson m + 1, r, rt << 1 | 1
 
void push_up(int l, int r, int rt) {
    if(cnt[rt]) sum[rt] = all[r + 1] - all[l];
    else if(l == r) sum[rt] = 0; //leaves have no sons
    else sum[rt] = sum[rt << 1] + sum[rt << 1 | 1];
}
 
void update(int L, int R, int v, int l, int r, int rt) {
    if(L <= l && r <= R) {
        cnt[rt] += v;
        push_up(l, r, rt);
        return;
    }
    int m = l + r >> 1;
    if(L <= m) update(L, R, v, lson);
    if(R > m) update(L, R, v, rson);
    push_up(l, r, rt);
}
 
int main() {
    int kase = 0;
    while(scanf("%d", &n) == 1 && n) {
        for(int i = 1; i <= n; ++i) {
            double x1, y1, x2, y2;
            scanf("%lf%lf%lf%lf", &x1, &y1, &x2, &y2);
            a[i] = Seg(x1, x2, y1, 1);
            a[i + n] = Seg(x1, x2, y2, -1);
            all[i] = x1; all[i + n] = x2;
        }
        n <<= 1;
        sort(a + 1, a + 1 + n);
        sort(all + 1, all + 1 + n);
        int m = unique(all + 1, all + 1 + n) - all - 1;
 
        memset(cnt, 0, sizeof cnt);
        memset(sum, 0, sizeof sum);
 
        double ans = 0;
        for(int i = 1; i < n; ++i) {
            int l = lower_bound(all + 1, all + 1 + m, a[i].l) - all;
            int r = lower_bound(all + 1, all + 1 + m, a[i].r) - all;
            if(l < r) update(l, r - 1, a[i].d, 1, m, 1);
            ans += sum[1] * (a[i + 1].h - a[i].h);
        }
        printf("Test case #%d\nTotal explored area: %.2f\n\n", ++kase, ans);
    }
    return 0;
}

基本算法

凸包

现在平面上有一堆点，让你找到一个点集，依次连接这些点，组成一个多边形，且所有的点都在这个多边形内部或边上，这个多边形我们称之为凸包。

就相当于平面上有很多钉子，你拿一个橡皮筋一围，围出来的多边形就是凸包。

大概张这个样子：

然后观察一下这个凸包，你会发现一下很好玩的性质。

上凸壳中的点的斜率是单调递减的。
下凸壳中的点的斜率的单点递增的。

这个性质有什么用呢？他主要是用来斜率优化的。

Graham 扫描算法

首先我们发现最下面的一个点肯定会在凸包上。所以一开始我们可以通过一遍扫描把这一个点找出来记做点 $A$ (如果有多个从坐标相同的点都在最下方，则取最左边的点)。

第二步：把所有的点按 $Ap_i$ 的极角大小进行排序，极角相同的按照到点 $A$ 的距离来排序。

第三步：维护一个栈，来保存当前凸包的形态。将每个点按顺序依次加入凸包中。如果当前这个点 $P$ 和栈顶的两个顶点 $B,C$ 。如果 $\vec{BC}$ 在 $\vec{PC}$ 的逆时针方向，显然 $B$ 点肯定不在凸包上，然后每加入一个点，都重复上述步骤，排除不符合的点即可。

顺时针的情况为什么不考虑呢？因为可能是在拐弯处。

口胡说的不太清楚，拿一些图演示一下。

（1）首先给你一堆点，最下面的点 $p_1$ 肯定是在凸包里面的，然后把所有的点按极角排序，编号。

（2）然后 $P_2$ 准备入栈，由于栈中元素过少，所以符合条件，可以直接进栈。

（3）之后因为 $\vec {P_2P_1}$ 在 $\vec{P_3P_1}$ 的顺时针方向，符合条件，暂时把他加进去

（4）之后 $\vec {P_3P_2}$ 在 $\vec{P_4P_2}$ 的逆时针方向，就把 $P_3$ 出栈，再去检查 $P_2$ 符合条件，于是 $P_4$ 入栈。

（5） $P_5$ 一切正常，入栈。

（6） $P_6$ 这里就麻烦多了，首先 $\vec{P_5P_4}$ 在 $\vec{P_6P_4}$ 的逆时针方向，把 $P_5$ 出栈

（7）又发现 $\vec{P_4P_2}$ 在 $\vec{P_6P_2} $ 的逆时针方向，把 $P_4$ 出栈。

（8） $P_7$ 一切正常，入栈

（9） $P_8$ 一切正常，入栈

(10) 最后就连到了起点 $P_1$ ，我们的凸包就求出来了。

#include<iostream>
#include<cstdio> 
#include<algorithm>
#include<cmath>
using namespace std;
const double eps = 1e-8;
int n,top;
double ans;
struct point
{
	double x,y;
	point(double X = 0, double Y = 0){x = X, y = Y;}
}p[100010],sta[100010];
point operator + (point a,point b){return point(a.x+b.x,a.y+b.y);}
point operator - (point a,point b){return point(a.x-b.x,a.y-b.y);}
point operator * (point a,double b){return point(a.x*b,a.y*b);}
double Dot(point a,point b){return a.x*b.x+a.y*b.y;}
double Cro(point a,point b){return a.x*b.y-a.y*b.x;}
int dcmp(double x)
{
	if(fabs(x) < eps) return 0;
	if(x > 0) return 1;
	return -1; 
}
double dis(point a,point b)
{
	return sqrt((a.x-b.x)*(a.x-b.x) + (a.y-b.y) * (a.y-b.y));
}
bool comp(point a,point b)
{
	double m = Cro(a-p[1],b-p[1]);//由叉积来判断极角大小，如果a的极角比b小，那么ap一定在bp的顺时针方向
	if(dcmp(m) > 0) return 1;
	if(dcmp(m) == 0 && dcmp(dis(a,p[1])-dis(b,p[1])) < 0) return 1;
	return 0;
}
int main()
{
	scanf("%d",&n);
	for(int i = 1; i <= n; i++)
	{
		scanf("%lf%lf",&p[i].x,&p[i].y);
		if(dcmp(p[i].y-p[1].y) < 0) swap(p[1],p[i]);
		else if(dcmp(p[i].y-p[1].y) == 0 && dcmp(p[i].x-p[1].x) < 0) swap(p[1],p[i]);
	}
	sort(p+2,p+n+1,comp);
	sta[++top] = p[1]; sta[++top] = p[2];
	for(int i = 3; i <= n; i++)
	{
		while(top >= 2 && Cro(sta[top]-sta[top-1],p[i]-sta[top-1]) <= 0) top--;
		sta[++top] = p[i];
	}
	for(int i = 1; i < top; i++) ans += dis(sta[i],sta[i+1]);
	ans += dis(sta[top],sta[1]);
	printf("%.2lf\n",ans);
	return 0;
}

例题：模板信用卡凸包

Graham 水平序扫描法

这个其实是上面那种算法的变种。一般来说这个是要比上面那个算法是快的，并且精度也比较高。

第一步：先把所有的点按横坐标为第一关键字纵坐标为第二关键字排序（这个和上面是不同的，这也是水平序比极角序要快的原因）。

第二步：从前往后扫一遍得到下凸壳。

第二步：从后往前倒着扫一遍得到上凸壳。

来张图模拟一下：

排完序之后的点集如图：

首先，我们先把前两个点入栈

因为 $p_1$ 的横坐标是最小的，所以他一定实在凸包上的
$P_2$ 无法确定，等之后在判断。

$ P_1P_2P_3$ 组成的是一个凸多边形，所以 $P_3$ 合法，入栈

$P_4$ 与 $P_3$ 的连线在 $P_2P_3$ 的逆时针方向，于是 $P_3$ 出栈，检查 $P_2$ 合法，之后 $P_4$ 入栈。

同理： $P_4$ 出栈， $P_5$ 入栈

$P_5$ 出栈， $P_6$ 出栈

$P_6$ 出栈， $P_7$ 入栈

$P_8$ 入栈

有木有发现一个很严重的问题，扫描到最后，我们只求出来个下凸壳。

那么上凸壳怎么求呢，很简单，倒着再跑一遍即可。

最后求出来的凸包长这样:

最后栈中的元素为： $P_1,P_2,P_7,P_8,P_8,P_5,P_3,P_1$

你会发现起点和终点在栈中出现了两次，但自己到自己的距离为零，所以不会对答案产生影响。

#include<iostream>
#include<cstdio>
#include<cmath>
#include<algorithm> 
using namespace std;
const int N = 1e5+10;
const double eps = 1e-8;
int n,top;
double ans;
struct point
{
	double x,y;
	point(){}
	point(double a,double b){x = a, y = b;}
}p[N],sta[N];
point operator + (point a,point b){return point(a.x+b.x,a.y+b.y);}
point operator - (point a,point b){return point(a.x-b.x,a.y-b.y);}
point operator * (point a,double k){return point(a.x*k,a.y*k);}
double Dot(point a,point b){return a.x*b.x+a.y*b.y;}
double Cro(point a,point b){return a.x*b.y-a.y*b.x;}
double dis(point a,point b)
{
	return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
int dcmp(double x)
{
	if(fabs(x) < eps) return 0;
	return x > 0 ? 1 : -1;
}
bool comp(point a,point b)
{
	if(dcmp(a.x-b.x) == 0) return a.y < b.y;
	return a.x < b.x;
}
int main()
{
	scanf("%d",&n);
	for(int i = 1; i <= n; i++) scanf("%lf%lf",&p[i].x,&p[i].y);
	sort(p+1,p+n+1,comp);
	sta[++top] = p[1]; sta[++top] = p[2];
	for(int i = 3; i <= n; i++) //下凸壳
	{
		while(top >= 2 && Cro(sta[top]-sta[top-1],p[i]-sta[top-1]) < 0) top--;
		sta[++top] = p[i];
	}
	int last = top;
	sta[++top] = p[n]; sta[++top] = p[n-1];//上凸壳
	for(int i = n-2; i >= 1; i--)
	{
		while(top - last >= 2 && Cro(sta[top]-sta[top-1],p[i]-sta[top-1]) < 0) top--;
		sta[++top] = p[i];
	}
	for(int i = 1; i < top; i++) ans += dis(sta[i],sta[i+1]);//起点和终点都在栈中出现了两次
	printf("%.2lf\n",ans);
	return 0;
}

动态凸包

先来道板子题：防线修建

我们考虑加入一个点之后，会对当前的凸包产生什么影响。

如图我们当前凸包长这样：

当我们加入一个点 $G$ 的时候，图中紫色的线就是新的凸包增加的那一部分。

根据静态凸包的构造方法（Anderw扫描法）：

将各点坐标以横坐标为第一关键字，纵坐标为第二关键字排序。
维护一下上凸壳和下凸壳。拿维护上凸壳举例，注意到如果之前有两点 $A$ 和 $B$ ，当插入新节点 $C$ 时，如果 $\overrightarrow{AB} \times \overrightarrow{AC} ≥0$ 时，这说明 $B$ 在线段 $AC$ 的下方，就可以将 $B$ 删除。我们可以看出此时凸包是一个类似单调栈的结构。

我们拿 $set$ 维护一下这个凸壳，当我们加入一个点的时候，检查一下周围的点是否会被删除即可。

显然每加入一个点，被删除的只能是他在 $set$ 中的前驱和后继。

这个题我们离线一下，倒着做，就只需要动态维护一个上凸壳，比较好写点.

感觉自己写的时候，尽量画个图参照着写，不然的话你就很容易搞混。

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cmath>
#include<set>
using namespace std;
const int N = 1e5+10;
const double eps = 1e-8;
int n,m,cntq;
double x,y,last;
bool vis[N];
struct node
{
	int opt,c;
	double ans;
}q[N];
struct point
{
	double x,y;
	double ang;
	point(){}
	point(double a,double b){x = a, y = b;}
	friend bool operator < (point a,point b)
	{
		if(a.x == b.x) return a.y < b.y;
		return a.x < b.x;
	}
}p[N],o;
int dcmp(double x)
{
	if(fabs(x) < eps) return 0;
	return x > 0 ? 1 : -1;
}
point operator + (point a,point b){return point(a.x+b.x,a.y+b.y);}
point operator - (point a,point b){return point(a.x-b.x,a.y-b.y);}
point operator * (point a,double k){return point(a.x*k,a.y*k);}
double Dot(point a,point b){return a.x*b.x+a.y*b.y;}
double Cro(point a,point b){return a.x*b.y-a.y*b.x;}
double dis(point a,point b){return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));}
set<point> s;
void insert(point now)
{
	set<point>::iterator itl,itr,tmp;
	itr = s.lower_bound(now);//itl 是他的前驱，itr是他的后继
	itl = itr; --itl;
	if(Cro(now-*itl,*itr-*itl) > 0) return;//如果这个点在凸包内的话就不用进行任何操作
	last -= dis(*itl,*itr);
	while(itr != s.end())//s.end返回的是set最后一个元素的下一个位置
	{
		tmp = itr; ++itr;//tmp指now的后继，itr为now的后继的后继
		if(itr == s.end() || Cro(now-*itr,*tmp-*itr) < 0) break;//可以参考上面的那个图理解一下
		last -= dis(*tmp,*itr);//把这条边从答案中减去
		s.erase(tmp);//tmp这个点不在凸包中了
	}
	while(itl != s.begin())
	{
		tmp = itl; --itl;//tmp为now的前驱，itl为now的前驱的前驱
		if(Cro(now-*itl,*tmp-*itl) > 0) break;
		last -= dis(*tmp,*itl);//把这条边的贡献从答案中减去
		s.erase(tmp);
	}
	s.insert(now);
	itl = itr = s.find(now);
	--itl; ++itr;//找now在凸包中的上一个点和下一个点
	last += dis(now,*itl) + dis(now,*itr);
}
int main()
{
	scanf("%d%lf%lf",&n,&x,&y);
	point a = point(0,0);
	point b = point(n,0); 
	point c = point(x,y);
	s.insert(a); s.insert(b); s.insert(c);
	last += dis(b,c) + dis(c,a);
	scanf("%d",&m);
	for(int i = 1; i <= m; i++) scanf("%lf%lf",&p[i].x,&p[i].y);
	scanf("%d",&cntq);
	for(int i = 1; i <= cntq; i++)
	{
		scanf("%d",&q[i].opt);
		if(q[i].opt == 1)
		{
			scanf("%d",&q[i].c);
			vis[q[i].c] = 1;
		}
	}
	for(int i = 1; i <= m; i++) if(!vis[i]) insert(p[i]);//先把所有没被消灭的点加入凸包中
	for(int i = cntq; i >= 1; i--)
	{
		if(q[i].opt == 1) insert(p[q[i].c]);
		else q[i].ans = last;
	}
	for(int i = 1; i <= cntq; i++) if(q[i].opt == 2) printf("%.2lf\n",q[i].ans);
	return 0;
}

旋转卡壳

这个他的原理是什么呢？主要是在一个凸包上，他的很多函数都是单峰的。比如固定一条多边形的边，其他的点和这条边组成的面积是单峰的。

求凸包直径

凸包的直径：定义凸包上任意两点之间距离的最大值。

怎么求？暴力枚举 $O(n^2) $

有一个结论：如果一个点到一条边的距离最远，那么点到线段两端点中一个的距离一定是最远的。

优化一下，首先我们可以固定一条边，然后点到这条边的距离是一个单峰的函数，显然我们是可以三分出来的，复杂度 $O(nlogn)$ 。

继续优化，旋转卡壳利用了一个更好的性质,假设说我们固定的这条边是不断滚动的，如图：

先固定 $BC$ 这条边，在固定 $CD$ , 之后接着固定 $DE$ ,之后…

你会发现，当你固定 $BC$ 这条边的时候，如果 $E$ 没有 $F$ 优，那么当你固定 $CD$ 的时候， $E$ 肯定不会比 $F$ 优。

也就是说：我们的边是顺时针的枚举，而此时距离他们最远的点也是顺时针的出现！！

然后我们就可以通过双指针扫一遍就可以得到这个凸包的直径。

基本上旋转卡壳就是这么一个流程。

一个动图演示一下：

double Get_zhijing()
{
	double res = 0;
	sta[top + 1] = p[1];//把 sta[top]sta[1] 这条边也要加进去 
	for(int i = 1, j = 2; i <= top; i++)
	{
		while(Cro(sta[i]-sta[j],sta[i+1]-sta[j]) < Cro(sta[i]-sta[j+1],sta[i+1]-sta[j+1])) j = j % top + 1;//找距离 sta[i]sta[i+1] 最远的点
		res = max(res,max(dis(sta[i],sta[j]),dis(sta[i+1],sta[j])));
	}
	return res;
}

例题：旋转卡壳

最小矩形覆盖

让你找一个矩形使得这个矩形能覆盖住点集中所有的点，且矩形的，面积最小。

怎么求？暴力枚举 $O(n^4)$ .

考虑到这个矩形他肯定有一条边是会和凸包的一条边重合的（这个并不好证，但反正他是正确的）

当我们固定一条边的时候，我们还需要知道，每个点到这条边投影的最大值和最小值（矩形的左右侧），和点到这条边距离的最大值（矩形的上底边）。

结合图理解一下：

假如说，我们要求这个多边形的最小矩形覆盖：

当我们固定 $DE$ 这条边的时候，我们最小的矩形长这样：

很容易发现，矩形的高就是点到 $ED$ 距离的最大值，矩形的宽就是点到 $ED$ 这条直线投影的最小值和最大值的差。

这三个东西，都可以拿旋转卡壳维护出来 (也可以在凸包上三分)，复杂度 $O(n)$

半平面交

如图：一条直线把平面分为了两部分，我们通常把直线的左侧平面称为他的半平面。

如下图所示，直线 $P_1P_2$ 的半平面就是黄色的部分。

注意直线的方向不同，所构成的半平面也是不同的，因此在做半平面交的时候要注意一下方向问题。

半平面交：多个半平面同时覆盖的区域叫做半平面交，如图，红色区域就是那三条直线的半平面交。

我们还可以观察到，第一个图中直线 $P_1P_2$ 的半平面区域 $ax+by+c > 0$ ,而直线 $P_2P_1$ 的半平面区域

则是 $ax+by+c < 0$ .

所以半平面交的本质就是方程组：$a_ix+b_iy+c > 0 $ 或 $a_ix+b_iy+c < 0$ 的合法解的区域（这好像是线性规划的柿子）。

一般半平面交的题都是让你列出等式，然后构造半平面交求解。

怎么求？

我们发现半平面交的区域是一个凸多边形（或者空集，线段）。

计算半平面交的一个好方法就是增量法，即初始答案为整个平面，然后注意加入各个平面，维护当前的半平面交。

一开始我们像凸包那样把所有的直线按极角排一下序，这样可以保证我们当前合法的直线的交就是我们的半平面交，然后考虑依次加入每一条线段。

我们开两个队列，一个维护合法的直线，另一个维护半平面交的交点，以便确定半平面交的形态。

如果前面两条合法直线 $A,B$ 的交点在这条直线 $C$ 的右侧，那么直线 $A$ 和直线 $B$ 组成的平面交会被直线 $C$ 覆盖掉一部分，那么直线 $B$ 就不合法了，把他从队列中弹出来即可，同时在第二个队列中把 $A$ 和 $B$ 的交点出队，重复上述过程，直到之前加入的直线都合法为止。最后在把直线 $C$ 和 $AC$ 的交点入队即可。

如图：

如果交点 $G$ 在直线左侧，说明我们队列中所有的直线都是合法的，这时候，我们把直线 $C$ 以及他与 $B$ 的交点入队即可。

因为最后我们要围成一个多边形，所以你现在加入的这一条直线可能会使对首的直线不合法，举个例子：

当我们加入 $DH$ 这条直线的时候，合法的半平面交的区域是图中的黄色部分，显然直线 $AB$ 不在这个区域内。

那么我们需要开个双端队列来维护一下，每次加入一条直线的时候，检查队头和队尾，去除不合法的元素即可。

当我们半平面交是无限大的时候怎么办？如图：

这时候我们在一开始可以人为地加入四条边界直线，这样半平面交就不会无限大啦。

愉快的代码时间：

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cmath>
using namespace std;
const double lim = 1000.0;
const double eps = 1e-8;
const int N = 1e5+10;
int n,m,l,r,cnt;
double area;
struct point
{
	double x,y;
	point(){}
	point(double a,double b){x = a, y = b;}
}sta[N],p[N],ans[N];
typedef point Vector;
struct line
{
	point x,y;
	double ang;
	line(){}
	line(point a,point b)
	{
		x = a; y = b;
		ang = atan2(b.y-a.y,b.x-a.x);
	}
}L[N],q[N];
point operator + (point a,point b){return point(a.x+b.x,a.y+b.y);}
point operator - (point a,point b){return point(a.x-b.x,a.y-b.y);}
point operator * (point a,double k){return point(a.x*k,a.y*k);}
double Dot(point a,point b){return a.x*b.x+a.y*b.y;}
double Cro(point a,point b){return a.x*b.y-a.y*b.x;}
double dis(point a,point b){return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));}
bool OnRight(point p,line s){return Cro(p-s.x,s.y-s.x) > 0;}
int dcmp(double x)
{
	if(fabs(x) < eps) return 0;
	return x > 0 ? 1 : -1;
}
bool comp(line a,line b)
{
	if(a.ang == b.ang) return OnRight(b.x,a);
	return a.ang < b.ang;
}
point Root_LL(line a,line b)
{
	Vector v = a.y-a.x, u = b.y-b.x, w = a.x-b.x;
	double t = Cro(w,u)/Cro(u,v);
	return a.x + v * t;
}
int HPI()
{
	int m = 0;
	sort(L+1,L+cnt+1,comp);
	l = 1, r = 1; q[1] = L[1];
	for(int i = 2; i <= cnt; i++)
	{
		if(dcmp(L[i].ang - L[i-1].ang) == 0) continue;
		while(l < r && OnRight(sta[r-1],L[i])) r--;
		while(l < r && OnRight(sta[l],L[i])) l++;
		q[++r] = L[i]; 
		if(l < r) sta[r-1] = Root_LL(q[r],q[r-1]);
	}
	while(l < r && OnRight(sta[r-1],q[l])) r--;//删去多余的点
	while(l < r && OnRight(sta[l],q[r])) l++;
	if(r - l + 1 <= 2) return 0;
	sta[r] = Root_LL(q[l],q[r]);
	for(int i = l; i <= r; i++) ans[++m] = sta[i];
	return m;
}
int main()
{
	scanf("%d",&n);
	for(int i = 1; i <= n; i++)
	{
		scanf("%d",&m);
		for(int j = 1; j <= m; j++) scanf("%lf%lf",&p[j].x,&p[j].y);
		for(int j = 1; j < m; j++) L[++cnt] = line(p[j],p[j+1]);
		L[++cnt] = line(p[m],p[1]);
	} 
	point a = point(-lim,lim);
	point b = point(-lim,-lim);
	point c = point(lim,-lim);
	point d = point(lim,lim);
	L[++cnt] = line(a,b); L[++cnt] = line(b,c); L[++cnt] = line(c,d); L[++cnt] = line(d,a);
	n = HPI();
	for(int i = 2; i < n; i++) area += Cro(ans[i]-ans[1],ans[i+1]-ans[1]);
	printf("%.3lf\n",area/2);
	return 0;
}

随机增量法

随机增量算法是计算几何的一个重要算法，它对理论知识要求不高，算法时间复杂度低，应用范围广大。(扯淡)

增量法 (Incremental Algorithm) 的思想与第一数学归纳法类似，它的本质是将一个问题化为规模刚好小一层的子问题。解决子问题后加入当前的对象。写成递归式是： $T(n) = T(n-1)+g(n)$ 。 --------wiki

直接看题吧：最小圆覆盖

怎么说呢？不会做，直接看题解。

首先，假设说我们已经求出来了前 $i-1$ 个点的最小圆覆盖，现在我们要新加入一个点 $i$ ,如果点 $i$ 在最小圆内，好继续往里面加点，否则的话第 $i$ 个点肯定在新的最小圆上。

这时候我们需要求出过点 $i$ 且覆盖前 $i-1$ 个点的最小圆。

由于三点确定一个圆，所以我们还需要确定一下其余的两个点。

第一步：不断枚举不在最小圆上的点 $j$ , 来扩大我们最小圆的面积,此时的圆心为 $i+j\over 2$ ,半径为 $dis(i,j)$ ，之后做过 $i,j$ 两个点且覆盖前 $j-1 $ 个点的最小圆。

第三步：枚举剩下的 $j-1$ 个点 $k$ ，如果 $k$ 不在最小圆中，过 $i,j,k$ 三个点的圆，这个圆就是过点 $i$ 且覆盖前 $k$ 个点的最小圆。

复杂度分析：

复杂度。PNG.PNG

上述的是期望复杂度，一般情况下出题人可能会造一些数据卡掉这种解法。

肿么办，random_shuffle 一下即可。 $n^3$ 过百万。

至于怎么求过三个点的最小圆，你可以尝试解方程（反正太复杂了，我没看懂，第二篇题解有证明），你也可以做中垂线，求两条直线的交点。

#include<bits/stdc++.h>
using namespace std;
const double eps = 1e-15;
int n;
double r;
struct point
{
	double x,y;
	point(){}
	point(double a,double b){x = a, y = b;}
}p[100010],o;
typedef point Vector;
struct line
{
	point x,y;
	line(){}
	line(point a,point b){x = a, y = b;}
};
point operator + (point a,point b){return point(a.x+b.x,a.y+b.y);}
point operator - (point a,point b){return point(a.x-b.x,a.y-b.y);}
point operator * (point a,double k){return point(a.x*k,a.y*k);}
double Dot(point a,point b){return a.x*b.x+a.y*b.y;}
double Cro(point a,point b){return a.x*b.y-a.y*b.x;}
point retate(point a){return point(a.y,-a.x);}//向量顺时针旋转90.
double dis(point a,point b)
{
	return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
int dcmp(double x)
{
	if(fabs(x) < eps) return 0;
	return x > 0 ? 1 : -1;
}
line ZhongChuiXian(point a,point b)
{
	point v = retate(b-a);
	point A = point ((a.x+b.x)/2,(a.y+b.y)/2);
	return line(A,v);
}
point Root_LL(line a,line b)
{
	Vector v = a.y, u = b.y, w = a.x-b.x;
	double t = Cro(w,u)/Cro(u,v);
	return a.x + v * t;
}
point get(point a,point b,point c)
{
	line s1 = ZhongChuiXian(a,b);
	line s2 = ZhongChuiXian(b,c);
	return Root_LL(s1,s2);	//两个点中垂线的交点
}
int main()
{
	srand(time(0)); rand();
	scanf("%d",&n);
	for(int i = 1; i <= n; i++) scanf("%lf%lf",&p[i].x,&p[i].y);
	for(int i = 1; i <= n; i++)
	{
		int l = rand() % n + 1;
		int r = rand() % n + 1;
		swap(p[l],p[r]);
	}
	o = p[1]; r = 0;
	for(int i = 2; i <= n; i++)//做覆盖前 i个点的最小圆
	{
		if(dcmp(dis(p[i],o)-r) <= 0) continue;
		o = p[i]; r = 0;//如果第i个点不在之前的最小圆上，那么点i一定在新的最小圆上
		for(int j = 1; j < i; j++)//做过点i，且覆盖前i-1个点的最小圆
		{
			if(dcmp(dis(p[j],o)-r) <= 0) continue;//不断往里面加点
			o.x = (p[i].x+p[j].x)/2;
			o.y = (p[i].y+p[j].y)/2;
			r = dis(p[i],o);
			for(int k = 1; k < j; k++)//做过点i和点j且覆盖前j-1个点的最小覆盖圆
			{
				if(dcmp(dis(p[k],o)-r) <= 0) continue;//不断的往里面加点
				o = get(p[i],p[j],p[k]);//三个点确定一个圆
				r = dis(p[i],o); 
			}
		}
	}
	printf("%.10lf\n",r);
	printf("%.10lf %.10lf\n",o.x,o.y);
	return 0;
}

例题：信号塔

三角剖分

太难了，弃坑。有兴趣的可以到 wiki 上了解一下。

闵可夫斯基和

就是说，平面上有两个凸包 $A$ 和 $B$ ,如图：

每个凸包可以看做一个点集。

定义两个凸包的闵可夫斯基和 $A+B = {\alpha + \beta }, \alpha\in A,\beta\in B $ 。

也就是说我们从 $A$ 和 $B$ 中任选两个点出来，在把这两个点的横纵坐标相加，把得到的新的点放到一个集合中，这个集合就叫做两个凸包的闵可夫斯基和。

如下图，粉色区域就是三角形和四边形的闵可夫斯基和。

结论1：闵可夫斯基和是凸性的。

证明：

一个多边形是凸性的，等价于 $\alpha ,\beta\in A, x\beta+(1-x)\alpha \in A$ 。

也就是说如果 $\alpha$ 和 $\beta$ 在这个多边形内，那么 $\alpha\beta$ 这一条线段上所有的点也是在这一个凸多边形上的。

有 $\alpha_1\alpha_2\in A,\beta_1\beta_2\in B$ ，则有 $\alpha_1+\beta_1,\alpha_2+\beta_2 \in {A+B}$

那么则有 $(\alpha_1+\beta_1)x+(1-x)\ (\alpha_2+\beta_2) = \alpha_1x+(1-x)\alpha_2 + \beta_1x+(1-x)\beta_2 $

其中 $\alpha_1x+(1-x)\alpha_2\in A$ , $\beta_1x+(1-x)\beta_2\in B$ ，那么就有 $(\alpha_1+\beta_1)x+(1-x)\ (\alpha_2+\beta_2)\in A+B$ ,所以闵可夫斯基和也是个凸集。

结论2：任意一条凸包 $A$ 或 $B$ 的边，那么一定也是 $A+B $ 的边。

证明：

一条线段在一个凸包上，等价于凸包上有一条线段和这一条线段的斜率相等，且线段外侧没有其他的点。

还是那刚才的图来举例：

比如说我们确定 $\alpha_1\alpha_2$ 这条边是否出现在 $A+B$ 上，我们取 $B$ 中离 $\alpha_1 \alpha_2$ 这条边最远的一个切点 $ \beta$ ，那么显然 $\alpha_1+\beta$ 和 $\alpha_2+\beta$ 都属于 $A+B$ , 且这一条线段的斜率等于 $\alpha_1\alpha_2$ 的斜率，在这一条线段之外没有其他的点（因为 $\alpha_1\alpha_2$ 是凸包 $ A$ 这个方向的最大值， $\beta$ 是凸包 $B$ 这个方向的最大值，两个最大值相加肯定也是最大的）。

闵可夫斯基和的求法:

假设我们有 $n$ 个点的凸包 $A$ 和 $m$ 个点的凸包 $B$ , 我们把凸包 $A$ 和 $B$ 上的边表示出来也就是 $\vec{A_iA_{i+1}}$ 和 $\vec{B_iB_{i+1}}$ 根据结论 $2$ 可得，这些边都是 $A$ 和 $B$ 的闵可夫斯基和上的边，现在我们就考虑怎么把这些边组成一个凸多边形。

首先，我们把每条边按斜率大小 $O(n)$ 的排序，然后我们把这些边按排完序之后的顺序依次接起来，就可以得到闵可夫斯基和的形态。然后把 $A$ 和 $B$ 中横纵坐标的最小值相加，就可以确定出闵可夫斯基和的横纵坐标的最小值，显然这个凸包的形态就是唯一的。

void Minkowski()
{
	for(ll i=1;i<n;i++) s1[i]=C1[i+1]-C1[i];s1[n]=C1[1]-C1[n];
	for(ll i=1;i<m;i++) s2[i]=C2[i+1]-C2[i];s2[m]=C2[1]-C2[m];
	A[tot=1]=C1[1]+C2[1];
	ll p1=1,p2=1;
	while(p1<=n&&p2<=m) ++tot,A[tot]=A[tot-1]+(s1[p1]*s2[p2]>=0?s1[p1++]:s2[p2++]);
	while(p1<=n) ++tot,A[tot]=A[tot-1]+s1[p1++];
	while(p2<=m) ++tot,A[tot]=A[tot-1]+s2[p2++];
}

例题： JSOI2018战争

这个题让我们判断 $A+\alpha$ 和 $B$ 是否重合。

对于 $A$ 中一个点 $x$ 和 $B$ 中一个点 $y$ , 那么 $\alpha = y-x$

则不合法的 $\alpha $ 的集合就是 $C= \{y-x\},y\in B,x\in A$ 。

我们把 $A$ 的坐标去相反数就变为了 $B$ 和 $A$ 的闵可夫斯基和。

每次询问判断一个点是否在这个凸包里面，二分即可。

pick定理

这个当公式记就好了，一般来说没什么用。

给定顶点均为整点的简单多边形，皮克指出了其面积 $A$ 与内部格点数 $i$ 和边上格点数 $b$ 的关系： $A= i + {b\over 2} - 1$

具体证明：论文

高维中，皮克定理和欧拉公式 $V - E + F = 2$ 是等价的. ( $V,E,F$ 表示简单几何体的顶点数，边数，面数)

例题：poj 1265

平面图转对偶图

简单来说，平面图就是点和点之间有连边的图，如图：

而他的对偶图呢，就是把每个面抠出来，标上号，如果两个面之间有交点，那么就在对偶图中连一条边。

最后的对偶图建出来长这样：

对偶图的性质：平面图的最小割等价于它对偶图的最短路。（~~这好像就不属于计算几何的内容了~~）。

至于平面图怎么转对偶图，啊我也不会，主要是网上的博客讲的太模糊了。

计算几何杂（水）题选讲

1.城墙

给你平面上一些点，让你求一个围墙能覆盖住这些点，且点到多边形的距离不能小于 $L$ .

求这个围墙的最小周长。 $n,L\leq 10^5$

首先，我们考虑一下当 $L = 0$ 的情况，很显然答案就是这个点集的凸包。

当 $L \neq 0$ 的时候，我们这时候只需要考虑凸包顶点上的点。

对于一个点距离小于 $L$ 的点，显然在以这个点为圆心，半径为 $L$ 的圆内。

我们把这个图画一下：

根据肉眼观察，要求的是最外面那一圈的周长，并且他可以分为两部分：原凸包上的边和若干个圆弧组成。

根据肉眼观察可得若干个圆弧可以拼接成一个圆。

证明： $\angle AKL = 180 - \angle CAI,\angle IMO = 180-\angle AIG,\angle PGN = 180-\angle FGI$ , $\angle SFQ= 180-\angle EFG, \angle UER = 180-\angle CEF,\angle JCT = 180-\angle ACE$

外角相加等于 $n*180 - (n-2)*180 = 360$

所以我们最后的答案就是：凸包的周长加上一个圆的周长。

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cmath>
using namespace std;
const int N = 1e5+10;
const double pi = acos(-1.0);
const double eps = 1e-8;
int n,r,top;
double ans;
struct point
{
	double x,y;
	point(){}
	point(double a,double b){x = a, y = b;}
}sta[N],p[N];
point operator + (point a,point b){return point(a.x+b.x,a.y+b.y);}
point operator - (point a,point b){return point(a.x-b.x,a.y-b.y);}
point operator * (point a,double k){return point(a.x*k,a.y*k);}
point retate(point a,double alpha){return point(a.x*cos(alpha)-a.y*sin(alpha),a.x*sin(alpha)+a.y*cos(alpha));}
double Dot(point a,point b){return a.x*b.x+a.y*b.y;}
double Cro(point a,point b){return a.x*b.y-a.y*b.x;}
double dis(point a,point b){return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));}
int dcmp(double x)
{
	if(fabs(x) < eps) return 0;
	return x > 0 ? 1 : -1;
}
bool comp(point a,point b)
{
	if(a.x == b.x) return a.y < b.y;
	return a.x < b.x;
}
int main()
{
	scanf("%d%d",&n,&r);
	for(int i = 1; i <= n; i++) scanf("%lf%lf",&p[i].x,&p[i].y);
	sort(p+1,p+n+1,comp);
	sta[++top] = p[1]; sta[++top] = p[2];
	for(int i = 3; i <= n; i++)
	{
		while(top >= 2 && Cro(sta[top]-sta[top-1],p[i]-sta[top-1]) < 0) top--;
		sta[++top] = p[i];
	}
	int tmp = top;
	sta[++top] = p[n]; sta[++top] = p[n-1];
	for(int i = n-2; i >= 1; i--)
	{
		while(top-tmp >= 2 && Cro(sta[top]-sta[top-1],p[i]-sta[top-1]) < 0) top--;
		sta[++top] = p[i];
	}
	for(int i = 1; i < top; i++) ans += dis(sta[i],sta[i+1]);
	printf("%.0lf\n",ans+2*pi*r);
	return 0;
}

2.信用卡凸包

给你一些奇怪的形状的图形，求他们并在一起的图形的周长.

和上个题一样，答案可以分为两部分：原凸包上的边和若干段圆弧。

最后的答案为凸包周长加上一个圆的周长。

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cmath>
using namespace std;
const int N = 1e5+10;
const double pi = acos(-1.0);
const double eps = 1e-12;
int n,cnt,top;
double x,y,theta,a,b,r,ans;
struct point
{
	double x,y;
	point (double a = 0, double b = 0){x = a, y = b;} 
}p[N],sta[N],w[5];
typedef point Vector;
point operator + (point a,point b){return point(a.x+b.x,a.y+b.y);}
point operator - (point a,point b){return point(a.x-b.x,a.y-b.y);}
point operator * (point a,double b){return point(a.x*b,a.y*b);}
double Dot(point a,point b){return a.x*b.x+a.y*b.y;}
double Cro(point a,point b){return a.x*b.y-a.y*b.x;}
double len(point a){return sqrt(a.x*a.x+a.y*a.y);}
double dis(point a,point b)
{
	return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
point retate(point a,double alpha)
{
	double coss = cos(alpha), sinn = sin(alpha);
	return point(a.x*coss-a.y*sinn,a.x*sinn+a.y*coss);
}
int dcmp(double x)
{
	if(fabs(x) < eps) return 0;
	if(x > 0) return 1;
	return -1;
}
bool comp(point a,point b)
{
	double m = Cro(a-p[1],b-p[1]);
	if(dcmp(m) > 0) return 1;
	if(dcmp(m) == 0 && dis(a,p[1]) < dis(b,p[1])) return 1;
	return 0;
}
int main()
{
	scanf("%d",&n);
	scanf("%lf%lf%lf",&b,&a,&r);//先输入的是b后输入的是a，不要搞反了。
	a = a * 0.5 - r; b = b * 0.5 - r;
	for(int i = 1; i <= n; i++)
	{
		scanf("%lf%lf%lf",&x,&y,&theta);//向量旋转
		w[1] = point(a,b);
		w[2] = point(a,-b);
		w[3] = point(-a,b);
		w[4] = point(-a,-b);
		p[++cnt] = retate(w[1],theta)+point(x,y);
		p[++cnt] = retate(w[2],theta)+point(x,y);
		p[++cnt] = retate(w[3],theta)+point(x,y);
		p[++cnt] = retate(w[4],theta)+point(x,y);
	}
	for(int i = 1; i <= cnt; i++)
	{
		if(dcmp(p[i].y-p[1].y) < 0) swap(p[i],p[1]);
		else if(dcmp(p[i].y-p[1].y) == 0 && dcmp(p[i].x-p[1].x) < 0) swap(p[i],p[1]);
	}
	sort(p+2,p+cnt+1,comp);
	sta[++top] = p[1]; sta[++top] = p[2];//求凸包
	for(int i = 3; i <= cnt; i++)
	{
		while(top >= 2 && Cro(sta[top]-sta[top-1],p[i]-sta[top-1]) <= 0) top--;
		sta[++top] = p[i];
	}
	for(int i = 1; i < top; i++) ans += dis(sta[i],sta[i+1]);
	ans += dis(sta[top],sta[1]);
	ans += 2 * pi * r;
	printf("%.2lf",ans);
	return 0;
}

3.冷东波

因为求的是最小的时间，考虑二分答案。

我们现在时间确定了，那么每个巫师能打死的小精灵的数量也是确定的。

考虑用网络流求解。

从源点向巫师连一条容量为巫师能打死的小精灵的数量，从每个巫师向他能打到的小精灵连一条容量为 $1$ 的边。

再由每个小精灵向汇点连一条容量为 $1$ 的边。

这样跑出来的最大流，就是我们最多能打死的精灵的数量。

现在想想，怎么求每个巫师能打死那个精灵。 $n^3$ 暴力枚举一下。

因为每棵树的范围是一个圆，所以就相当于问你这个圆是否和巫师和小精灵的连线有交。

可以求出圆心到这条线段的最小距离，在和半径比较即可。

此外，还应该判断小精灵是否在巫师的攻击范围内，即点是否在圆内。

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<queue>
#include<cmath>
#include<cstring>
using namespace std;
const double eps = 1e-8;
const int inf = 1e7+10;
const int N = 2010;
int n,m,k,u,v,s,t,res,tot = 1;
int head[N],dep[N],tim[N],r[N],R[N];
bool used[N][N];
struct node
{
	int to,net,w;
}e[100010];
struct point
{
	int x,y;
	point(){}
	point(int a,int b){x = a, y = b;}
}c[N],p[N],tr[N];
typedef point Vector;
point operator + (point a,point b){return point(a.x+b.x,a.y+b.y);}
point operator - (point a,point b){return point(a.x-b.x,a.y-b.y);}
point operator * (point a,double k){return point(a.x*k,a.y*k);}
double Dot(point a,point b){return a.x*b.x+a.y*b.y;}
double Cro(point a,point b){return a.x*b.y-a.y*b.x;}
double dis(point a,point b){return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));}
double dis_PS(point p,point A,point B)//点到直线的距离
{
	if(Dot(p-A,B-A) < 0) return dis(p,A);
	if(Dot(p-B,A-B) < 0) return dis(p,B);
	return abs(Cro(A-p,B-p))/dis(A,B);
}
void add(int x,int y,int w)
{
	e[++tot].to = y;
	e[tot].w = w;
	e[tot].net = head[x];
	head[x] = tot;
}
bool bfs()//网络流
{
	queue<int> q;
	for(int i = 0; i <= t; i++) dep[i] = 0;
	q.push(s); dep[s] = 1;
	while(!q.empty())
	{
		int x = q.front(); q.pop();
		for(int i = head[x]; i; i = e[i].net)
		{
			int to = e[i].to;
			if(e[i].w && !dep[to])
			{
				dep[to] = dep[x] + 1;
				q.push(to);
				if(to == t) return 1;
			}
		}
	}
	return 0;
}
int dinic(int x,int flow)
{
	if(x == t || !flow) return flow;
	int rest = flow, val = 0;
	for(int i = head[x]; i && rest; i = e[i].net)
	{
		int to = e[i].to;
		if(!e[i].w || dep[to] != dep[x] + 1) continue;
		val = dinic(to,min(rest,e[i].w));
		if(val == 0) dep[to] = 0;
		e[i].w -= val; e[i^1].w += val; rest -= val;
	}
	return flow - rest;
}
bool judge(int mid)//二分答案
{
	s = 0, t = n+m+1, tot = 1, res = 0;
	memset(head,0,sizeof(head));
	for(int i = 1; i <= n; i++) add(s,i,(mid/tim[i])+1), add(i,s,0);
	for(int i = 1; i <= m; i++) add(i+n,t,1), add(t,i+n,0);
	for(int i = 1; i <= n; i++)
	{
		for(int j = 1; j <= m; j++)
		{
			if(used[i][j]) add(i,j+n,1), add(j+n,i,0);
		}
	}
	int flow = 0;
	while(bfs())
	{
		while(flow = dinic(s,inf)) res += flow;
	}
	return res >= m;
}
int main()
{
	scanf("%d%d%d",&n,&m,&k);
	for(int i = 1; i <= n; i++) scanf("%d%d%d%d",&c[i].x,&c[i].y,&r[i],&tim[i]);
	for(int i = 1; i <= m; i++) scanf("%d%d",&p[i].x,&p[i].y);
	for(int i = 1; i <= k; i++) scanf("%d%d%d",&tr[i].x,&tr[i].y,&R[i]);
	for(int i = 1; i <= n; i++)
	{
		for(int j = 1; j <= m; j++)
		{
			int flag = 0;
			double tmp = dis(c[i],p[j]);
			if(tmp > r[i]) continue;
			for(int l = 1; l <= k; l++) 
			{
				double tmp = dis_PS(tr[l],c[i],p[j]);
				if(tmp <= R[l]) flag = 1;
			}
			if(flag == 0) used[i][j] = 1;
		}
	}
	int L = 0, R = 500000, ans = -1;
	while(L <= R)
	{
		int mid = (L + R)>>1;
		if(judge(mid)) 
		{
			R = mid - 1;
			ans = mid;
		}
		else L = mid + 1;
	}
	printf("%d\n",ans);
	return 0;
}

4.射箭

考虑二分答案，转化为判定问题。

设抛物线的方程为： $y = ax^2+bx$ 。

显然能通关一关当且仅当： $y_1\leq ax^2+bx \leq y_2$ 。

化简一下： ${y_1\over x}\leq ax+b \leq {y_2\over x} $。

即： $ax+b-{y_1\over x} > 0$ , $ax+b-{y_2\over x} < 0$ 。

把 $a,b$ 看成 $x,y$ 就是半平面交的形式，即符合条件的 $a,b$ 的区域就是两个一次函数的交。

每次判定的时候，把前 $mid$ 关的限制条件扔进去跑半平面交，最后判断是否有交区域即可

剩下的就是你卡精度的问题了

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cmath>
using namespace std;
#define double long double
const double lim = 1e14;
const double eps = 1e-14;
const int N = 3e5+10;
int n,cnt,ans;
double x,l,r;
struct point
{
	double x,y;
	point(){}
	point(double a,double b){x = a; y = b;}
}sta[N];
typedef point Vector;
struct line
{
	int id;
	point x,y;
	double ang;
	line(){}
	line(point a,point b,int i)
	{
		x = a; y = b;
		ang = atan2(b.y-a.y,b.x-a.x);
		id = i;
	}
}L[N],q[N];
inline int read()
{
	int s = 0,w = 1; char ch = getchar();
	while(ch < '0' || ch > '9'){if(ch == '-') w = -1; ch = getchar();}
	while(ch >= '0' && ch <= '9'){s = s * 10 + ch - '0'; ch = getchar();}
	return s * w;
}
point operator + (point a,point b){return point(a.x+b.x,a.y+b.y);}
point operator - (point a,point b){return point(a.x-b.x,a.y-b.y);}
point operator * (point a,double k){return point(a.x*k,a.y*k);}
double Dot(point a,point b){return a.x*b.x+a.y*b.y;}
double Cro(point a,point b){return a.x*b.y-a.y*b.x;}
double dis(point a,point b){return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));}
bool OnRight(point p,line s){return Cro(p-s.x,s.y-s.x) > 0;}
bool comp(line a,line b)
{
	if(a.ang == b.ang) return OnRight(b.x,a);
	return a.ang < b.ang; 
}
int dcmp(double x)
{
	if(fabs(x) < eps) return 0;
	return x > 0 ? 1 : -1;
}
point Root_LL(line a,line b)
{
	Vector v = a.y-a.x, u = b.y-b.x, w = a.x-b.x;
	double t = Cro(w,u)/Cro(u,v);
	return a.x + v * t;
}
bool HPI(int mid)
{
	int l = 1, r = 0;
	for(int i = 1; i <= cnt; i++)
	{
		if(dcmp(L[i].ang-L[i-1].ang) == 0 || L[i].id > mid) continue;
		while(l < r && OnRight(sta[r-1],L[i])) r--;
		while(l < r && OnRight(sta[l],L[i])) l++;
		q[++r] = L[i];
		if(l < r) sta[r-1] = Root_LL(q[r],q[r-1]); 
	}
	while(l < r && OnRight(sta[r-1],q[l])) r--;
	while(l < r && OnRight(sta[l],q[r])) l++;
	if(r-l+1 <= 2) return 0;
	return 1;
}
int main()
{
	n = read();
	point a = point(-lim,eps); point b = point(-eps,eps);
	point c = point(-eps,lim); point d = point(-lim,lim);
	L[++cnt] = line(a,b,0); L[++cnt] = line(b,c,0);
	L[++cnt] = line(c,d,0); L[++cnt] = line(d,a,0);
	for(int i = 1; i <= n; i++) 
	{
		cin>>x>>l>>r;
		L[++cnt] = line(point(0,l/x),point(1,(l/x)-x),i);
		L[++cnt] = line(point(1,(r/x)-x),point(0,r/x),i);//下降方向 
	}
	sort(L+1,L+cnt+1,comp);
	int lx = 1, rx = n;
	while(lx <= rx)
	{
		int mid = (lx + rx)>>1;
		if(HPI(mid)) 
		{
			lx = mid + 1;
			ans = mid;
		}
		else rx = mid - 1;
	}
	printf("%d\n",ans);
	return 0;
}

5.赛车

solution

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cmath>
using namespace std;
#define double long double
const int N = 1e5+10;
const long long lim = 1e18+10;
const double eps = 1e-18;
int n,m,cnt,num,l,r,ans[N];
double k[N],w[N];
struct point
{
	double x,y;
	point(){}
	point(double a,double b){x = a, y = b;}
}sta[N];
typedef point Vector;
struct line
{
	int id;
	point x,y;
	double ang;
	line(){}
	line(point a,point b,int i)
	{
		x = a; y = b; id = i;
		ang = atan2(b.y-a.y,b.x-a.x);
	}
}q[N],L[N];
int dcmp(double x)
{
	if(fabs(x) < eps) return 0;
	return x > 0 ? 1 : -1;
}
point operator + (point a,point b){return point(a.x+b.x,a.y+b.y);}
point operator - (point a,point b){return point(a.x-b.x,a.y-b.y);}
point operator * (point a,double k){return point(a.x*k,a.y*k);}
double Dot(point a,point b){return a.x*b.x+a.y*b.y;}
double Cro(point a,point b){return a.x*b.y-a.y*b.x;}
double dis(point a,point b){return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));}
bool OnRight(point p,line a){return dcmp(Cro(p-a.x,a.y-a.x)-eps) >= 0;}
point Root_LL(line a,line b)
{
	Vector v = a.y-a.x, u = b.y-b.x, w = a.x-b.x;
	double t = Cro(w,u)/Cro(u,v);
	return a.x + v * t;
}
bool comp(line a,line b)
{
	if(dcmp(a.ang-b.ang) == 0) return OnRight(b.x,a);
	return a.ang < b.ang;
}
int HPI()
{
	l = 1, r = 1; q[1] = L[1];
	for(int i = 2; i <= cnt; i++)
	{
		if(dcmp(L[i].ang-L[i-1].ang) == 0) continue;
		while(l < r && OnRight(sta[r-1],L[i])) r--;
		while(l < r && OnRight(sta[l],L[i])) l++;
		q[++r] = L[i];
		if(l < r) sta[r-1] = Root_LL(q[r],q[r-1]);
	}
	while(l < r && OnRight(sta[r-1],q[l])) r--;
	while(l < r && OnRight(sta[l],q[r])) l++;
	sta[r] = Root_LL(q[l],q[r]);
	if(r-l+1 == 2) return 0;
	return r-l+1;
}
int main()
{
	scanf("%d",&n);
	for(int i = 1; i <= n; i++) cin>>w[i];
	for(int i = 1; i <= n; i++) cin>>k[i];
	point a = point(0,0); point b = point(lim,0);
	point c = point(lim,lim); point d = point(0,lim);
	L[++cnt] = line(a,b,0); L[++cnt] = line(b,c,0);
	L[++cnt] = line(c,d,0); L[++cnt] = line(d,a,0);
	for(int i = 1; i <= n; i++) 
	{
		point a = point(0,w[i]); point b = point(1,k[i]+w[i]); 
		L[++cnt] = line(a,b,i);
	}
	sort(L+1,L+cnt+1,comp);
	m = HPI();
	for(int i = l; i <= r; i++)
	{
		for(int j = 1; j <= n; j++)
		{
			if(w[j] == w[q[i].id] && k[j] == k[q[i].id]) ans[++num] = j;
		}
	}
	sort(ans+1,ans+num+1);
	printf("%d\n",num);
	for(int i = 1; i <= num; i++) printf("%d ",ans[i]);
	printf("\n");
	return 0;
}

6.逃考

给你一个平面直角坐标系，每个亲戚有其监视范围，问你小凸从 $(x,y)$ 走到边界上，最少要经过多少个被人监视的区域。

不难想到：对于每个亲戚，向和他监视区域相临亲戚的连一条公共边，监视区域和边界相邻的向汇点连一条边，以一开始监视小凸的亲戚为源点，跑最短路，从源点到汇点的最短距离就是答案。

关键是怎么把图建出来。

对于两个点的情况，以两个点的中垂线为界，分成的两部分分别被两个人监控。如图：

黄色的区域就是 $A$ 的监视范围。

那么对于 $n$ 个点的情况呢？显然就是这几个半平面区域的交集。

所以，我们可以枚举每个亲戚，然后对他做一遍半平面交，最后队列中剩下的元素就是监视区域和他相连的亲戚。

在做半平面交的时候也要把四个边界加进去，如果队列中剩下的元素中有边界，那么就把这个点和汇点连一条边。

数据可能会有点不在这个矩形内部，特判掉即可。

最后每个亲戚的监视区域长成这个样子：

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cmath>
#include<queue>
#include<cstring>
using namespace std;
const double eps = 1e-8;
const int N = 1e5+10;
const int inf = 1e7+10;
int T,n,cnt,tot,s,l,r,d[N],head[N];
bool vis[N],used[N];
struct node
{
	int to,net,w;
}e[N];
struct point
{
	double x,y;
	point(){}
	point(double a,double b){x = a, y = b;}
}sta[N],p[N],st,en;
typedef point Vector;
struct line
{
	int id;
	double ang;
	point x; Vector v;
	line(){}
	line(point a,Vector b,int i)
	{
		x = a; v = b; id = i;
		ang = atan2(b.y,b.x);
	}
}L[N],q[N];
int dcmp(double x)
{
	if(fabs(x) < eps) return 0;
	return x > 0 ? 1 : -1;
}
point operator + (point a,point b){return point(a.x+b.x,a.y+b.y);}
point operator - (point a,point b){return point(a.x-b.x,a.y-b.y);}
point operator * (point a,double k){return point(a.x*k,a.y*k);}
point retate(point a){return point(-a.y,a.x);}//逆时针旋转90度 
double Dot(point a,point b){return a.x*b.x+a.y*b.y;}
double Cro(point a,point b){return a.x*b.y-a.y*b.x;}
double dis(point a,point b){return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));}
bool OnRight(point p,line a){return dcmp(Cro(p-a.x,a.v)) >= 0;}
point Root_LL(line a,line b)
{
	Vector v = a.v, u = b.v, w = a.x-b.x;
	double t = Cro(w,u)/Cro(u,v);
	return a.x + v * t;
}
bool comp(line a,line b)
{
	if(dcmp(a.ang-b.ang) == 0) return OnRight(b.x,a);
	return a.ang < b.ang;
}
void add(int x,int y,int w)
{
	e[++tot].to = y;
	e[tot].w = w;
	e[tot].net = head[x];
	head[x] = tot;
}
void HPI()
{
	sort(L+1,L+cnt+1,comp);
	l = 1, r = 1; q[1] = L[1];
	for(int i = 2; i <= cnt; i++)
	{
		if(dcmp(L[i].ang-L[i-1].ang) == 0) continue;
		while(l < r && OnRight(sta[r-1],L[i])) r--;
		while(l < r && OnRight(sta[l],L[i])) l++;
		q[++r] = L[i];
		if(l < r) sta[r-1] = Root_LL(q[r],q[r-1]);
	}
	while(l < r && OnRight(sta[r-1],q[l])) r--;
	while(l < r && OnRight(sta[l],q[r])) l++;
	sta[r] = Root_LL(q[l],q[r]);
}
void Dij()
{
	priority_queue<pair<int,int> > Q;
	for(int i = 0; i <= n+2; i++) d[i] = inf, vis[i] = 0;
	d[s] = 0; Q.push(make_pair(0,s));
	while(!Q.empty())
	{
		int t = Q.top().second; Q.pop();
		if(vis[t]) continue;
		vis[t] = 1;
		for(int i = head[t]; i; i = e[i].net)
		{
			int to = e[i].to;
			if(d[to] > d[t] + e[i].w)
			{
				d[to] = d[t] + e[i].w;
				Q.push(make_pair(-d[to],to));
			}
		}
	}
}
int main()
{
	scanf("%d",&T);
	while(T--)
	{
		scanf("%d",&n); tot = 0;
		memset(used,0,sizeof(used));
		memset(head,0,sizeof(head));
		scanf("%lf%lf%lf%lf",&en.x,&en.y,&st.x,&st.y);
		for(int i = 1; i <= n; i++) 
		{
			scanf("%lf%lf",&p[i].x,&p[i].y);
			if(p[i].x > en.x || p[i].y > en.y) used[i] = 1;
		}
		if(n == 0)
		{
			printf("%d\n",0);
			continue;
		}
		double res = inf;
		for(int i = 1; i <= n; i++)
		{
			if(res > dis(p[i],st)) s = i;
			res = min(res,dis(p[i],st));	
		}
		for(int i = 1; i <= n; i++)
		{
			cnt = 0;
			if(used[i]) continue;
			L[++cnt] = line(point(0,en.y),point(0,-1),n+1); 
			L[++cnt] = line(point(en.x,0),point(0,1),n+1);
			L[++cnt] = line(point(0,0),point(1,0),n+1); 
			L[++cnt] = line(en,point(-1,0),n+1);
			for(int j = 1; j <= n; j++)
			{
				if(j == i || used[j]) continue;
				Vector v = retate(p[j]-p[i]);//PiPj向量逆时针旋转90度 
				point o = point((p[i].x+p[j].x)/2,(p[i].y+p[j].y)/2);//注意方向问题
				L[++cnt] = line(o,v,j); 
			}
			HPI();
			for(int j = l; j <= r; j++) add(i,q[j].id,1);
		}
		Dij();
		printf("%d\n",d[n+1]);
	}
	return 0;
}

7.小凸想跑步

由题意可得： $S_{\Delta PA_0A_1} < S_{\Delta PA_iA_{i+1}}$

然后就开始暴力展开，推柿子。

$\vec {A_0A_1} \times \vec{A_0P} <\vec{A_iA_{i+1}} \times \vec {A_iP}$

设 $P = (x,y) $ , $A_i = (x_i,y_i)$ ,则有：

$(x_1-x_0,y_1-y_0) \times (x-x_0,y-y_0) < (x_{i+1}-x_i,y_{i+1}-y_i) \times (x-x_i,y-y_i)$

$(x_1-x_0)(y-y_0) - (x-x_0)(y_1-y_0) < (x_{i+1}-x_i)(y-y_i) - (x-x_i) (y_{i+1}-y_i)$

$x_1y-x_1y_0-x_0y+x_0y_0 - (xy_1-xy_0-x_0y_1+x_0y_0) < x_{i+1}y - x_{i+1}y_i- x_iy + x_iy_i - (xy_{i+1} - xy_i - x_{i}y_{i+1} + x_iy_i)$

$x_1y-x_1y_0-x_0y-xy_1+xy_0+x_0y_1 < x_{i+1}y-x_{i+1}y_i - x_iy - xy_{i+1} + xy_i + x_iy_{i+1} $

$(y_0-y_1) x + (x_1-x_0)y + x_0y_1 - x_1y_0 < (y_i-y_{i+1})x + (x_{i+1}-x_i)y + x_iy_{i+1} -x_{i+1}y_i $

$(y_0-y_1-y_i+y_{i+1})x + (x_1-x_0-x_{i+1}+x_i)y + x_0y_1-x_1y_0 + x_{i+1}y_{i} - x_iy_{i+1} < 0$

$(y_1-y_0-y_{i+1}+y_i) x + (x_0-x_1 - x_i+x_{i+1})y + x_1 y_0-x_0y_1 + x_{i}y_{i+1} - x_{i+1}y_{i} > 0$

那么合法的 $p$ 的区域就是半平面的交。

最后的答案就是多边形面积除以半平面交区域的面积。

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cmath>
using namespace std;
const int N = 1e5+10;
const double eps = 1e-8;
int n,l,r,cnt,top;
double ans1,ans2;
struct point
{
    double x,y;
    point(){}
    point(double a,double b){x = a, y = b;}
}p[N],sta[N],Ans[N];
typedef point Vector;
struct line
{
    point x,v;
    double ang;
    line(){}
    line(point a,point b)
    {
        x = a; v = b;
        ang = atan2(b.y,b.x);
    }
}L[N],q[N];
int dcmp(double x)
{
    if(fabs(x) < eps) return 0;
    return x > 0 ? 1 : -1;
}
point operator + (point a,point b){return point(a.x+b.x,a.y+b.y);}
point operator - (point a,point b){return point(a.x-b.x,a.y-b.y);}
point operator * (point a,double k){return point(a.x*k,a.y*k);}
double Dot(point a,point b){return a.x*b.x+a.y*b.y;}
double Cro(point a,point b){return a.x*b.y-a.y*b.x;}
double dis(point a,point b){return sqrt((a.x-b.x)*(a.x-b.x) + (a.y-b.y)*(a.y-b.y));}
bool OnRight(point p,line s){return dcmp(Cro(p-s.x,s.v)) >= 0;}
point Root_LL(line a,line b)
{
    Vector v = a.v, u = b.v, w = a.x-b.x;
    double t = Cro(w,u) / Cro(u,v);
    return a.x + v * t;
}
bool comp(line a,line b)
{
    if(dcmp(a.ang-b.ang) == 0) return OnRight(b.x,a);
    return a.ang < b.ang;
}
void HPI()
{
    sort(L+1,L+cnt+1,comp);
    l = 1, r = 1, q[1] = L[1];
    for(int i = 2; i <= cnt; i++)
    {
        if(dcmp(L[i].ang-L[i-1].ang) == 0) continue;
        while(l < r && OnRight(sta[r-1],L[i])) r--;
        while(l < r && OnRight(sta[l],L[i])) l++;
        q[++r] = L[i];
        if(l < r) sta[r-1] = Root_LL(q[r],q[r-1]);
    }
    while(l < r && OnRight(sta[r-1],q[l])) r--;
    while(l < r && OnRight(sta[l],q[r])) l++;
    sta[r] = Root_LL(q[l],q[r]);
    for(int i = l; i <= r; i++) Ans[++top] = sta[i];
}
int main()
{
    scanf("%d",&n);
    for(int i = 0; i < n; i++) scanf("%lf%lf",&p[i].x,&p[i].y);
    for(int i = 1; i < n-1; i++) ans1 += Cro(p[i]-p[0],p[(i+1) % n]-p[0]);
    for(int i = 0; i < n; i++) L[++cnt] = line(p[i],p[(i+1) % n]-p[i]);
    for(int i = 1; i < n; i++)
    {
        int t = (i+1) % n;
        double A = p[0].y - p[1].y - p[i].y + p[t].y;
        double B = p[1].x - p[0].x - p[t].x + p[i].x;
        double C = p[0].x * p[1].y - p[1].x * p[0].y + p[t].x * p[i].y - p[i].x * p[t].y;
        if(dcmp(B) == 0) L[++cnt] = line(point(-C/A,0),point(-B,A));
        else L[++cnt] = line(point(0,-C/B),point(-B,A));
    }
    HPI();
    for(int i = 2; i < top; i++) ans2 += Cro(Ans[i]-Ans[1],Ans[i+1]-Ans[1]);
    printf("%.4lf\n",ans2/ans1);
    return 0;
}

查看评论