寄算几何

Pre-words:

计算几何，即利用计算机建立数学模型解决几何问题。

非常开心能成为计算几何班级的课代表之一 虽然当时睡过了没能现场领头起立QWQ 。一只老🕊表示这次一定不咕咕，必定每次好好 做题 补题！

曾经浅尝过计算几何（大概是小破站上 jiry 和 dls 的老教程¹，但是当时也没学懂。感觉这个板块是可做的，因此报了这个班² （公号买的hhh） 这次一定！

概论

要解决的问题类型

1. 判定类 Determine
2. 计数类 Count
3. 列举类 + 构造类 Enumerate + Construct

题目特点

• 码量大、细节多、正确率低
• 团队核心¹

参考资料

• OI-WIKI/geometry
• 学堂在线 - 邓俊辉 计算几何
• Codeforces - catalog

工具

• cs academy - geometry_widget

btw，用来绘制树、图也很好用的 ： cs academy - graph_editor

• geogebra - calculator

前置知识 - 浮点数

精度

因为无限不循环小数无法准确表达，因而会出现精度问题。十进制照样出错，比如：

所以在什么位置上截断是重要的。C++ 中三种类型的浮点数中 float 的位数为 $3232$，(long) double 的位数为 $6464$。而其精度 随着名字的长度而 递增。

特殊值

$\pm 0.0\backslash pm0.0$	$\pm \inf \backslash pm\backslash inf$	$\mathrm{n}\mathrm{a}\mathrm{n}\backslash rm\; nan$ (not a number)
两者等值，但输出形式不同（符号位）	$1.0\mathrm{/}\pm 0.01.0\; /\; \backslash pm\; 0.0$ 可得到	$\mathrm{n}\mathrm{a}\mathrm{n}\mathrm{\ne }\mathrm{x}(\mathrm{x}\mathrm{不}\mathrm{为}\mathrm{n}\mathrm{a}\mathrm{n})\backslash rm\; nan\; \backslash ne\; x\backslash ,(x\; 不为\; nan)$ 为 true 其他比较都为 false

注意

• 若问题可以使用整数解决则不用浮点数。
• no float! (狗头) 视情况使用 long double
• 尽量减少数学函数的使用 (尤其是三角函数)
• 比较时加入 eps ( 即 $ϵ\backslash epsilon$ )

using db = double;
static const db eps = 1e-6;
int sign(db k) {
  return (k > eps) ? 1 : (k < -eps) ? -1 : 0;
}
int cmp(db k1, db k2) { return sign(k1 - k2); }

板子

不要盲目使用他人的板子

• https://xcpcio.com/code-library/code-library-list/ sto 葫芦 orz
• jiry 的板子（在一些XCPC群里搜索一个叫做 mygeo.cpp的程式）
• 构建自己的板子

第一讲：用代码表达(二维)几何

点 / 向量

• What is wasted?

using point = std::complex<db>;
using point = std::pair<db, db>;

• What is nice?

struct point {
  db x, y;
  // some functions ...
};

向量

• point 即可

点积

是个标量

几何意义：向量在另一向量上的投影有向长度。

用途：判垂直、求夹角、求向量的模。

db dot(point k1, point k2) { return k1.x * k2.x + k1.y * k2.y; }

叉积

是个向量。

几何意义：两向量围成的平行四边形的有向面积（不妨看作正弦定理的向量形式）。

用途：判平行、结合点积可判定三点的位置关系³。

db cross(point k1, point k2) { return k1.x * k2.y - k1.y * k2.x; }

• 点积的正负可以确定是在 y 轴左右（如果点积是正，那么在 y 轴右侧）
• 叉积的正负可以确定是在 x 轴上下（如果叉积是正（角度为正，即顺时针方向），那么在 x 轴上方）to-left-test

可以看到，单纯的某一运算都只能确定半平面。结合使用就可以确定具体在哪个象限。

// k1 k2 k3 逆时针 1 顺时针 -1 否则 0
// toLeft
int clockwise(point k1, point k2, point k3) {
	return sign(cross(k2 - k1, k3 - k1));
}

旋转 (逆时针)

• 请看这篇博文 vector-formula

