题目链接

题目描述

给定 \(n,p\)，求 \(\left(\sum_{i=1}^{n}{\sum_{j=1}^{n}{i\times j\times\gcd(i,j)}}\right)\bmod p\)。
数据范围：\(n\le {10}^{10}\)，\(5\times 10^8\le p\le 10^9\).
时间范围：\(4000\operatorname{ms}\)。

Solution

魔怔推式子：

\[\begin{aligned} &\sum_{i=1}^{n}{\sum_{j=1}^{n}{i\times j\times\gcd(i,j)}}\\ =&\sum_{d=1}^{n}{\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}{\sum_{j=1}^{\lfloor\frac{n}{d}\rfloor}{i\times j\times d^3\times [(i,j)=1]}}}\\ =&\sum_{d=1}^{n}{d^3\times\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}{i\times\sum_{j=1}^{\lfloor\frac{n}{d}\rfloor}{j\times[(i,j)=1]}}}\\ =&\sum_{d=1}^{n}{d^3\times\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}{i\times\sum_{j=1}^{\lfloor\frac{n}{d}\rfloor}{j\times\sum_{k\mid i\wedge k\mid j}{\mu(k)}}}}\\ =&\sum_{d=1}^{n}{d^3\times \sum_{k=1}^{\lfloor\frac{n}{d}\rfloor}{\mu(k)\times k^2\times\left(\sum_{i=1}^{\lfloor\frac{n}{dk}\rfloor}{i}\right)^2}}\\ =&\sum_{i=1}^{n}{s^2(\lfloor\frac{n}{i}\rfloor)}\sum_{d\mid i}{d^3\times(\frac{i}{d})^2\times \mu(\frac{i}{d})}\\ =&\sum_{i=1}^{n}{s^2(\lfloor\frac{n}{i}\rfloor)\times i^2\sum_{d\mid i}{d\times \mu(\frac{i}{d})}}\\ =&\sum_{i=1}^{n}{s^2(\lfloor\frac{n}{i}\rfloor)\times i^2\times \varphi(i)} \end{aligned}\]

其中 \(s(i)=\sum_{i=1}^{n}{i}\)。
考虑怎么快速筛出 \(f(i)=i^2\times\varphi(i)\)。
考虑杜教筛：由于有 \(i^2\) 项，令 \(g(i)=i^2\)，那么有

\[\begin{aligned} (g*f)(n)&=\sum_{d\mid n}{d^2\times\frac{n^2}{d^2}\times\varphi(\frac{n}{d})}\\ &=n^2\sum_{d\mid n}{\varphi(\frac{n}{d})}\\ &=n^3\\ \end{aligned}\]

\(\sum_{i=1}^{n}{i^3}\) 的前缀和为 \(s^2(n)\)，可以杜教筛 \(O(n^{\frac{2}{3}})\) 快速求解。
最后对原式套上数论分块即可求出最后答案。

Code

#include<cstdio>
#include<unordered_map>
using namespace std;
typedef long long LL;
const int maxn=10000010;
int pri[maxn],tot,phi[maxn],sum[maxn],MLY,inv2,inv6;
unordered_map<LL,int> ssum;
bool vis[maxn];
LL n;
inline int power(int a,int b){
	int ans=1;
	while(b){
		if(b&1)ans=1ll*ans*a%MLY;
		a=1ll*a*a%MLY;
		b>>=1;
	}
	return ans;
}
inline void Init(){
	inv2=power(2,MLY-2);inv6=power(6,MLY-2);
	phi[1]=1;
	for(int i=2;i<maxn;++i){
		if(!vis[i])pri[++tot]=i,phi[i]=i-1;
		for(int j=1;i*pri[j]<maxn;++j){
			vis[i*pri[j]]=1;
			if(i%pri[j])phi[i*pri[j]]=phi[i]*(pri[j]-1);
			else{
				phi[i*pri[j]]=phi[i]*pri[j];
				break;
			}
		}
	}
	for(int i=1;i<maxn;++i)
		sum[i]=(sum[i-1]+1ll*i*i%MLY*phi[i])%MLY;
}
inline LL get(LL n,LL l){
	return n/(n/l);
}
inline int S1(LL n){
	n%=MLY;
	return n%MLY*(n+1)%MLY*inv2%MLY;
}
inline int S2(LL n){
	n%=MLY;
	return n*(n+1)%MLY*(2*n+1)%MLY*inv6%MLY;
}
inline int S3(LL n){
	return 1ll*S1(n)*S1(n)%MLY;
}
inline int Sum(LL n){
	if(n<maxn)return sum[n];
	else if(ssum.count(n))return ssum[n];
	int ans=S3(n);
	for(LL l=2,r;r<n;l=r+1){
		r=get(n,l);
		ans=(ans-1ll*(S2(r)-S2(l-1))*Sum(n/l)%MLY+MLY)%MLY;
	}
	return ssum[n]=ans;
}
int main(){
	scanf("%d%lld",&MLY,&n);
	Init();
	int ans=0;
	for(LL l=1,r;r<n;l=r+1){
		r=get(n,l);
		ans=(ans+1ll*S1(n/l)*S1(n/l)%MLY*(Sum(r)-Sum(l-1))%MLY+MLY)%MLY;
	}
	printf("%d",ans);
	return 0;
}