BZOJ3994 莫比乌斯反演

SDOI2015题目描述

求∑ni=1∑mj=1d(ij)，其中d(x)表示x的约数个数和

首先我们证明：d(nm)=∑i|n∑j|m[gcd(i,j)==1]
易证：对于nm的每个质因数p，设n=n′∗pa，m=m′∗pb，那么p对d(nm)的贡献就是(a+b+1)。在等式右边，(i′∗pa,i′),(i′∗pa−1,j′),...,(i′,j′),...,(i′,j′∗pb)都是可行的数对，共a+b+1个，因此成立。
那么就可以推导：

∑i=1N∑j=1Md(ij)=∑i=1N∑j=1M∑k|i∑l|j[gcd(k.l==1)]=∑i=1N∑j=1M⌊Ni⌋⌊Mj⌋[gcd(i,j)==1]=∑i=1N∑j=1M⌊Ni⌋⌊Mj⌋∑d|gcd(i,j)μ(d)=∑d=1min(N,M)μ(d)∑i=1N∑j=1M⌊Ni⌋⌊Mj⌋=∑d=1min(N,M)μ(d)∑i=1Nd⌊Ni∗d⌋∑j=1Md⌊Md∗j⌋

有一个性质：

∑d|nμ(d)={1     n=10      n>1

令f(n)=∑ni=1⌊ni⌋，则瓦们要求的就是：

∑i=1min(N,M)μ(i)∗f(⌊Ni⌋)∗f(⌊Mi⌋)

Code如下

#include<algorithm>#include<cstring>#include<cstdio>using namespace std;typedef long long LL;const int N=50001;int tot,mu[N],is[N],p[N];LL f[N],Ans[1001][1001];template <class T> void read(T &x) {    x=0;int f=1;char ch=getchar();    for(;ch<'0'||ch>'9';) {if(ch=='-') f=-1;ch=getchar();}    for(;ch>='0'&&ch<='9';) x=x*10+ch-'0',ch=getchar();x*=f;}template <class T> void write(T x) {    int num=0,ch[20];    if(x<0) x=-x,putchar('-');    do ch[++num]=x%10+'0',x/=10; while(x);    for(;num;) putchar(ch[num--]);putchar('\n');}void prime(int n) {    int i,j;    mu[1]=1;    for(i=2;i<=n;++i) {        if(!is[i]) p[++tot]=i,mu[i]=-1;        for(j=1;j<=tot&&i*p[j]<=n;++j) {            is[i*p[j]]=1;            if(i%p[j]==0) {mu[i*p[j]]=0;break;}            mu[i*p[j]]=-mu[i];        }        mu[i]+=mu[i-1];    }    return;}int main() {    int n,m,i,j,nex,test;LL ans;    for(prime(50000),i=1;i<=50000;++i)        for(nex=0,j=1;j<=i;j=nex+1)            nex=i/(i/j),f[i]+=(i/j)*(nex-j+1);    for(read(test);test;test--) {        read(n),read(m);        if(n<=1000&&m<=1000&&Ans[n][m]) {write(Ans[n][m]);continue;}        if(n>m) swap(n,m);        for(ans=nex=0,i=1;i<=n;i=nex+1) {            nex=min(n/(n/i),m/(m/i));            ans+=(mu[nex]-mu[i-1])*f[n/i]*f[m/i];        }        if(n<=1000&&m<=1000) Ans[n][m]=Ans[m][n]=ans;        write(ans);    }    return 0;}