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)=sumlimits_{amid i}sumlimits_{bmid j}[gcd(a,b)=1])。证明的话就是考虑每一个质数的贡献,发现左右两边相等。有了这个结论就可以推式子了。

[ans=sumlimits_{i=1}^n sumlimits_{j=1}^n sigma_0(ij) ]

[ans=sumlimits_{i=1}^n sumlimits_{j=1}^nsumlimits_{amid i}sumlimits_{bmid j}[gcd(a,b)=1] ]

  莫比乌斯反演得:

[ans=sumlimits_{d=1}^n mu(d)sumlimits_{a=1}^{frac{n}{d}}sumlimits_{b=1}^{frac{n}{d}}sumlimits_{amid i}sumlimits_{bmid j}1 ]

  把后面的式子合并一下:
  $$ans=sumlimits_{d=1}nmu(d)(sumlimits_{i=1}{frac{n}{d}}frac{n}{id})^2$$

  发现这样就可以分块套分块去做了,(summu)用杜教筛即可。分块第一层可以看做枚举(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
",ans);
	return 0;
}
原文地址:https://www.cnblogs.com/sdfzsyq/p/10485671.html