「学习笔记」莫比乌斯反演

合集 - 学习笔记(3)

1.「学习笔记」二项式定理2023-10-09

2.「学习笔记」莫比乌斯反演2023-09-15

3.「学习笔记」初等数论2023-07-27

收起

更熟悉的阅读体验？

前置芝士：数论分块。（之后再说。QWQ）

积性函数

定义

一个数论函数 $f(n)$ 满足 $f(xy)=f(x)\times f(y)$（$gcd(x,y)=1$），则称 $f(n)$ 是积性函数。

莫比乌斯函数：

$\mu (n)=\{\begin{array}{ll}1& n=1\\ 0& n\text{含有平方因子}\\ (-1{)}^k& k\text{为}n\text{的本质不同质因子个数}\end{array}$

杜教筛

杜教筛代码真是好短啊。

一种时间复杂度为 $O(n^{\frac{2}{3}})$ 的筛，可以求出一些积性函数的前缀和。

具体来说，求 $S(n)=\sum _{i=1}^{n}f(i)$，$f$ 是积性函数。

构造函数 $h,g$ 满足 $h=f\ast g$。

那么：

$$\begin{array}{rl}\sum _{i=1}^{n}h(i)& =\sum _{i=1}^n\sum _{d|i}f(i)g(\frac{n}{i})\\ & =\sum _{d=1}^{n}g(d)\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }f(i)\\ & =\sum _{d=1}^{n}g(d)S(\lfloor \frac{n}{d}\rfloor )\end{array}$$

那么：

$$g(1)S(n)=\sum _{i=1}^{n}h(i)-\sum _{d=2}^{n}g(d)S(\lfloor \frac{n}{d}\rfloor )$$

（这里就是把 $d=1$ 的给提出来了）

$$S(n)=\frac{\sum _{i=1}^{n}h(i)-\sum _{d=2}^{n}g(d)S(\lfloor \frac{n}{d}\rfloor )}{g(1)}$$

对于 $\mu$ 函数，知道 $\mu \ast I=ϵ$，而 $ϵ$ 和 $\mu$ 的前缀和都很好求，所以可以用杜教筛。

而对于 $\phi$ 函数，知道 $\phi \ast I=id$，也可以用杜教筛。

P4213 【模板】杜教筛（Sum）

点击查看代码

#include<bits/stdc++.h>
#define XD 114514
#define MAXN 10000010
#define int long long
using namespace std;
int T,n;
int phi[MAXN],mu[MAXN],sumphi[MAXN],summu[MAXN];
bool vis[MAXN];
vector<int> prm;
map<int,int> mphi,mmu;
void sieve(){
	mu[1]=phi[1]=sumphi[1]=summu[1]=1;
	for(int i=2;i<=1e7;i++){
		if(!vis[i]){
			prm.push_back(i);
			phi[i]=i-1;mu[i]=-1;
		}
		sumphi[i]=sumphi[i-1]+phi[i];
		summu[i]=summu[i-1]+mu[i];
		for(auto j:prm){
			if(i*j>1e7) break;
			vis[i*j]=true;
			if(i%j==0){
				phi[i*j]=phi[i]*j;
				break;
			}
			phi[i*j]=phi[i]*phi[j];
			mu[i*j]=-mu[i];
		}
	}
}
int getmu(int x){
	if(x<=1e7) return summu[x];
	if(mmu[x]) return mmu[x];
	int nem=1,r;
	for(int l=2;l<=x;l=r+1){
		r=x/(x/l);
		nem-=(r-l+1)*getmu(x/l);
	}
	return mmu[x]=nem;
}
int getphi(int x){
	if(x<=1e7) return sumphi[x];
	if(mphi[x]) return mphi[x];
	int nem=x*(x+1)/2,r;
	for(int l=2;l<=x;l=r+1){
		r=x/(x/l);
		nem-=(r-l+1)*getphi(x/l);
	}
	return mphi[x]=nem;
}
signed main(){
	cin>>T;
	sieve();
	while(T--){
		cin>>n;
		cout<<getphi(n)<<" "<<getmu(n)<<"\n";
	}
	return 0;
}

莫比乌斯反演

P2398 GCD SUM

（GCDMAT - GCD OF MATRIX 双倍经验）

求：

$$\sum _{i=1}^n\sum _{j=1}^{n}gcd(i,j)$$

运用欧拉反演，可得：

$$\sum _{i=1}^n\sum _{j=1}^n\sum _{d|gcd(i,j)}\phi (d)$$

转换，得：

$$\sum _{i=1}^n\sum _{j=1}^n\sum _{d=1}^n\phi (d)[d|i][d|j]$$

$$\sum _{d=1}^n\phi (d)\sum _{i=1}^n[d|i]\sum _{j=1}^n[d|j]$$

$$\sum _{d=1}^n\phi (d)\lfloor \frac{n}{d}{\rfloor }^2$$

线性欧拉筛即可。（可以用数论分块，也可以不用，因为 $n$ 范围太小了）

不用数论分块：

点击查看代码

