【莫比乌斯反演】[BZOJ3994]约数个数和

设d(x)为x的约数个数，给定N、M，求

∑i=1n∑j=1md(i×j)

首先答案肯定是

Ans=∑i=1n∑j=1md(i×j)

发现

d(i)=∑i=1n⌊n/i⌋=f(i)

那么

Ans=∑i=1n∑j=1m[(i,j)==1]⌊n/i⌋⌊m/j⌋

为什么i 和 j要互质呢？
若i, j不互质 那么令 (i, j) = P
那么可以表示

i=a×P

j=b×P

那么该这次统计的个数表示的是含有因数

a×b×P2

可以发现这种次数在

a×b

的时候已经被统计过了， 那么就发生了重复计数的情况，所以必须要互质不然就会有重复
那么根据莫比乌斯函数的性质

∑d|nμ(d)=0(n>1)

所以就可以

Ans=∑i=1n∑j=1m∑d|i,d|jμ(d)∗⌊n/i⌋⌊m/j⌋

这个时候把

∑d|i,d|j

弄出去
就变成了

∑d=1min(n,m)μ(d)∑i=1n/d∑j=1m/d⌊nid⌋⌊mjd⌋

然后发现可以把后面的两个Σ分开，就变成了

∑d=1min(n,m)μ(d)∑i=1n/d⌊ndi⌋∑j=1m/d⌊mdj⌋

那么发现

∑i=1n/d⌊ndi⌋=f(⌊nd⌋)

那么就变成了

∑d=1min(n,m)μ(d)f(⌊nd⌋)f(⌊md⌋)

这个时候用筛法把f和d预处理出来就行了

下面看看筛法怎么做

void Init(int up){
    int tmp;
    d[1] = 1, mu[1] = 1;
    for(int i=2;i<=up;i++){
        if(!nprime[i]){f[i] = 1; mu[i] = -1; prime[++pcnt] = i; d[i] = 2;}
        for(int j=1;j<=pcnt&&(tmp = prime[j]*i)<=up;j++){
            nprime[tmp] = true;
            if(i % prime[j] == 0){
                mu[tmp] = 0;
                f[tmp] = f[i]+1;
                d[tmp] = d[i]/(f[i]+1)*(f[i]+2);
                break;
            }
            f[tmp] = 1;
            d[tmp] = d[i] * 2;
            mu[tmp] = -mu[i];
        }
    }
    for(int i=2;i<=up;i++)
        d[i] += d[i-1], mu[i] += mu[i-1];
}

其中mu数组就是莫比乌斯函数，然后d数组就是f，f数组存的是当前这个元素的最小质因子的次数是多少，根据这个我们可以发现/(当前最小质因子个数+1)*(当前质因子个数+2) 就等于下一种的d了，因为若当前的最小质因子为P 那么 P^a P参与组合的种数 就是a+1, 那么a+了1实际上就是多了一种P参与的方案（可能说的不太清楚。。。）

#include <bits/stdc++.h>
using namespace std;
const int MAXN = 50000;
int d[MAXN+10], f[MAXN+10], mu[MAXN+10], prime[MAXN/3], pcnt, ff, vv;
bool nprime[MAXN+10];
char ch;
int get()
{
    ff=0,vv=0;
    while(!isdigit(ch=getchar()))if(ch=='-')break;
    if(ch=='-')ff=1;else vv=ch-48;
    while(isdigit(ch=getchar()))vv=vv*10+ch-48;
    if(ff==1)return -vv;else return vv;
}
void Init(int up){
    int tmp;
    d[1] = 1, mu[1] = 1;
    for(int i=2;i<=up;i++){
        if(!nprime[i]){f[i] = 1; mu[i] = -1; prime[++pcnt] = i; d[i] = 2;}
        for(int j=1;j<=pcnt&&(tmp = prime[j]*i)<=up;j++){
            nprime[tmp] = true;
            if(i % prime[j] == 0){
                mu[tmp] = 0;
                f[tmp] = f[i]+1;
                d[tmp] = d[i]/(f[i]+1)*(f[i]+2);
                break;
            }
            f[tmp] = 1;
            d[tmp] = d[i] * 2;
            mu[tmp] = -mu[i];
        }
    }
    for(int i=2;i<=up;i++)
        d[i] += d[i-1], mu[i] += mu[i-1];
}
int main(){
    long long Ans = 0;
    Init(MAXN);
    int n, m, T;
    T = get();
    while(T--){
        Ans = 0;
        n = get(); m = get();
        if(n > m) swap(n, m);
        for(int i=1,nex;i<=n;i=nex+1){
            nex = min(n/(n/i), m/(m/i));
            Ans += 1LL * (mu[nex] - mu[i-1]) * d[n/i] * d[m/i];
        }
        printf("%lld\n", Ans);
    }

    return 0;
}

其中为什么c++ nex = min(n/(n/i), m/(m/i));
因为在(n/(n/i))和(m/(m/i))的时候如果不是i|n或者 i|m后边的f(n/i)*f(m/i)是不会发生变化的然后分配率。