一类积性函数的前缀和---刷题记录

题目来源于糖教主浅谈一类积性函数的前缀和...

51Nod 1244 莫比乌斯函数之和

考虑$mu(x)$的性质：$[n==1]=sum _{dmid n} mu(d)$

可以用上面哪个公式来推导：

$f(n)=sum _{i=1}^{n}$

$1=sum _{i=1}^{n} [i==1]$

$=sum _{i=1}^{n} sum _{dmid i} mu (d)$

$=sum _{frac{i}{d}=1}^{n} sum _{d=1}^{frac{n}{frac{i}{d}}} mu (d)$

$=sum _{i=1}^{n}sum _{d=1}^{frac{n}{i}} mu(d)$

$=sum _{i=1}^{n}f(frac{n}{i})$

$f(n)=1-sum _{i=2}^{n} f(frac{n}{i})$

然后，我们预处理出前$n^{frac{2}{3}}$个的$f(x)$，然后对于大于$n^{frac{2}{3}}$的数的答案，分块递归计算...

复杂度的证明请见糖教主的文章...

代码：

#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdio>
#include<map>
//by NeighThorn
using namespace std;

const int maxn=5000000+5;

int cnt,mu[maxn],pri[maxn],vis[maxn];
long long n,m,f[maxn];

map<long long,long long> mp;

inline void prework(void){
	mu[1]=1;
	for(int i=2;i<=5000000;i++){
		if(!vis[i])
			vis[i]=1,pri[++cnt]=i,mu[i]=-1;
		for(int j=1;j<=cnt&&1LL*i*pri[j]<=5000000;j++){
			vis[i*pri[j]]=1;
			if(i%pri[j]==0){
				mu[i*pri[j]]=0;
				break;
			}
			mu[i*pri[j]]=-mu[i];
		}
	}
	for(int i=1;i<=5000000;i++) f[i]=f[i-1]+mu[i];
}

inline long long calc(long long x){
	if(x<=5000000) return f[x];
	if(mp.find(x)!=mp.end()) return mp[x];
	long long ans=1;
	for(long long i=2,r;i<=x;i=r+1){
		r=x/(x/i);
		ans-=calc(x/i)*(r-i+1);
	}
	return mp[x]=ans;
}