#include<bits/stdc++.h>
#define XD 114514
#define MAXN 100010
#define int long long
using namespace std;
int ans;
int n,phi[MAXN];
vector<int> prm;
bool vis[MAXN];
void prime(){
	phi[1]=1;
	for(int i=2;i<=n;i++){
		if(!vis[i]){
			prm.push_back(i);
			phi[i]=i-1;
		}
		for(int j=0;j<prm.size();j++){
			if(i*prm[j]<=n){
				int nem=i*prm[j];
				vis[nem]=true;
				if(i%prm[j]==0){
					phi[nem]=phi[i]*prm[j];
					break;
				}else{
					phi[nem]=phi[i]*phi[prm[j]];
				}
			}else break;
		}
	}
}
signed main(){
	cin>>n;
	prime();
	for(int i=1;i<=n;i++){
		ans+=phi[i]*(n/i)*(n/i);
	}
	cout<<ans;
	return 0;
}

用数论分块：

点击查看代码

#include<bits/stdc++.h>
#define XD 114514
#define MAXN 100010
#define int long long
using namespace std;
int ans;
int n,sum[MAXN],phi[MAXN];
vector<int> prm;
bool vis[MAXN];
void prime(){
	sum[1]=phi[1]=1;
	for(int i=2;i<=n;i++){
		if(!vis[i]){
			prm.push_back(i);
			phi[i]=i-1;
		}
		for(int j=0;j<prm.size();j++){
			if(i*prm[j]<=n){
				int nem=i*prm[j];
				vis[nem]=true;
				if(i%prm[j]==0){
					phi[nem]=phi[i]*prm[j];
					break;
				}else{
					phi[nem]=phi[i]*phi[prm[j]];
				}
			}else break;
		}
		sum[i]=sum[i-1]+phi[i];
	}
}
signed main(){
	cin>>n;
	prime();
	int r;
	for(int l=1;l<=n;l=r+1){
		int nem=n/l;
		r=n/nem;
		ans+=(sum[r]-sum[l-1])*(n/l)*(n/l);
	}
	cout<<ans;
	return 0;
}

P2303 [SDOI2012] Longge 的问题

$$\sum _{i=1}^{n}gcd(i,n)$$

$$\sum _{k|n}k\sum _{i=1}^n[gcd(i,n)=k]$$

$$\sum _{k|n}k\sum _{i=1}^n[gcd(\frac{i}{k},\frac{n}{k})=1]$$

$$\sum _{k|n}k\sum _{i=1}^{\frac{n}{k}}[gcd(i,\frac{n}{k})=1]$$

$$\sum _{k|n}k\times \phi (\frac{n}{k})$$

单次求欧拉函数即可。

点击查看代码

#include<bits/stdc++.h>
#define XD 114514
#define int long long
using namespace std;
int n,ans;
int phi(int x){
	int num=x;
	for(int i=2;i*i<=x;i++){
		if(x%i==0) num=num/i*(i-1);
		while(x%i==0) x/=i;
	}
	if(x>1) num=num/x*(x-1);
	return num;
}
signed main(){
	cin>>n;
	for(int i=1;i*i<=n;i++){
		if(i*i==n) ans+=i*phi(n/i);
		else if(n%i==0) ans+=i*phi(n/i)+(n/i)*phi(i);
	}
	cout<<ans;
	return 0;
}

P1390 公约数的和

UVA11417 GCD / UVA11424 GCD - Extreme (I) /P1390 公约数的和/ SP3871 GCDEX - GCD Extreme /UVA11426 拿行李（极限版） GCD - Extreme (II) （5 倍经验，按数据范围从小到大排序）

求：

$$\sum _{i=1}^n\sum _{j=i+1}^{n}gcd(i,j)$$

转化为：

$$\sum _{i=1}^n\sum _{j=1}^{i-1}gcd(i,j)$$

后面的直接套用上一道题的式子，可以提前线性筛欧拉函数，然后提前把答案都求出来，这样就可以把前两倍经验给切了。QWQ

点击查看代码

#include<bits/stdc++.h>
#define XD 114514
#define MAXN 200010
using namespace std;
int n;
bool vis[MAXN];
int phi[MAXN];
vector<int> prm;
void sieve(){
	phi[1]=1;
	for(int i=2;i<=2e5+1;i++){
		if(!vis[i]){
			prm.push_back(i);
			phi[i]=i-1;
		}
		for(int j=0;j<prm.size() and i*prm[j]<=2e5+1;j++){
			vis[i*prm[j]]=true;
			if(i%prm[j]==0){
				phi[i*prm[j]]=phi[i]*prm[j];
				break;
			}
			phi[i*prm[j]]=phi[i]*phi[prm[j]];
		}
	}
}
long long ans[MAXN],num;
int main(){
	ios::sync_with_stdio(false);
	cin.tie(0);cout.tie(0);
	sieve();
	for(int i=1;i<=2e5+1;i++){
		for(int j=1;j*j<=i;j++){
			if(j*j==i) num+=j*phi[j];
			else if(i%j==0) num+=1ll*j*phi[i/j]+(i/j)*phi[j];
		}
		num-=i;
		ans[i]=num;
	}
	while(true){
		cin>>n;if(n==0) break;
		cout<<ans[n]<<"\n";
	}
	return 0;
}

时间复杂度预处理 $O(n\sqrt{n})$，查询 $O(1)$。

考虑优化。

$$\sum _{i=1}^n\sum _{j=1}^{i-1}gcd(i,j)$$

$$\sum _{i=1}^n\sum _{j=1}^{i}gcd(i,j)-i$$

