计算几何基础算法

计算几何常用函数

结构体与向量

struct Point {
    double x, y;
    Point(double x = 0, double y = 0) : x(x), y(y) {}
};

typedef Point Vector;

Vector operator+ (Vector A, Vector B) {
    return Vector(A.x+B.x, A.y+B.y);
}
Vector operator- (Vector A, Vector B) {
    return Vector(A.x-B.x, A.y-B.y);
}
Vector operator* (Vector A, double p) {
    return Vector(A.x*p, A.y*p);
}
Vector operator/ (Vector A, double p) {
    return Vector(A.x/p, A.y/p);
}
bool operator< (const Point &a, const Point &b) {
    return a.x < b.x || (a.x == b.x && a.y < b.y);
}

const double eps = 1e-8;
int dcmp(double x) {
    if (fabs(x) < eps) return 0;
    else return x < 0 ? -1 : 1;
}

bool operator== (const Point &a, const Point &b) {
    return dcmp(a.x-b.x) == 0 && dcmp(a.y-b.y) == 0;
}

向量运算：点积

double Dot(Vector A, Vector B) {
    return A.x*B.x + A.y*B.y;
}
double Length(Vector A) {
    return sqrt(Dot(A, A));
}
double Angle(Vector A, Vector B) {
    return acos(Dot(A, B) / Length(A) / Length(B));
}

其中角度计算是根据余弦定理，设 $\displaystyle\vec{a}, \displaystyle\vec{b}$ 的夹角为 $\theta$

$$\cos{\theta} = \frac{\vec{a} \cdot \vec{b}}{|\vec{a}|\cdot |\vec{b}|}$$

向量运算：叉积
$\bm{A} \times \bm{B}$ 的几何意义，是向量 $\bm{A}, \bm{B}$ 围成的平行四边形的有向面积
$\text{Cross}(\bm{v}, \bm{w})$，正方向定义为从 $\bm{v}$ 向 $\bm{w}$ 旋转

double Cross(Vector A, Vector B) {
    return A.x*B.y - A.y*B.x;
}
double Area2(Point A, Point B, Point C) {
    return Cross(B-A, C-A);
}

极角，极角排序

• 根据 STL 函数求极角，对于向量 $\bm{A}_i(x, y)$，极角为 $\text{atan2}(y, x)$，单位为弧度
  逆时针为正方向，值域为 $[-\pi, +\pi]$，保存在数组 $\text{ang}[i]$ 中
  特别地，很多问题需要对极角排序，然后用双指针计算向量 $\bm{A}_i, \bm{A}_j$ 的夹角 $\theta$
  由于极角往往是环形的，处理完 $\bm{A}(i \to j)$ 的夹角，往往还需要处理 $\bm{A}(j \to i)$ 的夹角
  比如 $i = n, j = 1$，求第 $n$ 个向量往第 $1$ 个向量的转角
  这里采用拆换成链并复制的方法，令 $\text{ang}[1\cdots n]$ 每个元素加上 $2\pi$ 之后
  放入 $\text{ang}[n+1 \cdots n+n]$ 中
• 此外还可以根据叉积来对极角排序
  $\text{Cross}(A, B) = 0$ 向量重合，$\text{Cross}(A, B) > 0$ 向量 $B$ 在 $A$ 上方
  $\text{Cross}(A, B) < 0$ 向量 $B$ 在 $A$ 下方

  double compare(Point A, Point B, Point C) {
      return Cross(B-A, C-A);
  }
  bool cmp(Point A, Point B) {
      Point O;
      if (compare(O, A, B) == 0) return A.x < B.x;
      else return compare(O, A, B) > 0;
  }

• 一般情况下，使用 $\text{atan2}(y, x)$ 库更快，其中 $y = \text{Vector}.y, \ x = \text{Vector}.x$

向量旋转

向量旋转的公式为 （转角为 $\alpha$）

$$\begin{cases} x' = x \cos \alpha - y \sin \alpha \\ y' = x \sin \alpha + y \cos \alpha \end{cases}$$

Vector Rotate(Vector A, double rad) {
    return Vector(A.x*cos(rad)-A.y*sin(rad), A.x*sin(rad)+A.y*cos(rad));
}

特殊情况，如果令 $\text{rad} = \displaystyle\frac{\pi}{2}$，并且 $A$ 不是 $\bm{0}$ 向量

可以求出法线，并且将长度归一化

Vector Normal(Vector A) {
    double L = Length(A);
    return Vector(-A.y/L, A.x/L);
}

点和直线的位置关系

直线的参数方程
直线上的点满足 $P = P_0 + t\bm{v}$，已知直线上两个点 $A, B$，那么
直线表示为 $A+ (B-A) \bm{t}$

