莫比乌斯反演:
莫比乌斯反演变形:
假设:n<m
:gcd(i,j)=d的倍数的数对的个数;
:gcd(i,j)=d的数对的个数;
参考blog:https://blog.sengxian.com/algorithms/mobius-inversion-formula
kuangbin的模板代码不错!
*****************************************************************************************************************************************************************
题意:给出两个区间x:(1,n),y:(1,m);求gcd(x,y)=d的数对个数;比如d=1、(4,5)(5,4)算一对;
思路:
根据数学知识易知gcd(x,y)=d的倍数的数对个数为:
利用莫比乌斯反演,使用变形公式:
则,
那么计算f(d)的复杂度为O(n),优化:可以对(n/dx)*(m/dx)进行分块:
找到某一段,这一段的(n/dx)都相等,(m/dx)都相等,求取μ的前缀和,批量计算这一段的答案。
复杂度变为O(sqrt(n)+sqrt(m));
分块举例:例如,n=32,m=40,d=2; x:1~min(n,m)/d
分块代码:
last=min(n/(n/i),m/(m/i))
代码:
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;
const int N=100100;
bool check[N+10];
int prime[N+10];
int mu[N+10];
void Moblus()
{
memset(check,false,sizeof(check));
mu[1]=1;
int tot=0;
for(int i=2;i<=N;i++)
{
if(!check[i])
{
prime[tot++]=i;
mu[i]=-1;
}
for(int j=0;j<tot;j++)
{
if(i*prime[j]>N) break;
check[i*prime[j]]=true;
if(i%prime[j]==0)
{
mu[i*prime[j]]=0;
break;
}
else
{
mu[i*prime[j]]=-mu[i];
}
}
}
}
int sum[N+10];
int main()
{
int t;
Moblus();
sum[0]=0;
for(int i=1;i<=N;i++)
sum[i]=sum[i-1]+mu[i];
// for(int i=1;i<=5;i++)
// printf("%d\n",sum[i]);
scanf("%d",&t);
int co=0;
while(t--)
{
int a,n,c,m,k;
scanf("%d%d%d%d%d",&a,&n,&c,&m,&k);
if(n>m) swap(n,m);
if(k==0)
{
printf("Case %d: 0\n",++co);
continue;
}
else
n/=k,m/=k;
ll ans1=0,ans2=0,ans=0;
for(int i=1,la=1;i<=n;i=la+1)
{
la=n/(n/i); //分块;
ans1+=(ll)(sum[la]-sum[i-1])*(n/i)*(n/i);
}
ans1/=2;
// printf("%d\n",ans1);
for(int i=1,la=1;i<=n;i=la+1)
{
la=min(n/(n/i),m/(m/i)); //分块;
ans2+=(ll)(sum[la]-sum[i-1])*(n/i)*(m/i);
}
// printf("%d\n",ans2);
ans=ans2-ans1;
printf("Case %d: %I64d\n",++co,ans);
}
return 0;
}