$$\sum _{i=1}^n\sum _{k|i}k\times \phi (\frac{i}{k})-\frac{(n+1)\times n}{2}$$

看左面的式子。

$$\sum _{i=1}^n\sum _{k=1}^{n}k\times \phi (\frac{i}{k})[k|i]$$

$$\sum _{k=1}^{n}k\sum _{i=1}^n\phi (\frac{i}{k})[k|i]$$

$$\sum _{k=1}^{n}k\sum _{i=1}^{\lfloor \frac{n}{k}\rfloor }\phi (i)$$

还是预处理线性筛欧拉函数，查询直接数论分块即可。

在代码里，为防止爆 long long，我把要减的 $\frac{(n+1)\times n}{2}$ 直接在数论分块中算了。QWQ

点击查看代码

#include<bits/stdc++.h>
#define XD 114514
#define MAXN 4000010
using namespace std;
int n;
bool vis[MAXN];
int phi[MAXN];
long long num[MAXN];
vector<int> prm;
void sieve(){
	phi[1]=1;num[1]=1;
	for(int i=2;i<=4e6;i++){
		if(!vis[i]){
			prm.push_back(i);
			phi[i]=i-1;
		}
		num[i]=num[i-1]+phi[i];
		for(int j=0;j<prm.size() and i*prm[j]<=4e6;j++){
			vis[i*prm[j]]=true;
			if(i%prm[j]==0){
				phi[i*prm[j]]=phi[i]*prm[j];
				break;
			}
			phi[i*prm[j]]=phi[i]*phi[prm[j]];
		}
	}
}
long long ans;
int main(){
	ios::sync_with_stdio(false);
	cin.tie(0);cout.tie(0);
	sieve();
	while(true){
		ans=0;
		cin>>n;
		if(n==0) break;
		int r;
		for(int l=1;l<=n;l=r+1){
			r=n/(n/l);
			ans+=(1ll*(r+1)*r/2-1ll*(l-1)*l/2)*(1ll*num[n/l]-1);
		}
		cout<<ans<<"\n";
	}
	
	return 0;
}

说句闲话，这篇题解的式子（是我上面推的式子的方法）有点问题。QWQ

膜运算

这是一个裸题，因为显然可以用欧几里得定理搞成原题。

而我才不会告诉你我想的有亿点复杂。

推式子。

P3455 [POI2007] ZAP-Queries

（P4450 双亲数 双倍经验）

求：

$$\sum _{i=1}^n\sum _{j=1}^m[gcd(i,j)=k]$$

经典把 $gcd(i,j)=k$ 转化为 $gcd(i,j)=1$。

$$\sum _{i=1}^n\sum _{j=1}^m[gcd(\frac{i}{k},\frac{j}{k})=1][i|k][j|k]$$

$$\sum _{i=1}^{\lfloor \frac{n}{k}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{k}\rfloor }[gcd(i,j)=1]$$

$$\sum _{i=1}^{\lfloor \frac{n}{k}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{k}\rfloor }\sum _{d|gcd(i,j)}\mu (d)$$

$$\sum _{i=1}^{\lfloor \frac{n}{k}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{k}\rfloor }\sum _{d|i,d|j}\mu (d)$$

$$\sum _{i=1}^{\lfloor \frac{n}{k}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{k}\rfloor }\sum _{d=1}^{min(\lfloor \frac{n}{k}\rfloor ,\lfloor \frac{m}{k}\rfloor )}\mu (d)[d|i][d|j]$$

$$\sum _{d=1}^{min(\lfloor \frac{n}{k}\rfloor ,\lfloor \frac{m}{k}\rfloor )}\mu (d)\sum _{i=1}^{\lfloor \frac{n}{k}\rfloor }[d|i]\sum _{j=1}^{\lfloor \frac{m}{k}\rfloor }[d|j]$$

$$\sum _{d=1}^{min(\lfloor \frac{n}{k}\rfloor ,\lfloor \frac{m}{k}\rfloor )}\mu (d)\lfloor \frac{\lfloor \frac{n}{k}\rfloor }{d}\rfloor \lfloor \frac{\lfloor \frac{m}{k}\rfloor }{d}\rfloor$$

线性筛一遍莫比乌斯函数再用数论分块即可。

点击查看代码

#include<bits/stdc++.h>
#define XD 114514
#define MAXN 50010
using namespace std;
int t,n,m,a,b,k,ans;
int mu[MAXN],num[MAXN];
bool vis[MAXN];
vector<int> prm;
void sieve(){
	mu[1]=1;num[1]=1;
	for(int i=2;i<=5e4;i++){
		if(!vis[i]){
			prm.push_back(i);
			mu[i]=-1;
		}
		num[i]=num[i-1]+mu[i];
		for(int j=0;j<prm.size() and i*prm[j]<=5e4;j++){
			vis[i*prm[j]]=true;
			if(i%prm[j]==0) break;
			mu[i*prm[j]]=-mu[i];
		}
	}
}
void calc(){
	int r;
	for(int l=1;l<=a;l=r+1){
		r=min(a/(a/l),b/(b/l));
		ans+=(num[r]-num[l-1])*(a/l)*(b/l);
	}
}
int main(){
	cin>>t;
	sieve();
	while(t--){
		ans=0;
		cin>>n>>m>>k;
		a=n/k;b=m/k;
		if(a>b) swap(a,b);
		calc();
		cout<<ans<<"\n";
	}
	return 0;
}

