luogu P3768 简单的数学题 杜教筛 + 欧拉反演 + 逆元

求 $sum_{i=1}^{n}sum_{j=1}^{n}ijgcd(i,j)$

 

考虑欧拉反演： $sum_{d|n}varphi(d)=n$

 

$Rightarrow sum_{i=1}^{n}sum_{j=1}^{n}ijsum_{d|gcd(i,j)}varphi(d)$

 

$Rightarrow sum_{i=1}^{n}sum_{j=1}^{n}ijsum_{d|i,d|j}varphi(d)$

 

$Rightarrow sum_{d=1}^{n}varphi(d)sum_{d|i}sum_{d|j}ij$

 

$Rightarrowsum_{d=1}^{n} varphi(d)d^2sum_{i=1}^{frac{n}{d}}isum_{j=1}^{frac{n}{d}}j$

 

对于 $sum_{i=1}^{frac{n}{d}}isum_{j=1}^{frac{n}{d}}j$，直接 $O(1)$求

 

令 $calc(n,m)=frac{n(n+1)}{2} imes frac{m(m+1)}{2}$

 

将 $frac{n}{d}$ 带入即可.

 

原式可化为 $sum_{d=1}^{n} varphi(d)d^2 imes calc(frac{n}{d},frac{n}{d})$

 

这个复杂度是 $O(sqrt n)$ 的.

 

然而，有一个问题：我们正常只能求出 $[1,1e7]$ 的欧拉函数值.

 

于是，要用杜教筛优化一下.

 

令 $f(x)=varphi(x)x^2$

 

搬出杜教是公式： $S(n)=frac{sum_{i=1}^{n}(f*g)(i)-sum_{i=2}^{n}g(i)S(frac{n}{i})}{g(1)}$

 

$f(x)$ 是固定的，现在只需选一个合适的 $g(x)$，使得可以快速算出 $sum_{i=1}^{n}(f*g)(i)$

 

$(f*g)(i)=sum_{d|i}f(d)g(frac{i}{d})$

 

$Rightarrow sum_{d|i}varphi(d)d^2g(frac{i}{d})$

 

$d^2$ 有点讨厌，要是能消掉就好了.

 

令 $g(x)=x^2$

 

$Rightarrow sum_{d|i}varphi(d)d^2(frac{i}{d})^2$

 

$Rightarrow sum_{d|i}varphi(d)i^2$

 

$Rightarrow i^2 imessum_{d|i}varphi(d)$ 

 

$Rightarrow i^3$    

 

将 $(f*g)(i)$ 的结果带入杜教筛公式中.

 

即当 $f(x)=varphi(x)x^2$, $g(x)=x^2$ 时，

 

$S(n)=sum_{i=1}^{n}i^3-sum_{i=2}^{n}i^2 imes S(frac{n}{i})$

 

搬出高中数学公式：

 

$sum_{i=1}^{n}i^2=frac{n(n+1)(2n+1)}{6}$

 

$sum_{i=1}^{n}i^3=frac{n^2(n+1)^2}{4}$

 

回到原式 $sum_{d=1}^{n} varphi(d)d^2 imes calc(frac{n}{d},frac{n}{d})$

 

直接用杜教筛算 $varphi(d)d^2$ 的前缀和并用整除分块算出结果即可. 

#include<bits/stdc++.h>
#define maxn 10200006
#define ll long long 
#define M 10000007
using namespace std; 
int cnt; 
ll sumv[maxn], rev4, rev6, mod, rev2; 
bool vis[maxn]; 
ll  phi[maxn], prime[maxn]; 
map<ll,ll>ansphi; 
void setIO(string s)
{
	string in=s+".in"; 
	freopen(in.c_str(),"r",stdin); 
}
ll qpow(ll base, ll k)
{
	ll tmp=1;
	while(k)
	{
		if(k&1) tmp=tmp*base%mod; 
		base=base*base%mod; 
		k>>=1; 
	}
	return tmp; 
}
void init()
{
	int i,j;                   
	rev4=qpow(4ll, mod-2), rev6=qpow(6ll, mod-2), rev2=qpow(2ll, mod-2); 
	phi[1]=1;             
	for(i=2;i<=M;++i) 
	{       
		if(!vis[i]) prime[++cnt]=i, phi[i]=i-1; 
		for(j=1;j<=cnt&&1ll*i*prime[j]<=M;++j) 
		{
			vis[i*prime[j]]=1; 
			if(i%prime[j]==0) 
			{                                  
				phi[i*prime[j]]=phi[i]*prime[j]; 
				break; 
			}
			phi[i*prime[j]]=phi[i]*(prime[j]-1); 
		}
	}
	for(i=1;i<=M;++i) sumv[i]=(sumv[i-1]+(1ll*phi[i]*i%mod*i%mod))%mod; 
}
// 平方 
ll cal1(ll i)
{
	i%=mod; 
	ll re=i%mod; 
	re=re*(i+1)%mod; 
	re=re*(i+i+1)%mod; 
	re=(re*rev6)%mod; 
	return re; 
}
// 立方 
ll cal2(ll i)
{
	i%=mod; 
	ll re=i%mod; 
	re=(re*i)%mod; 
	re=(re*(i+1))%mod; 
	re=re*(i+1)%mod; 
	re=(re*rev4)%mod; 
	return re;  
}  
ll get(ll n)
{
	if(n<=M) return sumv[n];     
	if(ansphi[n]) return ansphi[n];   
	ll i,j,re=cal2(n),tmp; 
	for(i=2;i<=n;i=j+1)
	{   
		j=n/(n/i);         
		tmp=(cal1(j)-cal1(i-1)+mod)%mod; 
		tmp=(tmp*get(n/i))%mod;  
		re=(re-tmp+mod)%mod; 
	}
	return ansphi[n]=re;   
}
ll calc(ll n)
{
	n%=mod; 
	return (((n*(n+1))%mod)*(rev2%mod))%mod ; 
}
int main()
{
	// setIO("input");         
	ll n,i,j,re=0,tmp=0; 
	scanf("%lld%lld",&mod,&n);        
	init(); 
	for(i=1;i<=n;i=j+1)
	{
		j=n/(n/i);                   
		tmp=(calc(n/i)*calc(n/i)%mod*(get(j)-get(i-1)+mod)%mod)%mod;   
		re=(re+tmp+mod)%mod;  
	}
	printf("%lld
",re); 
	return 0; 
}