洛谷 P3768 简单的数学题 杜教筛

P3768

虽然刷数论老是把智商刷掉线，但是比起上百行代码的数据结构，我还是更喜欢数论....

我们要求的是这个：

                                                

我们干脆枚举gcd的值吧：

                                                

枚举两数看两数gcd是否为d，要不我直接从d的倍数开始枚举这两个数吧：

                               

由于莫比乌斯函数性质，我们代换：

                                         

从枚举gcd(i,j)的因数换成枚举x，然后i,j分别表示x的倍数吧：

                                         

设sum(x)为1,2,3...到x的和：

                                       

要不我们先直接从1开始枚举dx的值，然后再枚举dx的因数吧，设T=dx：

                                      

前面那块好像可以用整除分块了，现在我们重点处理后面那一部分的前缀和，化简：

                                     

有一个公式：

                                                   

代换：

                                                    

这题需要的前缀和有1e10，因此我们只打表1e7的前缀和，其他的用杜教筛搞定，杜教筛：

                                       

设g(x)=x^2，则h=f 卷积 g：（用了欧拉函数性质：φ 卷积 I=id ）

                                       

设S2(n)=1^3+2^3+3^3+...+n^3，我们容易知道S2(n)=(1+2+3+..+n)^2，套一下杜教筛公式：

                                           

1^2+2^2+..+n^2=n*(n+1)*(2n+1)/6，那么我们就可以用杜教筛快速求S(T)了，配合整除分块，就搞定这题啦。

                                                

#include<bits/stdc++.h>
#include<tr1/unordered_map>
#define ll long long
using namespace std;
const int maxn=1e7+5,N=1e7;
tr1::unordered_map<ll,int>s;
int vis[maxn],pri[maxn],phi[maxn],cnt,inv6,inv2,mod;
void init(int n)
{
	phi[1]=1;
	for(int i=2;i<=n;i++)
	{
		if(!vis[i])
		phi[i]=i-1,pri[++cnt]=i;
		for(int j=1;j<=cnt&&pri[j]*i<=n;j++)
		{
			vis[pri[j]*i]=1;
			if(i%pri[j]==0)
			{
				phi[i*pri[j]]=1ll*phi[i]*pri[j]%mod;
				break;
			}
			phi[i*pri[j]]=1ll*phi[i]*(pri[j]-1)%mod;
		}
	}
	for(int i=1;i<=n;i++)
	phi[i]=(phi[i-1]+1ll*i*i%mod*phi[i])%mod;
}
ll ksm(ll x,int y)
{
	ll res=1;
	while(y)
	{
		if(y&1)res=res*x%mod;
		x=x*x%mod;
		y/=2;
	}
	return res;
}
ll cal(ll n)
{
	n%=mod;
	return n*(n+1)%mod*(2*n+1)%mod*inv6%mod;
}
ll cal2(ll n)
{
	n%=mod;
	return n*(n+1)%mod*inv2%mod;
}
ll dfs(ll n)
{
	if(n<=N)return phi[n];
	if(s[n])return s[n];
	ll ans=cal2(n)*cal2(n)%mod,res=0;
	for(ll l=2,r;l<=n;l=r+1)
	{
		r=n/(n/l);
		res+=(cal(r)-cal(l-1)+mod)*dfs(n/l)%mod;
	}
	s[n]=(ans-res%mod+mod)%mod;
	return s[n];
}
int main()
{
	ll n,ans=0;
	scanf("%d%lld",&mod,&n);
	inv6=ksm(6,mod-2);
	inv2=ksm(2,mod-2);
	init(N);
	for(ll l=1,r;l<=n;l=r+1)
	{
		r=n/(n/l);
		ans+=cal2(n/l)*cal2(n/l)%mod*(dfs(r)-dfs(l-1)+mod)%mod;
	}
	printf("%lld\n",ans%mod);
}