P2522 [HAOI2011] Problem b

就是用容斥原理搞一下，之后就和上面的一样了。

点击查看代码

#include<bits/stdc++.h>
#define XD 114514
#define MAXN 50010
using namespace std;
int t,n,m,a,b,c,d,k;
int mu[MAXN],num[MAXN];
bool vis[MAXN];
vector<int> prm;
void sieve(){
	mu[1]=1;num[1]=1;
	for(int i=2;i<=5e4;i++){
		if(!vis[i]){
			prm.push_back(i);
			mu[i]=-1;
		}
		num[i]=num[i-1]+mu[i];
		for(int j=0;j<prm.size() and i*prm[j]<=5e4;j++){
			vis[i*prm[j]]=true;
			if(i%prm[j]==0) break;
			mu[i*prm[j]]=-mu[i];
		}
	}
}
int calc(int x,int y){
	n=x/k;m=y/k;
	if(n>m) swap(n,m);
	int r,ans=0;
	for(int l=1;l<=n;l=r+1){
		r=min(n/(n/l),m/(m/l));
		ans+=(num[r]-num[l-1])*(n/l)*(m/l);
	}
	return ans;
}
int main(){
	sieve();
	cin>>t;
	while(t--){
		cin>>a>>b>>c>>d>>k;
		cout<<calc(b,d)-calc(a-1,d)-calc(b,c-1)+calc(a-1,c-1)<<"\n";
	}
	return 0;
}

P2257 YY的GCD 和 PGCD - Primes in GCD Table

下面正式推式子。

求：

$$\sum _{i=1}^n\sum _{j=1}^m[gcd(i,j)\in \mathrm{prime}]$$

经典把 $gcd(i,j)=k$ 转化为 $gcd(i,j)=1$。

$$\sum _{i=1}^n\sum _{j=1}^m\sum _{p\in \mathrm{prime},p|i,p|j}[gcd(\frac{i}{p},\frac{j}{p})=1]$$

$$\sum _{i=1}^n\sum _{j=1}^m\sum _{p\in \mathrm{prime}}[gcd(\frac{i}{p},\frac{j}{p})=1][p|i][p|j]$$

$$\sum _{p\in \mathrm{prime}}\sum _{i=1}^{\lfloor \frac{n}{p}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{p}\rfloor }[gcd(i,j)=1]$$

然后用莫比乌斯反演。

$$\sum _{p\in \mathrm{prime}}\sum _{i=1}^{\lfloor \frac{n}{p}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{p}\rfloor }\sum _{d|gcd(i,j)}\mu (d)$$

$$\sum _{p\in \mathrm{prime}}\sum _{i=1}^{\lfloor \frac{n}{p}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{p}\rfloor }\sum _{d|i,d|j}\mu (d)$$

$$\sum _{p\in \mathrm{prime}}\sum _{i=1}^{\lfloor \frac{n}{p}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{p}\rfloor }\sum _{d=1}^{min(\lfloor \frac{n}{p}\rfloor ,\lfloor \frac{m}{p}\rfloor )}\mu (d)[d|i][d|j]$$

这里就把 $n$ 当做 $min(n,m)$。

$$\sum _{p\in \mathrm{prime}}\sum _{i=1}^{\lfloor \frac{n}{p}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{p}\rfloor }\sum _{d=1}^{\lfloor \frac{n}{p}\rfloor }\mu (d)[d|i][d|j]$$

$$\sum _{p\in \mathrm{prime}}\sum _{d=1}^{\lfloor \frac{n}{p}\rfloor }\mu (d)\sum _{i=1}^{\lfloor \frac{n}{p}\rfloor }[d|i]\sum _{j=1}^{\lfloor \frac{m}{p}\rfloor }[d|j]$$

$$\sum _{p\in \mathrm{prime}}\sum _{d=1}^{\lfloor \frac{n}{p}\rfloor }\mu (d)\lfloor \frac{n}{pd}\rfloor \lfloor \frac{m}{pd}\rfloor$$

我们设 $T=pd$，并将式子转换成枚举 $T$。

$$\sum _{T=1}^n\lfloor \frac{n}{T}\rfloor \lfloor \frac{m}{T}\rfloor \sum _{p\in \mathrm{prime},p|T}\mu (\frac{T}{p})$$

可以发现右面的式子可以提前求出。左面的式子数论分块即可。

点击查看代码

#include<bits/stdc++.h>
#define XD 114514
#define MAXN 10000010
using namespace std;
int t,n,m;
int mu[MAXN],f[MAXN],sum[MAXN];
bool vis[MAXN];
vector<int> prm;
void sieve(){
	mu[1]=1;
	for(int i=2;i<=1e7;i++){
		if(!vis[i]){
			prm.push_back(i);
			mu[i]=-1;
		} 
		for(int j=0;j<prm.size();j++){
			if(i*prm[j]>1e7) break;
			vis[i*prm[j]]=true;
			if(i%prm[j]==0) break;
			mu[i*prm[j]]=-mu[i];
		}
	}
	for(int i=0;i<prm.size();i++){
		for(int j=1;j<=1e7;j++){
			if(prm[i]*j>1e7) break;
			f[prm[i]*j]+=mu[j];
		}
	}
	for(int i=1;i<=1e7;i++) sum[i]=sum[i-1]+f[i];
}
long long ans;
void solve(){
	int r=0;if(n>m) swap(n,m);
	for(int l=1;l<=n;l=r+1){
		int nem1=n/l,nem2=m/l;
		r=min(n/nem1,m/nem2);
		ans+=1ll*(sum[r]-sum[l-1])*(n/l)*(m/l);
	}
}
int main(){
	ios::sync_with_stdio(false);
	cin.tie(0);cout.tie(0);
	cin>>t;
	sieve();
	while(t--){
		ans=0;
		cin>>n>>m;
		solve();
		cout<<ans<<"\n";
	}
	return 0;
}

