PS 无限版

GitHub跳转原题关系图 Canvas 关系图返回列表

把平面操作统一为仿射变换，在线段树中维护区间的一次矩和二次矩并懒惰复合变换。

OJ: shumeng

题目 ID: CSP202206E

难度:提高+/省选-

标签:数据结构线段树仿射变换线性代数计算几何

日期: 2026-07-31 16:21

形式化题目

给定 $n$ 个平面点，支持对编号连续区间执行平移、旋转、放缩、关于直线对称、投影五种修改操作，并查询区间点集的重心，或到指定点 $(a,b)$ 的距离平方和。

坐标和操作参数都是实数，任意时刻点坐标绝对值不超过 $10^6$ ； $n\le 5\times 10^5$ ，操作数不超过 $2\times 10^4$ ，因此不能每次操作都逐点修改。

思路

先看逐点执行操作的朴素程序，它只在数据很小时可行，但能直观对应题目语义：

cpp

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/**
 * Author by Rainboy blog: https://rainboylv.com github: https://github.com/rainboylvx
 * rbook: -> https://rbook.roj.ac.cn  https://rbook2.roj.ac.cn
 * rainboy的学习导航网站: https://idx.roj.ac.cn
 * create_at: 2026-07-31 16:21
 * update_at: 2026-08-17 22:40
 */
// brute.cpp：小数据暴力解，逐点执行每次仿射变换，用来辅助对拍。
#include <bits/stdc++.h>
using namespace std;

struct Point {
    long double x, y;
};

int n, query_count;
vector<Point> point;

void translate(int left, int right, long double x, long double y) {
    for (int i = left; i <= right; i++) {
        point[i].x += x;
        point[i].y += y;
    }
}

void rotate_point(int left, int right, long double x, long double y, long double theta) {
    long double c = cos(theta), s = sin(theta);
    for (int i = left; i <= right; i++) {
        long double old_x = point[i].x - x;
        long double old_y = point[i].y - y;
        point[i].x = c * old_x - s * old_y + x;
        point[i].y = s * old_x + c * old_y + y;
    }
}

void scale_point(int left, int right, long double x, long double y, long double lambda) {
    for (int i = left; i <= right; i++) {
        point[i].x = x + lambda * (point[i].x - x);
        point[i].y = y + lambda * (point[i].y - y);
    }
}

void line_transform(int left, int right, long double theta, long double y0, bool projection) {
    long double c = cos(theta), s = sin(theta);
    long double m00, m01, m10, m11;
    if (projection) {
        m00 = c * c;
        m01 = c * s;
        m10 = c * s;
        m11 = s * s;
    } else {
        m00 = c * c - s * s;
        m01 = 2 * c * s;
        m10 = 2 * c * s;
        m11 = s * s - c * c;
    }
    long double vx = -m01 * y0;
    long double vy = y0 - m11 * y0;
    for (int i = left; i <= right; i++) {
        long double old_x = point[i].x;
        long double old_y = point[i].y;
        point[i].x = m00 * old_x + m01 * old_y + vx;
        point[i].y = m10 * old_x + m11 * old_y + vy;
    }
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    cin >> n >> query_count;
    point.resize(n + 1);
    for (int i = 1; i <= n; i++) cin >> point[i].x >> point[i].y;

    cout << fixed << setprecision(6);
    while (query_count--) {
        int type, left, right;
        cin >> type >> left >> right;
        if (type == 1) {
            long double a, b;
            cin >> a >> b;
            translate(left, right, a, b);
        } else if (type == 2) {
            long double a, b, theta;
            cin >> a >> b >> theta;
            rotate_point(left, right, a, b, theta);
        } else if (type == 3) {
            long double a, b, lambda;
            cin >> a >> b >> lambda;
            scale_point(left, right, a, b, lambda);
        } else if (type == 4 || type == 5) {
            long double theta, y0;
            cin >> theta >> y0;
            line_transform(left, right, theta, y0, type == 5);
        } else if (type == 6) {
            long double sx = 0, sy = 0;
            for (int i = left; i <= right; i++) {
                sx += point[i].x;
                sy += point[i].y;
            }
            long double count = right - left + 1;
            cout << sx / count << ' ' << sy / count << '\n';
        } else {
            long double a, b, answer = 0;
            cin >> a >> b;
            for (int i = left; i <= right; i++) {
                long double dx = point[i].x - a;
                long double dy = point[i].y - b;
                answer += dx * dx + dy * dy;
            }
            cout << answer << '\n';
        }
    }

    return 0;
}

正解的关键是把五种几何操作统一成同一个数学形式，再用线段树批量维护。

