卡常递推---题解 P4619 【[SDOI2018]旧试题】
Vocalise
·
·
题解
非常神仙的题目，三元环做法已经非常神仙，但是 @big_news 的做法让这题更加神仙。
这个做法大概没什么人提到吧。
首先有
d(ijk)=\sum_{x|i}\sum_{y|j}\sum_{z|k}[x\perp y][y\perp z][x\perp z]
非常显然有 d(ijk)=\sum_{x|ijk}1，分开 i,j,k 数会重复，于是我们让 ijk 相同的质因子在 x,y,z 中只出现一个。
P3327 [SDOI2015]约数个数和 中有介绍。
&\ \ \ \ \sum_{x=1}^a\sum_{y=1}^b\sum_{z=1}^cd(xyz) \\
&=\sum_{x=1}^a\sum_{y=1}^b\sum_{z=1}^c\sum_{i|x}\sum_{j|y}\sum_{k|z}[i\perp j][j\perp k][i\perp k] \\
&=\sum_{x=1}^a\sum_{y=1}^b\sum_{z=1}^c \lfloor\frac ax\rfloor \lfloor\frac by\rfloor \lfloor\frac cz\rfloor[x\perp y][y\perp z][x\perp z] \\
&=\sum_{x=1}^a\lfloor\frac ax\rfloor \sum_{y=1}^b\sum_{z=1}^c \lfloor\frac by\rfloor \lfloor\frac cz\rfloor[x\perp y]\sum_{d|y,d|z}\mu(d)[x\perp z] \\
&=\sum_{x=1}^a\lfloor\frac ax\rfloor\sum_{d=1}^{\min(b,c)}\mu(d)\sum_{y=1}^{b/d}\sum_{z=1}^{c/d}\lfloor\frac b{dy}\rfloor\lfloor\frac c{dz}\rfloor[x\perp dy][x\perp dz] \\
&=\sum_{x=1}^a\lfloor\frac ax\rfloor\sum_{d=1}^{\min(b,c)}\mu(d)[x\perp d]\sum_{y=1}^{b/d}\sum_{z=1}^{c/d}\lfloor\frac b{dy}\rfloor\lfloor\frac c{dz}\rfloor[x\perp y][x\perp z] \\
\end{aligned}
注意到我们在三个互质中只反演了一个。从这里开始就不同了。
后面式子的形式引出了两个函数的定义：
f(n,k)=\sum_{i=1}^n \lfloor\frac ni\rfloor[i\perp k]
g(n,k)=\sum_{i=1}^n\mu(i)[i\perp k]
则原式可以直接代入 f：
\sum_{x=1}^a \lfloor\frac ax\rfloor\sum_{d=1}^{\min(b,c)}\mu(d)[d\perp x]f(\lfloor\frac b{d} \rfloor,x)f(\lfloor\frac c{d} \rfloor,x)
则如果找出 f(n,k) 在整除点的值和 g(n,k) 在分块右端点前缀和的计算方法，式子后一部分可以整除分块计算。
忽略 $k$ 中平方因子，而 $p$ 是 $k$ 中一个质因子（同样参见链接），我们得到结论：
$$g(n,k)=g(\lfloor\frac np\rfloor,k)+g(n,\frac kp)$$
那我们尝试同样分析 $f(n,k)$。
$$\begin{aligned}
&\ \ \ \ \ f(n,k)=\sum_{i=1}^n\lfloor\frac ni\rfloor[i\perp k] \\
&=\sum_{i=1}^n \lfloor\frac ni\rfloor[i\perp\frac kp]-\sum_{i=1} ^n\lfloor\frac ni\rfloor[i\perp\frac kp][p|i] \\
&=f(n,\frac kp)-\sum_{i=1}^{n/p}\lfloor\frac n{ip}\rfloor[i\perp\frac kp] \\
&=f(n,\frac kp)-f(\lfloor\frac np \rfloor,\frac kp)
\end{aligned}$$
边界也显然。
$$f(0,k)=0,f(n,1)=\sum_{i=1}^n\lfloor\frac ni \rfloor=\sum_{i=1}^nd(i)$$
$$g(0,k)=0,g(n,1)=\sum_{i=1}^n\mu(d)$$
两者可以简单地线性筛出来。
---
$f(n,k),g(n,k)$ 的第一维对于一个 $n$，取遍 $b,c$ 的整除点有 $\sqrt n$ 个；
第二维对于一个 $k$，有取遍其素因子个，这个数低于 $\log k$。那 $[1,n]$ 中的无平方因子数有多少呢？最粗略的估计下，它有 $O(n)$ 个，因为即使有 $\log_2n$ 个 $2$ 作为素数进行乘积，也仅有 $2^{\log_2n}=n$ 个。
显然实际远低于 $O(n)$，但是下文仍写为 $O(n)$。
得到状态是 $O(n\sqrt n)$ 的，$O(1)$ 转移，加上分块，总复杂度是 $O(n\sqrt n)$ 的。
于是用 `unordered_map` 维护 $f,g$，枚举无平方因子的 $x$，分块时递归计算得到答案。
我们发现，根据 $k$ 的本质是一个素因子集合，对于该集合相同的 $x$，没有必要全部枚举。
记 $sv(x)$ 为 $x$ 因子去重得到的数，它可以线性筛得出。
我们统计 $sv(x)$ 相同的 $x$ 的 $\lfloor\dfrac ax\rfloor $ 之和，枚举存在的 $sv(x)$ 即可。
```cpp
#include <cstdio>
#include <iostream>
#include <cstring>
#include <cmath>
#include <algorithm>
#include <unordered_map>
#define umap std::unordered_map
typedef long long ll;
const int V = 100001;
const int p = 1000000007;
inline int read() {
int x = 0, f = 1; char ch = getchar();
while(ch > '9' || ch < '0') { if(ch == '-') f = -1; ch = getchar(); }
do x = x * 10 + ch - 48, ch = getchar(); while(ch >= '0' && ch <= '9');
return x * f;
}
int df[V],P[V],tt; ll mu[V],d[V],sv[V];
void Seive() {
mu[1] = d[1] = sv[1] = 1;
for(int i = 2;i < V;i++) {
if(!df[i]) P[++tt] = i, mu[i] = -1, d[i] = 2, sv[i] = i;
for(int j = 1;j <= tt && i * P[j] < V;j++) {
int t = i * P[j]; df[t] = 1;
if(i % P[j]) mu[t] = p - mu[i], d[t] = 2 * d[i] % p, sv[t] = sv[i] * P[j];
else mu[t] = 0, d[t] = (2 * d[i] % p - d[i / P[j]] + p) % p, sv[t] = sv[i];
}
}
for(int i = 1;i < V;i++) mu[i] = (mu[i] + mu[i - 1]) % p;
for(int i = 1;i < V;i++) d[i] = (d[i] + d[i - 1]) % p;
return;
}
ll pd[V],pdt,cl[V];
umap <int,umap <int,ll> > Fmap,Gmap;
ll F(int n,int k,int s) {
if(!n) return 0;
else if(k == 1) return d[n];
else if(Fmap[n].count(k)) return Fmap[n][k];
return Fmap[n][k] = (F(n,k / pd[s],s + 1) - F(n / pd[s],k / pd[s],s + 1) + p) % p;
}
ll G(int n,int k,int s) {
if(!n) return 0;
else if(k == 1) return mu[n];
else if(Gmap[n].count(k)) return Gmap[n][k];
return Gmap[n][k] = (G(n / pd[s],k,s) + G(n,k / pd[s],s + 1)) % p;
}
int main() {
Seive();
int tc = read();
while(tc--) {
int a = read(), b = read(), c = read(); ll ans = 0;
if(b > c) std::swap(b,c); if(a > b) std::swap(a,b);
for(int i = 1;i <= a;i++) cl[i] = 0;
for(int i = 1;i <= a;i++) cl[sv[i]] += a / i;
for(int x = 1;x <= a;x++) if(cl[x]) {
ll res = 0; pdt = 0;
for(int i = 1;i <= tt && P[i] <= x;i++)
if(!(x % P[i])) pd[++pdt] = P[i];
for(int i = 1;i <= b && i <= c;i++) {
int j = std::min(b / (b / i),c / (c / i));
res = (res + F(b / i,x,1) * F(c / i,x,1) % p * (G(j,x,1) - G(i - 1,x,1) + p) % p) % p;
i = j;
}
ans = (ans + cl[x] % p * res % p) % p;
}
std::printf("%lld\n",ans);
}
return 0;
}
```
---
结果 $a=b=c=10000$ 就跑了 $10s$，$a=b=c=100000$ 的根本出不来。
瓶颈在对 $f(n,k),g(n,k)$ 的维护上。考虑去掉二维的 `unordered_map`。
然后似乎无法可想？$f,g$ 的第一维是取满的 $O(\sqrt n)$，仍然考虑第二维。
如果对质因子集合的访问是随机的，就像我们初步得到的从小到大枚举的方法，无法解决问题。
那如果是顺序插入删除呢？
显然我们需要的是一个按 Dfs 序遍历质因子的顺序维护的 $f,g$ 的值。
那么可以把状态简单地顺序分层了。具体地，对于质数序列 $p_1,p_2,\cdots p_s$ 和 $k=\prod p_i$，将 $f(n,k)$ 和 $g(n,k)$ 记录在第 $s$ 层，转移时根据 $s-1$ 层的状态。
对于每一个 $k$ 统计答案。
$s$ 根据上面的分析，是 $O(\log_2n)$ 的，空间上完全可以承受。
注意转移中 $f$ 和 $g$ 用到的第一维永远只有 $O(\sqrt n)$ 个点，如上所述，$f$ 用到了整除点，$g$ 用到了整除点 $-1$ 的点。
实际复杂度比 $O(n\sqrt n)$ 要小不少，$100000$ 的样例大约 $4s$。
```cpp
#include <cstdio>
#include <iostream>
#include <cstring>
#include <cmath>
#include <algorithm>
typedef long long ll;
const int V = 100001;
const int p = 1000000007;
inline int read() {
int x = 0, f = 1; char ch = getchar();
while(ch > '9' || ch < '0') { if(ch == '-') f = -1; ch = getchar(); }
do x = x * 10 + ch - 48, ch = getchar(); while(ch >= '0' && ch <= '9');
return x * f;
}
int df[V],P[V],t;
ll mu[V],sv[V],d[V],cl[V];
void Seive() {
mu[1] = sv[1] = d[1] = 1;
for(int i = 2;i < V;i++) {
if(!df[i]) P[++t] = i, mu[i] = -1, sv[i] = i, d[i] = 2;
for(int j = 1;j <= t && i * P[j] < V;j++) {
int t = i * P[j]; df[t] = 1;
if(i % P[j]) mu[t] = -mu[i], sv[t] = sv[i] * P[j], d[t] = 2 * d[i] % p;
else mu[t] = 0, sv[t] = sv[i], d[t] = (2 * d[i] % p - d[i / P[j]] + p) % p;
}
}
for(int i = 1;i < V;i++) mu[i] = (mu[i - 1] + mu[i]) % p;
for(int i = 1;i < V;i++) d[i] = (d[i - 1] + d[i]) % p;
}
int a,b,c; ll ans;
ll f[V][20], g[V][20];
void Dfs(ll x,int s,int k) {
ll res = 0;
for(int i = 1;i <= b && i <= c;i++) {
int j = std::min(b / (b / i),c / (c / i));
res = (res + (g[j][k] - g[i - 1][k] + p) % p * f[b / i][k] % p * f[c / i][k] % p) % p;
i = j;
}
ans = (ans + cl[x] * res % p) % p;
for(int r = s;r <= t && x * P[r] <= a;r++) {
for(int i = 1;i <= b && i <= c;i++) {
int j = std::min(b / (b / i),c / (c / i));
f[b / i][k + 1] = (f[b / i][k] - f[b / i / P[r]][k] + p) % p;
f[c / i][k + 1] = (f[c / i][k] - f[c / i / P[r]][k] + p) % p;
g[j][k + 1] = (g[j / P[r]][k + 1] + g[j][k]) % p;
i = j;
}
Dfs(x * P[r],r + 1,k + 1);
}
return;
}
int main() {
Seive();
int t = read();
while(t--) {
a = read(), b = read(), c = read();
for(int i = 1;i <= a;i++) cl[i] = 0;
for(int i = 1;i <= a;i++) cl[sv[i]] += a / i;
for(int i = 1;i <= b && i <= c;i++) {
int j = std::min(b / (b / i),c / (c / i));
f[b / i][1] = d[b / i], f[c / i][1] = d[c / i];
g[j][1] = mu[j];
i = j;
}
ans = 0, Dfs(1,1,1);
std::printf("%lld\n",ans);
}
return 0;
}
```