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

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

i=1nj=1md(i×j)

首先答案肯定是
Ans=i=1nj=1md(i×j)

发现
d(i)=i=1nn/i=f(i)

那么
Ans=i=1nj=1m[(i,j)==1]n/im/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=1nj=1md|i,d|jμ(d)n/im/j

这个时候把
d|i,d|j
弄出去
就变成了
d=1min(n,m)μ(d)i=1n/dj=1m/dnidmjd

然后发现可以把后面的两个Σ分开,就变成了
d=1min(n,m)μ(d)i=1n/dndij=1m/dmdj

那么发现
i=1n/dndi=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)是不会发生变化的然后分配率。

posted on 2015-05-27 13:16  JeremyGuo  阅读(112)  评论(0编辑  收藏  举报

导航