「题解」旧试题

来自神仙 11Dimensions 的神仙做法……

题目

原题链接

求：
$$ \sum\limits_{i=1}^A \sum\limits_{j=1}^B \sum\limits_{k=1}^C d(ijk) $$
其中 $A,B,C ⩽ 200005$。

分析

记$x⊥y$表示$(x,y)=1$。
结论：

$$d(ijk)=\sum\limits_{x|i}\sum\limits_{y|i}\sum\limits_{z|i}[x⊥y][x⊥z][y⊥z]$$

其中[]表示艾弗森括号，当且仅当括号内命题为真时取值为1，否则为0。
在本篇题解中，这个括号的意义等价于单位函数，即$\epsilon((a,b))=[(a,b)=1]=[a⊥b]$。
那么：

$$\begin{aligned}& \sum\limits_{i=1}^A \sum\limits_{j=1}^B \sum\limits_{k=1}^C d(xyz) \newline
=& \sum\limits_{i=1}^A \sum\limits_{j=1}^B\sum\limits_{k=1}^C\sum\limits_{x|i}\sum\limits_{y|i}\sum\limits_{z|i}[x⊥y][x⊥z][y⊥z] \newline
=& \sum\limits_{x=1}^A \sum\limits_{y=1}^B \sum\limits_{z=1}^C[x⊥y][x⊥z][y⊥z] \lfloor\frac{A}{x}\rfloor \lfloor\frac{B}{y}\rfloor \lfloor\frac{C}{z}\rfloor \end{aligned}$$

从三个单位函数里面任选一个反演，利用$\epsilon=\mu*1$：

$$\begin{aligned}& \sum\limits_{x=1}^A \sum\limits_{y=1}^B \sum\limits_{z=1}^C[x⊥y][x⊥z][y⊥z] \lfloor\frac{A}{x}\rfloor \lfloor\frac{B}{y}\rfloor \lfloor\frac{C}{z}\rfloor \newline
=& \sum\limits_{x=1}^A \sum\limits_{y=1}^B \sum\limits_{z=1}^C (\sum\limits_{d|x,d|y}\mu(d)) [x⊥z][y⊥z] \lfloor\frac{A}{x}\rfloor \lfloor\frac{B}{y}\rfloor \lfloor\frac{C}{z}\rfloor \newline
=& \sum\limits_{z=1}^C\lfloor\frac{C}{z}\rfloor\sum\limits_{d=1}^{\min(A,B)}\mu(d)\sum\limits_{k_1=1}^{\lfloor\frac{A}{d}\rfloor}\sum\limits_{k_2=1}^{\lfloor\frac{B}{d}\rfloor}[k_1d⊥z][k_2d⊥z]
\lfloor\frac{A}{k_1d}\rfloor \lfloor\frac{B}{k_2d}\rfloor\end{aligned}$$

依据$[ab⊥c]\iff[a⊥c][b⊥c]$，整理得到：

$$\sum\limits_{z=1}^C\lfloor\frac{C}{z}\rfloor\sum\limits_{d=1}^{\min(A,B)}\mu(d)[d⊥z]
(\sum\limits_{x=1}^{\lfloor\frac{A}{d}\rfloor}[x⊥z]\lfloor\frac{A}{xd}\rfloor)
(\sum\limits_{y=1}^{\lfloor\frac{B}{d}\rfloor}
[y⊥z] \lfloor\frac{B}{yd}\rfloor)$$

记：

$$\begin{aligned} g(n,x) = \sum\limits_{i=1}^n\mu(i)[i⊥x] \newline
f(n,x) = \sum\limits_{i=1}^n\lfloor\frac{n}{i}\rfloor[i⊥x] \end{aligned}$$

答案变成：