直线的交点
两条直线分别为 $P+ t_1 \bm{v}$，$Q+ t_2 \bm{w}$，令 $\bm{u} = P-Q$

$$\begin{cases} u_x = -v_xt_1 + w_xt_2 \\ u_y = -v_yt_1 + w_yt_2 \end{cases} \Rightarrow$$

$$\begin{pmatrix} -v_x & w_x \\ -v_y & w_y \end{pmatrix} \begin{pmatrix} t_1 \\ t_2 \end{pmatrix} = \begin{pmatrix} u_x \\ u_y \end{pmatrix}$$

伴随矩阵定义为矩阵 $\bm{A}$ 的余子式矩阵的转置

$$K = \det \left| \begin{matrix} -v_x & w_x \\ -v_y & w_y \end{matrix}\right| \Rightarrow \begin{pmatrix} t_1 \\ t_2 \end{pmatrix} = \frac{\begin{pmatrix} w_y & -w_x \\ v_y & -v_x \end{pmatrix}}{K} \cdot \begin{pmatrix} u_x \\ u_y \end{pmatrix}$$

由此可以推出

$$\begin{pmatrix} t_1 \\ t_2 \end{pmatrix} = \frac{\begin{pmatrix} u_xw_y - u_yw_x \\ u_xv_y - u_yv_x \end{pmatrix}}{\bm{w} \times \bm{v}}$$

令 $\bm{u} = \overrightarrow{QP}$，直线分别为 $P+t\bm{v}$，$Q+t\bm{w}$
向量 $\bm{u}$ 从第二条直线指向第一条直线
两直线交点在第一条直线上的参数为 $t_1$，在第二条直线上的参数为 $t_2$

$$t_1 = \frac{\text{cross}(\bm{w}, \bm{u})}{\text{cross}(\bm{v}, \bm{w})}, \quad t_2 = \frac{\text{cross}(\bm{v}, \bm{u})}{\text{cross}(\bm{v}, \bm{w})}$$

Point lineIntersection(Point P, Vector v, Point Q, Vector w) {
    Vector u = P-Q;
    double t = Cross(w, u) / Cross(v, w);
    return P + v*t;
}

点到直线的距离
可以利用多边形面积公式，对于直线 $AB$，直线外的点为 $P$
先计算三角形 $PAB$ 的有向面积的 $2$ 倍，$S = \text{cross}(\overrightarrow{AB}, \overrightarrow{AP})$
$P$ 到 $AB$ 的距离为三角形的高，$D = \text{cross}(\overrightarrow{AB}, \overrightarrow{AP}) / \text{Length}(\overrightarrow{AB})$

double distToLine(Point P, Point A, Point B) {
    Vector v1 = B-A, v2 = P-A;
    return fabs(Cross(v1, v2) / Length(v1));
}

点到线段的距离
不妨设 $P$ 在线段 $AB$ 上的投影点为 $Q$，如果
$\langle \overrightarrow{AB}, \overrightarrow{AP} \rangle< 0$，$Q$ 在 $AB$ 左侧，此时长度为 $|\overrightarrow{AP}|$
$\langle \overrightarrow{AB}, \overrightarrow{BP} \rangle> 0$，$Q$ 在 $AB$ 右侧，此时长度为 $|\overrightarrow{BP}|$
其他情况，只要求出点到直线的距离即可

double distToSegment(Point P, Point A, Point B) {
    if (A == B) return Length(P-A);
    Vector v1 = B-A, v2 = P-A, v3 = P-B;
    if (dcmp(Dot(v1, v2)) < 0) return Length(v2);
    else if (dcmp(Dot(v1, v3)) > 0) return Length(v3);
    else return fabs( Cross(v1, v2) / Length(v1) );
}

点在直线上的投影
直线为 $\overrightarrow{AB} = \bm{v}$，参数方程为 $A+t\bm{v}$
假设投影点为 $Q$，$Q$ 对应的参数为 $t_0$，那么 $\overrightarrow{Q}=A+t_0\bm{v}$
从而有 $\overrightarrow{QP} = P-(A+t_0\bm{v})$
$\langle \bm{v}, P-(A+t_0 \bm{v}) \rangle = 0$，即 $\langle \bm{v}, P-A \rangle - t_0 \langle \bm{v}, \bm{v} \rangle = 0$

计算出 $t_0 = \displaystyle \frac{\text{Dot}(\bm{v}, P-A)}{\text{Dot}(\bm{v}, \bm{v})}$，从而 $A + \bm{v} \cdot t_0$ 就是 $Q$

Point getLineProjection(Point P, Point A, Point B) {
    Vector v = B-A;
    return A + v * (Dot(v, P-A) / Dot(v, v));
}

线段树维护矩阵变换

Gentle Jena II

$x, y$ 绕原点逆时针旋转，可以通过矩阵变换描述

