Jzoj5425 数论

聪明的0v0正在学习莫比乌斯反演。
她看到了这样的一道题:有n*m个人站成了一个n*m的方阵……

剩下的题面,聪明的0v0不记得了。

但是,她通过自己高超的数论技巧,给出了一个转化后的模型:给出n和m,求

ΣΣmin(n/i,m/j)*[gcd(i,j)=1]{1<=i<=n,1<=j<=m}

聪明的0v0当然知道怎么做了,但是她想考考你。

是个反演题也,可以用杜教筛也可以分块

不会繁衍反演怎么办?

我们发现,gcd(i,j)=1可以视为代表一种斜率

而min(n/i,m/j)是一条直线在二维区间[1,n][1,m]的整点个数

所以答案就是n*m

code就贴一个大佬写得反演吧,code fromBAJim_H

#include <cstdio>
#include <algorithm>
#include <cstring>
#include <iostream>
#include <cstdlib>
#include <cmath>
#define fo(i,a,b) for(int i=a;i<=b;i++)
#define fod(i,a,b) for(int i=a;i>=b;i--)
#define N 10000005
#define LL long long
using namespace std;
LL mo,pr[N],su[N],s[N];
int n,m,l;
bool bz[N];
void prp()
{
	su[1]=1;
	fo(i,2,n)
	{
		if(!bz[i]) pr[++l]=i,su[i]=-1;
		for(int j=1;j<=l&&pr[j]*i<=n;j++)
		{
			bz[i*pr[j]]=1;
			if(i%pr[j]==0) break;
			su[i*pr[j]]=-su[i];
		}
		(su[i]+=su[i-1])%=mo;
	}
}
LL get(LL n,LL m)
{
	if(n==0) return 0;
	LL d1=m/(m/n);
	return (s[d1]-(d1-n)*(m/n)%mo)%mo;
	
}
void fd(LL m)
{
	s[0]=0;
	LL d=1;
	while(d<=m)
	{
		LL d1=m/(m/d);
		s[d1]=(s[d-1]+(d1-d+1)*(m/d)%mo)%mo;
		d=d1+1;
	}
}
int main()
{
	freopen("math.in","r",stdin);
	freopen("math.out","w",stdout);
	cin>>n>>m>>mo;
	if(n<m) swap(n,m);
	prp();
	LL d=1,ans=0;
	while(d<=m)
	{
		LL d1=min(n/(n/d),m/(m/d)),i=1,s1=0;
		fd(m/d);
		while(i<=n/d)
		{
			LL i1=(n/d)/(n/d/i),j=(m/d)/(n/d/i);
			(s1+=(LL)(i1-i+1)%mo*(j*(n/d/i)%mo+get(m/d,m/d)-get(j,m/d)))%=mo;
			i=i1+1;
		}
		(ans+=s1*(su[d1]-su[d-1])%mo)%=mo;
		d=d1+1;
	}
	printf("%lld
",(ans+mo)%mo);
}

原文地址:https://www.cnblogs.com/Extended-Ash/p/7846037.html