P3327 [SDOI2015] 约数个数和

求：

$$\sum _{i=1}^n\sum _{j=1}^{m}d(ij)$$

首先，要知道 $d$ 函数的一个性质：

$$d(ij)=\sum _{x|i}\sum _{y|j}[gcd(x,y)=1]$$

所以原式就变成了：

$$\sum _{i=1}^n\sum _{j=1}^m\sum _{x|i}\sum _{y|j}[gcd(x,y)=1]$$

$$\sum _{x=1}^n\lfloor \frac{n}{x}\rfloor \sum _{y=1}^m\lfloor \frac{m}{y}\rfloor [gcd(x,y)=1]$$

$$\sum _{x=1}^n\lfloor \frac{n}{x}\rfloor \sum _{y=1}^m\lfloor \frac{m}{y}\rfloor \sum _{d=1}^{min(n,m)}\mu (d)[d|x][d|y]$$

$$\sum _{d=1}^{min(n,m)}\mu (d)\sum _{x=1}^{\lfloor \frac{n}{d}\rfloor }\lfloor \frac{n}{xd}\rfloor \sum _{y=1}^{\lfloor \frac{m}{d}\rfloor }\lfloor \frac{m}{yd}\rfloor$$

设 $p=\lfloor \frac{n}{d}\rfloor$，$q=\lfloor \frac{m}{d}\rfloor$，则:

$$\sum _{d=1}^{min(n,m)}\mu (d)\sum _{x=1}^p\lfloor \frac{p}{x}\rfloor \sum _{y=1}^q\lfloor \frac{q}{y}\rfloor$$

可以发现 $\sum _{x=1}^p\lfloor \frac{p}{x}\rfloor$ 和 $\sum _{y=1}^q\lfloor \frac{q}{y}\rfloor$ 可以预处理时存在数组里，我们设 $f(x)=\sum _{i=1}^x\lfloor \frac{x}{i}\rfloor$。

$$\sum _{d=1}^{min(n,m)}\mu (d)f(p)f(q)$$

$$\sum _{d=1}^{min(n,m)}\mu (d)f(\lfloor \frac{n}{d}\rfloor )f(\lfloor \frac{m}{d}\rfloor )$$

于是查询时数论分块即可。（相当于数论分块套数论分块 QWQ）

点击查看代码

#include<bits/stdc++.h>
#define XD 114514
#define MAXN 50010
#define int long long
using namespace std;
int T,n,m,ans;
bool vis[MAXN];
int mu[MAXN],sum[MAXN],num[MAXN];
vector<int> prm;
void sieve(){
	mu[1]=1;sum[1]=1;
	for(int i=2;i<=5e4;i++){
		if(!vis[i]){
			prm.push_back(i);
			mu[i]=-1;
		}
		for(int j=0;j<prm.size() and prm[j]*i<=5e4;j++){
			vis[i*prm[j]]=true;
			if(i%prm[j]==0) break;
			mu[i*prm[j]]=-mu[i];
		}
		sum[i]=sum[i-1]+mu[i];
	}
	for(int i=1;i<=5e4;i++){
		int r,nem=0;
		for(int l=1;l<=i;l=r+1){
			r=i/(i/l);
			nem+=(r-l+1)*(i/l);
		}
		num[i]=nem;
	}
}
void calc(){
	int r;
	for(int l=1;l<=n;l=r+1){
		r=min(n/(n/l),m/(m/l));
		ans+=(sum[r]-sum[l-1])*num[n/l]*num[m/l];
	}
}
signed main(){
	cin>>T;
	sieve();
	while(T--){
		ans=0;
		cin>>n>>m;
		if(n>m) swap(n,m);
		calc();
		cout<<ans<<"\n";
	}
	return 0;
}

P3768 简单的数学题

求：

$$\sum _{i=1}^n\sum _{j=1}^{n}ijgcd(i,j)$$

$$\sum _{i=1}^n\sum _{j=1}^{n}ij\sum _{d=1}^n\phi (d)[d|i][d|j]$$

$$\sum _{d=1}^n\phi (d)\sum _{i=1}^{n}i[d|i]\sum _{j=1}^{n}j[d|j]$$

$$\sum _{d=1}^n\phi (d)\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }i\times d\sum _{j=1}^{\lfloor \frac{n}{d}\rfloor }j\times d$$

为方便，设 $m=\lfloor \frac{n}{d}\rfloor$。

$$\sum _{d=1}^n\phi (d)d^2\times \frac{m^2(m+1{)}^2}{4}$$

$\phi (n)n^2$ 这个东西要用杜教筛。

设 $f=\phi \times id\times id$，$g=id\times id$，则：

