旧试题

Posted on 2020-08-23 Edited on 2026-03-02 In OI Word count in article: 1.5k Reading time ≈ 5 mins.

莫比乌斯反演好题

P4619 [SDOI2018]旧试题

问题描述

这里没有奇奇怪怪的描述

这道题说白了就是求： \[ \sum_{i=1}^A\sum_{j=1}^B\sum_{k=1}^C\text{d}(ijk)\\ 其中\text{d}(x)表示x的约数个数 \] 这道题和P3327 [SDOI2015]约数个数和比较类似，其实就是多了一个$j$

问题化简

我们可以知道： \[ \text{d}(ijk)=\sum_{x|i}\sum_{y|j}\sum_{z|k}[x\perp y][y\perp z][z\perp x] \]

即，原式可以化为： \[ \begin{align*} \sum_{i=1}^A\sum_{j=1}^B\sum_{k=1}^C\text{d}(ijk)&=\sum_{i=1}^A\sum_{j=1}^B\sum_{k=1}^C\sum_{x|i}\sum_{y|j}\sum_{z|k}[x\perp y][y\perp z][x\perp z]\\ &=\sum_{x=1}^A\sum_{y=1}^B\sum_{z=1}^C[x\perp y][y\perp z][x\perp z][\frac Ax][\frac By][\frac Cz]\\ &\text{选择x$\perp$y进行反演，即：}\\ &=\sum_{x=1}^A\sum_{y=1}^B\sum_{z=1}^C(\sum_{d|\gcd(x,y)}\mu(d))[y\perp z][x\perp z][\frac Ax][\frac By][\frac Cz]\\ &=\sum_{d=1}^{\min(A,B)}\mu(d)\sum_{x=1}^{\lfloor\frac Ad\rfloor}\sum_{y=1}^{\lfloor\frac Bd\rfloor}\sum_{z=1}^C[yd\perp z][xd\perp z][\frac A{xd}][\frac B{yd}][\frac Cz]\\ &\text{整理一下可以得到：}\\ &=\sum_{z=1}^C[\frac Cz]\sum_{d=1}^{\min(A,B)}\mu(d)[d\perp z](\sum_{x=1}^{\frac Ad}[x\perp z][\frac A{dx}])(\sum_{y=1}^{\frac Bd}[y\perp z][\frac B{dx}]) \end{align*} \] 这样一来，我们就成功了一大半了

现在，我们发现有两串看似比较相似的运算，我们可以用一个函数来代替它: \[ f(n,k)=\sum_{i=1}^n[i\perp k][\frac ni]\\ g(n,k)=\sum_{i=1}^n\mu(i)[i\perp k] \] 那么： \[ \begin{align*} Ans&=\sum_{z=1}^C[\frac Cz]\sum_{d=1}^{\min(A,B)}\mu(d)[d\perp z]f(\frac Ad,z)f(\frac Bd,z)\\ &=\sum_{z=1}^C[\frac Cz]\sum (g(r,z)-g(l-1,z))f(\frac Al,z)f(\frac Bl,z)\\ \end{align*} \] 这里，我们发现$\sum (g(r,z)-g(l-1,z))f(\frac Al,z)f(\frac Bl,z)$其实与$z$本身无关

简言之，我们只需要考虑$z$的无平方因子数

即： \[ lw(z)=\prod_{i=1} P_i^1 \]

但是如何快速求解这两个函数呢？

按照老套路，我们可以用$f(n,k/x)\;or\;g(n,k/x)$减去这其中不可行的得到

那么就有： \[ \begin{align*} f(n,k)&=\sum_{i=1}^n[i\perp k][\frac ni]\\ &=\sum_{i=1}^n[i\perp \frac kx][\frac ni]-\sum_{i=1}^n[i\perp \frac kx][\frac ni][x|i]\\ &=f(n,k/x)-\sum_{d=1}^{[\frac nx]}[dx\perp\frac kx][\frac n{dx}][x|kx]\\ &=f(n,k/x)-\sum_{d=1}^{[\frac nx]}[dx\perp\frac kx][\frac n{dx}]\\ &=f(n,k/x)-[x\perp \frac kx]\sum_{d=1}^{[\frac nx]}[d\perp\frac kx][\frac n{dx}]\\ &=f(n,k/x)-f([\frac nx],k/x) \end{align*} \] 那么还有： \[ \begin{align*} g(n,k)&=\sum_{i=1}^n\mu(i)[i\perp k] \\&=\sum_{i=1}^n\mu(i)[i\perp \frac kx]-\sum_{i=1}^n\mu(i)[i\perp\frac kx][x|i] \\&=g(n,k/x)-\sum_{d=1}^{[\frac nx]}\mu(dx)[dx\perp\frac kx][x|dx] \\&=g(n, k/x)-\sum_{d=1}^{[\frac nx]}\mu(d)\mu(x)[d\perp x][dx\perp \frac kx][x|dx] \\&=g(n,k/x)-\mu(x)[x\perp\frac kx]\sum_{d=1}^{[\frac nx]}\mu(d)[d\perp x][d\perp\frac kx] \\&=g(n,k/x)-\mu(x)[x\perp\frac kx]\sum_{d=1}^{[\frac nx]}\mu(d)[d\perp (x*\frac kx)] \\&=g(n,k/x)+\sum_{d=1}^{[\frac nx]}\mu(d)[d\perp k] \\&=g(n,k/x)+g([\frac nx],k) \end{align*} \] 根据上文，我们知道所有的$k$都一定是我平方因字数，所以这个结论很好得到

