莫比乌斯反演经典例题(1)

链接:P2257 YY的GCD - 洛谷 | 计算机科学教育新生态 (luogu.com.cn)
题意:给定 n,m,求 ∑ i = 1 n ∑ j = 1 m [ g c d ( i , j ) = = p r i m e ] \sum_{i=1}^{n}\sum_{j=1}^{m}[gcd(i,j)==prime] i=1nj=1m[gcd(i,j)==prime]
题解:

  1. 首先枚举质数可化为 ∑ d ∈ p r i m e m i n ( n , m ) ∑ i = 1 n / d ∑ j = 1 m / d [ g c d ( i , j ) = = 1 ] \sum_{d \in prime }^{min(n,m)}\sum_{i=1}^{n/d}\sum_{j=1}^{m/d}[gcd(i,j)==1] dprimemin(n,m)i=1n/dj=1m/d[gcd(i,j)==1]
  2. 对于后面的式子先来求其函数的反演,
  3. f ( n , m , d ) = ∑ i = 1 n ∑ j = 1 m [ g c d ( i , j ) = = d ] f(n,m,d)=\sum_{i=1}^{n}\sum_{j=1}^{m}[gcd(i,j)==d] f(n,m,d)=i=1nj=1m[gcd(i,j)==d]
  4. F ( n , m , d ) = ∑ i = 1 n ∑ j = 1 m [ g c d ( i , j ) = = d ∗ k ] , k 为 任 意 数 F(n,m,d)=\sum_{i=1}^{n}\sum_{j=1}^{m}[gcd(i,j)==d*k],k为任意数 F(n,m,d)=i=1nj=1m[gcd(i,j)==dk],k,可以看出 F 为求 g c d ( i , j ) = = d gcd(i,j)==d gcd(i,j)==d 倍数的个数,那么 F ( n , m , d ) = ⌊ n d ⌋ ⌊ m d ⌋ F(n,m,d)=\lfloor \frac{n}{d} \rfloor \lfloor \frac{m}{d} \rfloor F(n,m,d)=dndm 可以 O ( 1 ) O(1) O(1) 取得。
  5. 考虑 f ( a , b , d ) f(a,b,d) f(a,b,d) F ( a , b , d ) F(a,b,d) F(a,b,d) 的关系。实际上, F ( a , b , d ) = ∑ d ∣ k m i n ( a , b ) f ( a , b , k ) F(a,b,d)=\sum_{d|k}^{min(a,b)}f(a,b,k) F(a,b,d)=dkmin(a,b)f(a,b,k),也就是 g c d gcd gcd 为 d 的倍数的个数都可以取。
  6. 对 F 进行反演有 f ( a , b , n ) = ∑ n ∣ d m i n ( a , b ) μ ( d n ) F ( a , b , d ) f(a,b,n)=\sum_{n|d}^{min(a,b)}\mu(\frac{d}{n})F(a,b,d) f(a,b,n)=ndmin(a,b)μ(nd)F(a,b,d)。(反演公式 f ( n ) = ∑ n ∣ d μ ( d n ) F ( d ) f(n)=\sum_{n|d}\mu(\frac{d}{n})F(d) f(n)=ndμ(nd)F(d),不熟的点这里,借用一下别人的博客)
  7. 回到 1 ,发现 ∑ i = 1 n / d ∑ j = 1 m / d [ g c d ( i , j ) = = 1 ] \sum_{i=1}^{n/d}\sum_{j=1}^{m/d}[gcd(i,j)==1] i=1n/dj=1m/d[gcd(i,j)==1] 的值为 f ( n d , m d , 1 ) = ∑ k = 1 m i n ( n d , m d ) μ ( k ) F ( n d , m d , k ) = ∑ k = 1 m i n ( n d , m d ) μ ( k ) ⌊ n k ∗ d ⌋ ⌊ m k ∗ d ⌋ f(\frac{n}{d},\frac{m}{d},1)=\sum_{k=1}^{min(\frac{n}{d},\frac{m}{d})}\mu(k)F(\frac{n}{d},\frac{m}{d},k)=\sum_{k=1}^{min(\frac{n}{d},\frac{m}{d})}\mu(k)\lfloor \frac{n}{k*d} \rfloor \lfloor \frac{m}{k*d} \rfloor f(dn,dm,1)=k=1min(dn,dm)μ(k)F(dn,dm,k)=k=1min(dn,dm)μ(k)kdnkdm (代入 4)。
  8. 答案就是 ∑ d ∈ p r i m e m i n ( n , m ) ∑ k = 1 m i n ( n d , m d ) μ ( k ) ⌊ n k ∗ d ⌋ ⌊ m k ∗ d ⌋ \sum_{d \in prime }^{min(n,m)}\sum_{k=1}^{min(\frac{n}{d},\frac{m}{d})}\mu(k)\lfloor \frac{n}{k*d} \rfloor \lfloor \frac{m}{k*d} \rfloor dprimemin(n,m)k=1min(dn,dm)μ(k)kdnkdm
  9. 很明显枚举质数这一位不削掉是不行的,令 T = k ∗ d T=k*d T=kd,那么原式子化为 ∑ d ∈ p r i m e m i n ( n , m ) ∑ T = 1 m i n ( n , m ) μ ( T d ) ⌊ n T ⌋ ⌊ m T ⌋ \sum_{d \in prime }^{min(n,m)}\sum_{T=1}^{min(n,m)}\mu(\frac{T}{d})\lfloor \frac{n}{T} \rfloor \lfloor \frac{m}{T} \rfloor dprimemin(n,m)T=1min(n,m)μ(dT)TnTm,把第二个累加往外提有 ∑ T = 1 m i n ( n , m ) ∑ d ∈ p r i m e m i n ( n , m ) μ ( T d ) ⌊ n T ⌋ ⌊ m T ⌋ \sum_{T=1}^{min(n,m)}\sum_{d \in prime }^{min(n,m)}\mu(\frac{T}{d})\lfloor \frac{n}{T} \rfloor \lfloor \frac{m}{T} \rfloor T=1min(n,m)dprimemin(n,m)μ(dT)TnTm,发现后面的两个除法与 d 无关,往外提为 ∑ T = 1 m i n ( n , m ) ⌊ n T ⌋ ⌊ m T ⌋ ∑ d ∈ p r i m e m i n ( n , m ) μ ( T d ) \sum_{T=1}^{min(n,m)}\lfloor \frac{n}{T} \rfloor \lfloor \frac{m}{T} \rfloor\sum_{d \in prime }^{min(n,m)}\mu(\frac{T}{d}) T=1min(n,m)TnTmdprimemin(n,m)μ(dT),观察 ∑ d ∈ p r i m e m i n ( n , m ) μ ( T d ) \sum_{d \in prime }^{min(n,m)}\mu(\frac{T}{d}) dprimemin(n,m)μ(dT),是可以预处理出前缀和的,具体为一个数的该累加值为其质数的倍数的 μ \mu μ 值,按质数的倍数去赋值可以做到 n l o g ( n ) nlog(n) nlog(n) 的预处理,然后在求其前缀和 t。
  10. 答案很明显了,通过整除分块求 ⌊ n T ⌋ ⌊ m T ⌋ \lfloor \frac{n}{T} \rfloor \lfloor \frac{m}{T} \rfloor TnTm 的同时乘上那一段的前缀和 t r − t l − 1 t_{r}-t_{l-1} trtl1 。总体复杂度为 O ( n l o g ( n ) + T m i n ( n , m ) ) O(nlog(n)+T\sqrt{min(n,m)}) O(nlog(n)+Tmin(n,m) ) 。需要注意的是整除分块时两个变量的块并不同步,需要取块长较短的那一个作为右端点去求。
