bzoj3994/洛谷P3327 莫比乌斯反演

前言

话说这是让我来学莫比乌斯反演的入坑题呢，然而学了莫比乌斯反演还是不会做=_=，连题解都看不懂QAQ
注意事项：
1.数恐症患者慎入！
2.所有除号如果没有作特殊说明，都是向下取整。
3.写的又臭又长，每一步都有一堆啰嗦的解释，因为本蒟蒻每一步都理解了半天（看神犇题解，一个等号看半个小时系列）
4.本蒟蒻发现自己学了假的懵逼钨丝繁衍，所以我重构了一遍这个博客

题目分析

第一步：预备的结论

首先，我们要证明一条结论：

$$d(nm)=\sum_{i|n}\sum_{j|m}[gcd(i,j)==1]$$
怎么证明呢？对于 $nm$的每一个质因数 $p_i$，我们设 $n=p_i^{a_i}*n'$, $m=p_i^{b_i}*m'$,那么 $$d(nm)=\prod_i(a_i+b_i+1)$$
( $\prod$是累乘的意思）
这条结论意会一下就可以了吧？对了，我们记 $p_i$对 $d(nm)$的贡献为 $a_i+b_i+1$,“贡献”这个词第三步要用
然后，我们看着要证的那个结论的等号右边部分，设 $l$和$t$互质且都不含质因子 $p_i$，那么 $(l*p_i^{a_i},t),(l*p_i^{a_i-1},t)...(l,t)...(l,t*p_i^{b_i-1}),(l,t*p_i^{b_i})$都是互质的数对，而 $l$和$t$的取法正好有 $\prod_{k \neq i}(a_k+b_k+1)$种。
啊哈，得证了！

第二步：公式要变形

有了上面的结论，我们就可以玩公式变形了！

$$\sum_{i=1}^n\sum_{j=1}^m d(ij)=\sum_{i=1}^n\sum_{j=1}^m \sum_{k|i}\sum_{l|j}[gcd(k,l)==1]$$
现在我们枚举每一个数对 $(k,l)$,那么它们会被算多少遍呢？每一个 $k$的倍数都会算一次$k$，每一个 $l$的倍数又会算一次$l$,而用 $\frac{n}{k}$来表示 $n$以内$k$的倍数有多少个是显然的，所以：
$$\sum_{k=1}^n\sum_{l=1}^m\frac{n}{k}*\frac{m}{l}[gcd(k,l)==1]$$
然后大家应该知道莫比乌斯函数有一个神奇的性质（详见： 莫比乌斯反演的讲解）：
$$\sum_{d|n}\mu(d)=1(n=1)$$
$$\sum_{d|n}\mu(d)=0(n>1)$$
所以呢：
$$\sum_{k=1}^n\sum_{l=1}^m\frac{n}{k}*\frac{m}{l}\sum_{d|gcd(k,l)}\mu(d)$$
这样子如果 $gcd(i,j)=1$的话后面那团就是1，否则为0.
然后这一步变形非常显然对吧：
$$\sum_{k=1}^n\sum_{l=1}^m\frac{n}{k}*\frac{m}{l}\sum_{d|k且d|l}\mu(d)$$
现在我们先枚举 $\mu(d)$，而哪些数会和 $\mu(d)$一起乘呢？根据上一个式子，显然是 $d$的倍数啊，那么是几倍呢？枚举啊！
$$\sum_{d=1}^{min(n,m)}\mu(d)\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{m}{d}}\frac{n}{di}\frac{m}{dj}$$

第三步：转化成代码

已经变形的差不多了，再小小变一下得到最终公式：

$$\sum_{d=1}^{min(n,m)}\mu(d)(\sum_{k=1}^{\frac{n}{d}}\frac{n}{dk})(\sum_{l=1}^{\frac{m}{d}}\frac{m}{dl})$$
现在我们令 $g(x)=\sum_{i=1}^i\frac{x}{i}$，原式就是:
$$\sum_{d=1}^{min(n,m)}\mu(d)*g(\frac{n}{d})*g(\frac{m}{d})$$
那么我们只要预处理 $g(x)$就可以了！
现在观察一下 $g(x)$，就会发现其实式子可以这么理解： $x$以内有约数$i$的数的个数和，也就是可看作原题中 $d(x)$的前缀和
要求 $d(x)$的前缀和，可以求出每一个 $d(x)$,也就是用线性筛求啦，要记录一下 $c(x)$:将 $x$分解质因数后，最小质因数的指数。
那么参见第一步，$x$的最小质因数对 $d(x)$的贡献为 $c(x)+1$，如果还不懂怎么筛求看代码吧。

代码

经过刚才含辛茹苦地淬炼，蒟蒻终于写出一段精简的代码：

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
#define ll long long
int read(){
    int q=0;char ch=' ';
    while(ch<'0'||ch>'9')ch=getchar();
    while(ch>='0'&&ch<='9')q=q*10+ch-'0',ch=getchar();
    return q;
}
const int maxn=50001;
int tot,lim=50000,T,n,m;
int pri[30001],mu[maxn],c[maxn];
ll d[maxn];bool is[maxn];
void init(){
    int i,j,k;
    mu[1]=c[1]=d[1]=1;
    for(i=2;i<=lim;i++){
        if(!is[i]){
            pri[++tot]=i;mu[i]=-1;
            c[i]=1;d[i]=2;
        }
        for(j=1;j<=tot;j++){
            k=pri[j]*i;if(k>lim)break;
            is[k]=1;
            if(i%pri[j]){ d[k]=d[i]*d[pri[j]];c[k]=1;mu[k]=-mu[i];}//有新的最小质因子了
            else {d[k]=d[i]/(c[i]+1)*(c[i]+2);c[k]=c[i]+1;break;}
            //d:删去原来的贡献，添加新的贡献
        }
    }
    for(i=2;i<=lim;i++)d[i]=d[i]+d[i-1],mu[i]=mu[i]+mu[i-1];//前缀和
}
int main()
{
    int i,j;
    init();T=read();
    while(T--){
        n=read();m=read();
        if(n>m)swap(n,m);ll ans=0;
        for(i=1;i<=n;i=j+1){
            j=min(n/(n/i),m/(m/i));//不懂这句请看我的莫比乌斯反演讲解
            ans+=(mu[j]-mu[i-1])*d[n/i]*d[m/i];
        }
        printf("%lld\n",ans);
    }
    return 0;
}