即： \[ \begin{align*} f(n,k)&=f(n,k/x)-f(\lfloor\frac nx\rfloor,k/x)\\ g(n,k)&=g(n,k/x)+g(\lfloor\frac nx\rfloor,k) \end{align*} \] 现在，最后的问题就是关于这两个函数的边界的问题了 \[ \begin{align*} g(n,1)&=\sum_{i=1}^n\mu(i)\qquad\text{这个可以直接前缀和预处理} \\f(n,1)&=\sum_{i=1}^n[\frac ni]\qquad\,\text{这个用数论分块吗？} \end{align*} \] 下面有一个非常巧妙的，线性求接的办法

inline void Prework(int N)
{
	D[1] = 1;
	for (int i = 2; i <= N; ++i)
	{
		if (!Ip[i]) P[++P[0]] = i, D[i] = 2;
		for (int j = 1; j <= P[0] && i * P[j] <= N; ++j)
		{
			Ip[i * P[j]] = 1;
			if (i % P[j]) D[i * P[j]] = D[i] << 1;
			else
			{
				D[i * P[j]] = (D[i] << 1) - D[i / P[j]];
				break;
			}
		}
	}
	for (int i = 1; i <= N; ++i) D[i] += D[i - 1];
}

关于这个的证明我先鸽了，不会

Code

#include <bits/stdc++.h>

using namespace std;

typedef long long ll;
const int Maxn = 2e5 + 10, Lim = 2e5, Mod = 1e9 + 7;

ll D[Maxn], F[Maxn][10], G[Maxn][10], S[Maxn], Ans;
int Prim[Maxn], Mu[Maxn], Lw[Maxn];
bool IP[Maxn];
int T, A, B, C;

inline int Read()
{
    int X(0);
    char O(getchar());
    while (O < '0' || O > '9') O = getchar ();
    for (; O >= '0' && O <= '9'; O = getchar ()) X = (X << 1) + (X << 3) + (O ^ 48);
    return X; 
}

inline void Prework()
{
    Mu[1] = D[1] = Lw[1] = 1;
    for (int i = 2; i <= Lim; ++i)
    {
        if (!IP[i]) Prim[++Prim[0]] = i, Mu[i] = -1, D[i] = 2, Lw[i] = i;
        for (int j = 1; j <= Prim[0] && i * Prim[j] <= Lim; ++j)
        {
            const int New = i * Prim[j];
            IP[New] = 1;
            if (i % Prim[j]) D[New] = D[i] << 1, Lw[New] = Lw[i] * Prim[j], Mu[New] = -Mu[i];
            else 
            {
                D[New] = (D[i] << 1) - D[i / Prim[j]];
                Lw[New] = Lw[i];
                break;
            }
        }
    }
    for (int i = 1; i <= Lim; ++i) D[i] += D[i - 1], Mu[i] += Mu[i - 1];
} 

inline void Init()
{
    for (int l = 1, r; l <= A; l = r + 1)
    {
        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];

    }
    memset (S, 0, sizeof S);
    for (int z = 1; z <= C; ++z) S[Lw[z]] += 1ll * C / z;
    Ans = 0;
}

inline void UpDate(int x, int k)
{
    for (int l = 1, r; l <= A; l = r + 1)
    {
        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]; 
    }
}//x就是上文中推公式用的x，k表示的是第k个素数，因为空间的问题，就写成这个鬼样子了

void DFS(int z, int u, int k)
{
    ll Tmp(0);
    for (int l = 1, r; l <= A; l = r + 1)
    {
        r = min(A / (A / l), B / (B / l));
        Tmp += (G[r][k] - G[l - 1][k]) * F[A / l][k] * F[B / l][k];
        Tmp = (Tmp % Mod + Mod) % Mod;
    }
    Ans = (Ans + S[z] * Tmp % Mod) % Mod;
    for (int v = u; 1ll * Prim[v] * z <= C; v++) UpDate(Prim[v], k + 1), DFS(z * Prim[v], v + 1, k + 1);
}

int main()
{
    T = Read();
    Prework();
    while (T--)
    {
        A = Read(), B = Read(), C = Read();
        if (A > B) swap (A, B);
        Init();
        DFS (1, 1, 1);
        printf ("%lld\n", Ans);
    }
    system ("pause");
}