题解 P4619 【[SDOI2018]旧试题】

· · 题解

来自神仙 \text{1}\color{red}\text{1Dimensions} 的神仙做法......

在blog中查看？

分析

记x\perp y表示(x,y)=1。 结论：

d(ijk)=\sum\limits_{x|i}\sum\limits_{y|j}\sum\limits_{z|k}[x\perp y][x\perp z][y\perp z]

那么：

=& \sum\limits_{i=1}^A \sum\limits_{j=1}^B\sum\limits_{k=1}^C\sum\limits_{x|i}\sum\limits_{y|j}\sum\limits_{z|k}[x\perp y][x\perp z][y\perp z] \\ =& \sum\limits_{x=1}^A \sum\limits_{y=1}^B \sum\limits_{z=1}^C[x\perp y][x\perp z][y\perp z] \lfloor\frac{A}{x}\rfloor \lfloor\frac{B}{y}\rfloor \lfloor\frac{C}{z}\rfloor \end{aligned}

从三个单位函数里面任选一个反演，按老套路来：

=& \sum\limits_{x=1}^A \sum\limits_{y=1}^B \sum\limits_{z=1}^C (\sum\limits_{d|x,d|y}\mu(d)) [x\perp z][y\perp z] \lfloor\frac{A}{x}\rfloor \lfloor\frac{B}{y}\rfloor \lfloor\frac{C}{z}\rfloor \\ =& \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\perp z][k_2d\perp z] \lfloor\frac{A}{k_1d}\rfloor \lfloor\frac{B}{k_2d}\rfloor\end{aligned}

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

(\sum\limits_{x=1}^{\lfloor\frac{A}{d}\rfloor}[x\perp z]\lfloor\frac{A}{xd}\rfloor) (\sum\limits_{y=1}^{\lfloor\frac{B}{d}\rfloor} [y\perp z] \lfloor\frac{B}{yd}\rfloor)

记：

f(n,x) = \sum\limits_{i=1}^n\lfloor\frac{n}{i}\rfloor[i\perp x] \end{aligned}

答案变成：

=& \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)，将出现何种变化？应当有一部分多加了，要减去：

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

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

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

即：

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})的。如果常数写得好那么它可以跑的飞快（预处理整除分段的端点），但是本人代码常数没那么小，最慢一个点大约4s？ 比三元环是好得多了。

Orz \text{1}\color{red}\text{1Dimensions}

代码

#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);
    }
}