point turn(db k1) {
  return (point){x * std::cos(k1) - y * std::sin(k1),
                 x * std::sin(k1) + y * std::cos(k1)};
}
point turn90() { return (point){-y, x}; }

线段

一个点 + 一个向量就够了。

struct line {
  // p[0]->p[1]
  point p[2];
  // some functions ...
};

判点在线段上

$PP$ 为点， $\overrightarrow{\mathrm{A}\mathrm{B}}\backslash rm\backslash vec\{AB\}$ 为线段。

• $\overrightarrow{\mathrm{P}\mathrm{A}}\times \overrightarrow{\mathrm{P}\mathrm{B}}=\vec{0}\backslash rm\; \backslash vec\{PA\}\backslash times\backslash vec\{PB\}\; =\; \backslash vec0$ 这是在判平行。
• $\overrightarrow{\mathrm{P}\mathrm{A}}\cdot \overrightarrow{\mathrm{P}\mathrm{B}}\le 0\backslash rm\; \backslash vec\{PA\}\; \backslash cdot\; \backslash vec\{PB\}\; \backslash le\; 0$ 从 $PP$ 向端点引向量，方向相反。也可以判定横纵坐标是否夹在 $ABAB$ 之间。

bool onS(point a1, point a2, point p) {
  if (sign(p.x - a1.x) == 0 && sign(p.y - a1.y) == 0) return true;
  if (sign(p.x - a2.x) == 0 && sign(p.y - a2.y) == 0) return true;
  return sign(dot(a1 - p, a2 - p)) <= 0;
}

两线段判交⁴

注意不是直线判交。

先检查一下横纵坐标是否交叉出现（防止误判三四点共线）。

随后跨立实验：$ABAB$ 分居 $CDCD$ 两侧，反之亦然。

int intersect(db l1, db r1, db l2, db r2) {
  if (l1 > r1) std::swap(l1, r1);
  if (l2 > r2) std::swap(l2, r2);
  return ~ cmp(r1, l2) && ~ cmp(r2, l1);
}
int checkSS(point k1, point k2, point k3, point k4) {
  return intersect(k1.x, k2.x, k3.x, k4.x) &&
         intersect(k1.y, k2.y, k3.y, k4.y) &&
         sign(cross(k3 - k1, k4 - k1)) * sign(cross(k3 - k2, k4 - k2)) <= 0 &&
         sign(cross(k1 - k3, k2 - k3)) * sign(cross(k1 - k4, k2 - k4)) <= 0;
}

射线交线段

这算是... 课堂小作业？

和上面的思路相同，来想一想如果射线与线段有交，会满足哪些条件？

• 交点应当与这条射线同侧。
• 交点应当在该线段的两端点之间。

这样的话，大概就是：

// 1 表示相交, -1 表示重合(或反向), 0 表示平行
int checkSH(point k1, point k2, point k3, point k4) {
  point p1 = k2 - k1, p2 = k3 - k4, p3 = k3 - k1;
  db d1 = cross(p1, p2), d2 = cross(p3, p2), d3 = cross(p1, p3);
  if (cmp(d1, 0) > 0) {
    // return k1 + p1 * (d2 / d1);
    return 1;
  } else if (!cmp(d2, 0) && !cmp(d3, 0)) {
    return -1;
  } else return 0;
}

直线

依然可以使用：

struct line {
  // p[0]->p[1]
  point p[2];
  // some functions ...
};

点到直线的距离

$PP$ 为点， $\overrightarrow{\mathrm{A}\mathrm{B}}\backslash rm\backslash vec\{AB\}$ 为直线。从 $PP$ 向 $\overrightarrow{AB}\backslash vec\{AB\}$ 引出两条边，以此构建平行四边形。这个四边形的面积为$PA\mathrm{、}PBPA、PB$ 的叉积的绝对值，同时也是距离 $\times AB\backslash times\; AB$ 的模。因此：

