BZOJ 4176: Lucas的数论(莫比乌斯反演+杜教筛)

Description

去年的Lucas非常喜欢数论题，但是一年以后的Lucas却不那么喜欢了。
在整理以前的试题时，发现了这样一道题目“求Sigma(f(i)),其中1<=i<=N”，其中 表示i的约数个数。他现在长大了，题目也变难了。
求如下表达式的值：

其中 表示ij的约数个数。
他发现答案有点大，只需要输出模1000000007的值。
Input

第一行一个整数n。
Output

一行一个整数ans，表示答案模1000000007的值。
Sample Input
2
Sample Output
8
HINT

对于100%的数据n <= 10^9。

解题思路

　　好神仙的一道题。首先有一个结论就是\(\sigma_0(ij)=\sum\limits_{a\mid i}\sum\limits_{b\mid j}[gcd(a,b)=1]\)。证明的话就是考虑每一个质数的贡献，发现左右两边相等。有了这个结论就可以推式子了。

\[ans=\sum\limits_{i=1}^n \sum\limits_{j=1}^n \sigma_0(ij) \]

\[ans=\sum\limits_{i=1}^n \sum\limits_{j=1}^n\sum\limits_{a\mid i}\sum\limits_{b\mid j}[gcd(a,b)=1] \]

　　莫比乌斯反演得：

\[ans=\sum\limits_{d=1}^n \mu(d)\sum\limits_{a=1}^{\frac{n}{d}}\sum\limits_{b=1}^{\frac{n}{d}}\sum\limits_{a\mid i}\sum\limits_{b\mid j}1 \]

　　把后面的式子合并一下：
　　$$ans=\sum\limits_{d=1}^{n\mu(d)(\sum\limits_{i=1}}{d}}\frac{n}{id})^2$$

　　发现这样就可以分块套分块去做了，\(\sum\mu\)用杜教筛即可。分块第一层可以看做枚举\(n/d\)，第二层就是里面的\(n/(id)\)。

代码

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<map>

using namespace std;
const int N=4000005;
const int MOD=1e9+7;
typedef long long LL;

int n,sum[N],prime[N],cnt,miu[N],ans;
bool vis[N];
map<int,int> mp;

inline int Add(int x,int y){
	x+=y; return (x>=MOD)?(x-MOD):x;
}
inline int Sub(int x,int y){
	x-=y; return x<0?x+MOD:x;
}

inline void prework(int n){
	miu[1]=1;
	for(int i=2;i<=n;i++){
		if(!vis[i]) prime[++cnt]=i,miu[i]=-1;
		for(int j=1;j<=cnt && 1ll*prime[j]*i<=n;j++){
			vis[i*prime[j]]=1;
			if(!(i%prime[j])) break;
			miu[i*prime[j]]=-miu[i];
		}
	}
	for(int i=1;i<=n;i++) sum[i]=sum[i-1]+miu[i];
}

int calc(int x){
	if(x<=N-5) return sum[x];
	if(mp.count(x)) return mp[x];
	int ret=1;
	for(int l=2,r;l<=x;l=r+1){
		r=x/(x/l);
		ret-=(r-l+1)*calc(x/l);
	} 
	return mp[x]=ret;
}

int calc2(int x){
	int ret=0;
	for(int l=1,r;l<=x;l=r+1){
		r=x/(x/l);
		ret=Add(ret,1ll*(r-l+1)*(x/l)%MOD);
	}
	return 1ll*ret*ret%MOD;
}

int main(){
	scanf("%d",&n); prework(4000000);
	for(int l=1,r;l<=n;l=r+1){
		r=n/(n/l);
		ans=Add(ans,1ll*Sub(calc(r),calc(l-1))*calc2(n/l)%MOD);
	}
	printf("%d\n",ans);
	return 0;
}