$$\begin{pmatrix} \cos t & -\sin t \\ \sin t & \cos t \end{pmatrix} \begin{pmatrix} x \\ y \end{pmatrix}$$

如果绕特定的点 $(X, Y)$ 旋转角度 $t$，可以写成
先平移 $(-a, -b)$，再绕原点旋转，最后再平移 $(a, b)$，构造仿射变换如下

$$\begin{pmatrix} 1 & 0 & X \\ 0 & 1 & Y \\ 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} \cos t & - \sin t & 0 \\ \sin t & \cos t & 0 \\ 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} 1 & 0 & -X \\ 0 & 1 & -Y \\ 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} x \\ y \\ 1 \end{pmatrix}$$

另外注意矩阵变换是满足分配律的，比如对于当前星星 $\bm{A}_i = \displaystyle \binom{x_i}{y_i}$，对于变换 $\bm{P}$

$\bm{P}(\bm{A}_1 + \bm{A}_2 + \cdots \bm{A}_n) = \bm{PA}_1 + \bm{PA}_2 + \cdots \bm{PA}_n = \bm{A}'_1 + \bm{A}'_2 + \cdots \bm{A}'_n$

所以只需要用线段树维护 $sx = \sum x, \ sy = \sum y$，分别表示区间中所有星星 $x, y$ 坐标的和

构造仿射变换 $\bm{P}$

• 初始化，对于每一个星星 $i$，$sx = x_i, sy = y_i$，建立线段树
• 对于操作 $1$，即执行仿射变换 $\bm{P}$，只需要 $\bm{P}\cdot (sx, sy)^{T}$，难点在于 $\text{push}$ 和 $\text{pull}$ 怎么写，以及懒标记处理
  $(sx', sy')^{T} \leftarrow \bm{P}\cdot (sx, sy)^{T}$
• 每次修改，难点在于 $\text{lazy}$ 如何处理
  $\textbf{lazy}$ 初始化为单位矩阵 $\bm{I}$，每一次对线段树节点 $p$ 执行变换 $\bm{C}$ 时
  $\bm{C} \cdot (t_p.sx, t_p.sy)^{T}$ 执行矩阵乘法，并且下传延迟标记
  递归到的每一个子区间都被修改了，后面的修改是基于前面修改基础上的
  对于子区间，执行 $t(p \ll 1).\textbf{lazy} \leftarrow \bm{C} \cdot t(p \ll 1).\textbf{lazy}$，右子树同理
  $\text{pull}$ 操作只需要执行，$t(p).sx = t(p \ll 1).sx + t(p\ll 1 \mid 1).sx$
• 对于操作 $2$ 的问询，只需要注意每次递归查询时候，下传延迟标记，最后找到表示区间 $[li, ri]$ 的节点 $u$
  $t(u).sx, t(u).sy$ 就是答案，注意在修改操作中更新 $sx, sy$ 需要取 $\text{mod}$

坑点提示

$$\sum \left(\bm{PA}_i \right) = \bm{P} \left(\sum \bm{A}_i \right)$$

所以维护的矩阵变换应该是

$$\bm{P} \cdot \begin{pmatrix} sx \\ sy \\ \text{len}(p) \end{pmatrix}$$

其中 $\text{len}(p)$ 表示线段树 $p$ 节点的区间长度

const int maxn = 1e5 + 500;
const int mod = 998244353;
int n, m;
typedef vector<vector<ll> > Matrix;
Matrix I(3, vector<ll>(3, 0));
 
void init() {
    for (int i = 0; i < 3; i++) I[i][i] = 1;
}
 
struct Point {
    int x, y;
} a[maxn];
 
 
Matrix trans(const Matrix &A, const Matrix &B) {
    int _n = A.size(), _m = B[0].size();
    Matrix res(_n, vector<ll> (_m, 0));
 
    for (int i = 0; i < _n; i++) {
        for (int j = 0; j < _m; j++)
            for (int k = 0; k < (int)A[0].size(); k++)
                res[i][j] = res[i][j] + A[i][k] * B[k][j];
    }
    return res;
}
 
Matrix mul(const Matrix &A, const Matrix &B) {
    int _n = A.size(), _m = B[0].size();
    Matrix res(_n, vector<ll> (_m, 0));
 
    for (int i = 0; i < _n; i++) {
        for (int j = 0; j < _m; j++)
            for (int k = 0; k < (int)A[0].size(); k++)
                res[i][j] = (res[i][j] + A[i][k] * B[k][j] + mod) % mod;
    }
    return res;
}
 