signed main(void){
	prework();scanf("%lld%lld",&n,&m);
	printf("%lld
",calc(m)-calc(n-1));
	return 0;
}

　　

51Nod 1239 欧拉函数之和 

和上面的题目差不多...

这次利用的是$phi(x)$的这个性质：$sum _{dmid n} phi(d)=n$

$phi(n)=n-sum _{dmid n  d<n}phi(d)$

$f(n)=sum _{i=1}^{n} (i-sum _{dmid i  d<i} phi(d))$

$=frac{n(n+1)}{2}-sum _{i=2}^{n} sum _{dmid i  d<i} phi(d)$

$=frac{n(n+1)}{2}-sum _{i=2}^{n}sum _{d=1}^{frac{n}{i}} phi(d)$

$=frac{n(n+1)}{2}-sum _{i=2}^{n} f(frac{n}{i})$

代码：

#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdio>
#include<map>
//by NeighThorn
using namespace std;

const int maxn=5000000+5,mod=1e9+7;

int cnt,f[maxn],pri[maxn],phi[maxn],vis[maxn];
long long n;

map<long long,int> mp;

inline int mul(long long x,long long y){
	int res=0;x%=mod;
	while(y){
		if(y&1) res=(res+x)%mod;
		x=(x+x)%mod,y>>=1;
	}
	return res;
}

inline int power(int x,int y){
	int res=1;
	while(y){
		if(y&1) res=1LL*res*x%mod;
		x=1LL*x*x%mod,y>>=1;
	}
	return res;
}

inline void prework(void){
	phi[1]=1;
	for(int i=2;i<=5000000;i++){
		if(!vis[i])
			pri[++cnt]=i,vis[i]=1,phi[i]=i-1;
		for(int j=1;j<=cnt&&1LL*i*pri[j]<=5000000;j++){
			vis[i*pri[j]]=1;
			if(i%pri[j]==0){
				phi[i*pri[j]]=phi[i]*pri[j];
				break;
			}
			phi[i*pri[j]]=phi[i]*(pri[j]-1);
		}
	}
	for(int i=1;i<=5000000;i++) f[i]=(f[i-1]+phi[i])%mod; 
}

inline int calc(long long n){
	if(n<=5000000) return f[n];
	if(mp.find(n)!=mp.end()) return mp[n];
	int ans=mul(mul(n,n+1),power(2,mod-2));
	for(long long i=2,r;i<=n;i=r+1){
		r=n/(n/i);
		ans=(ans-1LL*calc(n/i)*((r-i+1)%mod)%mod+mod)%mod;
	}
	return mp[n]=ans;
}

signed main(void){
	prework();scanf("%lld",&n);
	printf("%lld
",calc(n));
	return 0;
}

　　

51Nod 1190 最小公倍数之和 V2

 推了好久的式子发现本来还挺对的，后来越来越麻烦...一定是化简方向错了...所以看了题解...发现果然如此...

懒得写一遍式子了...一起%题解吧...

需要注意的是，对于枚举判断一个数的质因子的时候，记得加上一个小的剪枝，否则就是把9000+的素数全部枚举一边，绝逼GG啊...

if(1LL*pri[i]*pri[i]>tmp) break;

 代码：

#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdio>
//by NeighThorn
using namespace std;

const int maxn=100000+5,mod=1e9+7;

int n,m,ans,cas,cnt,inv,tot,tot_pri,pri[maxn],exp[maxn],vis[maxn],fact[maxn],fac_pri[maxn];

inline int power(int x,int y){
	int res=1;
	while(y){
		if(y&1)
			res=1LL*res*x%mod;
		x=1LL*x*x%mod,y>>=1;
	}
	return res;
}

inline void prework(void){
	for(int i=2;i<=100000;i++){
		if(!vis[i])
			vis[i]=1,pri[++cnt]=i;
		for(int j=1;j<=cnt&&1LL*i*pri[j]<=100000;j++){
			vis[i*pri[j]]=1;
			if(i%pri[j]==0) break;
		}
	}
}

inline void dfs(int id,int fac){
	if(id>tot_pri){
		fact[++tot]=fac;return;
	}
	dfs(id+1,fac);
	for(int i=1;i<=exp[id];i++)
		fac*=fac_pri[id],dfs(id+1,fac);
}

inline void getfact(void){
	int tmp=m;
	for(int i=1;i<=cnt;i++){
		if(1LL*pri[i]*pri[i]>tmp) break;
		if(tmp%pri[i]==0){
			fac_pri[++tot_pri]=pri[i],exp[tot_pri]=0;
			while(tmp%pri[i]==0)
				exp[tot_pri]++,tmp/=pri[i];
		}
	}
	if(tmp>1) fac_pri[++tot_pri]=tmp,exp[tot_pri]=1;
	dfs(1,1);
}

signed main(void){
//	freopen("51nod_Problem_1190_Test_16_In.txt","r",stdin);
//	freopen("out.txt","w",stdout);
	scanf("%d",&cas);prework();inv=power(2,mod-2);
	while(cas--){
		scanf("%d%d",&n,&m);
		tot_pri=tot=ans=0;getfact();
		for(int i=1,tmp,ttmp,lala;i<=tot;i++){
			tmp=1;ttmp=n+fact[i]-1;lala=fact[i];
			for(int j=1;j<=tot_pri;j++)
				if(fact[i]%fac_pri[j]==0)
					tmp=1LL*tmp*(1-fac_pri[j]+mod)%mod;
			tmp=1LL*tmp*inv%mod;
			tmp=1LL*tmp*((ttmp/lala+m/lala)%mod)%mod;
			tmp=1LL*tmp*(m/lala-ttmp/lala+1)%mod;
			ans=(ans+tmp)%mod;
		}
		ans=1LL*ans*m%mod;
		printf("%d
",ans);
	}
	return 0;
}

　　

 BZOJ 3944: Sum

就是第一道题和第二道题的结合...

然后BZOJ就tm有人卡OJ了...

#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdio>
#include<map>
//by NeighThorn
using namespace std;

const int maxn=5000000+5;

bool vis[maxn];
int n,cas,cnt,pri[maxn/5];
long long f1[maxn],f2[maxn];

map<int,long long> mp1,mp2;

inline void prework(void){
	f2[1]=f1[1]=1;
	for(int i=2;i<=5000000;i++){
		if(!vis[i])
			vis[i]=1,pri[++cnt]=i,f2[i]=-1,f1[i]=i-1;
		for(int j=1;j<=cnt&&1LL*i*pri[j]<=5000000;j++){
			vis[i*pri[j]]=1;
			if(i%pri[j]==0){
				f2[i*pri[j]]=0;
				f1[i*pri[j]]=f1[i]*pri[j];
				break;
			}
			f1[i*pri[j]]=f1[i]*(pri[j]-1),f2[i*pri[j]]=-f2[i];
		}
	}
	for(int i=1;i<=5000000;i++) f1[i]+=f1[i-1],f2[i]+=f2[i-1];
}

inline long long calc1(long long n){
	if(n<=5000000) return f1[n];
	if(mp1.find(n)!=mp1.end()) return mp1[n];
	long long ans=1LL*n*(n+1)/2;
	for(long long i=2,r;i<=n;i=r+1){
		r=n/(n/i);
		ans-=(r-i+1)*calc1(n/i);
	}
	return mp1[n]=ans;
}

inline long long calc2(long long n){
	if(n<=5000000) return f2[n];
	if(mp2.find(n)!=mp2.end()) return mp2[n];
	long long ans=1;
	for(long long i=2,r;i<=n;i=r+1){
		r=n/(n/i);
		ans-=(r-i+1)*calc2(n/i);
	}
	return mp2[n]=ans;
}

signed main(void){
	scanf("%d",&cas);prework();
	while(cas--){
		scanf("%d",&n);
		printf("%lld %lld
",calc1(n),calc2(n));
	}
	return 0;
}

　　

 HDU 5608 Function

看到题目的第一眼觉得是莫比乌斯反演...

然后看了看数据范围，貌似是杜教筛...

考虑$g(n)=sum _{i=1}^{n} f(i)$

$sum _{i=1}^{n} sum _{dmid i} f(d)$

$=sum _{i=1}^{n} f(i)left lfloor frac{n}{i} ight floor $

$=sum _{i=1}^{n} g(left lfloor frac{n}{i} ight floor)$

又因为$sum _{i=1}^{n} sum _{dmid i} f(d)=sum _{i=1}^{n} (i-1)*(i-2)$

所以说$g(n)=sum _{i=1}^{n} (i-1)(i-2)-sum _{i=2}^{n} g(left lfloor frac{n}{i} ight floor)$

至于预处理我是用莫比乌斯反演写的...

$sum _{dmid n}f(d)=F(n)$

$f(n)=sum _{dmid n} mu(d) F(frac{n}{d})$

代码：

#include<algorithm> 
#include<iostream>
#include<cstring>
#include<cstdio>
#include<map>
//by NeighThorn
using namespace std;

const int mod=1e9+7,maxn=1000000+5;

int n,cas,cnt,inv,f[maxn],mu[maxn],pri[maxn],vis[maxn];
map<int,int> mp;

inline int power(int x,int y){
	int res=1;
	while(y){
		if(y&1)
			res=1LL*res*x%mod;
			x=1LL*x*x%mod,y>>=1;
	}
	return res;
}

inline int F(int x){
	return 1LL*(x-1)*(x-2)%mod;
}

inline void prework(void){
	mu[1]=1;
	for(int i=2;i<=1000000;i++){
		if(!vis[i])
			vis[i]=1,pri[++cnt]=i,mu[i]=-1;
		for(int j=1;j<=cnt&&1LL*i*pri[j]<=1000000;j++){
			vis[i*pri[j]]=1;
			if(i%pri[j]==0){
				mu[i*pri[j]]=0;
				break;
			}
			mu[i*pri[j]]=-mu[i];
		}
	}
	for(int i=1;i<=1000000;i++)
		for(int j=i;j<=1000000;j+=i){
			f[j]+=mu[i]*F(j/i);
			if(f[j]<0) f[j]+=mod;
			f[j]%=mod;
		}
	for(int i=1;i<=1000000;i++) f[i]=(f[i]+f[i-1])%mod;
}

inline int calc(int n){
	if(n<=1000000) return f[n];
	if(mp.find(n)!=mp.end()) return mp[n];
	int ans=1LL*n*(n-1)%mod*(n-2)%mod*inv%mod;
	for(int i=2,r;i<=n;i=r+1){
		r=n/(n/i);
		ans=ans-1LL*(r-i+1)*calc(n/i)%mod;
		if(ans<0) ans+=mod;
		ans%=mod;
	}
	return mp[n]=ans;
}

signed main(void){
#ifndef ONLINE_JUDGE
	freopen("in.txt","r",stdin);
	freopen("out.txt","w",stdout);
#endif
	scanf("%d",&cas);inv=power(3,mod-2);prework();
	while(cas--){
		scanf("%d",&n);
		printf("%d
",calc(n));
	}
	return 0;
}

　　

 51Nod 1238 最小公倍数之和 V3

$sum _{i=1}^{n} sum _{j=1}^{m} lcm(i,j)$

$=sum _{i=1}^{n} sum _{j=1}^{m} frac {i*j}{gcd(i,j)}$

$=sum _{d=1}^{n} d sum _{i=1}^{frac{n}{d}} sum _{j=1}^{frac{m}{d}} i*j*[gcd(i,j)==1]$

$=sum _{d=1}^{n} d sum _{i=1}^{frac{n}{d}} i sum _{j=1}^{frac{m}{d}} j*[gcd(i,j)==1]$

因为有以下的公式：

$sum _{i=1}^{n} i*[gcd(i,n)==1]=frac{nphi(n)}{2}$

证明请见HDU 3501

于是可以接着化简得到如下式子：

$sum _{d=1}^{n} sum _{i=1}^{frac{n}{d}} i*i*phi(i)$

那么我们的问题就剩下了如何快速求$Sum(n)=sum _{i=1}^{n}f(i)=sum _{i=1}^{n} i*i*phi(i)$

考虑$sum _{i=1}^{n} sum _{dmid i} phi(d) =frac{n(n+1)}{2}$

$sum _{i=1}^{n} sum _{dmid i} phi(d) i^2=frac{n^2(n+1)^2}{4}$

$=sum _{i=1}^{n} sum _{dmid i} f(d)(frac{i}{d})^2$

$=sum _{i=1}^{n} sum _{d=1}^{frac{n}{i}} f(d)i^2$

$Sum(n)=frac{n^2(n+1)^2}{4}-sum _{i=2}^{n} i^2Sum(left lfloor frac{n}{i} ight floor$

没有推出来的原因：不知道$sum _{i=1}^{n} i*[gcd(i,n)==1]=frac{nphi(n)}{2}$这个公式，推了很久发现自己的方向是错误的...

对于杜教筛的化简一定要考虑一个积性函数前缀和的形式...

代码：

#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdio>
#include<map>
//by NeighThorn
using namespace std;

const int mod=1e9+7,maxn=5000000+5;

int cnt,inv,inv2,f[maxn],phi[maxn],pri[maxn],vis[maxn];
long long n; 

map<long long,int> mp;

inline int power(int x,int y){
	int res=1;
	while(y){
		if(y&1) res=1LL*res*x%mod;
		x=1LL*x*x%mod,y>>=1;
	}
	return res;
}

inline void prework(void){
	phi[1]=1;
	for(int i=2;i<=5000000;i++){
		if(!vis[i])
			vis[i]=1,pri[++cnt]=i,phi[i]=i-1;
		for(int j=1;j<=cnt&&i*pri[j]<=5000000;j++){
			vis[i*pri[j]]=1;
			if(i%pri[j]==0){
				phi[i*pri[j]]=phi[i]*pri[j];
				break;
			}
			phi[i*pri[j]]=phi[i]*(pri[j]-1);
		}
	}
	for(int i=1;i<=5000000;i++) f[i]=1LL*i*i%mod*phi[i]%mod;
	for(int i=1;i<=5000000;i++) f[i]=(f[i]+f[i-1])%mod;
}

inline int calc(long long n){
	if(n<=5000000) return f[n];
	if(mp.find(n)!=mp.end()) return mp[n];
	int ans,tmp;tmp=n%mod;ans=1LL*tmp*(tmp+1)%mod*inv2%mod;ans=1LL*ans*ans%mod;
	for(long long i=2,r,x,y;i<=n;i=r+1){
		r=n/(n/i);x=(i-1)%mod,y=r%mod;
		tmp=1LL*y*(y+1)%mod*(y<<1|1)%mod*inv%mod;
		tmp-=1LL*x*(x+1)%mod*(x<<1|1)%mod*inv%mod;
		if(tmp<0) tmp+=mod;
		ans=(ans-1LL*tmp*calc(n/i)%mod);
		if(ans<0) ans+=mod;
	}
	return mp[n]=ans;
}

signed main(void){
	scanf("%lld",&n);prework();
	int ans=0;inv=power(6,mod-2),inv2=power(2,mod-2);
	for(long long i=1,r,x;i<=n;i=r+1){
		r=n/(n/i);x=(n/i)%mod;
		ans=(ans+1LL*x*(x+1)%mod*inv2%mod*(calc(r)-calc(i-1)+mod)%mod)%mod;
	}
	printf("%d
",ans);
	return 0;
}

　　

---

By NeighThorn