【整除分块&莫比乌斯反演】洛谷P2522 Problem b

搂鱼114514

已于 2024-03-10 10:35:41 修改

阅读量158

点赞数 1

文章标签：算法

于 2024-03-09 14:11:46 首次发布

本文链接：https://blog.csdn.net/m0_60506105/article/details/136581761

版权

题目

给定a,b,c,d,k，求：

$\sum_{i=a}^{b}\sum_{i=c}^{d}[gcd(i,j)=k]$

思路

显然两段区间所产生的的答案可以用容斥原理“归一”处理。

将两段区间用容斥原理处理，那么仅需考虑：

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

同时除以k：

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

这时候大名鼎鼎的莫比乌斯反演出场了，运用狄利克雷卷积可以推导出反演结论：

$[gcd(i,j)=1]=\sum_{d|gcd(i,j)}\mu (d)$

带入反演结论得：

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

交换求和顺序得：

$\sum_{d=1}\mu (d)\sum_{i=1}^{\left \lfloor \frac{n}{k} \right \rfloor}[d|i]\sum_{j=1}^{\left \lfloor \frac{m}{k} \right \rfloor}[d|j]$

因为 $[1,\left \lfloor \frac{n}{k} \right \rfloor]$ 中有 $\left \lfloor \frac{n}{kd} \right \rfloor$ 个d的倍数， $[1,\left \lfloor \frac{m}{k} \right \rfloor]$ 中有 $\left \lfloor \frac{m}{kd} \right \rfloor$ 个d的倍数，所以可以把后面两个大 $\sum$ 省掉：

$\sum_{d=1}^{min(\left \lfloor \frac{n}{k} \right \rfloor,\left \lfloor \frac{m}{k} \right \rfloor)}\mu (d)\left \lfloor \frac{n}{kd} \right \rfloor\left \lfloor \frac{m}{kd} \right \rfloor$

对于 $\left \lfloor \frac{n}{kd} \right \rfloor$ 和 $\left \lfloor \frac{m}{kd} \right \rfloor$ ，可以使用整除分块处理，详参上一篇题解。

代码

#pragma GCC optimise(2)
#include<bits/stdc++.h>
using namespace std;
#define ll long long
const ll z=1e7+9,N=1e7+9;
ll mu[z],p[z],np,f[z],xb[z];
bool isp[z];
ll n,a,b,c,d,k,ans;
void getmu()
{
	memset(isp,1,sizeof(isp));
	isp[1]=0;
	mu[1]=1;
	for(ll i=2;i<=N;i++)
	{
		if(isp[i])
		{
			p[++np]=i;
			mu[i]=-1;
		}
		for(ll j=1;j<=np&&i*p[j]<=N;j++)
		{
			ll x=i*p[j];
			isp[x]=0;
			if(i%p[j]!=0)mu[x]=-mu[i];
			else 
			{
				mu[x]=0;
				break;
			}
		}
	}
	for(ll i=1;i<=np;i++)
	for(ll j=1;j*p[i]<=N;j++)
	f[j*p[i]]+=mu[j];//考虑每一个质数k，对于k的倍数kk，将其加上mu[kk/k]
	for(ll i=1;i<=N;i++)
	xb[i]=xb[i-1]+f[i];
}
ll cal(ll n,ll m)
{
	ll res=0,rr=min(n,m);
	for(ll l=1,r;l<=rr;l=r+1)
	{
		r=min(n/(n/l),m/(m/l));
		res+=(xb[r]-xb[l-1])*(n/l)*(m/l);
	}
	return res;
}
int main()
{
	scanf("%lld",&n);
	getmu();
	for(int i=1;i<=n;i++)
	{
		scanf("%lld%lld",&b,&d);
		if(b>d)swap(b,d);
		ans=cal(b,d);
		printf("%lld\n",ans);
	}
	return 0;
}