统一表示平面操作

把点写成列向量 $p=(x,y)^T$ 。五种修改都可以写成仿射变换：

p' = Mp+v,

其中 $M$ 是 $2\times2$ 矩阵， $v$ 是平移向量。例如：

平移 $(a,b)$ ： $M=I$ ， $v=(a,b)^T$ ；
绕 $(a,b)$ 逆时针旋转 $\theta$ ： $M=\begin{bmatrix}\cos\theta&-\sin\theta\\\sin\theta&\cos\theta\end{bmatrix}$ ， $v=(a,b)^T-M(a,b)^T$ ；
绕 $(a,b)$ 放缩 $\lambda$ ： $M=\lambda I$ ， $v=(1-\lambda)(a,b)^T$ 。

令直线方向单位向量为 $u=(\cos\theta,\sin\theta)^T$ ，直线上点 $p_0=(0,y_0)^T$ 。投影的线性部分是 $uu^T$ ，对称的线性部分是 $2uu^T-I$ ，两者的平移部分都为 $p_0-Mp_0$ 。

若连续执行两个变换 $T_1(p)=M_1p+v_1$ 和 $T_2(p)=M_2p+v_2$ ，合成结果仍是仿射变换：

T_1(T_2(p))=(M_1M_2)p+(M_1v_2+v_1)。

因此线段树的懒标记只需保存一个仿射变换。

维护区间矩

对每个线段树节点维护一次矩与二次矩：

(N,S_x,S_y,S_{xx},S_{yy},S_{xy}) =\left(m,\sum x_i,\sum y_i,\sum x_i^2,\sum y_i^2,\sum x_iy_i\right)。

设一次仿射变换为 $x'=ax+by+u$ 、 $y'=cx+dy+v$ ，那么整段覆盖时矩可以整体更新：

\begin{aligned} S_x'&=aS_x+bS_y+Nu,\\ S_y'&=cS_x+dS_y+Nv,\\ S_{xx}'&=a^2S_{xx}+b^2S_{yy}+2abS_{xy}+2u(aS_x+bS_y)+Nu^2,\\ S_{yy}'&=c^2S_{xx}+d^2S_{yy}+2cdS_{xy}+2v(cS_x+dS_y)+Nv^2,\\ S_{xy}'&=acS_{xx}+(ad+bc)S_{xy}+bdS_{yy}\\ &\quad+u(cS_x+dS_y)+v(aS_x+bS_y)+Nuv。 \end{aligned}

这样整段覆盖时只更新节点统计量和懒标记，不需要访问区间内的每个点。

回答查询

重心： $\left(\dfrac{S_x}{N},\dfrac{S_y}{N}\right)$ 。
到 $(a,b)$ 的距离平方和： $S_{xx}+S_{yy}-2aS_x-2bS_y+N(a^2+b^2)$ 。

区间查询用线段树合并矩；下传懒标记前先把同一个仿射变换应用到左右儿子。

代码

cpp

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
175/**
 * Author by Rainboy blog: https://rainboylv.com github: https://github.com/rainboylvx
 * rbook: -> https://rbook.roj.ac.cn  https://rbook2.roj.ac.cn
 * rainboy的学习导航网站: https://idx.roj.ac.cn
 * create_at: 2026-07-31 16:21
 * update_at: 2026-08-17 22:40
 */
#include <bits/stdc++.h>
using namespace std;

struct Moment {
    // 区间内点的数量、一次矩和二次矩。
    long double count, sx, sy, sxx, syy, sxy;
};

struct Affine {
    // (x,y) -> M(x,y) + v。
    long double m00, m01, m10, m11, vx, vy;
};

struct Node {
    Moment data;
    Affine lazy;
};

int n, query_count, base;
vector<Node> tree;

Affine identity() {
    return {1, 0, 0, 1, 0, 0};
}

Affine compose(const Affine &first, const Affine &second) {
    // 先执行 second，再执行 first。
    return {
        first.m00 * second.m00 + first.m01 * second.m10,
        first.m00 * second.m01 + first.m01 * second.m11,
        first.m10 * second.m00 + first.m11 * second.m10,
        first.m10 * second.m01 + first.m11 * second.m11,
        first.m00 * second.vx + first.m01 * second.vy + first.vx,
        first.m10 * second.vx + first.m11 * second.vy + first.vy,
    };
}

Moment merge_moment(const Moment &left, const Moment &right) {
    return {left.count + right.count, left.sx + right.sx, left.sy + right.sy,
            left.sxx + right.sxx, left.syy + right.syy, left.sxy + right.sxy};
}

