Problem A. balloon
将两个序列从小到大排序，从跳跃高度较低的人考虑到跳跃高度较高的人，使用一个变量 ( ext{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 << "
";
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 << "
";
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)，最后的答案便为
(minlimits_{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 << "
";
return ;
}
signed main() {
int t = 1;
while(t--)
solve();
return 0;
}
Problem G. road
注意到异或不进位，每位之间是相互独立的，于是我们可以考虑计算一下每一位的贡献。
比如说我们现在要计算一下第 ( ext{index}) 位的贡献。
设原图中的邻接矩阵为 (A)，将第 ( ext{index}) 位提取出来可以得到一个 (0 / 1) 矩阵 (A')。我们让点在新的邻接矩阵上走，注意到按位异或的性质，一条非简单路径只有经过了奇数个边，才会对答案有 (2^{ ext{index}}) 的贡献。
考虑 dp。
设 (f(i, x, 0 / 1)) 表示：有多少条 (S) 到 (x) 的边数为 (i) 的路径，且经过了 偶数 / 奇数 个 (1)。
则有状态转移方程：
(
f(i, v, 0) = sumlimits_{(u, v) in E, a'_{u, v} = 0} f(i - 1, u, 0) + sumlimits_{(u, v) in E, a'_{u, v} = 1} f(i - 1, u, 1) \
f(i, v, 1) = sumlimits_{(u, v) in E, a'_{u, v} = 0} f(i - 1, u, 1) + sumlimits_{(u, v) in E, a'_{u, v} = 1} f(i - 1, u, 0) \
)
则第 ( ext{index}) 位对答案的贡献即为 (f_{k, T, 1} ast 2^{ ext{index}})。
直接做的话显然会 TLE，但是这个转移方程对于所有的阶段 (i) 转移都一样，于是矩阵快速幂加速即可。
单次回答询问的时间复杂度是 (mathcal{O}(n^3 log^2 ext{size})) 的，其中 ( ext{size}) 表示值域大小。
还不够优秀。
因为回答的询问只关注单源信息，注意到 (1 imes n) 的矩阵和 (n imes n) 的矩阵相乘的复杂度是 (mathcal{O}(n^2)) 的，(n imes n) 的矩阵和 (n imes 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 ext{size} + Q ast n^2 log^2 ext{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
", 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 =sumlimits_{x=0}^{29}2^xa_{i, x})，(a_j =sumlimits_{y=0}^{29}2^ya_{j, y})
(a_i ⊕ a_j=sumlimits_{k=0}^{29} 2^k[a_{i, k} ot=a_{j, k}])
将原式展开，
(sumlimits_{i=1}^{n}sumlimits_{j=1}^{n} (a_i ⊕a_j)^2)
(=sumlimits_{i=1}^{n}sumlimits_{j=1}^{n}(sumlimits_{x=0}^{29} 2^x[a_{i, x} ot=a_{j, x}])(sumlimits_{y=0}^{29} 2^y[a_{i, y} ot=a_{j, y}]))
(=sumlimits_{i=1}^{n}sumlimits_{j=1}^{n}sumlimits_{x=0}^{29}sumlimits_{y=0}^{29} 2^{x+y}[a_{i, x} ot=a_{j, x}][a_{i, y} ot=a_{j, y}])
(=sumlimits_{x=0}^{29}sumlimits_{y=0}^{29} 2^{x+y}sumlimits_{i=1}^{n}sumlimits_{j=1}^{n}[a_{i, x} ot=a_{j, x}][a_{i, y} ot=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 << '
';
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
直接上莫比乌斯反演：
(
egin{aligned}
sumlimits_{i = 1}^n sumlimits_{j = 1}^m varphi(gcd(i, j))
& = sumlimits_{d = 1} varphi(d) cdot sumlimits_{i = 1}^n sumlimits_{j = 1}^m [gcd(i, j) = d]\
& = sumlimits_{d = 1} varphi(d) cdot sumlimits_{i = 1}^{leftlfloorfrac{n}{d}
ight
floor} sumlimits_{j = 1}^{leftlfloorfrac{m}{d}
ight
floor} [gcd(i, j) = 1]\
& = sumlimits_{d = 1} varphi(d) cdot sumlimits_{i = 1}^{leftlfloorfrac{n}{d}
ight
floor} sumlimits_{j = 1}^{leftlfloorfrac{m}{d}
ight
floor} sumlimits_{p mid gcd(i, j)} mu(p) \
& = sumlimits_{d = 1} varphi(d) cdot sumlimits_{p = 1} mu(p) cdot leftlfloorfrac{n}{pd}
ight
floor cdot leftlfloorfrac{m}{pd}
ight
floor
end{aligned}
)
记 (T = pd)，则式子化简为：
(
sumlimits_{T = 1}^{min(n, m)} leftlfloorfrac{n}{T}
ight
floor cdot leftlfloorfrac{m}{T}
ight
floor cdot left(sumlimits_{d mid T} varphi(d) cdot muleft(frac{T}{d}
ight)
ight)
)
设 (g = varphi ast mu)（这里的 (ast) 是狄利克雷卷积，下文同），则式子化简为：
(
sumlimits_{T = 1}^{min(n, m)} leftlfloorfrac{n}{T}
ight
floor cdot leftlfloorfrac{m}{T}
ight
floor cdot g(T)
)
注意到函数 (g) 显然为积性函数，可以使用线性筛将 (g) 的前缀和筛出。
配合数论分块即可做到 (mathcal{O}(n)) 预处理、(mathcal{O}(sqrt{n})) 回答询问。但还不够优秀。
可以考虑杜教筛，设：
(
S(n) = sumlimits_{i = 1}^n g(i)
)
记 (d = 1 ast 1)，其中函数 (1(n) = 1)，显然函数 (d) 为约数个数函数。
那么可以得到 (S(n)) 有关 (Sleft(leftlfloorfrac{n}{d}
ight
floor
ight)) 的递推式：
(
d(1)S(n) = sumlimits_{i = 1}^n (g ast d)(i) - sumlimits_{i = 2}^n d(i) cdot Sleft(leftlfloorfrac{n}{i}
ight
floor
ight) \
)
化简：
(
d(1)S(n) = sumlimits_{i = 1}^n (varphi ast 1 ast mu ast 1)(i) - sumlimits_{i = 2}^n d(i) cdot Sleft(leftlfloorfrac{n}{i}
ight
floor
ight) \
d(1)S(n) = frac{n cdot (n + 1)}{2} - sumlimits_{i = 2}^n d(i) cdot Sleft(leftlfloorfrac{n}{i}
ight
floor
ight)
)
该式显然需要数论分块来做。取一个阈值 (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;
}