2021 ccpc 新疆省赛
Problem A. balloon
将两个序列从小到大排序，从跳跃高度较低的人考虑到跳跃高度较高的人，使用一个变量 \(\text{now}\) 维护当前取到的气球，若气球能取就取，不能取就换成下一个人。
时间复杂度的瓶颈在于排序。
#include <iostream>
#include <algorithm>
#include <cstring>
#include <cstdio>
using namespace std;
const int N = 2e5 + 10;
int n, m;
struct Node {
int h;
int id;
} e[N];
int ans[N];
int h[N];
bool cmp(Node a, Node b) {
return a.h < b.h;
}
void solve() {
cin >> n >> m;
for(int i = 1; i <= n ; i ++) {
cin >> e[i].h;
e[i].id = i;
}
sort(e + 1, e + 1 + n, cmp);
for(int i = 1; i <= m; i ++)
cin >> h[i];
int now = 1;
sort(h + 1, h + 1 + m);
for(int i = 1; i <= n; i ++) {
int num = 0;
while(now <= m && h[now] <= e[i].h)
num ++, now ++;
ans[e[i].id] = num ;
}
for(int i = 1; i <= n; i ++)
cout << ans[i] << endl;
return ;
}
signed main() {
int t = 1;
while(t--)
solve();
return 0;
}
Problem B. sophistry
可以从后往前倒序 dp。
设 \(f_i\) 表示：考虑到了第 \(i\) 天到第 \(n\) 天时，造成的最大伤害。
- 当 \(a_i \leq m\) 时，有转移：
- 当 \(a_i > m\) 时，有转移：
最后答案即为 \(f_1\)，时间复杂度 \(\mathcal{O}(n)\)。
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxn=2e5+10;
ll dp[maxn];
ll a[maxn];
int main()
{
int n,d,m;
scanf("%d%d%d",&n,&d,&m);
for(int i=1;i<=n;i++) scanf("%lld",&a[i]);
for(int i=n;i>=1;i--)
{
if(a[i]>m)
{
dp[i]=max(dp[i+d+1]+a[i],dp[i+1]);
}
else
{
dp[i]=dp[i+1]+a[i];
}
}
cout<<dp[1]<<endl;
}
Problem C. bomb
先考虑一下没有环的情况。
结论：答案为有向图的最长链长度。
证明：
-
由于最长链上的点两两不能在同一轮染色。所以最优解 \(\ge\) 最长链长度。
-
令 \(dist_x\) 表示从 \(x\) 开始的最长链长度，可以在第 \(i\) 轮中对所有 \(dist_x = i\) 的点 \(x\) 进行染色，恰好可以得到一个轮数为最长链长度的方案。所以最优解 \(\leq\) 最长链长度。
故答案为有向图的最长链长度。
再考虑一下有环的情况。
考虑使用 tarjan 将这张有向图中的所有强连通分量求出，那么对于任意一个强连通分量，如果在一轮染色中对这个强连通分量中的其中一个点进行染色，那么这个强连通分量中的其他点、以及它所能到达的强连通分量里的所有点都不能染色。
我们发现这与题目的定义类似，对于一个大小为 \(s\) 的强连通分量，可以将其缩为一个需要被染色 \(s\) 次的点。这样就得到了一张新图，答案即为这张新图中的最长链，使用拓扑排序即可解决。
时间复杂度 \(\mathcal{O}(n + m)\)。
#include <iostream>
#include <algorithm>
#include <cstring>
#include <cstdio>
#include <map>
#include <queue>
using namespace std;
typedef long long ll;
const int N = 2e6 + 10;
int e[N], h[N], ne[N], idx;
int cnt;
int stk[N];
int top;
int in_stk[N];
int dfn[N];
int low[N];
int scc_cnt;
int in_scc_num[N];
int scc_size[N];
int n, m;
void add(int a, int b) {
e[idx] = b;
ne[idx] = h[a];
h[a] = idx ++;
}
void tarjan(int u) {
dfn[u] = low[u] = ++ cnt;
stk[++ top] = u;
in_stk[u] = 1;
for(int i = h[u]; i != -1; i = ne[i]) {
int v = e[i];
if(!dfn[v]) {
tarjan(v);
low[u] = min(low[u], low[v]);
} else if(in_stk[v]) {
low[u] = min(low[u], dfn[v]);
}
}
if(dfn[u] == low[u]) {
scc_cnt ++;
int x;
do {
x = stk[top --];
in_stk[x] = 0;
in_scc_num[x] = scc_cnt;
scc_size[scc_cnt] ++;
} while(x != u);
}
}
int h1[N], e1[N], ne1[N], idx1;
int du[N];
void add1(int a, int b) {
e1[idx1] = b;
ne1[idx1] = h1[a];
h1[a] = idx1 ++;
du[b] ++;
}
ll dist[N];
void topsort() {
queue<int>q;
for(int i = 1; i <= scc_cnt; i ++)
if(du[i] == 0)
q.push(i), dist[i] = scc_size[i];
while(q.size()) {
int u = q.front();
q.pop();
for(int i = h1[u]; i != -1; i = ne1[i]) {
int v = e1[i];
dist[v] = max(dist[v], dist[u] + scc_size[v]);
if(--du[v] == 0)
q.push(v);
}
}
}
void solve() {
cin >> n >> m;
idx = idx1 = 0;
memset(h, -1, sizeof h);
memset(h1, -1, sizeof h1);
for(int i = 1; i <= m; i ++) {
int x, y;
cin >> x >> y;
add(x, y);
}
for(int i = 1; i <= n ; i++) {
if(!dfn[i])
tarjan(i);
}
for(int u = 1; u <= n; u ++) {
for(int i = h[u]; i != -1; i = ne[i]) {
int v = e[i];
if(in_scc_num[u] == in_scc_num[v])
continue;
add1(in_scc_num[u], in_scc_num[v]);
}
}
topsort();
ll ans = 0;
for(int i = 1; i <= scc_cnt; i ++)
ans = max(ans, dist[i]);
cout << ans << "\n";
return ;
}
signed main() {
int t = 1;
while(t--)
solve();
return 0;
}
Problem D. maxsum
题目要求的就是前 \(w\) 大子段和。
因为 \(a_i\) 为非负整数，所以如果当前最大的是 \([l,r]\) 子段，那么易得 \([l,r+1]\) 子段和 \([l-1,r]\) 子段一定之前就已经取出。
最大的子段一定为 \([1,n]\)，所以首先将 \([1,n]\) 加入堆。设堆顶为 \([l,r]\)，则 \([l+1,r]\) 和 \([l,r-1]\) 可以成为备选答案，加入堆并用 map 判重即可。
时间复杂度 \(O(w \log n)\)。
#include <iostream>
#include <algorithm>
#include <cstring>
#include <cstdio>
#include <map>
#include <queue>
using namespace std;
typedef long long ll;
const int N = 1e5 + 10;
typedef pair<int, int> pii;
int n, w;
ll a[N];
ll sum;
map<pii, bool>mp;
struct node {
int l, r;
ll sum;
friend bool operator < (node x, node y) {
return x.sum < y.sum;
}
};
priority_queue<node>q;
void solve() {
cin >> n >> w;
for(int i = 1; i <= n; i++)
cin >> a[i], sum += a[i];
q.push({1, n, sum});
mp[ {1, n}] = 1;
while(w) {
auto now = q.top();
q.pop();
w--;
cout << now.sum << " ";
if(now.l != now.r) {
if(!mp[ {now.l + 1, now.r}])
q.push({now.l + 1, now.r, now.sum - a[now.l]});
if(!mp[ {now.l, now.r - 1}])
q.push({now.l, now.r - 1, now.sum - a[now.r]});
mp[ {now.l + 1, now.r}] = mp[ {now.l, now.r - 1}] = 1;
}
}
return ;
}
signed main() {
int t = 1;
while(t--)
solve();
return 0;
}
Problem E. array
令 \(A, B\) 分别为两个数组的最大值，则 \(c_i ≥ \max(A, B)\)。
另一方面，如果 \(c_i > \max(A, B)\)，则对应的 \(a_i\) 和 \(b_j\) 都非零。
根据题意，两个数组最多只有 \(5000\) 个数非零，暴力计算即可。
时间复杂度 \(\mathcal O(n + 5000^2)\)。
#include <iostream>
#include <algorithm>
#include <cstring>
#include <cstdio>
#include <map>
#include <queue>
using namespace std;
typedef long long ll;
const int N = 2e5 + 10;
const int mod = 1e9 + 7;
int a[N];
int b[N];
int c[N];
int n;
int nexta[N];
int nextb[N];
int maxv = 0;
void solve() {
cin >> n;
for(int i = 0; i <= n - 1; i ++)
cin >> a[i], maxv = max(maxv, a[i]);
for(int i = 0; i <= n - 1; i ++)
cin >> b[i], maxv = max(maxv, b[i]);
for(int i = 0; i < n; i ++)
c[i] = maxv;
nexta[n - 1] = nextb[n - 1] = n;
for(int i = n - 1; i >= 1; i --) {
if(a[i])
nexta[i - 1] = i;
else
nexta[i - 1] = nexta[i];
}
for(int i = n - 1; i >= 1; i --) {
if(b[i])
nextb[i - 1] = i;
else
nextb[i - 1] = nextb[i];
}
for(int i = 0 ; i < n; i = nexta[i])
for(int j = 0; j < n; j = nextb[j])
c[(i + j) % n] = max(c[(i + j) % n], a[i] + b[j]);
for(int i = 0; i < n ; i++)
cout << c[i] << " ";
cout << "\n";
return ;
}
signed main() {
int t = 1;
while(t--)
solve();
return 0;
}
Problem F. fare
设 \(u\) 是树的一个结点，设 \(num1(u),num2(u),dp1(u)\) 分别表示以 \(u\) 为根的子树内所有参加活动的选手到 \(u\) 的距离的零次方和，一次方和，二次方和，再设 \(tmp1(u),tmp2(u),dp2(u)\) 分别表示除去以 \(u\) 为根的子树内所有参加活动的选手到 \(u\) 的距离的零次方和，一次方和，二次方和。
遍历整棵树就可以求出所有结点的 \(num1,num2,dp1\)，再次遍历整棵树，我们使用换根 \(dp\) 的方法，就可以求出所有结点的 \(tmp1,tmp2,dp2\)，最后的答案便为
\(\min\limits_{1 \leq u \leq n} \{ dp1(u)+dp2(u) \}\)
时间复杂度 \(O(n)\)。
#include <iostream>
#include <algorithm>
#include <cstring>
#include <cstdio>
#include <map>
#include <queue>
using namespace std;
typedef long long ll;
const int N = 4e5 + 10;
const int mod = 1e9 + 7;
int h[N], e[N], ne[N], idx, w[N];
ll c[N];
ll num1[N];// num
ll num2[N];// edge*num
ll dp1[N];// 1 as root
ll tmp1[N];
ll tmp2[N];
ll dp2[N];//change root
ll n;
void add(int a, int b, int c) {
e[idx] = b;
w[idx] = c;
ne[idx] = h[a];
h[a] = idx ++;
}
void dfs1(int u, int fa) {
num1[u] = c[u];
for(int i = h[u]; i != -1; i = ne[i]) {
int v = e[i];
if(v == fa)
continue;
dfs1(v, u);
num1[u] += num1[v];
num2[u] += num2[v] + w[i] * num1[v];
dp1[u] += dp1[v] + 2 * w[i] * num2[v] + w[i] * w[i] * num1[v];
}
}
void dfs2(int u, int fa) {
for(int i = h[u]; i != -1; i = ne[i]) {
int v = e[i];
if(v == fa)
continue;
tmp1[v] = num1[u] - num1[v] + tmp1[u];
tmp2[v] = num2[u] - num1[v] * w[i] - num2[v] + tmp2[u] + tmp1[v] * w[i];
dp2[v] = dp1[u] - dp1[v] - 2 * w[i] * num2[v] - w[i] * w[i] * num1[v] + dp2[u] +
2 * w[i] * (num2[u] - num1[v] * w[i] - num2[v] + tmp2[u] )
+ w[i] * w[i] * tmp1[v];
dfs2(v, u);
}
}
void solve() {
cin >> n;
for(int i = 1; i <= n; i ++)
cin >> c[i];
idx = 0;
memset(h, -1 , sizeof h);
for(int i = 1; i < n; i ++) {
int a, b, c;
cin >> a >> b >> c;
add(a, b, c);
add(b, a, c);
}
dfs1(1, 0);
// cout << dp1[1] << " ";
dfs2(1, 0);
ll ans = 1e18;
for(int i = 1; i <= n; i++)
ans = min(ans , dp1[i] + dp2[i]);
cout << ans << "\n";
return ;
}
signed main() {
int t = 1;
while(t--)
solve();
return 0;
}
Problem G. road
注意到异或不进位，每位之间是相互独立的，于是我们可以考虑计算一下每一位的贡献。
比如说我们现在要计算一下第 \(\text{index}\) 位的贡献。
设原图中的邻接矩阵为 \(A\)，将第 \(\text{index}\) 位提取出来可以得到一个 \(0 / 1\) 矩阵 \(A'\)。我们让点在新的邻接矩阵上走，注意到按位异或的性质，一条非简单路径只有经过了奇数个边，才会对答案有 \(2^{\text{index}}\) 的贡献。
考虑 dp。
设 \(f(i, x, 0 / 1)\) 表示：有多少条 \(S\) 到 \(x\) 的边数为 \(i\) 的路径，且经过了 偶数 / 奇数 个 \(1\)。
则有状态转移方程：
\(
f(i, v, 0) = \sum\limits_{(u, v) \in E, a'_{u, v} = 0} f(i - 1, u, 0) + \sum\limits_{(u, v) \in E, a'_{u, v} = 1} f(i - 1, u, 1) \\
f(i, v, 1) = \sum\limits_{(u, v) \in E, a'_{u, v} = 0} f(i - 1, u, 1) + \sum\limits_{(u, v) \in E, a'_{u, v} = 1} f(i - 1, u, 0) \\
\)
则第 \(\text{index}\) 位对答案的贡献即为 \(f_{k, T, 1} \ast 2^{\text{index}}\)。
直接做的话显然会 TLE，但是这个转移方程对于所有的阶段 \(i\) 转移都一样，于是矩阵快速幂加速即可。
单次回答询问的时间复杂度是 \(\mathcal{O}(n^3 \log^2 \text{size})\) 的，其中 \(\text{size}\) 表示值域大小。
还不够优秀。
因为回答的询问只关注单源信息，注意到 \(1 \times n\) 的矩阵和 \(n \times n\) 的矩阵相乘的复杂度是 \(\mathcal{O}(n^2)\) 的，\(n \times n\) 的矩阵和 \(n \times n\) 的矩阵相乘的复杂度是 \(\mathcal{O}(n^3)\) 的。
于是我们可以先把 \(A^1, A^2, A^4, A^8, ...\) 都预处理出来。
然后将 \(k\) 二进制分解，若 \(k\) 的第 \(i\) 位上为 \(1\)，则令 \(f = f \ast A^{2^i}\)。
这样做的话，单次询问的时间复杂度就成功降到了 \(\mathcal{O}(n^2 \log^2 n)\)。
总时间复杂度 \(\mathcal{O}(n^3 \log^2 \text{size} + Q \ast n^2 \log^2 \text{size})\)。
#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
inline int read() {
int x = 0, f = 1; char s = getchar();
while (s < '0' || s > '9') { if (s == '-') f = -f; s = getchar(); }
while (s >= '0' && s <= '9') { x = x * 10 + s - '0'; s = getchar(); }
return x * f;
}
const int N = 45;
const int mod = 1e9 + 7;
int n, m, Q;
int a[N][N];
int f[31][31][N][N][2];
void mul(int f[N][N][2], int a[N][N][2]) {
static int c[N][N][2]; memset(c, 0, sizeof(c));
for (int i = 1; i <= n; i ++)
for (int j = 1; j <= n; j ++)
for (int k = 1; k <= n; k ++) {
c[i][j][0] = (c[i][j][0] + 1ll * a[i][k][0] * a[k][j][0]) % mod;
c[i][j][0] = (c[i][j][0] + 1ll * a[i][k][1] * a[k][j][1]) % mod;
c[i][j][1] = (c[i][j][1] + 1ll * a[i][k][0] * a[k][j][1]) % mod;
c[i][j][1] = (c[i][j][1] + 1ll * a[i][k][1] * a[k][j][0]) % mod;
}
memcpy(f, c, sizeof(c));
}
void mulstar(int a[N][2], int b[N][N][2]) {
static int c[N][2]; memset(c, 0, sizeof(c));
for (int j = 1; j <= n; j ++)
for (int k = 1; k <= n; k ++) {
c[j][0] = (c[j][0] + 1ll * a[k][0] * b[k][j][0]) % mod;
c[j][0] = (c[j][0] + 1ll * a[k][1] * b[k][j][1]) % mod;
c[j][1] = (c[j][1] + 1ll * a[k][0] * b[k][j][1]) % mod;
c[j][1] = (c[j][1] + 1ll * a[k][1] * b[k][j][0]) % mod;
}
memcpy(a, c, sizeof(c));
}
void calc(int index) {
for (int u = 1; u <= n; u ++)
for (int v = 1; v <= n; v ++) {
if (a[u][v] == 0x3f3f3f3f) continue;
int w = a[u][v] >> index & 1;
f[index][0][u][v][w] = 1;
}
for (int i = 1; i <= 30; i ++)
mul(f[index][i], f[index][i - 1]);
}
void work() {
int S = read(), T = read(), k = read() - 1;
int ans = 0;
for (int index = 0; index <= 30; index ++) {
int d[N][2]; memset(d, 0, sizeof(d));
for (int i = 1; i <= n; i ++) {
if (a[S][i] == 0x3f3f3f3f) continue;
int w = a[S][i] >> index & 1;
d[i][w] = 1;
}
for (int i = 0; i <= 30; i ++)
if (k >> i & 1) mulstar(d, f[index][i]);
ans = (ans + 1ll * (1 << index) * d[T][1]) % mod;
}
printf("%d\n", ans);
}
int main() {
n = read(), m = read();
memset(a, 0x3f, sizeof(a));
for (int i = 1; i <= m; i ++) {
int u = read(), v = read(), w = read();
a[u][v] = a[v][u] = w;
}
for (int index = 0; index <= 30; index ++)
calc(index);
Q = read();
while (Q --) work();
return 0;
}
Problem H. xor
将 \(a_i\) 和 \(a_j\) 按二进制每一位展开：
\(a_i =\sum\limits_{x=0}^{29}2^xa_{i, x}\)，\(a_j =\sum\limits_{y=0}^{29}2^ya_{j, y}\)
\(a_i ⊕ a_j=\sum\limits_{k=0}^{29} 2^k[a_{i, k} \not=a_{j, k}]\)
将原式展开，
\(\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n} (a_i ⊕a_j)^2\)
\(=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}(\sum\limits_{x=0}^{29} 2^x[a_{i, x}\not=a_{j, x}])(\sum\limits_{y=0}^{29} 2^y[a_{i, y}\not=a_{j, y}])\)
\(=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}\sum\limits_{x=0}^{29}\sum\limits_{y=0}^{29} 2^{x+y}[a_{i, x}\not=a_{j, x}][a_{i, y}\not=a_{j, y}]\)
\(=\sum\limits_{x=0}^{29}\sum\limits_{y=0}^{29} 2^{x+y}\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}[a_{i, x}\not=a_{j, x}][a_{i, y}\not=a_{j, y}]\)
枚举 \(x, y, a_{i,x}, a_{j,x}, a_{i,y}, a_{j,y}\) 的取值，乘以对应的方案数（即有多少个 \(i\) 满足第 \(x\) 位为 \(a_{i,x}\) 且第 \(y\) 位为 \(a_{i,y}\)）即可。
预处理出 \(f_{i,j},x,y\) 表示有多少个数满足第 \(i\) 位为 \(x\) 且第 \(j\) 位为 \(y\)，则方案数可以查表 \(O(1)\) 得到。
时间复杂度 \(O(n \log^2 a)\)。
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
const int N = 50010;
const int mod = 1e9 + 7;
typedef long long ll;
int n;
ll a[N];
ll inv[N];
ll f[40][40][2][2];
signed main() {
/*
5
12133
544123
47778
69885
11234
*/
cin >> n;
for(int i = 1; i <= n; i ++)
cin >> a[i];
// ll res = 0 ;
// for(int i = 1; i <= n; i ++)
// for(int j = 1; j <= n; j ++)
// res = (res + (a[i] ^ a[j]) * (a[i] ^ a[j]) % mod) % mod;
// cout << res << endl;
inv[0] = 1;
for(int i = 1; i <= 63; i ++)
inv[i] = 1ll * 2 * inv[i - 1] % mod;
for(int i = 1; i <= n; i ++)
for(int j = 0; j <= 31; j ++)
for(int k = 0; k <= 31; k ++)
f[j][k][a[i] >> j & 1][a[i] >> k & 1] ++;
ll ans = 0;
for(int i = 0; i <= 31; i ++)
for(int j = 0; j <= 31; j ++)
for(int k = 0; k <= 1; k ++)
for(int z = 0; z <= 1; z ++)
ans = (ans + 1ll * f[i][j][k][z] * f[i][j][k ^ 1][z ^ 1] % mod * inv[i + j] % mod) % mod;
cout << ans << '\n';
return 0;
}
Problem I. Fibonacci sequence
存在递推式
\(\sqrt{3 + f_n f_{n - 1}} = \sqrt{3 + f_{n - 1} f_{n - 2}} + f_{n - 1}\)
证明如下
\((\sqrt{3 + f_{n - 1} f_{n - 2}} + f_{n - 1})^2\)
\(=3 + f_{n - 1} f_{n - 2} + 2 f_{n - 1} \sqrt{3 + f_{n - 1} f_{n - 2}} + f_{n - 1}^2\)
\(=3 + (f_{n - 1} + f_{n - 2} + 2 \sqrt{3 + f_{n - 1}f_{n - 2}}) f_{n - 1}\)
\(=3 + f_n f_{n - 1}\)
于是可以构建矩阵来优化递推，但是考虑到有 \(T(T \leq 10^6)\) 组，我们可以考虑光速幂。
具体来说，我们想要知道 \(G^t\)（\(G\) 为转移矩阵），我们预处理出矩阵 \(G\) 的 \(G^0, G^1, G^2 \cdots, G^s, G^{2s},G^{3s}, G^{⌈\frac{MAX}{s}⌉s}\)（其中，\(s=\sqrt{MAX}\)）这样，我们就可以 \(O(3^3)\)（即矩阵乘法的复杂度） 求出一个 \(G^t (t = ks + b)\)，进而优化复杂度。
#include <bits/stdc++.h>
#define int long long
const int N = 3 + 10;
const int RN = 1e6 + 10;
using namespace std;
struct matrix
{
int l, r;
long long val[N][N];
} a, G, t[100000], t2[100000];
int x[RN];
int G0, G1, M, n, ans = 1;
inline matrix mul(matrix A, matrix B)
{
matrix C; C.l = A.l, C.r = B.r;
memset(C.val, 0, sizeof(C.val));
for (int i = 1; i <= C.l; ++ i)
for (int j = 1; j <= C.r; ++ j)
for (int k = 1; k <= A.r; ++ k)
C.val[i][j] = (C.val[i][j] + A.val[i][k] * B.val[k][j]) % M;
return C;
}
inline matrix pow(matrix A, int b)
{
matrix res = A; b --;
for (; b; b >>= 1)
{
if (b & 1) res = mul(res, A);
A = mul(A, A);
}
return res;
}
inline matrix solve(int b)
{
a.l = 3, a.r = 3;
a.val[1][1] = 1, a.val[1][2] = 1, a.val[1][3] = 1;
a.val[2][1] = 1, a.val[2][2] = 0, a.val[3][3] = 0;
a.val[3][1] = 2, a.val[3][2] = 0, a.val[3][3] = 1;
return pow(a, b);
}
int Max;
signed main()
{
cin >> G0 >> G1 >> M >> n;
G.l = 1, G.r = 3;
G.val[1][1] = G1, G.val[1][2] = G0, G.val[1][3] = (int)sqrt(3 + G1 * G0) % M;
a.l = 3, a.r = 3;
a.val[1][1] = 1, a.val[1][2] = 1, a.val[1][3] = 1;
a.val[2][1] = 1, a.val[2][2] = 0, a.val[3][3] = 0;
a.val[3][1] = 2, a.val[3][2] = 0, a.val[3][3] = 1;
for (int i = 1; i <= n; ++ i)
{
cin >> x[i];
Max = max(Max, x[i]);
}
Max = sqrt(Max) + 1;
t[1] = solve(1);
for (int i = 2; i <= Max; ++ i)
t[i] = mul(t[i - 1], t[1]);
t2[1] = solve(Max);
for (int i = 2; i <= Max; ++ i)
t2[i] = mul(t2[i - 1], t2[1]);
for (int i = 1; i <= n; ++ i)
{
if (x[i] == 0)
ans = ans * G0 % M;
else if (x[i] == 1)
ans = ans * G1 % M;
else
{
x[i] -= 1;
int s = x[i] / Max;
matrix c1 = t2[s], c2 = t[x[i] - Max * s];
matrix res = mul(c1, c2);
if (s == 0) res = c2;
if (x[i] - Max * s == 0) res = c1;
ans = ans * mul(G, res).val[1][1] % M;
}
}
cout << ans << endl;
return 0;
}
Problem J. mesh
直接上莫比乌斯反演：
\(
\begin{aligned}
\sum\limits_{i = 1}^n \sum\limits_{j = 1}^m \varphi(\gcd(i, j))
& = \sum\limits_{d = 1} \varphi(d) \cdot \sum\limits_{i = 1}^n \sum\limits_{j = 1}^m [\gcd(i, j) = d]\\
& = \sum\limits_{d = 1} \varphi(d) \cdot \sum\limits_{i = 1}^{\left\lfloor\frac{n}{d}\right\rfloor} \sum\limits_{j = 1}^{\left\lfloor\frac{m}{d}\right\rfloor} [\gcd(i, j) = 1]\\
& = \sum\limits_{d = 1} \varphi(d) \cdot \sum\limits_{i = 1}^{\left\lfloor\frac{n}{d}\right\rfloor} \sum\limits_{j = 1}^{\left\lfloor\frac{m}{d}\right\rfloor} \sum\limits_{p \mid \gcd(i, j)} \mu(p) \\
& = \sum\limits_{d = 1} \varphi(d) \cdot \sum\limits_{p = 1} \mu(p) \cdot \left\lfloor\frac{n}{pd}\right\rfloor \cdot \left\lfloor\frac{m}{pd}\right\rfloor
\end{aligned}
\)
记 \(T = pd\)，则式子化简为：
\(
\sum\limits_{T = 1}^{\min(n, m)} \left\lfloor\frac{n}{T}\right\rfloor \cdot \left\lfloor\frac{m}{T}\right\rfloor \cdot \left(\sum\limits_{d \mid T} \varphi(d) \cdot \mu\left(\frac{T}{d}\right) \right)
\)
设 \(g = \varphi \ast \mu\)（这里的 \(\ast\) 是狄利克雷卷积，下文同），则式子化简为：
\(
\sum\limits_{T = 1}^{\min(n, m)} \left\lfloor\frac{n}{T}\right\rfloor \cdot \left\lfloor\frac{m}{T}\right\rfloor \cdot g(T)
\)
注意到函数 \(g\) 显然为积性函数，可以使用线性筛将 \(g\) 的前缀和筛出。
配合数论分块即可做到 \(\mathcal{O}(n)\) 预处理、\(\mathcal{O}(\sqrt{n})\) 回答询问。但还不够优秀。
可以考虑杜教筛，设：
\(
S(n) = \sum\limits_{i = 1}^n g(i)
\)
记 \(d = 1 \ast 1\)，其中函数 \(1(n) = 1\)，显然函数 \(d\) 为约数个数函数。
那么可以得到 \(S(n)\) 有关 \(S\left(\left\lfloor\frac{n}{d}\right\rfloor\right)\) 的递推式：
\(
d(1)S(n) = \sum\limits_{i = 1}^n (g \ast d)(i) - \sum\limits_{i = 2}^n d(i) \cdot S\left(\left\lfloor\frac{n}{i}\right\rfloor\right) \\
\)
化简：
\(
d(1)S(n) = \sum\limits_{i = 1}^n (\varphi \ast 1 \ast \mu \ast 1)(i) - \sum\limits_{i = 2}^n d(i) \cdot S\left(\left\lfloor\frac{n}{i}\right\rfloor\right) \\
d(1)S(n) = \frac{n \cdot (n + 1)}{2} - \sum\limits_{i = 2}^n d(i) \cdot S\left(\left\lfloor\frac{n}{i}\right\rfloor\right)
\)
该式显然需要数论分块来做。取一个阈值 \(K\)，分段处理：
- 对于 \(n \leq K\) 的部分，使用线性筛求出每个 \(S(n), \sum_{i = 1}^n d(i)\)。
- 对于 \(n > K\) 的部分，使用数论分块求解 \(\sum_{i = 1}^n d(i)\)，使用上述递推式求解 \(S(i)\)。
可以使用积分近似证明当 \(K = \mathcal{O}(n^{\frac{2}{3}})\) 时，取到杜教筛的理论最优复杂度 \(\mathcal{O}(n^{\frac{2}{3}})\)。
#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <map>
#include <ctime>
using namespace std;
const int N = 1000100, SIZE = 1e6;
const int mod = 1e9 + 7;
int n, m;
int k, prime[N], low[N];
int h[N], d[N];
int Sd[N], Sh[N];
void GetPrimes(int N) {
h[1] = 1, d[1] = 1;
for (int i = 2; i <= N; i ++) {
if (!low[i]) {
prime[++ k] = i;
low[i] = i;
d[i] = 2;
h[i] = i - 2;
if (1ll * i * i <= N) {
low[i * i] = i * i;
d[i * i] = 3;
h[i * i] = (i * i - i) - (i - 1);
long long v = 1ll * i * i * i;
int c = 3;
while (v <= N) {
low[v] = v;
d[v] = c + 1;
h[v] = h[v / i] * i;
v *= i, c ++;
}
}
}
for (int j = 1; j <= k; j ++) {
if (prime[j] > N / i) break;
if (i % prime[j] == 0)
low[i * prime[j]] = low[i] * prime[j];
else
low[i * prime[j]] = prime[j];
d[i * prime[j]] = d[(i * prime[j]) / low[i * prime[j]]] * d[low[i * prime[j]]];
h[i * prime[j]] = h[(i * prime[j]) / low[i * prime[j]]] * h[low[i * prime[j]]];
if (i % prime[j] == 0) break;
}
}
for (int i = 1; i <= N; i ++)
Sd[i] = (Sd[i - 1] + d[i]) % mod,
Sh[i] = (Sh[i - 1] + h[i]) % mod;
}
map<int, int> Gd, Gh;
int sig(int n) {
if (n <= SIZE) return Sd[n];
if (Gd[n]) return Gd[n];
int ans = 0;
for (int x = 1, Nx; x <= n; x = Nx + 1) {
Nx = n / (n / x);
ans = (ans + 1ll * (Nx - x + 1) * (n / x)) % mod;
}
return Gd[n] = ans;
}
int S(int n) {
if (n <= SIZE) return Sh[n];
if (Gh[n]) return Gh[n];
int ans = (1ll * n * (n + 1) / 2) % mod;
int sub = 0;
for (int x = 2, Nx, cur; x <= n; x = Nx + 1) {
Nx = n / (n / x);
cur = ((sig(Nx) - sig(x - 1)) % mod + mod) % mod;
cur = 1ll * cur * S(n / x) % mod;
sub = (sub + cur) % mod;
}
ans = ((ans - sub) % mod + mod) % mod;
return Gh[n] = ans;
}
int main() {
cin >> n >> m;
GetPrimes(SIZE);
int ans = 0;
for (int x = 1, Nx, cur; x <= min(n, m); x = Nx + 1) {
Nx = min(n / (n / x), m / (m / x));
cur = ((S(Nx) - S(x - 1)) % mod + mod) % mod;
cur = 1ll * cur * (n / x) % mod;
cur = 1ll * cur * (m / x) % mod;
ans = (ans + cur) % mod;
}
cout << ans << endl;
return 0;
}
