[SDOI2018]旧试题 题解

Description

求 \(\sum_{i=1}^A\sum_{j=1}^B\sum_{k=1}^Cd(ijk)\)。
\(A,B,C\le10^5\)

Sol

直接莫反，枚举 \(i,j,k\) 的所有约数来得到 d(ijk) 的值：

\[\sum_{i=1}^A\sum_{j=1}^B\sum_{k=1}^Cd(ijk) \]

\[=\sum_{i=1}^A\sum_{j=1}^B\sum_{k=1}^C\sum_{d_1|i}\sum_{d_2|j}\sum_{d_3|k}[\gcd(d_1,d_2)=1][\gcd(d_2,d_3)=1][\gcd(d_1,d_3)=1] \]

枚举 \(d_1,d_2,d_3\)

\[=\sum_{d_1=1}^A\sum_{d_2=1}^B\sum_{d_3=1}^C\lfloor\frac A{d_1}\rfloor\lfloor\frac B{d_2}\rfloor\lfloor\frac C{d_3}\rfloor[\gcd(d_1,d_2)=1][\gcd(d_2,d_3)=1][\gcd(d_1,d_3)=1] \]

莫反

\[=\sum_{d_1=1}^A\sum_{d_2=1}^B\sum_{d_3=1}^C\lfloor\frac A{d_1}\rfloor\lfloor\frac B{d_2}\rfloor\lfloor\frac C{d_3}\rfloor\sum_{x|\gcd(d_1,d_2)}\sum_{y|\gcd(d_2,d_3)}\sum_{z|\gcd(d_1,d_3)}\mu(x)\mu(y)\mu(z) \]

枚举 x,y,z

\[=\sum_{x=1}^{\min\{A,B\}}\sum_{y=1}^{\min\{B,C\}}\sum_{z=1}^{\min\{A,C\}}\mu(x)\mu(y)\mu(z)\sum_{lcm(x,y)|d_1}\sum_{lcm(y,z)|d_2}\sum_{lcm(x,z)|d_3}\lfloor\frac A{d_1}\rfloor\lfloor\frac B{d_2}\rfloor\lfloor\frac C{d_3}\rfloor \]

我们设\(f(n,x)=\sum_{i=1}^{\lfloor\frac nx\rfloor}\lfloor\frac n{ix}\rfloor=\sum_{i=1}^{\lfloor\frac nx \rfloor}\lfloor\frac{\lfloor\frac nx\rfloor}{i}\rfloor\)，\(g(x)=\sum_{i=1}^x\lfloor\frac xi\rfloor \)
那么

\[\sum_{lcm(x,y)|d_1}\lfloor\frac A{d_1}\rfloor=\sum_{\frac{d_1}{lca(x,y)}=1}^{\frac{A}{lca(x,y)}}\lfloor\frac{A}{d_1\times lcm(x,y)}\rfloor \]

代入回原式

\[=\sum_{x=1}^{\min\{A,B\}}\sum_{y=1}^{\min\{B,C\}}\sum_{z=1}^{\min\{A,C\}}\mu(x)\mu(y)\mu(z)g(\lfloor\frac{A}{lcm(x,y)}\rfloor)g(\lfloor\frac{B}{lcm(y,z)}\rfloor)g(\lfloor\frac{C}{lcm(x,z)}\rfloor) \]

莫反完毕后，我们尝试建立图论模型，如果我们将所有点连成一个完全图的话，那么我们的答案就是所有三元环的权值之和。

观察到很多三元环的贡献为 \(0\)，考虑三元环 \((x,y,z)\) 的贡献 \(\mu(x)\mu(y)\mu(z)g(\lfloor\frac{A}{lcm(x,y)}\rfloor)g(\lfloor\frac{B}{lcm(y,z)}\rfloor)g(\lfloor\frac{C}{lcm(x,z)}\rfloor)\)，当且仅当 \(\mu(x),\mu(y),\mu(z)\) 的其中一个为 \(0\) 或 \(lcm(x,y),lcm(y,z),lcm(x,z)\) 其中一个为 \(0\) 时才没有贡献。

考虑先删掉所有 \(\mu\) 值为 \(0\) 的点，然后枚举公因数 \(k\) ，连接所有 \((ik,jk)\) 使得满足 \(i \neq j,\mu(ik)\neq0,\mu(jk)\neq0\)，然后 \(O(m\sqrt m)\) 找三元环即可。

Code