$$\begin{array}{rl}(h\ast g)(n)& =\sum _{d|n}f(d)g(\frac{n}{d})\\ & =\sum _{d|n}\phi (d)\times d^2\times (\frac{n}{d}{)}^2\\ & =n^2\times \sum _{d|n}\phi (d)\\ & =n^3\end{array}$$

所以用杜教筛后数论分块即可。

点击查看代码

#include<bits/stdc++.h>
#define XD 114514
#define MAXN 5000010 
#define int long long
using namespace std;
int p,n,ans;
int inv4,inv6;
int phi[MAXN],sumphi[MAXN];
bool vis[MAXN];
vector<int> prm;
void sieve(){
	phi[1]=sumphi[1]=1;
	for(int i=2;i<=5e6;i++){
		if(!vis[i]){
			prm.push_back(i);
			phi[i]=i-1;
		}
		sumphi[i]=(sumphi[i-1]+phi[i]*i%p*i%p)%p;
		for(auto j:prm){
			if(i*j>5e6) break;
			vis[i*j]=true;
			if(i%j==0){
				phi[i*j]=phi[i]*j;
				break;
			}
			phi[i*j]=phi[i]*phi[j];
		}
	}
}
int power(int x,int y){
	int nem=1;
	while(x){
		if(x&1) nem=(nem*y)%p;
		y=(y*y)%p;x>>=1;
	}
	return nem;
}
int pow2(int x){
	x%=p;return x*(x+1)%p*(2*x+1)%p*inv6%p;
}
int pow3(int x){
	x%=p;return x*x%p*(x+1)%p*(x+1)%p*inv4%p;
}
map<int,int> mphi;
int getphi(int x){
	if(x<=5e6) return sumphi[x];
	if(mphi[x]) return mphi[x];
	int nem=pow3(x),r;
	for(int l=2;l<=x;l=r+1){
		r=x/(x/l);
		nem=(nem-((pow2(r)-pow2(l-1)+p)%p*getphi(x/l)%p)%p+p)%p;
	}
	return mphi[x]=(nem%p+p)%p;
}
signed main(){
	cin>>p>>n;
	sieve();
	inv4=power(p-2,4);
	inv6=power(p-2,6);
	int r;
	for(int l=1;l<=n;l=r+1){
		r=n/(n/l);
		ans=(ans+((getphi(r)-getphi(l-1))%p+p)%p*pow3(n/l))%p;
	}
	cout<<(ans%p+p)%p;
	return 0;
}

P1891 疯狂 LCM 和 LCMSUM - LCM Sum

$$\sum _{i=1}^n\mathrm{lcm}(n,i)$$

$$\sum _{i=1}^n\sum _{d|n}\frac{n\times i}{d}[gcd(n,i)=d]$$

$$n\sum _{d|n}\sum _{i=1}^{\frac{n}{d}}i[gcd(i,\frac{n}{d})=1]$$

有一个式子： $\sum _{i=1}^{n}i[gcd(n,i)=1]=\frac{\phi (n)n}{2}$。

$$n\sum _{d|n}\frac{\phi (d)d}{2}$$

提前求出 $\sum _{d|n}\frac{\phi (d)d}{2}$ 的值即可。（时间复杂度 $O(n\ln n)$）

点击查看代码

#include<bits/stdc++.h>
#define XD 114514
#define MAXN 1000010
#define ll long long
using namespace std;
int t,n;
int phi[MAXN];
ll ans[MAXN];
vector<int> prm;bool vis[MAXN];
void sieve(){
	phi[1]=1;
	for(int i=2;i<=1e6;i++){
		if(!vis[i]) prm.push_back(i),phi[i]=i-1;
		for(auto j:prm){
			if(i*j>1e6) break;
			vis[i*j]=true;
			if(i%j==0){
				phi[i*j]=phi[i]*j;
				break;
			}
			phi[i*j]=phi[i]*phi[j];
		}
	}
	for(int d=1;d<=1e6;d++){
		for(int i=d;i<=1e6;i+=d){
			if(d==1){ans[i]++;continue;}
			ans[i]+=1ll*phi[d]*d/2;
		}
	}
}
int main(){
	ios::sync_with_stdio(false);
	cin.tie(0);cout.tie(0);
	sieve();
	cin>>t;
	while(t--){
		cin>>n;
		cout<<1ll*n*ans[n]<<"\n";
	}
	return 0;
}

P1829 [国家集训队] Crash的数字表格 / JZPTAB

设 $n<m$。求：

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

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

为方便，设 $x=\lfloor \frac{n}{d}\rfloor$，$y=\lfloor \frac{m}{d}\rfloor$。

$$\sum _{d=1}^n\sum _{i=1}^x\sum _{j=1}^{y}d\times i\times j[gcd(i,j)=1]$$

$$\sum _{d=1}^{n}d\sum _{i=1}^x\sum _{j=1}^{y}i\times j\sum _{k|i,k|j}\mu (k)$$

$$\sum _{d=1}^{n}d\sum _{k=1}^x\mu (k)\times k\times k\sum _{i=1}^{\lfloor \frac{x}{d}\rfloor }i\sum _{j=1}^{\lfloor \frac{y}{d}\rfloor }j$$

设 $\mathrm{sum}(x)=\sum _{i=1}^{x}i$，$f(x,y)=\sum _{k=1}^x\mu (k)\times k\times k\times \mathrm{sum}(\lfloor \frac{x}{d}\rfloor )\times \mathrm{sum}(\lfloor \frac{y}{d}\rfloor )$。

