【莫比乌斯反演】[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)是不会发生变化的然后分配率。