Project Euler 790

Project Euler 790

题目

Clock Grid

There is a grid of length and width \(50515093\) points. A clock is placed on each grid point. The clocks are all analogue showing a single hour hand initially pointing at \(12\).

A sequence \(S_t\) is created where:

\[ \begin{aligned} S_0 &= 290797\\ S_t &= S_{t-1}^2 \bmod 50515093 &t>0 \end{aligned} \]

The four numbers \(N_t = (S_{4t-4}, S_{4t-3}, S_{4t-2}, S_{4t-1})\) represent a range within the grid, with the first pair of numbers representing the x-bounds and the second pair representing the y-bounds. For example, if \(N_t = (3,9,47,20)\), the range would be \(3\le x\le 9\) and \(20\le y\le47\), and would include \(196\) clocks.

For each \(t\) \((t>0)\), the clocks within the range represented by \(N_t\) are moved to the next hour \(12\rightarrow 1\rightarrow 2\rightarrow \cdots\).

We define \(C(t)\) to be the sum of the hours that the clock hands are pointing to after timestep \(t\).

You are given \(C(0) = 30621295449583788\), \(C(1) = 30613048345941659\), \(C(10) = 21808930308198471\) and \(C(100) = 16190667393984172\).

Find \(C(10^5)\).

解决方案

本题将使用线段树和扫描线完成。令\(N=10^5\)

对每个点 \((x,y)\),记它被拨动的总次数为 \(k(x,y)\)。由于每次拨动相当于“加 1 小时”,并且周期为 \(12\),因此只需要关心\(r(x,y)=k(x,y)\bmod 12\in\{0,1,\dots,11\}.\)

钟面显示的小时数 \(H(r)\) 为:如果\(r=0\) 表示仍在 \(12\),如果\(r=1,\dots,11\) 表示显示 \(1,\dots,11\)。因此\(H(r)=\begin{cases}12,& r=0,\\r,& 1\le r\le 11.\end{cases}\)。于是\(\displaystyle{C(t)=\sum_{x=0}^{M-1}\sum_{y=0}^{M-1} H(r(x,y)).}\)

题目中,每步给的是包含端点的矩形:\(x_{\min}\le x\le x_{\max}, y_{\min}\le y\le y_{\max}.\)

为了让扫描线更自然,将其转为半开形式:\([x_0,x_1)\times [y_0,y_1),x_0=x_{\min},x_1=x_{\max}+1,y_0=y_{\min},y_1=y_{\max}+1.\)这样矩形内的点数恰好是\((x_1-x_0)(y_1-y_0).\)每次操作就是对该半开矩形内所有点的 \(r(x,y)\)\(1\pmod{12}\)

虽然 \(M\) 很大,但只有 \(N\) 个矩形,每个矩形只引入两个 \(x\) 边界和两个 \(y\) 边界。接下来我们把所有坐标进行离散化。把所有矩形的 \(x_0,x_1\) 全部收集,再加入 \(0\)\(M\),排序去重得到:\(X_0<X_1<\cdots <X_{p-1}.\)同理得到:\(Y_0<Y_1<\cdots <Y_{q-1}.\)

这些边界把整个点阵划分成许多压缩块:\(x\) 方向的第 \(i\) 段长度:\(\Delta x_i=X_{i+1}-X_i\)\(y\) 方向的第 \(j\) 段长度:\(\Delta y_j=Y_{j+1}-Y_j\)

在任意一个块 \([X_i,X_{i+1})\times [Y_j,Y_{j+1})\) 内,所有点被哪些矩形覆盖是完全相同的,因此 \(r(x,y)\) 在该块上是常数。所以最终答案可以写成\(\displaystyle{C=\sum_{i=0}^{p-2}\sum_{j=0}^{q-2}\Delta x_i\Delta y_j\cdot H(r_{i,j}),}\)

接下来使用线段树和扫描线进行加速。对一个矩形 \([x_0,x_1)\times [y_0,y_1)\),当扫描线位于 \(y\) 时,若 \(y\in[y_0,y_1)\),它对区间 \([x_0,x_1)\) 的所有点产生 \(+1(mod12)\) 的影响,否则无影响。

因此它在 \(y=y_0\) 产生加入的事件,在 \(y=y_1\) 产生移除的事件:

  • \(y=y_0\):对 \([x_0,x_1)\)\(+1\)
  • \(y=y_1\):对 \([x_0,x_1)\)\(-1\)

把所有事件按 \(y\) 排序。扫描过程中,在相邻两个 \(y\) 边界之间 \([Y_j,Y_{j+1})\),活动矩形集合不变,于是每个 \(x\) 段的余数 \(r\) 不变。

设当前 \(y\) 水平带内,\(x\)\(i\) 段的余数为 \(r_i\),那么该水平带内每一行(每个 \(y\) 点)对应的 \(x\) 总和为:\(\displaystyle{S=\sum_{i=0}^{p-2}\Delta x_i\cdot H(r_i).}\)该水平带高度为 \(\Delta y_j\),对总答案贡献:\(\Delta y_j\cdot S.\)

因此问题变成:

  • 动态维护一维 \(x\) 段数组,每段有长度 \(\Delta x_i\) 和状态 \(r_i\in[0,11]\)
  • 支持区间加 \(\pm 1\pmod{12}\)
  • 随时能求 \(S=\sum \Delta x_i H(r_i)\)

这正是线段树 + 懒标记的典型应用。