void apply_node(int node, const Affine &transform) {
    Moment old = tree[node].data;
    long double nx = transform.m00 * old.sx + transform.m01 * old.sy + old.count * transform.vx;
    long double ny = transform.m10 * old.sx + transform.m11 * old.sy + old.count * transform.vy;
    long double nxx = transform.m00 * transform.m00 * old.sxx
        + transform.m01 * transform.m01 * old.syy
        + 2 * transform.m00 * transform.m01 * old.sxy
        + 2 * transform.vx * (transform.m00 * old.sx + transform.m01 * old.sy)
        + old.count * transform.vx * transform.vx;
    long double nyy = transform.m10 * transform.m10 * old.sxx
        + transform.m11 * transform.m11 * old.syy
        + 2 * transform.m10 * transform.m11 * old.sxy
        + 2 * transform.vy * (transform.m10 * old.sx + transform.m11 * old.sy)
        + old.count * transform.vy * transform.vy;
    long double nxy = transform.m00 * transform.m10 * old.sxx
        + (transform.m00 * transform.m11 + transform.m01 * transform.m10) * old.sxy
        + transform.m01 * transform.m11 * old.syy
        + transform.vx * (transform.m10 * old.sx + transform.m11 * old.sy)
        + transform.vy * (transform.m00 * old.sx + transform.m01 * old.sy)
        + old.count * transform.vx * transform.vy;
    tree[node].data.sx = nx;
    tree[node].data.sy = ny;
    tree[node].data.sxx = nxx;
    tree[node].data.syy = nyy;
    tree[node].data.sxy = nxy;
    tree[node].lazy = compose(transform, tree[node].lazy);
}

void pull(int node) {
    tree[node].data = merge_moment(tree[node << 1].data, tree[node << 1 | 1].data);
}

void push(int node) {
    Affine lazy = tree[node].lazy;
    if (lazy.m00 == 1 && lazy.m01 == 0 && lazy.m10 == 0 && lazy.m11 == 1
            && lazy.vx == 0 && lazy.vy == 0) return;
    apply_node(node << 1, lazy);
    apply_node(node << 1 | 1, lazy);
    tree[node].lazy = identity();
}

void update(int node, int left, int right, int ql, int qr, const Affine &transform) {
    if (ql <= left && right <= qr) {
        apply_node(node, transform);
        return;
    }
    push(node);
    int middle = (left + right) >> 1;
    if (ql <= middle) update(node << 1, left, middle, ql, qr, transform);
    if (qr > middle) update(node << 1 | 1, middle + 1, right, ql, qr, transform);
    pull(node);
}

Moment query(int node, int left, int right, int ql, int qr) {
    if (ql <= left && right <= qr) return tree[node].data;
    push(node);
    int middle = (left + right) >> 1;
    Moment result = {};
    if (ql <= middle) result = merge_moment(result, query(node << 1, left, middle, ql, qr));
    if (qr > middle) result = merge_moment(result, query(node << 1 | 1, middle + 1, right, ql, qr));
    return result;
}

Affine transform_line(long double theta, long double y0, bool projection) {
    long double c = cos(theta), s = sin(theta);
    Affine result;
    if (projection) result = {c * c, c * s, c * s, s * s, 0, 0};
    else result = {c * c - s * s, 2 * c * s, 2 * c * s, s * s - c * c, 0, 0};
    result.vx = -result.m01 * y0;
    result.vy = y0 - result.m11 * y0;
    return result;
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    cin >> n >> query_count;
    base = 1;
    while (base < n) base <<= 1;
    tree.resize(base << 1);
    for (int i = 1; i < (base << 1); i++) tree[i].lazy = identity();
    for (int i = 1; i <= n; i++) {
        long double x, y;
        cin >> x >> y;
        tree[base + i - 1].data = {1, x, y, x * x, y * y, x * y};
    }
    for (int i = base - 1; i >= 1; i--) pull(i);

    cout << fixed << setprecision(6);
    while (query_count--) {
        int type, left, right;
        cin >> type >> left >> right;
        if (type <= 5) {
            long double a, b, theta, lambda;
            Affine transform;
            if (type == 1) {
                cin >> a >> b;
                transform = {1, 0, 0, 1, a, b};
            } else if (type == 2) {
                cin >> a >> b >> theta;
                long double c = cos(theta), s = sin(theta);
                transform = {c, -s, s, c, a - c * a + s * b, b - s * a - c * b};
            } else if (type == 3) {
                cin >> a >> b >> lambda;
                transform = {lambda, 0, 0, lambda, a * (1 - lambda), b * (1 - lambda)};
            } else {
                cin >> theta >> b;
                transform = transform_line(theta, b, type == 5);
            }
            update(1, 1, base, left, right, transform);
        } else if (type == 6) {
            Moment result = query(1, 1, base, left, right);
            cout << result.sx / result.count << ' ' << result.sy / result.count << '\n';
        } else {
            long double a, b;
            cin >> a >> b;
            Moment result = query(1, 1, base, left, right);
            long double answer = result.sxx + result.syy - 2 * a * result.sx - 2 * b * result.sy
                + result.count * (a * a + b * b);
            cout << answer << '\n';
        }
    }

    return 0;
}

