hdu4746-莫比乌斯反演+交换求和顺序+预处理

https://vjudge.net/contest/386727#problem/D

如果你在csdn上查本题的其他题解，你会很恼火，因为你不知道他们的F函数是怎么来的。我这篇题解将会把F函数的来历，这个唯一的关键点讲清楚！

不妨设n<=m。做这一假设是因为gcd值的范围为1~n。

考虑到只有gcd(i,j)==val这种形式的式子才能用莫比乌斯反演处理，我们不妨就加多一层求和符号，并增加限制。
$$ans=\sum _{1<=x<=n}[pf[x]<=k]\sum _{1<=i<=n}\sum _{1<=j<=m}[gcd(i,j)==x]$$
里面的两层求和符号是模板题，直接摆结论了。
$$ans=\sum _{1<=x<=n}[pf[x]<=k]\sum _{1<=d<=n/x}\mu (d)\frac{n}{x\ast d}\frac{m}{x\ast d}$$
预处理前缀和（注：k范围是1~18，因为全部取2就是上限。）+整除分块，直接求这个式子的复杂度应该不到O(n)（因为我不是很清楚），但也会TLE。

考虑交换求和顺序继续寻找预处理机会。设T=x*d。
$$ans=\sum _{1<=x<=n}[pf[x]<=k]\sum _{1<=T<=n}\mu (\frac{T}{x})\ast \frac{n}{T}\ast \frac{m}{T}\ast [x|T]$$
T摆脱了束缚，可以放最外层。
$$ans=\sum _{1<=T<=n}\frac{n}{T}\ast \frac{m}{T}\sum _{x|T}[pf[x]<=k]\ast \mu (\frac{T}{x})$$
记
$$F(k,T)=\sum _{x|T}[pf[x]<=k]\ast \mu (\frac{T}{x})$$

预处理所有的 F(k,T) 即可（代码里的数组sum）。为了方便，我们先固定pf[x]==k，求完以后最后遍历一遍数组求前缀和。我们枚举每个k（代码里的变量x），则k固定后，可以枚举倍数来求贡献。

re_(x,0,D){
    rep(i,1,lim){
        if(pf[i] != x) continue;
        for(int j = i;j <= lim;j += i)
            sum[x][j] += mu[j/i];
    }
}

完整代码

#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
#define rep(i,a,b) for(int i = (a);i <= (b);++i)
#define re_(i,a,b) for(int i = (a);i < (b);++i)
#define dwn(i,a,b) for(int i = (a);i >= (b);--i)

const int N = 5e5 + 5,D = 19;

int n,m,k,pf[N],sum[D][N];
bool vis[N];vector<int> primes;int mu[N];

template<typename Type>inline void read(Type &xx){
    Type f = 1;char ch;xx = 0;
    for(ch = getchar();ch < '0' || ch > '9';ch = getchar()) if(ch == '-') f = -1;
    for(;ch >= '0' && ch <= '9';ch = getchar()) xx = xx * 10 + ch - '0';
    xx *= f;
}

void init_primes(int lim){
    vis[1] = true;pf[1] = 0;mu[1] = 1;
    rep(i,2,lim){
        if(!vis[i]) primes.push_back(i),pf[i] = 1,mu[i] = -1;
        for(int &p: primes){
            if(i > lim / p) break;
            vis[i * p] = true;
            pf[i * p] = pf[i] + 1;
            if(i % p == 0) break;
            mu[i * p] = -mu[i];
        }
    }
    re_(x,0,D){
        rep(i,1,lim){
            if(pf[i] != x) continue;
            for(int j = i;j <= lim;j += i)
                sum[x][j] += mu[j/i];
        }
    }
    re_(x,0,D){
        rep(i,1,lim) sum[x][i] += sum[x][i-1];
    }
    rep(i,1,lim){
        re_(x,1,D) sum[x][i] += sum[x-1][i];
    }
}

int main(int argc, char** argv) {
    init_primes(N-5);
    int T;read(T);
    while(T--){
        read(n);read(m);read(k);
        k = min(k,D-1);
        int lim = min(n,m);
        LL ans = 0;
        for(int L = 1;L <= lim;){
            int R = min(n/(n/L),m/(m/L));
            ans += 1LL*(n/L)*(m/L)*(sum[k][R]-sum[k][L-1]);
            L = R+1;
        }
        printf("%lld\n",ans);
    }
    return 0;
}