$$\begin{aligned}& \sum\limits_{z=1}^C\lfloor\frac{C}{z}\rfloor\sum\limits_{d=1}^{\min(A,B)}\mu(d)[d⊥z]f(\lfloor\frac{A}{d}\rfloor, z)f(\lfloor\frac{B}{d}\rfloor, z) \newline
=& \sum\lfloor\frac{C}{z}\rfloor\sum(g(r,z) - g(l - 1, z))f(\lfloor\frac{A}{d}\rfloor, z)f(\lfloor\frac{B}{d}\rfloor, z)\end{aligned}$$

即对$\mu()$做前缀和然后对后面的$f()$分段，其中$[l,r]$表示整除分段的一段区间。

假设$O(n)$枚举$z$，那么求出后面的$\sum$的值是$O(\sqrt{n})$的。但是发现$z$是变化的，当$z$变化时暴力维护$f(),g()$是$O(n^2)$的，这成为了代码复杂度的瓶颈。如何解决？
能不能减少更新$f(),g()$的次数？不难发现$z$在函数中发挥作用的地方是判断一个数与其互质。考虑将$z$质因数分解，容易看出一个数与其互质仅和$z$的质因子的种类有关，而与质因子的幂次无关。
记：

$$lw(z) = \prod\limits_{i=1}^np_i, \text{where }z = \prod\limits_{i=1}^np_i^{\alpha_i}$$

那么我们只需要考虑$z\in [1,C]$的所有$lw(z)$值即可（即所有无平方因子的数），它们共用一套$f(),g()$的函数值。通过dfs暴力生成无平方因子数，我们可以把它们一并更新。

但是这远远不够，考虑通过递推来维护$f(),g()$。这个套路参考NOI2016 循环之美。
考虑在$z$中删去其一个质因子$x$，即$f(n,z)\to f(n,z/x)$，将出现何种变化？应当有一部分多加了，要减去：

$$\begin{aligned}f(n,z) &= \sum\limits_{i=1}^n\lfloor\frac{n}{i}\rfloor[i⊥z/x] -
\sum\limits_{i=1}^n\lfloor\frac{n}{i}\rfloor[i⊥z/x][x|i] \newline
&= f(n, z / x) - \sum\limits_{k=1}^{\lfloor\frac{n}{x}\rfloor}[kx⊥z/x]\lfloor\frac{n}{kx}\rfloor \newline
&= f(n, z / x) - [x⊥z/x]\sum\limits_{k=1}^{\lfloor\frac{n}{x}\rfloor}[k⊥z/x]\lfloor\frac{n}{kx}\rfloor\newline
&= f(n, z / x) - f(\lfloor\frac{n}{x}\rfloor, z/x)\end{aligned}$$

对$g()$的推导也同理，利用$\mu(ab)=\mu(a)\mu(b)[a⊥b]$，可以得出：

$$\begin{aligned}g(n,z) &= \sum\limits_{i=1}^n \mu(i) [i⊥z/x] -
\sum\limits_{i=1}^n \mu(i) [i⊥z/x][x|i] \newline
&= g(n, z / x) - \sum\limits_{k=1}^{\lfloor\frac{n}{x}\rfloor}[kx⊥z/x] \mu(kx) \newline
&= g(n, z / x) - \mu(x)\sum\limits_{k=1}^{\lfloor\frac{n}{x}\rfloor}\mu(k)[k⊥x][k⊥z/x] \newline
&= g(n, z / x) + \sum\limits_{k=1}^{\lfloor\frac{n}{x}\rfloor}\mu(k)[k⊥x(z/x)] \newline
&= g(n, z / x) + g(\lfloor\frac{n}{x}\rfloor, z)\end{aligned}$$

即：

$$\begin{aligned} g(n,z) = g(n, z / x) + g(\lfloor\frac{n}{x}\rfloor, z) \newline
f(n,z) = f(n, z / x) - f(\lfloor\frac{n}{x}\rfloor, z/x) \end{aligned}$$

那么就可以递推了，但是直接更新是$O(n)$的。容易发现整除分段并不会用到所用的$f(),g()$值，所以我们边做分段边更新就好了，这是$O(\sqrt{n})$的。
空间复杂度呢？f[2e5][2e5],g[2e5][2e5]看似存不下来，但是容易发现后面一维是不连续使用的，那么将其离散化，考虑在dfs构造无平方因子数的时候，下层状态的转移依赖于上层状态，而dfs树的深度是$O(\log_2n)$的，所以开f[2e5][10]即可。