复杂度

建树时间为 $O(n)$ 。每次区间修改或查询访问 $O(\log n)$ 个线段树节点，每个节点的矩更新和合并都是 $O(1)$ ，总时间复杂度为 $O((n+q)\log n)$ ，空间复杂度为 $O(n)$ 。

总结

这道题的关键是把“对很多点做同一个几何操作”转成“对点集的低阶矩做一次更新”。仿射变换保证二阶矩仍能用有限个统计量表示，线段树负责组织区间覆盖与查询。使用 long double 可以减少大数二次矩在展开式中相减时的精度损失。

题目来自第 26 次 CSP 认证 T5，评测使用自造高质量复刻数据。我们承认原始题面与大样例版权均归中国计算机学会（CCF）所有，因此题面与评测服务均免费对外开放。如果您认为我们侵犯了您的权益，可联系我们。

时间限制： ~~10.0 秒~~ 5.0 秒

空间限制： 2048 MB

题目背景

众所周知，PS 是一款图片编辑软件，编辑图片的本质是操作各像素。

但是，传统的图片编辑只能对有限个像素进行操作，而这对于一名数学系学生是不可忍受的——竟然不能把有限的、离散的问题推广到无穷的、连续的问题，这真是不可忍受。

正如在线性代数的理论我们为了将有限维线性空间推广到无穷维线性空间所做的那样，现在我们可以假定一张图片是一个无穷大的二维平面（方便起见，我们假定它是一个平面直角坐标系），其上的每个像素可以用 $(a,b)$ 表示（注意， $a,b$ 是实数）。类似于线性代数无穷维线性空间关于基的讨论，我们实际上不关心所有的像素，而只关注于其中的有限个像素，通过对每一组有限大小的像素集的刻画来描述图片整体的编辑情况。

当然，尽管在原理上成功的把 PS 升级成了 PSI(PS Infinite)，但就结论而言，我们应当讨论传统 PS 中的各种操作在 PSI 上的推广和实现。出于简单起见，我们只考虑平移、旋转、放缩、对称和投影这些基本的编辑操作。

题目描述

给定正整数 $n$ 与平面上的一些点 $(x_i,y_i)_{i=1}^n\subset \Re^2$ ，支持以下操作:

$1~l~r~a~b$ ：将编号在 $[l,r]$ 中的点平移 $\vec v=(a,b)$ 。
- 即沿 $\vec v$ 方向平移 $|\vec v|$ 的距离。
$2~l~r~a~b~\theta$ ：将编号在 $[l,r]$ 中的点以 $(a,b)$ 为中心逆时针旋转 $\theta$ 。
- 保证 $\theta \in (-\pi ,\pi)$ ，以弧度制给出。
$3~l~r~a~b~\lambda$ ：将编号在 $[l,r]$ 中的点以 $(a,b)$ 为中心放缩 $|\lambda|$ 倍。
- 即在指向 $(a,b)$ 的方向所在直线上移动，距离缩小 $(|\lambda|<1)$ 或变大 $(|\lambda|>1)$ 。
- 例如 $\lambda = 0$ 即变为 $(a,b)$ ， $\lambda<0$ 则其相对于 $(a,b)$ 的方向会相反。
$4~l~r~\theta~y_0$ ：将编号在 $[l,r]$ 中的点以 $(\tan\theta)x+y_0$ 为对称轴做对称变换。
- 保证 $\theta \in(-\dfrac{\pi}{2},\dfrac{\pi}{2})$ ，以弧度制给出。
- 例如， $\theta = 0,y_0=0$ 即沿 $x$ 轴对称。
$5~l~r~\theta~y_0$ ：将编号在 $[l,r]$ 中的点投影到 $(\tan\theta)x+y_0$ 。
- 保证 $\theta \in(-\dfrac{\pi}{2},\dfrac{\pi}{2})$ ，以弧度制给出。
- 例如， $\theta = 0,y_0=0$ 即投影到 $x$ 轴上。
$6~l~r$ ：求编号在 $[l,r]$ 中的点的重心。
- 点集 $\{(a_i,b_i)\}_{i=1}^m$ 的重心定义为 $(\sum\limits_{i=1}^m \dfrac{a_i}{m},\sum\limits_{i=1}^m\dfrac{b_i}{m})$ 。
$7~l~r~a~b$ ：求编号在 $[l,r]$ 中的点到 $(a,b)$ 的距离的平方的和（注意，不是距离的和的平方）。
- 点集 $\{(a_i,b_i)\}_{i=1}^m$ 到 $(a,b)$ 的距离的平方的和即 $\sum\limits_{i=1}^m (a_i-a)^2+(b_i-b)^2$ 。