$$\sum _{d=1}^{n}d\times f(x,y)$$

数论分块套数论分块即可。

点击查看代码

#include<bits/stdc++.h>
#define XD 114514
#define MAXN 10000010
#define int long long
using namespace std;
const int mod=20101009;
int t,n,m,ans;
int mu[MAXN],num[MAXN];
bool vis[MAXN];
vector<int> prm;
void sieve(){
	mu[1]=1;num[1]=1;
	for(int i=2;i<=1e7;i++){
		if(!vis[i]){
			prm.push_back(i);
			mu[i]=-1;
		}
		num[i]=(num[i-1]+mu[i]*i%mod*i%mod)%mod;
		for(int j=0;j<prm.size() and i*prm[j]<=1e7;j++){
			vis[i*prm[j]]=true;
			if(i%prm[j]==0) break;
			mu[i*prm[j]]=-mu[i];
		}
	}
}
int getnum(int l,int r){return ((r+1)*r/2%mod-(l-1)*l/2%mod+mod)%mod;}
int f(int x,int y){
	int r,nem=0;
	for(int l=1;l<=x;l=r+1){
		r=min(x/(x/l),y/(y/l));
		nem=(nem+(num[r]-num[l-1])%mod*getnum(1,x/l)%mod*getnum(1,y/l))%mod;
	}
	return nem%mod;
}
void calc(){
	int r;
	for(int l=1;l<=n;l=r+1){
		r=min(n/(n/l),m/(m/l));
		ans=(ans+getnum(l,r)*f(n/l,m/l))%mod;
	}
}
signed main(){
	sieve();
	cin>>n>>m;
	if(n>m) swap(n,m);
	calc();
	cout<<(ans%mod+mod)%mod;
	return 0;
}

P4917 天守阁的地板 和 P5221 Product

P4917 天守阁的地板

Product 凉心出题人又卡时间又卡空间

P3172 [CQOI2015] 选数

$$\sum _{a_1=L}^R\sum _{a_2=L}^R...\sum _{a_n=L}^R[\underset{i=1}{\overset{n}{gcd}}{a}_i=k]$$

设 $l=\lceil \frac{L}{k}\rceil$，$r=\lfloor \frac{R}{k}\rfloor$。

$$\sum _{a_1=l}^r\sum _{a_2=l}^r...\sum _{a_n=l}^r[\underset{i=1}{\overset{n}{gcd}}{a}_i=1]$$

$$\sum _{d=1}^r\mu (d)\sum _{a_1=l}^r[d|a_1]\sum _{a_2=l}^r[d|a_2]...\sum _{a_n=l}^r[d|a_n]$$

$$\sum _{d=1}^r\mu (d)(\lfloor \frac{r}{d}\rfloor -\lfloor \frac{l-1}{d}\rfloor {)}^n$$

点击查看代码

#include<bits/stdc++.h>
#define XD 114514
#define MAXN 10000010
#define int long long
using namespace std;
const int mod=1e9+7;
int n,k,L,R,ans;
int mu[MAXN],musum[MAXN];bool vis[MAXN];
vector<int> prm;
map<int,int> mmu;
void sieve(){
	mu[1]=musum[1]=1;
	for(int i=2;i<=5e6;i++){
		if(!vis[i]) prm.push_back(i),mu[i]=-1;
		musum[i]=musum[i-1]+mu[i];
		for(auto j:prm){
			if(i*j>5e6) break;
			vis[i*j]=true;
			if(i%j==0) break;
			mu[i*j]=-mu[i];
		}
	}
}
int getmu(int x){
	if(x<=5e6) return musum[x];
	if(mmu.count(x)) return mmu[x];
	int nem=1,r;
	for(int l=2;l<=x;l=r+1){
		r=x/(x/l);
		nem-=getmu(x/l)*(r-l+1);
	}
	return mmu[x]=nem;
}
int power(int x,int y){
	int nem=1;
	while(y){
		if(y&1) nem=nem*x%mod;
		x=x*x%mod;y>>=1;
	}
	return nem;
}
void calc(){
	int r;L--;
	for(int l=1;l<=R;l=r+1){
		if(!(L/l)) r=R/(R/l);
		else r=min(R/(R/l),L/(L/l));
		ans=(ans+power(R/l-L/l,n)*(getmu(r)-getmu(l-1))%mod)%mod;
	}
}
signed main(){
	sieve();
	cin>>n>>k>>L>>R;
	L=L/k+(L%k?1:0);R=R/k;
	calc();
	cout<<(ans%mod+mod)%mod;
	return 0;
}

P6055 [RC-02] GCD

$$\sum _{i=1}^n\sum _{j=1}^n\sum _{p=1}^{\lfloor \frac{n}{j}\rfloor }\sum _{q=1}^{\lfloor \frac{n}{j}\rfloor }[gcd(i,j)=1][gcd(p,q)=1]$$

$$\sum _{i=1}^n\sum _{j=1}^n\sum _{p=1}^j\sum _{q=1}^j[gcd(i,j)=1][gcd(p,q)=j]$$

$$\sum _{i=1}^n\sum _{p=1}^n\sum _{q=1}^n[gcd(i,p,q)=1]$$

$$\sum _{d=1}^n\mu (d){\lfloor \frac{n}{d}\rfloor }^3$$