总复杂度？粗略估计，爆搜的复杂度是$O(2^{\log_2n})$即$O(n)$的，后面的求和通过整除分段可以$O(\sqrt{n})$得出，那么总复杂度是$O(n\sqrt{n})$的。但实际上复杂度要小很多，也就是说如果常数写得好那么它可以跑的飞快（预处理整除分段的端点、cache-friendly之类的），但是本人代码常数没那么小，最慢一个点大约4s？
比三元环是好得多了。

$\text{1}{\color{Red}{1Dimensions}}$ Orz %%%%

代码

#include<iostream>
#include<cstdio>
#include<cstring>
#include<string>
using namespace std;

#define LL long long

const int CN = 2e5+5;
const LL P = 1e9+7;

LL read(){    
    LL s = 0,ne = 1; char c = getchar();
    for(;c < '0' || c > '9'; c = getchar()) if(c == '-') ne = -1;
    for(;c >= '0' && c <= '9';c = getchar()) s = (s << 1) + (s << 3) + c - '0';
    return s * ne;
}

int p[CN],mu[CN],lw[CN]; LL d[CN]; bool np[CN];
void sieve(int n){
    np[1] = true; mu[1] = 1, d[1] = 1, lw[1] = 1;
    for(int i = 2;i <= n; i++){
        if(!np[i]) p[ ++p[0] ] = i, d[i] = 2, lw[i] = i, mu[i] = -1;
        for(int j = 1;j <= p[0] && i * p[j] <= n; j++){
            int x = i * p[j]; np[x] = true;
            if(i % p[j]) d[x] = d[i] << 1, lw[x] = lw[i] * p[j], mu[x] = -mu[i];
            else {d[x] = (d[i] << 1) - d[i / p[j]], lw[x] = lw[i]; break;}
        }
    }
    for(int i = 1;i <= n;i++) d[i] += d[i - 1],mu[i] += mu[i - 1];
}

// 
int t,A,B,C;
LL f[CN][10],g[CN][10],s[CN],ans;
void init(){
    g[0][1] = f[0][1] = 0;
    for(int l = 1;l <= A; l++){
        int r = min(A / (A / l), B / (B / l));
        g[r][1] = mu[r]; f[A / l][1] = d[A / l]; f[B / l][1] = d[B / l];
        l = r;
    }

    memset(s, 0 , sizeof(s));
    for(int z = 1;z <= C; z++) s[ lw[z] ] += 1ll * C / z;

    ans = 0;
}
void upd(int x,int k){
    for(int l = 1;l <= A; l++){
        int r = min(A / (A / l), B / (B / l));
        g[r][k] = g[r][k - 1] + g[r / x][k];
        f[A / l][k] = f[A / l][k - 1] - f[(A / l) / x][k - 1];
        f[B / l][k] = f[B / l][k - 1] - f[(B / l) / x][k - 1];
        l = r;
    }
}
void dfs(int z0, int u, int k){
    LL cur = 0;
    for(int l = 1;l <= A; l++){
        int r = min(A / (A / l), B / (B / l));
        cur += (g[r][k] - g[l - 1][k]) * f[A / l][k] * f[B / l][k];
        cur = (cur % P + P) % P;
        l = r;
    }
    ans += (s[z0] * cur) % P; ans %= P;
    for(int v = u; 1ll * p[v] * z0 <= C; v++) upd(p[v], k + 1), dfs(z0 * p[v], v + 1, k + 1);
}

int main()
{
    freopen("_in.in", "r", stdin);
    
    t = read(); sieve(2e5); 
    while(t--){
        A = read(), B = read(), C = read();
        if(A > B) swap(A, B);

        init(); dfs(1, 1, 1);
        printf("%lld\n", ans);
    }
}