#pragma GCC optimize("Ofast")
#pragma GCC optimize("unroll-loops")
#include<iostream>
#include<algorithm>
#include<vector>
#include<cstring>
#include<functional>
#include<queue>
#include<unordered_map>
#include<map>
#include<set>

using namespace std;
using ll=long long;
using P=pair<int,int>;
const ll inf=1e18;

struct MobiusInversion{
	static constexpr int maxn=1e7+5;
	int cnt;
	vector<int>pri,vs,mu,t;
	MobiusInversion(int x=maxn):pri(x+5),vs(x+5),mu(x+5),t(x+5){
		const int N=x-5;
		cnt=0;
		vs[0]=vs[1]=1,mu[1]=1;
		for(int i=2;i<=N;i++)
		{
			if(!vs[i])
			{
				pri[++cnt]=i,mu[i]=-1;
			}
			for(int j=1;j<=cnt&&i*pri[j]<=N;j++)
			{
				vs[i*pri[j]]=1;
				if(i%pri[j]==0)break;
				mu[i*pri[j]]=-mu[i];
			}
		}
		for(int i=1;i<=cnt;i++)
		{
			for(int j=1;pri[i]*j<=N;j++)
			{
				t[pri[i]*j]+=mu[j];
			}
		}
		for(int i=1;i<=N;i++)t[i]+=t[i-1];
	}

	ll getans(ll a,ll b){
		ll p=min(a,b);
		ll ans=0;
		for(int l=1,r;l<=p;l=r+1)
		{
			r=min(a/(a/l),b/(b/l));
			ans+=(a/l)*(b/l)*(t[r]-t[l-1]);
		}
		return ans;
	};
}tt;

void solve()
{
	ll n,m; cin>>n>>m;
	cout<<tt.getans(n,m)<<"\n";
}

int main()
{
	ios::sync_with_stdio(false);
	cin.tie(0); cout.tie(0);
	int t; cin>>t;
	while(t--)solve();
	return 0;
} 
  • 1
    点赞
  • 3
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值