bzoj3739: DZY loves math VIII 莫比乌斯反演

bzoj3739: DZY loves math VIII

Description

在XYZ的dzy loves math6问世后，dzy一直觉得这道题答案太大，一点都不优美，于是他随手在外面套上一个μ。同时，他又觉得输入两个数实在太麻烦，于是题目变成了，你能解决这个问题吗？

Input

第一行一个整数T表示询问组数，接下来T行每行一个整数n。

Output

对于每一个询问输出一行表示答案

Sample Input

1
2

Sample Output

0

HINT

T<=10^3 n<=10^7

分析

首先要注意到当且仅当$gcd(i,j)==1$才有贡献。
于是变成
$ans=\sum_i^n\sum_j^i \mu(ij) [gcd(i,j)==1]$
$=\sum_i^n\sum_j^i \mu(i)\mu(j) \sum_{k|i,k|j}\mu(k)$
注意到这里不好直接枚举$k$，因为如果枚举$k$
$=\sum_k^n\mu(k)\sum_i^{\frac{n}{k}}\sum_j^i \mu(i*k)\mu(j*k)$
$=\sum_k^n\mu(k)(\sum_i^{\frac{n}{k}}\mu(i*k) \sum_j^{\frac{n}{k}} \mu(j*k) + \sum_x^{\frac{n}{k} }\mu^2(x*k))/2$
那个$k$很不好处理，虽然可以得到一个$O(Tnlogn)$的算法(暴力处理$\sum\mu(i*k)$,$\sum\mu^2(i*k)$)，但是复杂度不可以接受。
这个时候我们注意到$T$比较大，所以需要设计 一个类似递推的方法。
故先枚举$i$，再枚举$k$
$ans=\sum_i^n\mu(i)\sum_{k|i} \mu(k) \sum_j^{\frac{i}{k}}\mu(j*k)$
这样的好处是，后面的$\sum_{k|i} \mu(k) \sum_j^{\frac{i}{k}}\mu(j*k)$只和$i$有关，每次只要在之前答案的基础上计算$\sum_{k|i}\mu(k) \sum_j^{\frac{i}{k}}\mu(j*k)$
还有一个优秀的性质是，当且仅当$i$时square-free(无平方因子)，才会对答案有贡献，否则有$ans[n]=ans[n-1](\mu(n)==0)$
而$10^7<2*3*5*7*11*13*15*17*19*23$ 这意味着每个数最多只有9个不同的素因子，所以枚举$i$的复杂度变成了$O(2^9)$
现在的问题变成了如何快速求出$\sum_{k|i} \mu(k) \sum_j^{\frac{i}{k}}\mu(j*k)$
枚举$k$的复杂度是$O(\sqrt i)$的。
剩下处理$f(i,k)=\sum_j^{\frac{i}{k}}\mu(j*k)(k|i)$
注意到$(k|i)$，故对于每个$k$只需要考虑$f(i+k,k)$与$f(i,k)$的关系。
显然有$f(i+k,k)=\sum_j^{\frac{i}{k}+1}\mu(j*k)=f(i,k)+\mu(i+k)$
于是我们仅仅记录$f(k)$的值，随$i$的变化边做边更新即可。
复杂度$O(2^9 \sqrt n + T)$，优秀的算法！

分析

#include<cstdio>
const int N = 1e7 + 10;
int f[N], u[N], pr[N], p[N], s[N], q[N], mn[N], tp;
void Pre(int N) {
    u[1] = 1;
    for(int i = 2;i <= N; ++i) {
        if(!mn[i]) pr[++tp] = i, mn[i] = i, u[i] = -1;
        for(int j = 1;j <= tp && i * pr[j] <= N; ++j) {
            mn[i * pr[j]] = pr[j];
            if(i % pr[j]) u[i * pr[j]] = -u[i];
            else break;
        }
    }
}
int Dfs(int k, int x, int v) {
    if(k > tp) return f[x] += v, u[x] * f[x];
    return Dfs(k + 1, x * p[k], v) + Dfs(k + 1, x, v);
}
int main() {
    int T; scanf("%d", &T); int n = 0;
    for(int i = 1;i <= T; ++i)
        scanf("%d", &q[i]), n < q[i] ? n = q[i] : 0;
    Pre(n);
    for(int i = 1;i <= n; ++i) {
        s[i] = s[i - 1];
        if(u[i]) {
            tp = 0; for(int x = i; x != 1; x /= mn[x]) p[++tp] = mn[x];
            s[i] += u[i] * Dfs(1, 1, u[i]);
        }
    }
    for(int i = 1;i <= T; ++i) printf("%d\n", s[q[i]]);
    return 0;
}