输入格式

从标准输入读入数据。

第一行一个整数 $n,q$ 表示点数和操作数。

接下来 $n$ 行，每行两个实数表示 $(x_i,y_i)$ 。

接下来 $q$ 行，每行若干个实数表示一次操作，保证格式同题面。

输出格式

输出到标准输出。

若干行，每行依次对 $6$ 和 $7$ 操作输出两个或一个实数，表示所求的重心坐标或距离平方和。

input1

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
3110 20
26.389153 -31.339463
-98.664509 -58.061567
16.023894 14.489272
-67.840842 -74.793309
19.790708 -87.062719
31.541964 88.441505
-75.918013 24.526470
57.288832 -39.033977
38.274184 -67.446883
-90.906424 -73.528612
3 4 4 32.938694 -6.774595 1.000221
1 2 6 69.965610 -39.563795
4 3 10 -1.399075 38.282976
4 6 7 -1.016301 61.080461
7 9 10 76.549276 22.856189
7 3 7 -96.501727 5.585970
6 8 9
4 2 8 1.215917 -90.918350
7 4 8 55.948842 38.373278
1 5 9 -83.845362 -6.619437
5 6 9 -1.202044 -90.146760
7 1 4 -81.574047 -56.555229
3 1 5 75.690820 60.620104 0.980271
4 5 9 1.512746 89.531420
5 2 5 0.071305 79.784122
6 2 4
1 3 6 90.288492 72.829660
6 4 4
7 1 10 -51.991614 -6.732535
5 5 6 0.087950 71.164056

output1

1
2
3
4
5
6
7
821029.678359
120220.146461
-14.172376 -63.985055
95006.134951
52111.910474
2.849235 79.987632
35.040886 148.667661
302347.683678

样例 2

见题目文件区的 2.in 和 2.ans。

样例 2 解释

该样例中仅有 $1,3,6$ 操作，且 $n,q\le 2000$ 。

样例 3

见题目文件区的 3.in 和 3.ans。

样例 3 解释

该样例中仅有 $1,3,6,7$ 操作，且 $n,q\le 2000$ 。

样例 4

见题目文件区的 4.in 和 4.ans。

样例 4 解释

该样例中仅有 $1,2,3,6,7$ 操作，且 $n,q\le 2000$ 。

样例 5

见题目文件区的 5.in 和 5.ans。

样例 5 解释

该样例中 $n,q\le 2000$ ，包含全部操作。

样例 6

见题目文件区的 6.in 和 6.ans。

样例 6 解释

该样例中 $n\le 5\times 10^5,q\le 2\times 10^4$ ，包含全部操作。

子任务

对于所有数据，保证任意时刻任意一点的横纵坐标的绝对值均不超过 $10^6$ 。

本题采用捆绑测试，你只有通过一个子任务中的所有测试点才能得到该子任务的分数。

子任务编号	$n\le$	$q\le$	可能出现的操作编号	分值
$1$	$2000$	$2000$	$1,3,6,7$	$10$
$2$	$2000$	$2000$	$1,2,3,6,7$	$10$
$3$	$2000$	$2000$	$1,2,3,4,5,6,7$	$10$
$4$	$5\times 10^5$	$2\times 10^4$	$1,3,6$	$20$
$5$	$5\times 10^5$	$2\times 10^4$	$1,3,6,7$	$20$
$6$	$5\times 10^5$	$2\times 10^4$	$1,2,3,6,7$	$20$
$7$	$5\times 10^5$	$2\times 10^4$	$1,2,3,4,5,6,7$	$10$

提示

为了避免精度误差，评测时选手的输出与标准程序的输出相对或绝对误差不超过 $10^{-3}$ 即算通过。

其中，实数 $a,b$ 的绝对误差即 $|a-b|$ ，相对误差即 $\dfrac{|a-b|}{\max(|a|,|b|)}$ 。

PS 无限版

目录

形式化题目

思路

统一表示平面操作

维护区间矩

回答查询

代码

复杂度

总结