计算几何常用函数

结构体与向量

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
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;
}

向量运算:点积

1
2
3
4
5
6
7
8
9
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));
}

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

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

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

1
2
3
4
5
6
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 函数求极角,对于向量Ai(x,y)\bm{A}_i(x, y),极角为atan2(y,x)\text{atan2}(y, x),单位为弧度
    逆时针为正方向,值域为[π,+π][-\pi, +\pi],保存在数组ang[i]\text{ang}[i]
    特别地,很多问题需要对极角排序,然后用双指针计算向量Ai,Aj\bm{A}_i, \bm{A}_j 的夹角θ\theta
    由于极角往往是环形的,处理完A(ij)\bm{A}(i \to j) 的夹角,往往还需要处理A(ji)\bm{A}(j \to i) 的夹角
    比如i=n,j=1i = n, j = 1,求第nn 个向量往第11 个向量的转角
    这里采用拆换成链并复制的方法,令ang[1n]\text{ang}[1\cdots n] 每个元素加上2π2\pi 之后
    放入ang[n+1n+n]\text{ang}[n+1 \cdots n+n]
  • 此外还可以根据叉积来对极角排序
    Cross(A,B)=0\text{Cross}(A, B) = 0 向量重合,Cross(A,B) >0\text{Cross}(A, B) > 0 向量BBAA 上方
    Cross(A,B)<0\text{Cross}(A, B) < 0 向量BBAA 下方
    1
    2
    3
    4
    5
    6
    7
    8
    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;
    }
  • 一般情况下,使用atan2(y,x)\text{atan2}(y, x) 库更快,其中y=Vector.y, x=Vector.xy = \text{Vector}.y, \ x = \text{Vector}.x

向量旋转
01
向量旋转的公式为 (转角为α\alpha

{x=xcosαysinαy=xsinα+ycosα\begin{cases} x' = x \cos \alpha - y \sin \alpha \\ y' = x \sin \alpha + y \cos \alpha \end{cases}

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

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

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

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

点和直线的位置关系

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

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

{ux=vxt1+wxt2uy=vyt1+wyt2\begin{cases} u_x = -v_xt_1 + w_xt_2 \\ u_y = -v_yt_1 + w_yt_2 \end{cases} \Rightarrow

(vxwxvywy)(t1t2)=(uxuy)\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}

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

K=detvxwxvywy(t1t2)=(wywxvyvx)K(uxuy)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}

由此可以推出

(t1t2)=(uxwyuywxuxvyuyvx)w×v\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}}

u=QP\bm{u} = \overrightarrow{QP},直线分别为P+tvP+t\bm{v}Q+twQ+t\bm{w}
向量u\bm{u} 从第二条直线指向第一条直线
两直线交点在第一条直线上的参数为t1t_1,在第二条直线上的参数为t2t_2

t1=cross(w,u)cross(v,w),t2=cross(v,u)cross(v,w)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})}

1
2
3
4
5
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;
}

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

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

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

1
2
3
4
5
6
7
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) );
}

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

计算出t0=Dot(v,PA)Dot(v,v)t_0 = \displaystyle \frac{\text{Dot}(\bm{v}, P-A)}{\text{Dot}(\bm{v}, \bm{v})},从而A+vt0A + \bm{v} \cdot t_0 就是QQ

1
2
3
4
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,yx, y 绕原点逆时针旋转,可以通过矩阵变换描述

(costsintsintcost)(xy)\begin{pmatrix} \cos t & -\sin t \\ \sin t & \cos t \end{pmatrix} \begin{pmatrix} x \\ y \end{pmatrix}

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

(10X01Y001)(costsint0sintcost0001)(10X01Y001)(xy1)\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}

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

P(A1+A2+An)=PA1+PA2+PAn=A1+A2+An\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=x, sy=ysx = \sum x, \ sy = \sum y,分别表示区间中所有星星x,yx, y 坐标的和

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

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

坑点提示

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

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

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

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
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);
}
}
}

行列式计算

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
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;
}