$\mathrm{d}\mathrm{i}\mathrm{s}(\mathrm{P},\mathrm{A}\mathrm{B})={\displaystyle \frac{{\mathrm{S}}_{\mathrm{P}\mathrm{A}\mathrm{B}{\mathrm{P}}^{\mathrm{\prime }}}}{\mathrm{\mid }\overrightarrow{\mathrm{A}\mathrm{B}}\mathrm{\mid }}}={\displaystyle \frac{\mathrm{\mid }\overrightarrow{\mathrm{P}\mathrm{A}}\times \overrightarrow{\mathrm{P}\mathrm{B}}\mathrm{\mid }}{\mathrm{\mid }\overrightarrow{\mathrm{A}\mathrm{B}}\mathrm{\mid }}}\backslash rm\; dis(P,\; AB)\; =\; \backslash dfrac\{S\_\{PABP\text{'}\}\}\{|\backslash vec\{AB\}|\}\; =\; \backslash dfrac\{|\backslash vec\{PA\}\backslash times\; \backslash vec\{PB\}|\}\{|\backslash vec\{AB\}|\}$

db disLP(point k1, point k2, point p) {
  return std::abs(cross(p - k1, p - k2)) / point(k2 - k1).abs();
}

这段代码中含有一些依赖：

point point::operator- (const point &k1) const {
  return (point){x - k1.x, y - k1.y};
}
db point::abs() { return sqrt(x * x + y * y); }

点线投影⁵

而 $\mathrm{\mid }{K}_1{t}^{\mathrm{\prime }}\mathrm{\mid }|K\_1t\text{'}|$ 可以使用点积 ${\displaystyle \frac{\overrightarrow{K_1{K}_2}\cdot \overrightarrow{K_{1}t}}{\mathrm{\mid }\overrightarrow{K_1{K}_2}\mathrm{\mid }}}\backslash dfrac\{\backslash vec\{K\_1K\_2\}\backslash cdot\; \backslash vec\{K\_1t\}\}\{|\backslash vec\{K\_1K\_2\}|\}$ 求得。而对于单位化，可以使用向量除去其模长。

因此最终的公式就是：$\overrightarrow{O{t}^{\mathrm{\prime }}}=K_1+{\displaystyle \frac{\overrightarrow{K_1{K}_2}}{\mathrm{\mid }\overrightarrow{K_1{K}_2}\mathrm{\mid }}}\times {\displaystyle \frac{\overrightarrow{K_1{K}_2}\cdot \overrightarrow{K_{1}t}}{\mathrm{\mid }\overrightarrow{K_1{K}_2}\mathrm{\mid }}}\backslash vec\{Ot\text{'}\}\; =\; K\_1\; +\; \backslash dfrac\{\backslash vec\{K\_1K\_2\}\}\{|\backslash vec\{K\_1K\_2\}|\}\backslash times\backslash dfrac\{\backslash vec\{K\_1K\_2\}\backslash cdot\; \backslash vec\{K\_1t\}\}\{|\backslash vec\{K\_1K\_2\}|\}$

point proj(point k1, point k2, point q) {
	// q 到直线 k1, k2 的投影
	point k = k2 - k1;
	return k1 + k * (dot(q - k1, k) / k.abs2());
}

这段代码中有一些依赖：

point point::operator+ (const point &k1) const {
  return (point){k1.x + x, k1.y + y};
}
db point::abs2() { return x * x + y * y; }

线段与直线之交

point getLS(point k1, point k2, point k3, point k4) {
  return k1 + k2 * (cross(k4, k1 - k3) / cross(k2, k4));
}

两线之交⁶

画的好像立体图形啊 草

* 注：JJ 使用正弦定理推导，详情可以看视频，这里我用相似来推：

设欲求的交点为 $XX$，$OO$ 为原点。由相似关系，推出：

因而

point getLL(point k1, point k2, point k3, point k4) {
	db w1 = cross(k4 - k3, k1 - k3), w2 = cross(k2 - k3, k4 - k3);
	return (k1 * w2 + k2 * w1) / (w1 + w2);
}

这段代码中含有一些依赖：

point point::operator* (db k1) const {
  return (point){x * k1, y * k1};
}
point point::operator/ (db k1) const {
  return (point){x / k1, y / k1};
}

多边形

• 逆时针存顶点 $\mathrm{n}\mathrm{e}\mathrm{x}\mathrm{t}(\mathrm{A}\mathrm{.}\mathrm{b}\mathrm{a}\mathrm{c}\mathrm{k}())=\mathrm{A}\mathrm{.}\mathrm{f}\mathrm{r}\mathrm{o}\mathrm{n}\mathrm{t}()\backslash rm\; next\backslash big(A.back()\backslash big)\; =\; A.front()$，因而不妨 a.push_back(a.front())
• 不一定凸。

using polygon = std::vector<point>;

简单多边形的面积

是那道 HDU 典

实际上是一堆有向三角形之和，用叉积可以得到一堆有向平行四边形之和。

需要以逆时针序进行计算，否则将会得到负号。

db area(polygon A) {
  db ans = 0;
  for (int i = 0; i < A.size(); i++)
    ans += cross(A[i], A[(i + 1) % A.size()]);
  return ans / 2;
}

点与简单多边形的关系

算法一:光线投射算法 Crossing Number

从该点向 $xx$ 轴引出一条射线，看撞击的次数。若为偶数，则在多边形外，反之则反。

其实很直观：每碰撞一次简单多边形的边，内外状态就会变化一次。

• 问题1：射线或与边重合。
• 问题2：射线或碰到顶点。

解决方案：仅关注纵坐标的变化，同时将相邻两点标记为左闭右开（撞到左边才记作有效），这样就可以避免上面两情形影响算法。

// 2 内部 1 边界 0 外部
int contain(polygon A, point q) {
  int pd = 0;
  A.push_back(A.front());
  for (int i = 1; i < A.size(); i++) {
    point u = A[i - 1], v = A[i];
    if (onS(u, v, q)) return 1;
    if (cmp(u.y, v.y) > 0) std::swap(u, v);
    if (cmp(u.y, q.y) >= 0 || cmp(v.y, q.y) < 0) continue;
    if (sign(cross(u - v, q - v)) < 0) pd ^= 1;
  }
  return pd << 1;
}

这段代码中含有一些依赖：

// k3 在 [k1, k2] 内
int inmid(db k1, db k2, db k3) {
  return sign(k1 - k3) * sign(k2 - k3) <= 0;
}
int inmid(point k1, point k2, point k3) {
  return inmid(k1.x, k2.x, k3.x) && inmid(k1.y, k2.y, k3.y);
}
int onS(point k1, point k2, point q) {
  return inmid(k1, k2, q) && sign(cross(k1 - q, k2 - k1)) == 0;
}

算法二:回转数法 + Sunday's algo Winding Number

面内闭合曲线沿逆时针绕过该点的总次数称作回转数。若回转数为 $00$​，则点在曲线外，反之则反。

多边形显然是闭合曲线 （不过这种东西顶点处的曲率是怎么定义的）。

实现方式：计算相邻两边夹角之和，若为 $2\pi 2\backslash pi$ 则在曲线内，反之则反。

• 问题：使用 arcsin(theta) 会使得精度忍耐能力较差。

解决方案：是一种称作 Sunday's algo 的技巧：不过话说我并没有找到这个叫做Sunday的算法

如果沿着逆时针经过了这条边，windings ++，否则 windings --。最后判定 windings ：

还在调试中 qwq，下面是JJls的代码。

std::pair<bool, int> winding(const Point &a) const {
  int cnt = 0;
  for (size_t i = 0; i < p.size(); i++) {
    Point u = p[i], v = p[nxt(i)];
    if (std::abs((a - u) ^ (a - v)) <= eps && (a - u) * (a - v) <= eps)
      return {true, 0};
    if (abs(u.y - v.y) <= eps) continue;
    Line uv = {u, v - u};
    if (u.y < v.y - eps && uv.toleft(a) <= 0) continue;
    if (u.y > v.y + eps && uv.toleft(a) >= 0) continue;
    if (u.y < a.y - eps && v.y >= a.y - eps) cnt++;
    if (u.y >= a.y - eps && v.y < a.y - eps) cnt--;
  }
  return {false, cnt};
}

THE - END

题单链接

• 牛客竞赛计算几何专题班二维基础

解题笔记

禁止禁止咕咕咕