点击查看代码

#include<bits/stdc++.h>
#define XD 114514
#define MAXN 5000010
#define int long long
using namespace std;
const int mod=998244353;
int n;
int mu[MAXN],musum[MAXN];bool vis[MAXN];
vector<int> prm;
map<int,int> mmu;
void sieve(){
	mu[1]=musum[1]=1;
	for(int i=2;i<=5e6;i++){
		if(!vis[i]) prm.push_back(i),mu[i]=-1;
		musum[i]=musum[i-1]+mu[i];
		for(auto j:prm){
			if(i*j>5e6) break;
			vis[i*j]=true;
			if(i%j==0) break;
			mu[i*j]=-mu[i];
		}
	}
}
int getmu(int x){
	if(x<=5e6) return musum[x];
	if(mmu.count(x)) return mmu[x];
	int nem=1,r;
	for(int l=2;l<=x;l=r+1){
		r=x/(x/l);
		nem-=(r-l+1)*getmu(x/l);
	}
	return mmu[x]=nem;
}
int ans;
signed main(){
	sieve();
	cin>>n;
	int r;
	for(int l=1;l<=n;l=r+1){
		r=n/(n/l);
		ans=(ans+(getmu(r)-getmu(l-1))%mod*(n/l)%mod*(n/l)%mod*(n/l)%mod)%mod;
	}
	cout<<(ans%mod+mod)%mod;
	return 0;
}

P4449 于神之怒加强版

假设 $n<m$。

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

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

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

$$\sum _{d=1}^n{d}^k\sum _{p=1}^{\lfloor \frac{n}{d}\rfloor }\mu (p)\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }[p|i]\sum _{j=1}^{\lfloor \frac{m}{d}\rfloor }[p|j]$$

$$\sum _{d=1}^n{d}^k\sum _{p=1}^{\lfloor \frac{n}{d}\rfloor }\mu (p)\lfloor \frac{n}{dp}\rfloor \lfloor \frac{m}{dp}\rfloor$$

设 $T=dp$。

$$\sum _{T=1}^n\lfloor \frac{n}{T}\rfloor \lfloor \frac{m}{T}\rfloor \sum _{d|T}{d}^k\mu (\frac{T}{d})$$

不难发现 $\sum _{d|T}{d}^k\mu (\frac{T}{d})$ 是个积性函数，所以可以像求欧拉函数一样用线性筛求出。

具体来说，设 $f(T)=\sum _{d|T}{d}^k\mu (\frac{T}{d})$，则：

$$f(ab)=\{\begin{array}{l}f(a)\times f(b)gcd(a,b)=1\\ f(a)\times b^{k}gcd(a,b)\ne 1,b\in \mathbb{P}\end{array}$$

然后知道 $f(1)=1$，$f(p)=p^k-1$（$p\in \mathbb{P}$），所以就可以提前线性筛出来了。

点击查看代码

#include<bits/stdc++.h>
#define XD 114514
#define MAXN 5000010
#define int long long
using namespace std;
const int mod=1e9+7;
int t,n,m,k,ans;
vector<int> prm;bool vis[MAXN];
int sum[MAXN],num[MAXN];
int power(int x,int y){
	int nem=1;
	while(y){
		if(y&1) nem=nem*x%mod;
		x=x*x%mod;y>>=1;
	}
	return nem;
}
void sieve(){
	num[1]=sum[1]=1;
	for(int i=2;i<=5e6;i++){
		if(!vis[i]){
			prm.push_back(i);
			num[i]=power(i,k)-1;
		}
		sum[i]=(sum[i-1]+num[i])%mod;
		for(auto j:prm){
			if(i*j>5e6) break;
			vis[i*j]=true;
			if(i%j==0){
				num[i*j]=num[i]*(num[j]+1)%mod;
				break;
			} 
			num[i*j]=num[i]*num[j]%mod;
		}
	}
}
void calc(){
	int r;ans=0;
	for(int l=1;l<=n;l=r+1){
		r=min(n/(n/l),m/(m/l));
		ans=(ans+(sum[r]-sum[l-1])*(n/l)%mod*(m/l)%mod+mod)%mod;
	}
}
signed main(){
	ios::sync_with_stdio(false);
	cin.tie(0);cout.tie(0);
	cin>>t>>k;
	sieve();
	while(t--){
		cin>>n>>m;
		if(n>m) swap(n,m);
		calc();
		cout<<(ans+mod)%mod<<"\n";
	}
	return 0;
}

学习建议

建议看以下的博客，讲的非常好，反正我的也看不懂。QWQ

浅谈莫反

狄利克雷卷积和莫比乌斯反演

（我本人最推荐这两个。QWQ）

莫比乌斯反演-让我们从基础开始（洛谷日报上的）

「笔记」积性函数

莫比乌斯反演

peng-ym 的杜教筛

算法学习笔记(35): 狄利克雷卷积

__EOF__

本文作者：wdgm4
本文链接：https://www.cnblogs.com/wdgm4/p/17704213.html
关于博主：评论和私信会在第一时间回复。或者直接私信我。
版权声明：本博客所有文章除特别声明外，均采用 BY-NC-SA 许可协议。转载请注明出处！
声援博主：如果您觉得文章对您有帮助，可以点击文章右下角【推荐】一下。您的鼓励是博主的最大动力！