我们对线段树的每个节点(覆盖一段连续的 \(x\) 段),维护一个长度为 \(12\) 的数组:\(L_u\),含义是该节点覆盖范围内:余数为 \(r\)\(x\) 总长度之和为 \(L_u[r]\)

那么该节点的小时和就是:\(\displaystyle{12\cdot L_u[0]+\sum_{r=1}^{11} r\cdot L_u[r].}\)

区间加 \(d\pmod{12}\) 的效果是把所有余数整体平移:\(L_u'[ (r+d)\bmod 12 ]=L_u[r].\)这是一个纯旋转操作,因此可以在 \(O(1)\) 内完成,并用懒标记记录旋转量即可。

最终,在扫描线每次跨过 \(\Delta y\) 的高度时,根节点的小时和就是当前 \(S\),贡献 \(\Delta y\cdot S\)

代码

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
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
#include <bits/stdc++.h>
using namespace std;

typedef long long ll;

int MOD = 50515093;
const int N = 100000;
struct Rect
{
int x0, x1, y0, y1;
};

struct Event
{
int l, r, delta;
};

struct SegTree
{
int n;
vector<int> freq;
vector<int> lazy;

SegTree() : n(0) {}

SegTree(const vector<int> &w)
{
n = (int)w.size();
freq.assign((4 * n + 5) * 12, 0);
lazy.assign(4 * n + 5, 0);
build(1, 0, n, w);
}

inline void rotate_node(int u, int sh)
{
sh %= 12;
if (sh < 0)
sh += 12;
if (!sh)
return;
int tmp[12];
int *p = &freq[u * 12];
for (int i = 0; i < 12; i++)
tmp[(i + sh) % 12] = p[i];
for (int i = 0; i < 12; i++)
p[i] = tmp[i];
lazy[u] = (lazy[u] + sh) % 12;
}

inline void push(int u)
{
if (!lazy[u])
return;
int sh = lazy[u];
rotate_node(u << 1, sh);
rotate_node(u << 1 | 1, sh);
lazy[u] = 0;
}

inline void pull(int u)
{
int *p = &freq[u * 12];
int *L = &freq[(u << 1) * 12];
int *R = &freq[(u << 1 | 1) * 12];
for (int i = 0; i < 12; i++)
p[i] = L[i] + R[i];
}

void build(int u, int l, int r, const vector<int> &w)
{
if (l + 1 == r)
{
freq[u * 12 + 0] = (int)w[l];
return;
}
int m = (l + r) >> 1;
build(u << 1, l, m, w);
build(u << 1 | 1, m, r, w);
pull(u);
}

void update(int ql, int qr, int delta)
{
update(1, 0, n, ql, qr, delta);
}

void update(int u, int l, int r, int ql, int qr, int delta)
{
if (qr <= l || r <= ql)
return;
if (ql <= l && r <= qr)
{
rotate_node(u, delta);
return;
}
push(u);
int m = (l + r) >> 1;
update(u << 1, l, m, ql, qr, delta);
update(u << 1 | 1, m, r, ql, qr, delta);
pull(u);
}

inline ll total() const
{
const int *p = &freq[1 * 12];
ll res = 12LL * p[0];
for (int i = 1; i < 12; i++)
res += 1LL * i * p[i];
return res;
}
};

ll solve(int N)
{
ll s = 290797;
vector<Rect> rects;
rects.reserve(N);

vector<int> xs, ys;
xs.reserve(2 * N + 2);
ys.reserve(2 * N + 2);
xs.push_back(0);
xs.push_back(MOD);
ys.push_back(0);
ys.push_back(MOD);

for (int i = 0; i < N; i++)
{
int a = s;
s = s * s % MOD;
int b = s;
s = s * s % MOD;
int c = s;
s = s * s % MOD;
int d = s;
s = s * s % MOD;

int x0 = min(a, b);
int x1 = max(a, b) + 1;
int y0 = min(c, d);
int y1 = max(c, d) + 1;

rects.push_back({x0, x1, y0, y1});
xs.push_back(x0);
xs.push_back(x1);
ys.push_back(y0);
ys.push_back(y1);
}

sort(xs.begin(), xs.end());
xs.erase(unique(xs.begin(), xs.end()), xs.end());
sort(ys.begin(), ys.end());
ys.erase(unique(ys.begin(), ys.end()), ys.end());

int nx = (int)xs.size();
int ny = (int)ys.size();

vector<int> w(nx - 1);
for (int i = 0; i + 1 < nx; i++)
w[i] = xs[i + 1] - xs[i];

SegTree seg(w);

vector<vector<Event>> events(ny);
for (const auto &rc : rects)
{
int xl = lower_bound(xs.begin(), xs.end(), rc.x0) - xs.begin();
int xr = lower_bound(xs.begin(), xs.end(), rc.x1) - xs.begin();
int yl = lower_bound(ys.begin(), ys.end(), rc.y0) - ys.begin();
int yr = lower_bound(ys.begin(), ys.end(), rc.y1) - ys.begin();
events[yl].push_back({xl, xr, +1});
events[yr].push_back({xl, xr, -1});
}

ll ans = 0;
int yPrev = 0;
for (int i = 0; i < ny; i++)
{
int y = ys[i];
ll dy = (ll)y - yPrev;
ans += dy * seg.total();
for (auto &ev : events[i])
seg.update(ev.l, ev.r, ev.delta);
yPrev = y;
}
return ans;
}

int main()
{
cout << solve(N) << "\n";
return 0;
}