Matrix get_matrix(int X, int Y) {
    Matrix A(3, vector<ll>(3, 0));
    for (int i = 0; i < 3; i++) A[i][i] = 1;
    A[0][2] = X, A[1][2] = Y;
 
    Matrix B(3, vector<ll>(3, 0));
    B[0][1] = -1, B[1][0] = 1, B[2][2] = 1;
 
    Matrix C(3, vector<ll>(3, 0));
    for (int i = 0; i < 3; i++) C[i][i] = 1;
    C[0][2] = -X, C[1][2] = -Y;
 
    Matrix res = A;
    res = trans(res, B);
    res = trans(res, C);
    return res;
}
 
 
namespace segTree {
    struct node {
        int l, r, tag;
        ll sx, sy;
        Matrix lazy;
    } t[maxn << 2];
 
    void pull(int p) {
        t[p].sx = (t[p<<1].sx + t[p<<1|1].sx + mod) % mod;
        t[p].sy = (t[p<<1].sy + t[p<<1|1].sy + mod) % mod;
    }
 
    void apply(int p, const Matrix& C) {
        Matrix cur(3, vector<ll>(1, 0));
        cur[0][0] = t[p].sx, cur[1][0] = t[p].sy, cur[2][0] = t[p].r-t[p].l+1;
        Matrix nxt = trans(C, cur);
 
        t[p].sx = nxt[0][0], t[p].sy = nxt[1][0];
    }
 
    void push(int p) {
        if (!t[p].tag) return;
        if (t[p].l == t[p].r) return;
 
        t[p].tag = 0;
        t[p<<1].tag = t[p<<1|1].tag = 1;
 
        Matrix C = t[p].lazy;
        t[p<<1].lazy = trans(C, t[p<<1].lazy);
        t[p<<1|1].lazy = trans(C, t[p<<1|1].lazy);
        t[p].lazy = I;
 
        apply(p<<1, C), apply(p<<1|1, C);
    }
 
    void build(int p, int l, int r) {
        t[p].l = l, t[p].r = r;
        t[p].lazy = I;
        if (l >= r) {
            t[p].sx = a[l].x, t[p].sy = a[l].y;
            return;
        }
        int mid = (l+r) >> 1;
        build(p<<1, l, mid);
        build(p<<1|1, mid+1, r);
 
        pull(p);
    }
 
    void change(int p, const int l, const int r, const Matrix &P) {
        if (l <= t[p].l && t[p].r <= r) {
            t[p].tag = 1;
            apply(p, P);
            t[p].lazy = trans(P, t[p].lazy);
            return;
        }
 
        push(p);
        t[p].sx = t[p].sy = 0;
        int mid = (t[p].l + t[p].r) >> 1;
        if (l <= mid) change(p<<1, l, r, P);
        if (r > mid) change(p<<1|1, l, r, P);
 
        pull(p);
    }
 
    ll query(int p, const int l, const int r, bool fl) {
        if (l <= t[p].l && t[p].r <= r) {
            if (fl == 0) return (t[p].sx + mod) % mod;
            else return (t[p].sy + mod) % mod;
        }
        if (fl == 0) push(p);
 
        ll ans = 0;
        int mid = (t[p].l + t[p].r) >> 1;
        if (l <= mid) ans = (ans + query(p<<1, l, r, fl) + mod) % mod;
        if (r > mid) ans = (ans + query(p<<1|1, l, r, fl) + mod) % mod;
 
        return ans;
    }
}
 
 
int main() {
    //freopen("input.txt", "r", stdin);
    cin >> n >> m;
    init();
     
    using namespace segTree;
    // get data
    for (int i = 1; i <= n; i++) scanf("%d%d", &a[i].x, &a[i].y);
    build(1, 1, n);
 
    // get change
    while (m--) {
        int op, li, ri;
        scanf("%d%d%d", &op, &li, &ri);
        //debug(op), debug(li), debug(ri);
        if (op == 1) {
            int X, Y;
            scanf("%d%d", &X, &Y);
            Matrix P = get_matrix(X, Y);
            //dbg(P);
 
            change(1, li, ri, P);
        }
        if (op == 2) {
            ll sumx = query(1, li, ri, 0);
            ll sumy = query(1, li, ri, 1);
            //debug(sumy);
            sumx = (sumx + mod) % mod, sumy = (sumy + mod) % mod;
            ll ans = (sumx * sumy + mod) % mod;
            printf("%lld\n", ans);
        }
    }
}

行列式计算

直接从定义出发，递归计算

typedef vector<vector<double> > Matrix;

double calc(Matrix det) {
    int n = det.size();
    if (n == 1) return det[0][0];

    double detval = 0;
    Matrix tmp = vector<vector<double> >(n-1, vector<double> (n-1, 0));
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n-1; j++) for (int k = 0; k < n-1; k++) {
            if (k < i) tmp[j][k] = det[j+1][k];
            else tmp[j][k] = det[j+1][k+1];
        }
        detval += det[0][i] * pow(-1.0, i) * calc(tmp);
    }
    return detval;
}