文章较长，若数学公式渲染失败，请刷新界面。
还可以到团队作业页面，下载markdown源文件，在本地渲染，以获得更好的阅读体验。
向量
平面向量
向量： 既有大小又有方向的量称为向量，记作 $\vec {a}$ 或 $\boldsymbol a$。
有向线段： 带方向的线段。用有向线段来直观地表示向量。起点为 $A$ 终点为 $B$ 的有向线段表示的向量，用符号简记为 $\overrightarrow{AB}$.
向量的模：向量的大小（或长度），用 $|\overrightarrow{AB}|$ 或 $|\boldsymbol{a}|$ 表示。
零向量：模为 $0$ 的向量。零向量的方向任意。记为：$\vec 0$ 或 $\boldsymbol{0}$。
单位向量：模为 $1$ 的向量称为该方向上的单位向量。
平行向量：方向相同或相反的两个 非零 向量。规定 $\boldsymbol 0$ 与任意向量平行。$\boldsymbol{a}$ 与 $\boldsymbol{b}$ 平行，记作：$\boldsymbol{a}\parallel \boldsymbol{b}$ 。
共线向量：与平行向量的定义相同。任一组平行向量都可以平移到同一直线上。
向量的夹角：已知两个非零向量 $\boldsymbol a,\boldsymbol b$，作 $\overrightarrow{OA}=\boldsymbol a,\overrightarrow{OB}=\boldsymbol b$，那么 $\theta=\angle AOB$ 就是向量 $\boldsymbol a$ 与向量 $\boldsymbol b$ 的夹角。记作：$\langle \boldsymbol a,\boldsymbol b\rangle$。当 $\theta=\frac{\pi}{2}$ 时，称这两个向量垂直，记作 $\boldsymbol a\perp \boldsymbol b$。规定 $\theta \in [0,\pi]$。
向量的线性运算
向量的加法
向量加法的三角形法则
对于平面上的任意两个向量 $\boldsymbol{a}$ 和 $\boldsymbol{b}$，在平面内任取一点 $A$，作 $\overrightarrow{AB}=\boldsymbol{a}$，$\overrightarrow{BC}=\boldsymbol{b}$，作向量 $\overrightarrow{AC}$.
称向量 $\overrightarrow{AC}$ 为向量 $\boldsymbol{a}$ 和 $\boldsymbol{b}$ 的 和向量，$\overrightarrow{AB}+\overrightarrow{BC}=\overrightarrow{AC}.$
如图1，把向量首尾顺次相连，向量的和为第一个向量的起点指向最后一个向量的终点；
向量加法的平行四边形法则
若要求和的两个向量 共起点，那么它们的和向量为以这两个向量为邻边的平行四边形的对角线，起点为两个向量共有的起点，方向沿平行四边形对角线方向。
向量的减法
减法可以写成加上相反数的形式，即：$a-b=a+(-b)$，如图1，$b=c-a$，$a=c-b$.
向量的数乘
给定一个实数 $\lambda$ 和一个向量 $\boldsymbol a$，规定其乘积为一个向量，记作 $\lambda \boldsymbol a$，其模与方向定义如下：
- $|\lambda \boldsymbol a|=|\lambda| |\boldsymbol a|$；
- 当 $\lambda >0$ 时，$\lambda\boldsymbol a$ 与 $\boldsymbol a$ 同向，当 $\lambda =0$ 时，$\lambda \boldsymbol a=\boldsymbol 0$，当 $\lambda<0$ 时，$\lambda \boldsymbol a$ 与 $\boldsymbol a$ 方向相反。
这种运算是向量的数乘运算。
坐标表示
\boldsymbol a + \boldsymbol b=(\boldsymbol a_x+\boldsymbol b_x,\;\boldsymbol a_y+\boldsymbol b_y)\\
\boldsymbol a – \boldsymbol b=(\boldsymbol a_x-\boldsymbol b_x,\;\boldsymbol a_y-\boldsymbol b_y)\\
\lambda \boldsymbol a=(\lambda\boldsymbol a_x,\;\lambda\boldsymbol a_y)
$$
向量的数量积
已知两个向量 $\boldsymbol a,\boldsymbol b$，它们的夹角为 $\theta$，那么：
\boldsymbol a \cdot \boldsymbol b=|\boldsymbol a||\boldsymbol b|\cos \theta=\boldsymbol a_x\boldsymbol b_x+\boldsymbol a_y\boldsymbol b_y
$$
就是这两个向量的 数量积，也叫 点积 或 内积。其中称 $|\boldsymbol a|\cos \theta$ 为 $\boldsymbol a$ 在 $\boldsymbol b$ 方向上的投影。数量积的几何意义即为：数量积 $\boldsymbol a \cdot \boldsymbol b$ 等于 $\boldsymbol a$ 的模与 $\boldsymbol b$ 在 $\boldsymbol a$ 方向上的投影的乘积。
这种运算得到的结果是一个实数，为标量。
可以方便的计算出 $\cos \theta$，于是有如下应用：
- $\boldsymbol a \perp \boldsymbol b$ $\iff$ $\boldsymbol a\cdot \boldsymbol b=0$
- $\boldsymbol a = \lambda \boldsymbol b$ $\iff$ $|\boldsymbol a\cdot \boldsymbol b|=|\boldsymbol a||\boldsymbol b|$
- $|\boldsymbol a|=\sqrt {m^2+n^2}$
- $\cos \theta=\cfrac{\boldsymbol a\cdot\boldsymbol b}{|\boldsymbol a||\boldsymbol b|}$
向量的向量积
给定两个向量 $\boldsymbol a,\boldsymbol b$，规定其向量积为一个向量，记作 $\boldsymbol a\times \boldsymbol b$，其模与方向定义如下：
- $|\boldsymbol a\times \boldsymbol b|=|\boldsymbol a||\boldsymbol b|\sin \langle \boldsymbol a,\boldsymbol b\rangle$；
- $\boldsymbol a\times \boldsymbol b$ 与 $\boldsymbol a,\boldsymbol b$ 都垂直，且 $\boldsymbol a,\boldsymbol b,\boldsymbol a\times \boldsymbol b$ 符合右手法则。
向量积也叫外积，其几何意义是：$|\boldsymbol a\times \boldsymbol b|$ 是以 $\boldsymbol a,\boldsymbol b$ 为邻边的平行四边形的面积。
向量积与 $\boldsymbol a,\boldsymbol b$ 所在平面垂直，其竖坐标为 $\boldsymbol a_x\boldsymbol b_y-\boldsymbol a_y\boldsymbol b_x$.
我们根据右手法则可以推断出 $\boldsymbol b$ 相对于 $\boldsymbol a$ 的方向，逆时针方向竖坐标为正值，反之为负值。
坐标旋转公式
若将向量 $\boldsymbol a=(x,\;y)$ 逆时针旋转 $\alpha$，得到向量 $\boldsymbol b$，则有：
\boldsymbol b=(x\cos \alpha-y\sin \alpha,\;y\cos \alpha+x\sin\alpha)
$$
可以通过三角恒等变换证明。
参考代码
const double pi = 3.1415926535898;
const double eps = 1e-8;
inline int sgn(double x) {
if(fabs(x) < eps) return 0;
else return x > 0 ? 1 : -1;
}
inline int fcmp(double x, double y) {
if(fabs(x - y) < eps) return 0;
else return x > y ? 1 : -1;
}
struct Point{
double x, y;
Point(){};
Point(double a, double b): x(a), y(b) {}
Point(Point a, Point b): x(b.x - a.x), y(b.y - a.y) {}
friend Point operator + (const Point &a, const Point &b) {
return Point(a.x + b.x, a.y + b.y);
}
friend Point operator - (const Point &a, const Point &b) {
return Point(a.x - b.x, a.y - b.y);
}
friend bool operator == (const Point &a, const Point &b) {
return fcmp(a.x, b.x) == 0 && fcmp(a.y, b.y == 0);
}
friend Point operator * (const double &a, const Point &b) {
return Point(a * b.x, a * b.y);
}
friend Point operator * (const Point &a, const double &b) {
return Point(a.x * b, a.y * b);
}
friend double operator * (const Point &a, const Point &b) {
return a.x * b.y - a.y * b.x;
}
friend double operator & (const Point &a, const Point &b) {
return a.x * b.x + a.y * b.y;
}
friend Point operator / (const Point &a, const double &b) {
return Point(a.x / b, a.y / b);
}
friend bool operator < (const Point &a, const Point &b) {
return (a.x < b.x) || (a.x == b.x && a.y < b.y);
}
inline double len() {
return sqrt(1.0 * x * x + 1.0 * y * y);
}
friend double area(Point &a, Point &b, Point &c) {
return (b - a) * (c - a) / 2.0;
}
};
typedef Point Vec;
inline double dis(Point &a, Point &b) {
return (a - b).len();
}
inline double angle(Point &a, Point &b) {
return acos((a & b) / a.len() / b.len());
}
inline Vec rotate(Vec &a, double k) {
return Vec(a.x * cos(k) - a.y * sin(k), a.x * sin(k) - a.y * cos(k));
}
平面凸包
相关概念
凸多边形：所有内角大小都在 $[0,\pi]$ 范围内的 简单多边形。
平面凸包：平面上的一个子集 $S$ 被称为凸的，当且仅当对于任意两点 $p, q\in S$，线段 $\overline{pq}$ 都完全属于 $S$.集合 $S$ 的凸包 $\mathcal{CH}(S)$，是包含 $S$ 的最小凸集，也就是包含 $S$ 的所有凸集的交。
如上图，凸包还可以理解为平面上有若干柱子，用橡皮筋套住所有柱子，绷紧后形成的多边形即为凸包。
所以有更友好的定义（不一定准确）。
凸包：在平面上能包含所有给定点的最小凸多边形叫做凸包。
暴力求解
用二维坐标 $(x_i, y_i)$ 的形式给定点集 $P$，考虑如何暴力求解。
注意到，若线段 $\overline{pq}$ 在凸包上，则 $P$ 中的点均位于直线 $pq$ 的同一侧。若我们钦定 $p\rightarrow q$ 按顺时针方向，则有更强的限制，需要 $P$ 中的点都在直线的右侧。
于是可以枚举有序点对 $(p, q)\in P\times P$，若 $P$ 中的点都在有向线段 $\overline{pq}$ 的右侧，则 $\overline{pq}$ 是 $\mathcal{CH}(P)$ 中的一条边。
需要用到向量的叉积，点 $t$ 在 $\overrightarrow {pq}$ 右侧 $\Longleftrightarrow$ $\overrightarrow{pt}\times \overrightarrow{pq}>0$.
这样的复杂度是 $\mathcal{O}(n^3)$ 的，有很多可以优化的地方。
Andrew算法
Andrew 算法是一种递增式算法，流程如下。
递增式算法 (incremental algorithm)，在计算几何中常见。算法思想：逐一加入 $P$ 中的点，每增加一个点，都更新一次目前的解，加入最后一个点后，即可得到答案。
首先把所有点 排序 ，以横坐标为第一关键字，纵坐标为第二关键字。
排序后，第一个点和末尾的点，一定在凸包上，容易通过反证法证明。
从左往右看，上下凸壳斜率的 单调性相反，即所旋转的方向不同，所以要分开求。
我们 升序枚举 求出下凸壳，然后 降序枚举 求出上凸壳，这样凸包的每条边都是向 逆时针方向 旋转的。
设当前枚举到点 $P$，即将把其加入凸包；当前栈顶的点为 $S_1$，栈中第二个点为 $S_2$.
求凸包时，若 $P$ 与 $S_1$ 构成的新线段是顺时针旋转的，即叉积满足：$\overrightarrow{S_2S_1}\times \overrightarrow{S_1P}<0$，则弹出栈顶，继续检查，直到 $\overrightarrow{S_2S_1}\times \overrightarrow{S_1P}\ge 0$ 或者栈内仅剩一个元素为止。
上图是一个弹栈的例子，$p_i$ 是新加入的点，细线是加入 $p_i$ 之前的凸包状态。
记 $n=|P|$，则时间复杂度为 $\mathcal{O}(n\log n)$，瓶颈在排序部分。
如上图，若将弹出栈顶元素的条件改为 $\overrightarrow{S_2S_1}\times \overrightarrow{S_1P}\leq 0$，同时停止条件改为 $\overrightarrow{S_2S_1}\times \overrightarrow{S_1P}> 0$，则求出的凸包中不存在三点共线。可视情况更改。
下面是参考代码。函数返回值为凸包的点数，Point ret[]的下标从 $0$ 开始。
inline bool check(Point s1, Point s2, Point p) {
return Vec(s2, s1) * Vec(s1, p) > 0;
}
int Convex_hull_2d(int n, Point *p, Point *ret) {
sort(p, p + n, cmp1);
int top = -1;
for (int i = 0; i < n; i++) {
while (top > 0 && !check(ret[top], ret[top - 1], p[i]))
top--;
ret[++top] = p[i];
}
int k = top;
for (int i = n - 2; i >= 0; i--) {
while (top > k && !check(ret[top], ret[top - 1], p[i]))
top--;
ret[++top] = p[i];
}
return top;
}
Graham算法
Andrew 算法是 Graham 算法的改进版本。在 Graham 算法中，点集按照极角序排序。
极角：任取一个顶点 $O$ 作为极点，作射线 $OX$，称为极轴。平面上一点 $p$ 的极角，即为向量 $\overrightarrow{Op}$ 与极轴 $OX$ 的夹角。一般地，取 $x$ 轴作为极轴，以逆时针方向为正。
可以利用 atan2(double y, double x) 进行极角排序。函数返回值为 $(x, y)$ 与 $x$ 轴的极角，数值 $\in(-\pi, \pi]$.
bool cmp(Point a, Point b) {
if(atan2(a.y, a.x) - atan2(b.y, b.x) == 0)
return a.x < b.x;
return atan2(a.y, a.x) < atan2(b.y, b.x);
}
另外一种方式是利用叉积排序。
bool cmp(Point a, Point b) {
return a * b > 0;
}
注意极角排序时，无论用 atan2 还是叉积，精度上都会出现不少问题，尽量避免使用这种方法。
平面凸包的周长与面积
先求出按照顺时针排序的，构成凸包的点集 $p$，记 $n=|p|$.
求周长：把相邻两点组成的向量的模长求和，即：
l=\sum_{i=1}^n|\overrightarrow{p_{i}p_{i+1}}|+|\overrightarrow{p_{1}p_{n}}|
$$
double dis(Point a, Point b) {
return (a - b).len();
}
double Convex_hull_2d_L(int n, Point *p) {
Point convex[N];
int siz = Convex_hull_2d(n, p, convex);
double ans = dis(convex[0], convex[siz - 1]);
for (int i = 1; i < siz; i++)
ans += dis(convex[i - 1], convex[i]);
return ans;
}
求面积：任取凸包内一点（一般取 $p_1$），则有：
s=\sum_{i=2}^{n-1} \operatorname{area}(p_1, p_i, p_i+1)=\sum_{i=2}^{n-1}\frac{|(p_i-p_1)\times (p_{i+1}-p_1)|}{2}
$$
double area(Point a, Point b, Point c) {
return (b - a) * (c - a) / 2.0;
}
double Convex_hull_2d_S(int n, Point *p) {
Point convex[N];
int siz = Convex_hull_2d(n, p, convex);
double ans = 0;
for (int i = 2; i < siz; i++)
ans += area(convex[0], convex[i - 1], convex[i]);
return ans;
}
动态凸包（CF70D）
维护一个点集 $S$ 的凸包，需要支持如下操作：
- 询问点 $p$ 是否在当前的凸包中，
- 向 $S$ 中添加点 $p$，
保证坐标均为整数。
分析
和 Andrew 算法一样，这里的算法按照坐标字典序排序。相对于极角排序，能够减小精度误差。
用两个 std::map<int, int>，用 $\operatorname{top}$ 记录上凸包，$\operatorname{down}$ 记录下凸包。
存储方法：若上凸包中存在横坐标为 $x$ 的点，则这个点的纵坐标为 $\operatorname{top}[x]$，$\operatorname{down}$ 同理。
询问操作
只需满足：在上凸包之下且在下凸包之上。
以上凸包为例。$i$ 在上凸包之下，当且仅当 $|\overrightarrow{ik}\times \overrightarrow{ij}|\geq 0$.
bool check_top(int x, int y) { //是否在上凸包下面
auto k = top.lower_bound(x);
if(k == top.end())
return false;
if(k -> first == x)
return y <= k->second;
if(k == top.begin()) return false;
auto j = k; j--;
return Point(k->first - x, k->second - y) *
Point(j->first - x, j->second - y) >= 0;
}
下凸包同理，$i$ 在下凸包之上，当且仅当 $|\overrightarrow{ik}\times \overrightarrow{ij}|\leq 0$
bool check_down(int x, int y) { //是否在下凸包上面
auto k = down.lower_bound(x);
if(k == down.end())
return false;
if(k -> first == x)
return y >= k->second;
if(k == down.begin()) return false;
auto j = k; j--;
return Point(k->first - x, k->second - y) *
Point(j->first - x, j->second - y) <= 0;
}
插入操作
把 $p$ 点加入凸包，上下凸包都要尝试。把加入 $p$ 点后，删掉不满足凸性的点。
如上图，这些点一定是分布在 $p_x$ 左右的连续段。
因此找到 $p$ 点在上/下凸壳中的位置，向左右分别删点，直到满足凸性。
注意迭代器的边界问题。如果已经删没了，要及时退出循环，否则会 RE。
void insert_top(int x, int y) {
if(check_top(x, y)) return;
top[x] = y;
auto it = top.find(x);
auto jt = it;
if(it != top.begin()) { //remove left
jt--;
while(remove_top(jt++)) jt--;
}
if(++jt != top.end()) { //remove right
while(remove_top(jt--)) jt++;
}
}
void insert_down(int x, int y) {
if(check_down(x, y)) return;
down[x] = y;
auto it = down.find(x);
auto jt = it;
if(it != down.begin()) { //remove left
jt--;
while(remove_down(jt++)) jt--;
}
if(++jt != down.end()) { //remove right
while(remove_down(jt--)) jt++;
}
}
下面的函数用于：判断能否删除当前点，若能删，则执行删除操作。
bool remove_top(map<int, int>::iterator it) {
if(it == top.begin()) return false; //到边界就不删了
if(++it == top.end()) return false; it--;
auto jt = it, kt = it;
jt--; kt++;
if(Point(it -> first - jt -> first, it->second - jt->second) *
Point(it -> first - kt -> first, it->second - kt->second) <= 0) {
top.erase(it);
return true;
}
return false;
}
bool remove_down(map<int, int>::iterator it) {
if(it == down.begin()) return false;
if(++it == down.end()) return false; it--;
auto jt = it, kt = it;
--jt; ++kt;
if(Point(it -> first - jt -> first, it->second - jt->second) *
Point(it -> first - kt -> first, it->second - kt->second) >= 0) {
down.erase(it);
return true;
}
return false;
}
三维凸包
三维向量类
存储：用结构体记录三维坐标，方便重载运算符。
struct Point3 {
double x, y, z;
Point3(){};
Point3(double a, double b, double c) : x(a), y(b), z(c) {}
};
加法、减法和点乘等操作：与二维向量类似。
Point3 operator + (const Point3 &b) {
return Point3(x + b.x, y + b.y, z + b.z);
}
Point3 operator - (const Point3 &b) {
return Point3(x - b.x, y - b.y, z - b.z);
}
Point3 operator * (const Point3 &b) {
return Point3(y*b.z - z*b.y, z*b.x - x*b.z, x*b.y - y*b.x);
}
bool operator == (const Point3 &b) {
return fcmp(x, b.x) == 0 && fcmp(y, b.y) == 0 && fcmp(z, b.z) == 0;
}
double len() {
return sqrt(x * x + y * y + z * z);
}
叉乘：$\boldsymbol a\times \boldsymbol b$ 的结果为一个三维向量 $\boldsymbol c$，$\boldsymbol c\perp \boldsymbol a$ 且 $\boldsymbol c\perp \boldsymbol b$，结果向量的模长为 $|\boldsymbol a||\boldsymbol b|\sin \langle \boldsymbol a,\boldsymbol b\rangle$，代表以 $\boldsymbol a$、$\boldsymbol b$ 为两边的平行四边形的面积。
在三维向量体系中，我们需要用坐标表示结果向量 $\boldsymbol{a}$，推导过程如下。（来源：叉积 – 维基百科）
右手坐标系中，基向量 $\boldsymbol{i},\boldsymbol{j},\boldsymbol{k}$ 满足以下等式：
\begin{array}{c}
\boldsymbol{i}\times\boldsymbol{j} = \boldsymbol{k}\\
\boldsymbol{j\times i} = -\boldsymbol{k}\\
\boldsymbol{j}\times\boldsymbol{k} = \boldsymbol{i}\\
\boldsymbol{k\times j} = -\boldsymbol{i}\\
\boldsymbol{k}\times\boldsymbol{i} = \boldsymbol{j}\\
\boldsymbol{i\times k} = -\boldsymbol{j}
\end{array}
$$
根据外积的定义可以得出：$\boldsymbol{i}\times\boldsymbol{i} = \boldsymbol{j}\times\boldsymbol{j} = \boldsymbol{k}\times\boldsymbol{k} = \boldsymbol{0}$.
根据以上等式，结合外积的分配律，就可以确定任意向量的外积。
任取向量 $\boldsymbol{u}=u_1\boldsymbol{i} + u_2\boldsymbol{j} + u_3\boldsymbol{k} $ 和 $\boldsymbol{v}=v_1\boldsymbol{i} + v_2\boldsymbol{j} + v_3\boldsymbol{k}$，两者的外积 $\boldsymbol{u} \times \boldsymbol{v}$ 可以根据分配率展开：
\begin{align}
\boldsymbol{u}\times\boldsymbol{v} = &(u_1\boldsymbol{i} + u_2\boldsymbol{j} + u_3\boldsymbol{k}) \times (v_1\boldsymbol{i} + v_2\boldsymbol{j} + v_3\boldsymbol{k})\\
= &u_1v_1(\boldsymbol{i} \times \boldsymbol{i}) + u_1v_2(\boldsymbol{i} \times \boldsymbol{j}) + u_1v_3(\boldsymbol{i} \times \boldsymbol{k}) + \\
&u_2v_1(\boldsymbol{j} \times \boldsymbol{i}) + u_2v_2(\boldsymbol{j} \times \boldsymbol{j}) + u_2v_3(\boldsymbol{j} \times \boldsymbol{k}) + \\
&u_3v_1(\boldsymbol{k} \times \boldsymbol{i}) + u_3v_2(\boldsymbol{k} \times \boldsymbol{j}) + u_3v_3(\boldsymbol{k} \times \boldsymbol{k})\\
\end{align}
$$
把前面的 $6$ 个等式代入，则有：
\begin{align}
\boldsymbol{u}\times\boldsymbol{v} = {} &- u_1v_1\boldsymbol{0} + u_1v_2\boldsymbol{k} – u_1v_3\boldsymbol{j} \\
&- u_2v_1\boldsymbol{k} – u_2v_2\boldsymbol{0} + u_2v_3\boldsymbol{i} \\
&+ u_3v_1\boldsymbol{j} – u_3v_2\boldsymbol{i} – u_3v_3\boldsymbol{0} \\
= {} &(u_2v_3 – u_3v_2)\boldsymbol{i} + (u_3v_1 – u_1v_3)\boldsymbol{j} + (u_1v_2 – u_2v_1)\boldsymbol{k}\\
\end{align}
$$
因此结果向量 $\boldsymbol{s} = \boldsymbol{u} \times \boldsymbol{v} = s_1\boldsymbol{i} + s_2\boldsymbol{j} + s_3\boldsymbol{k}$ 的三维坐标为：
\begin{align}
s_1 &= u_2v_3-u_3v_2\\
s_2 &= u_3v_1-u_1v_3\\
s_3 &= u_1v_2-u_2v_1
\end{align}
$$
Point3 operator * (const Point3 &b) {
return Point3(y*b.z - z*b.y, z*b.x - x*b.z, x*b.y - y*b.x);
}
平面类
用三个向量表示一个三角形的平面。一个多面体可以通过三角剖分，用若干个三角形表示。
为了节省空间，用 point3 p[N] 存储所有可能出现的向量，结构体plane只记录向量在 p[] 中的下标。
记录的三个向量按逆时针首尾相接，这样在判断方向时比较方便。
struct plane{
int v[3]; //逆时针
plane(){};
plane(int a, int b, int c) { v[0] = a, v[1] = b, v[2] = c; }
};
平面的法向量：是指垂直于该平面的三维向量。一个平面具有无限个法向量，这些法向量有两个方向。
根据叉积的性质，将三角形的两条邻边叉乘，得到的向量即为法向量。
Point3 normal() {
return (p[v[1]] - p[v[0]]) * (p[v[2]] - p[v[0]]);
}
利用法向量的模长，也可以算出三角形的面积。
double area() {
return normal().len() / 2.0;
}
三维凸包及性质
由 $n$ 个点构成的凸多面体。
性质：根据欧拉公式，任意包含 $n$ 个顶点的凸多面体，所含的边不会超过 $3n-6$ 条，所含的小平面不会超过 $2n-4$ 张。
随机增量法
算法思想
和 Andrew 算法类似，考虑每次把 $p_r$ 加入到前 $r-1$ 个点的凸包中，也就是将 $\mathcal{CH}(P_{r-1})$ 转化为 $\mathcal{CH}(P_{r})$.
第一种情况：$p_r$ 在 $\mathcal{CH}(P_{r-1})$ 的内部或边界上，则 $\mathcal{CH}(P_{r-1})\rightarrow \mathcal{CH}(P_{r})$.
第二种情况：$p_r$ 在 $\mathcal{CH}(P_{r-1})$ 外部。设想你站在 $p_r$ 所在的位置，看向 $\mathcal{CH}(P_{r-1})$. $\mathcal{CH}(P_{r-1})$ 中的某些小平面会被看到，其余在背面的平面不会被看到。如下图，从 $p_r$ 可见的平面构成了一片连通的区域。
这片区域由一条封闭折线围成，称这条线为 $p_r$ 在 $\mathcal{CH}(P_{r-1})$ 上的边界线 $(\rm{horizon})$。
根据这条地平线，我们可以判断出，在原先 $\mathcal{CH}(P_{r-1})$ 表面上的哪些部分需要被保留，哪些需要被替换。
显然，不可见的平面在 $\mathcal{CH}(P_{r})$ 中被保留，并且我们用 $p_r$ 与地平线之间连接出新的小平面，来替换所有可见的小平面，如下图。
判断平面对点的可见性
如何用几何语言表达：一个平面对 $p_r$ 是可见的？
对于一个凸包上的小平面，它能将空间分为两半，一侧是凸包外部，一侧是凸包内部。容易发现，如果点 $p$ 位于小平面的外侧，那么这个平面对于 $p$ 点就是可见的，因为凸包上的其它平面都不会有遮挡。
形式化地，记 $a,b,c$ 为平面三角形的三个顶点，从凸包外部看，三点按照逆时针排列。
利用叉乘的性质，记 $\boldsymbol{s}=\overrightarrow{ab}\times \overrightarrow{ac}$，则结果向量 $\boldsymbol{s}$ 是一个平面的法向量，且指向凸包外部。
对于空间内任意一点 $p$，若 $\overrightarrow{ap}\times \boldsymbol{s}>0$，则这个平面对点 $p$ 是可见的。
下面是定义在结构体 plane 中的函数，用于判断点 $A$ 是否位于平面的外侧。
bool is_above(Point3 A) {
return (normal() & (A - p[v[0]])) >= 0;
}
求出边界线
要想把凸包从 $\mathcal{CH}(P_{r-1})$ 转化为 $\mathcal{CH}(P_{r})$，我们需要准确地求出凸包上的哪些边在边界线上。求出边界线之后，才能用 $p$ 与边界线构成的小平面替换需要被删掉的小平面。
定义 bool g[N][N]，$g[i][j]$ 表示 $\overrightarrow{p_ip_j}$ 所在的平面是否可见。
若规定平面 $(a,b,c)$ 只包含 $\overrightarrow{ab},\overrightarrow{bc},\overrightarrow{ca}$，则对于任意有序数对 $(i, j)$，向量 $\overrightarrow{p_ip_j}$ 最多被包含在一个平面内。
注意到，位于边界线上的向量 $\overrightarrow{p_ip_j}$ 一定满足 $g[i][j]=1$ 且 $g[j][i]=0$。所以我们只需对每个平面判断其可见性，并更新在 g[][] 中对应的数值，即可求出边界线。
利用边界线更新凸包
在上一步遍历小平面时，若遇到的小平面是不可见的，则把它加入新的凸包中；若可见，则单独记录。
之后遍历所有可见的小平面，若 $\overrightarrow{p_ip_j}$ 在边界线上，则把 $(p_i, p_j, p_r)$ 加入凸包中。$p_r$ 是新加入凸包的点。这样加入后的三点也满足逆时针排列。
参考代码
函数返回值为三维凸包的平面数，plane ret[]的下标从 $0$ 开始。
int Convex_hull_3d(int n, plane *ret) {
plane tmp[N];
bool g[N][N];
for (int i = 0; i < n; i++) p[i].shake();
int top = -1;
ret[++top] = plane(0, 1, 2);
ret[++top] = plane(0, 2, 1);
for (int i = 3; i < n; i++) {
int cnt = -1;
for (int j = 0; j <= top; j++) {
bool flag = ret[j].is_above(p[i]);
if (!flag)
tmp[++cnt] = ret[j];
for (int k = 0; k < 3; k++)
g[ret[j].v[k]][ret[j].v[(k + 1) % 3]] = flag;
}
for (int j = 0; j <= top; j++) {
for (int k = 0; k < 3; k++) {
int a = ret[j].v[k], b = ret[j].v[(k + 1) % 3];
if (g[a][b] && !g[b][a])
tmp[++cnt] = plane(a, b, i);
}
}
for (int j = 0; j <= cnt; j++) ret[j] = tmp[j];
top = cnt;
}
return (top + 1);
}
旋转卡壳
概述
旋转卡壳算法用于：在线性时间内，求凸包直径、最小矩形覆盖等于凸包性质相关的问题。线性时间是指求出凸包之后的算法时间复杂度。
求凸包直径
给定平面上的 $n$ 个点，求所有点对之间的最长距离。$2\leq n\leq 50000, |x|, |y|\leq 10^4$.
首先，求出这 $n$ 个点的凸包，复杂度可以做到为 $\mathcal{O}(n\log n)$，如何求直径？
暴力做法：可以遍历每个点对，求出最大距离，复杂度为 $\mathcal{O}(n^2)$.
算法流程
可以遍历凸包上的边，对每条边 $(a, b)$，去找距离这条边最远的点 $p$。对于 $p$ 点，距离它最远的点，一定是 $a,b$ 中的一个。
我们发现，若逆时针遍历凸包上的边，那么随着边的转动，对应的最远点也在逆时针旋转。
因此，我们可以在逆时针枚举边的同时，实时维护最远点，利用单调性，复杂度为 $\mathcal{O}(n)$.
算法实现
求凸包时，若使用 Andrew 算法，则凸包上的点已经按照逆时针排序了。问题在如何判断下一个点到当前边的距离是否更大。
一种方法是用点到直线距离公式，但利用叉积的精度更高，也更方便。
struct Point{
int x, y;
//......
int sqr_len() { return x * x + y * y; }
};
inline int sqr_dis(Point a, Point b) { return (a - b).sqr_len(); }
int Get_Max(int n, Point *ch) {//传入convex-hull
int ret = 0;
ch[n] = ch[0];
int j = 1;
for(int i = 0; i < n; i++) {
while((ch[i] - ch[j+1]) * (ch[i+1] - ch[j+1]) >
(ch[i] - ch[j]) * (ch[i+1] - ch[j]))
j = (j + 1) % n;
ret = max(ret, max(sqr_dis(ch[i], ch[j]), sqr_dis(ch[i+1], ch[j])));
}
return ret;
}
最小矩形覆盖
题意
给定一些点的坐标，求能够覆盖所有点的最小面积的矩形。
求出矩形面积、顶点坐标。$3\leq n\leq 50000$.
分析
先求出凸包，之后用旋转卡壳维护三个边界点。
利用叉积和点积，可以求出矩形面积及四个顶点的坐标。
实现
下面是一种可行的表示方法。
设 $H=|P_0P_3|=\frac{|\overrightarrow{AB}\times \overrightarrow{BU}|}{|AB|}$，$L=\frac{|\overrightarrow{BA}\cdot \overrightarrow{AL}|}{|BA|}$，$R=\frac{|\overrightarrow{AB}\cdot \overrightarrow{BR}|}{|AB|}$.
则矩形面积为：
S=H\times (L+|\overrightarrow{AB}|+R)
$$
顶点坐标为：
\begin{array}{l}
&\overrightarrow{P_0}&=&\overrightarrow A+L\times \frac{\overrightarrow{BA}}{|\overrightarrow{BA}|}\\
&\overrightarrow{P_1}&=&\overrightarrow B+R\times \frac{\overrightarrow{AB}}{|\overrightarrow{AB}|}\\
&\overrightarrow{P_2}&=&\overrightarrow P_1+H\times \frac{\overrightarrow{P_1R}}{|\overrightarrow{P_1R}|}\\
&\overrightarrow{P_3}&=&\overrightarrow P_0+H\times \frac{\overrightarrow{P_0L}}{|\overrightarrow{P_0L}|}
\end{array}
$$
若不需要求出顶点坐标，也可以这样表示矩形面积：
S=\frac{|\overrightarrow{AB}\times \overrightarrow{BU}|\times (|\overrightarrow{AD}\cdot\overrightarrow{AB}|+|\overrightarrow{BC}\cdot\overrightarrow{BA}|-|\overrightarrow{AB}\cdot\overrightarrow{BA}|)}{|\overrightarrow{AB}\cdot\overrightarrow{BA}|}
$$
设 $|\overrightarrow{P_0B}|=m$，$|\overrightarrow{AP_1}|=n$ 后易证。
代码
double Get_Max(int n, Point *ch) {
ch[n] = ch[0];
int u = 2, l, r = 2;
//u是距离AB最远的点；在AB为底时，l和r是两个最靠边的点
double ret = 1e100, H, L, R, S;
for(int i = 0; i < n; i++) {
Point A = ch[i], B = ch[i+1]; Vec AB = B - A, BA = A - B;
while((AB * Vec(B, ch[u+1])) >= (AB * Vec(B, ch[u])))
u = (u + 1) % n;
while((AB & Vec(B, ch[r+1])) >= (AB & Vec(B, ch[r])))
r = (r + 1) % n;
if(i == 0) l = r;
while((AB & Vec(B, ch[l+1])) <= (AB & Vec(B, ch[l])))
l = (l + 1) % n;
H = (AB * Vec(B, ch[u])) / AB.len(); //以AB所在直线为底边，矩形的高
L = (BA & Vec(A, ch[l])) / BA.len(); //A距离左侧顶点的距离
R = (AB & Vec(B, ch[r])) / AB.len(); //B距离右侧顶点的距离
S = H * (L + AB.len() + R); //矩形面积
if(S < ret) { //求矩形顶点坐标
ret = S;
P[0] = A + L * BA.unit();
P[1] = B + R * AB.unit();
P[2] = P[1] + H * (ch[r]-P[1]).unit();
P[3] = P[0] + H * (ch[l]-P[0]).unit();
}
}
return ret;
}
POJ2079 – Triangle
题意
给定平面上的 $n$ 个点，从其中任选三个点作为顶点，求能构成的最大三角形面积。
$1\leq n\leq 50000$.
分析
显然，三角形的顶点一定都在这 $n$ 个点的凸包上，所以先求出凸包。
考虑旋转卡壳。
这题中不能固定一条边，再枚举另外两个点，因为三角形的边不一定在凸包上。
因此，先逆时针枚举一个固定点 $i$，再逆时针旋转另外两个顶点 $j$ 和 $k$.·
由于凸包上的一些单峰性，我们旋转 $j$ 时停止的条件是 $S_{\triangle ijk}$ 最大，停止后固定 $j$，以相同停止条件旋转 $k$.
$S_{\triangle ijk}$ 最大，是指 $S_{\triangle i,j\operatorname{\_last},k}<s_{\triangle ijk}$=”” 且=”” $s_{\triangle=”” ijk}=””>S_{\triangle i,j\operatorname{\_next},k}$.</s_{\triangle>
代码
注意特判 $n\leq 2$ 以及凸包大小 $\operatorname{siz}\leq 2$ 的情况。用叉积求面积。
double Get_Max(int n, Point *ch) {
double ret = 0;
ch[n] = ch[0];
int j = 1, k = 2;
for(int i = 0; i < n; i++) {
while(Vec(ch[i], ch[j]) * Vec(ch[i], ch[k]) <
Vec(ch[i], ch[j]) * Vec(ch[i], ch[k+1]))
k = (k + 1) % n;
while(Vec(ch[i], ch[j]) * Vec(ch[i], ch[k]) <
Vec(ch[i], ch[j+1]) * Vec(ch[i], ch[k]))
j = (j + 1) % n;
ret = max(ret, Vec(ch[i], ch[j]) * Vec(ch[i], ch[k]));
}
return ret;
}
半平面交
半平面
半平面是一条直线和直线的一侧，是一个点集。
当点集包含直线时，称为闭半平面；当不包含直线时，称为开半平面。
直线类
在计算几何中，用 $(\boldsymbol s, \boldsymbol t)$ 表示直线 ${st}$，一般统一保留向量的左侧半平面。
求两直线的交
若交点存在，只需求出 $t=\dfrac{|A_sP|}{|A_sA_t|}$，则 $\overrightarrow P=\overrightarrow{A_s}+t\cdot \overrightarrow{A_sA_t}$.
注意到 $t=\dfrac{|A_sP|}{|A_sA_t|}=\dfrac{|\overrightarrow{B_sB_t}\times \overrightarrow{B_sA_s}|}{|\overrightarrow{B_sB_t}\times \overrightarrow{A_sA_t}|}$，证明也不难。
如上图，把 $\overrightarrow{A_sA_t}$ 平移到 $\overrightarrow{B_sC}$ 的位置，则可以把叉乘的模看作平行四边形面积。
记 $S_1=|\overrightarrow{B_sB_t}\times \overrightarrow{B_sA_s}|$，$S_2=|\overrightarrow{B_sB_t}\times \overrightarrow{A_sA_t}|$，由于两个平行四边形同底，所以有 $\dfrac{S_1}{S_2}=\dfrac{|A_sM|}{|CN|}$.
又因为 $\triangle A_sMP\sim \triangle CNB_s$，所以 $\dfrac{S_1}{S_2}=\dfrac{|A_sM|}{|CN|}=\dfrac{|A_sP|}{|CB_s|}=\dfrac{|A_sP|}{|A_sA_t|}$.
参考代码
struct Line{
Point s, t;
Line() {};
Line(Point a, Point b) : s(a), t(b) {}
double ang() { return atan2((t - s).y, (t - s).x); };
Line(double a, double b, double c) { //ax + by + c = 0
if(sgn(a) == 0) s = Point(0, -c/b), t = Point(1, -c/b);
else if(sgn(b) == 0) s = Point(-c/a, 0), t = Point(-c/a, 1);
else s = Point(0, -c/b), t = Point(1, (-c-a)/b);
}
friend bool parallel(const Line &A, const Line &B) {
return sgn((A.s - A.t) * (B.s - B.t)) == 0;
}
friend bool Calc_intersection(const Line &A, const Line &B, Point &res) {
if(parallel(A, B)) return false;
double s1 = (B.t - B.s) * (B.s - A.s);
double s2 = (B.t - B.s) * (A.t - A.s);
res = A.s + (A.t - A.s) * (s1 / s2);
return true;
}
};
半平面交
半平面交是多个半平面的交集。
半平面交是一个点集，并且是一个凸集。在直角坐标系上是一个区域。
半平面交在代数意义下，是若干个线性约束条件，每个约束条件形如：
a_ix+b_iy\leq c_i
$$
其中 $a_i, b_i, c_i$ 为常数，且 $a_i$ 和 $b_i$ 不都为零。
半平面交有下面五种可能情况，每个半平面位于边界直线的阴影一侧。
$\rm{(iii)}$ 和 $\rm{(v)}$ 比较特殊，我们一般假设产生的交集总是有界或空的。
性质
注意到，位于半平面交边界上的任何一个点，必定来自某张半平面的边界。
并且，因为半平面交是凸集，所以每条边界线最多只能为边界贡献一条边。
因此：由 $n$ 张半平面相交而成的凸多边形，其边界最多由 $n$ 条边围成。
Sort-and-Incremental Algorithm
S&I 算法，排序增量法。
算法思想
S&I 算法利用了性质：半平面交是凸集。
因为半平面交是凸集，所以我们维护凸壳。
若我们把半平面按极角排序，那么在过程中，只可能删除队首或队尾的元素，因此使用双端队列维护。
下面是一个例子。
TODO.
算法流程
- 把半平面按照极角序排序，需要 $\mathcal{O}(n\log n)$ 的时间。
- 对每个半平面，执行一次增量过程，每次根据需要弹出双端队列的头部或尾部元素。这是线性的，因为每个半平面只会被增加或删除一次。
- 最后，通过双端队列中的半平面，在线性时间内求出半平面交（一个凸多边形，用若干顶点描述）。
这样，我们得到了一个时间复杂度为 $\mathcal{O}(n\log n)$ 的算法，瓶颈在于排序。因此，若题目给定的半平面满足极角序，则我们可以在线性的时间内求出半平面交。
极角排序
注意判断不要写成 <= 或 >=。
inline bool cmp(Line A, Line B) {
//极角相等时，位置靠右的排在前面
if(!sgn(A.ang - B.ang)) return (A.t - A.s) * (B.t - A.s) > 0;
return A.ang < B.ang;
}
求半平面交
bool Halfplane_intersection(int n, Line *hp, Point *p) {
if(n < 3) return false;
sort(hp, hp + n, cmp);
Halfplane_unique(n, hp);
st = 0; ed = 1;
que[0] = 0; que[1] = 1;
if(parallel(hp[0], hp[1])) return false;
Calc_intersection(hp[0], hp[1], p[1]);
for(int i = 2; i < n; i++) {
while(st < ed &&
sgn((hp[i].t - hp[i].s) * (p[ed] - hp[i].s)) < 0)
ed--;
while(st < ed &&
sgn((hp[i].t - hp[i].s) * (p[st + 1] - hp[i].s)) < 0)
st++;
que[++ed] = i;
assert(ed >= 1);
if(parallel(hp[i], hp[que[ed - 1]])) return false;
Calc_intersection(hp[i], hp[que[ed - 1]], p[ed]);
}
while(st < ed &&
sgn((hp[que[st]].t - hp[que[st]].s) * (p[ed] - hp[que[st]].s)) < 0)
ed--;
while(st < ed &&
sgn((hp[que[ed]].t - hp[que[ed]].s) * (p[st + 1] - hp[que[ed]].s)) < 0)
st++;
if(st + 1 >= ed) return false;
return true;
}
求凸多边形的顶点
int Get_convex_hull(Line *hp, Point *p, Point *ch) {
Calc_intersection(hp[que[st]], hp[que[ed]], p[st]);
for(int i = 0, j = st; j <= ed; i++, j++) ch[i] = p[j];
return ed - st + 1;
}
计算面积
double Calc_area(int n, Point *ch) {
double ans = 0;
for (int i = 2; i < n; i++)
ans += area(ch[0], ch[i - 1], ch[i]);
return ans;
}
P3256 [JLOI2013]赛车
题目描述
赛场上一共有 $n$ 辆车，分别称为 $g_1,g_2,…,g_n$.
赛道是一条无限长的直线。最初，$g_i$ 位于距离起跑线前进 $x_i$ 的位置。比赛开始后，车辆 $g_i$ 将会以 $v_i$ 单位每秒的恒定速度行驶。
过程中，如果一辆赛车曾经处于领跑位置的话，即没有其他的赛车跑在他的前面，这辆赛车最后就可以得奖。
求哪些赛车将会得奖。输出赛车编号。
$1\leq n\leq 10^4$，$0\leq x_i\leq 10^9$，$0\leq v_i\leq 10^9$.
分析
可以想到，每辆车的 $x-t$ 图像都是一条直线。在平面直角坐标系上画出这些直线。
若取每条直线的左侧半平面，再与 $x$ 轴上侧平面、$y$ 轴右侧平面一起求交，则所有凸多边形上的直线可以获奖。
直接跑排序增量法即可，时间复杂度 $\mathcal{O}(n\log n)$.
注意到这题 $x_i$ 和 $v_i$ 的值域均为 $10^9$，所以运算过程中的数值会达到 $10^9\times 10^9=10^{18}$，
因此需要使用 long double，并且需要 const long double eps = 1e-18.
代码
这道题不能直接去重，因为需要输出获奖车辆的编号，所以用 std::map 维护。
同时，还需要在 Line 类中存储所有重合直线的编号，用 std::vector 维护。
struct Line{
Point s, t;
vector<int> id;
Line() {};
Line(Point a, Point b) : s(a), t(b) {}
Line(Point a, Point b, vector<int> c) : s(a), t(b), id(c) {}
//...
};
//...
signed main() {
int n, x[MAXN], v[MAXN], N = 0;
map<pair<int, int>, vector<int> > mp;
cin >> n;
for(int i = 1; i <= n; i++) cin >> x[i];
for(int i = 1; i <= n; i++) cin >> v[i];
for(int i = 1; i <= n; i++) mp[make_pair(x[i], v[i])].push_back(i);
L[N++] = Line(Point(0, 1), Point(0, 0));
L[N++] = Line(Point(0, 0), Point(1, 0));
for(auto it = mp.begin(); it != mp.end(); it++) {
pair<int, int> tmp = it -> first;
L[N++] = Line(Point(0, tmp.first),
Point(1, tmp.first+tmp.second), it -> second);
}
Halfplane_intersection(N, L, p);
int ans[MAXN], cnt_ans = 0;
for(int i = st + 1; i <= ed; i++) {
for(int j = 0; j < L[que[i]].id.size(); j++) {
ans[cnt_ans++] = L[que[i]].id[j];
}
}
sort(ans, ans + cnt_ans);
cout << cnt_ans << endl;
for(int i = 0; i < cnt_ans; i++) cout << ans[i] << " "; cout << endl;
return 0;
}
P4250 [SCOI2015]小凸想跑步
题目描述
小凸晚上喜欢到操场跑步，今天他跑完两圈之后，他玩起了这样一个游戏。
操场是个凸 $n$ 边形， $n$ 个顶点按照逆时针从 $0$ ∼ $n – 1$ 编号。
现在小凸随机站在操场中的某个位置，标记为 $p$ 点。将 $p$ 点与 $n$ 个顶点各连一条边，形成 $n$ 个三角形。
如果这时 $p$ 点， $0$ 号点， $1$ 号点形成的三角形的面积是 $n$ 个三角形中最小的一个，则认为这是一次正确站位。
现在小凸想知道他一次站位正确的概率是多少。
分析
推一点式子。
设 $\overrightarrow A=(x_a,y_a)$，$\overrightarrow B=(x_b, y_b)$，$\overrightarrow C=(x_c, y_c)$，$\overrightarrow D=(x_d, y_d)$，$\overrightarrow P=(x, y)$.
由题意，
\begin{array}{rl}
|\overrightarrow {AB}\times \overrightarrow{AP}|&<&|\overrightarrow {CD}\times \overrightarrow{CP}|\\
(x_b-x_a)(y-y_a) – (y_b-y_a)(x-x_a) &<& (x_d-x_c)(y-y_c)-(y_d-y_c)(x-x_c)\\
\end{array}
$$
展开得，
\begin{array}{rl}
\text{左式}&=&(x_b-x_a)y-x_by_a+x_ay_a-(y_b-y_a)x+x_ay_b-x_ay_a\\
\text{右式}&=&(x_d-x_c)y-x_dy_c+x_cy_c-(y_d-y_c)x+x_cy_d-x_cy_c
\end{array}
$$
合并同类项，
\begin{array}{l}
&(x_b – x_a)(y – y_a) – (y_b – y_a)(x – x_a) < (x_d – x_c)(y – y_c) – (y_d – y_c)(x – x_c) \\
\Rightarrow & (x_b – x_a + x_c – x_d)y + (y_a – y_b + y_d – y_c)x + (y_b x_a – x_b y_a – y_d x_c + x_d y_c) < 0
\end{array}
$$
注意到，这是直线解析式的形式，于是转化为了半平面交问题。
参考资料
https://cp-algorithms.com/geometry/halfplane-intersection.html