#include<bits/stdc++.h>
#define int long long
#define Mod 1000000007
using namespace std;
int Read() {
	int x = 0, f = 1; char ch = getchar();
	while(!isdigit(ch)) {if(ch == '-')  f = -1; ch = getchar();}
	while(isdigit(ch)) {x = (x << 3) + (x << 1) + ch - '0'; ch = getchar();}
	return x * f;
}
struct line {int x, y, w;}l[1000005];
int first[2000005], nxt[2000005], to[2000005], w[2000005], tot = 0;
void Add(int x, int y, int z) {nxt[++tot] = first[x]; first[x] = tot; to[tot] = y; w[tot] = z;}
int n, A, B, C, prime[100005], isp[100005], g[100005], mn[100005], mu[100005], id[100005], val[100005], cnt, num, lnum;
void prework() {
    mu[1] = 1, isp[1] = 1, g[1] = 1;
    for(int i = 1; i <= 100000; i++) {
        if(!isp[i])  prime[++cnt] = i, mu[i] = -1, g[i] = 2, mn[i] = 1;
        if(mu[i])  val[++num] = i, id[i] = num;
        for(int j = 1; j <= cnt && i * prime[j] <= 100000; j++) {
            isp[i * prime[j]] = 1;
            if(i % prime[j]) {
                mu[i * prime[j]] = -mu[i];
                mn[i * prime[j]] = mn[i];
                g[i * prime[j]] = g[i] * 2;
            }
            else {
                mu[i * prime[j]] = 0;
                g[i * prime[j]] = 2 * g[i] - g[i / prime[j]];
                mn[i * prime[j]] = mn[i] + 1;
                break;
            }
        }
    }
    for(int i = 1; i <= 100000; i++)  g[i]= (g[i - 1] + g[i]) % Mod;
}
int gcd(int a, int b) {return (a % b == 0) ? b : gcd(b, a % b);}
int rd[100005], ans, vis[100005];
int Calc(int x, int y, int z, int u, int v, int w) {
    return mu[x] * mu[y] * mu[z] * g[A / u] * g[B / v] * g[C / w];
}
signed main() {
    prework();
    int T = Read();
    while(T--) {
        ans = 0;
        memset(first, 0, sizeof(first)); tot = 0; lnum = 0;
        memset(rd, 0, sizeof(rd));
        A = Read(), B = Read(), C = Read();
        n = max(A, max(B, C));
        for(int i = 1; i <= num && val[i] <= n; i++) {
            for(int j = 1; j <= num && val[i] * val[j] <= n; j++) {
                if(mu[val[i] * val[j]] == 0)  continue;
                for(int k = j + 1; k <= num && val[i] * val[j] * val[k] <= n; k++) {
                    if(mu[val[i] * val[k]] == 0)  continue;
                    if(gcd(val[k], val[j]) != 1)  continue;
                    int u = id[val[i] * val[j]], v = id[val[i] * val[k]], w = val[i] * val[j] * val[k];
                    l[++lnum] = (line){u, v, w};
                    ++rd[u]; ++rd[v];
                    Add(u, v, w); Add(v, u, w);
                }
            }
        }
        n = min(A, B);
        for(int i = 1; i <= num && val[i] <= n; i++) {
            for(int e = first[i]; e; e = nxt[e]) {
                int y = to[e], z = w[e];
                ans += Calc(val[i], val[i], val[y], val[i], z, z);
                ans += Calc(val[i], val[y], val[i], z, val[i], z);
                ans += Calc(val[i], val[y], val[y], z, z, val[y]);
            }
        }
        n = min(min(A, B), C);
        for(int i = 1; i <= num && val[i] <= n; i++)
            ans += Calc(val[i], val[i], val[i], val[i], val[i], val[i]);
        memset(first, 0, sizeof(first)); tot = 0;
        for(int i = 1; i <= lnum; i++) {
            int x = l[i].x, y = l[i].y, z = l[i].w;
            if(rd[x] > rd[y])  Add(x, y, z);
            else if(rd[y] > rd[x])  Add(y, x, z);
            else if(x < y)  Add(x, y, z);
            else  Add(y, x, z);
        }
        n = max(max(A, B), C);
        for(int x = 1; x <= num && val[x] <= n; x++) {
            for(int y = first[x]; y; y = nxt[y]) {
                int z = to[y], v = w[y];
                vis[z] = v;
            }
            for(int y = first[x]; y; y = nxt[y]) {
                int u = to[y];
                for(int z = first[u]; z; z = nxt[z]) {
                    int v = to[z];
                    if(vis[v]) {
                        int a = val[x], b = val[u], c = val[v], d = w[y], e = w[z], f = vis[v];
                        ans += Calc(a, b, c, f, d, e); ans += Calc(a, c, b, d, f, e);
                        ans += Calc(b, a, c, e, d, f); ans += Calc(b, c, a, d, e, f);
                        ans += Calc(c, a, b, e, f, d); ans += Calc(c, b, a, f, e, d);
                    }
                }
            }
            for(int y = first[x]; y; y = nxt[y]) {
                vis[to[y]] = 0;
            }
        }
        cout << ans % Mod << endl;
    }
	return 0;
}