P3768 简单的数学题 题解

求：

\[\left(\sum_{i=1}^n\sum_{j=1}^n ij \gcd(i,j)\right) \bmod p \]

对于100%的数据，\(n \leq 10^{10}\)，\(5 \times 10^8 \leq p \leq 1.1 \times 10^9\) 且 \(p\) 为质数。

开始推式子：

\[\sum_{i=1}^n\sum_{j=1}^n ij \gcd(i,j)\\ =\sum_{g=1}^ng\sum_{i=1}^n\sum_{j=1}^n ij \left[\gcd(i,j)=g\right]\\ =\sum_{g=1}^ng^3\sum_{i=1}^{\lfloor\frac{n}{g}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{g}\rfloor} ij \left[\gcd(i,j)=1\right]\\ =\sum_{g=1}^ng^3\sum_{i=1}^{\lfloor\frac{n}{g}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{g}\rfloor} ij \sum_{d|g}\mu(d)\\ =\sum_{g=1}^ng^3\sum_{d=1}^{\lfloor\frac{n}{g}\rfloor}\mu(d)\sum_{i=1}^{\lfloor\frac{n}{g}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{g}\rfloor} ij \left[d|g\right]\\ =\sum_{g=1}^ng^3\sum_{d=1}^{\lfloor\frac{n}{g}\rfloor}\mu(d)\sum_{i=1}^{\lfloor\frac{n}{g}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{g}\rfloor} ij \left[d|i\right] \left[d|j\right]\\ =\sum_{g=1}^ng^3\sum_{d=1}^{\lfloor\frac{n}{g}\rfloor}\mu(d)d^2\sum_{i=1}^{\lfloor\frac{n}{gd}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{gd}\rfloor} ij\\ =\sum_{g=1}^ng^3\sum_{d=1}^{\lfloor\frac{n}{g}\rfloor}\mu(d)d^2\sum_{i=1}^{\lfloor\frac{n}{gd}\rfloor}\left(i\sum_{j=1}^{\lfloor\frac{n}{gd}\rfloor} j\right)\\ =\sum_{g=1}^ng^3\sum_{d=1}^{\lfloor\frac{n}{g}\rfloor}\mu(d)d^2\left(\sum_{i=1}^{\lfloor\frac{n}{gd}\rfloor}i\right)\left(\sum_{j=1}^{\lfloor\frac{n}{gd}\rfloor} j\right)\\ \]

记：

\[\operatorname{sum}(n)=\sum_{i=1}^ni=\dfrac{n*(n+1)}{2} \]

则：

\[Ans =\sum_{g=1}^ng^3\sum_{d=1}^{\lfloor\frac{n}{g}\rfloor}\mu(d)d^2 \operatorname{sum}^2(\lfloor\frac{n}{gd}\rfloor) \]

再记 \(T=gd\)：

\[Ans =\sum_{g=1}^ng^3\sum_{d=1}^{\lfloor\frac{n}{g}\rfloor}\mu(d)d^2 \operatorname{sum}^2(\lfloor\frac{n}{T}\rfloor)\\ =\sum_{T=1}^n\operatorname{sum}^2(\lfloor\frac{n}{T}\rfloor)\sum_{g|T}g^3\mu(\frac{T}{g}) \left(\dfrac{T}{g}\right)^2\\ =\sum_{T=1}^n\operatorname{sum}^2(\lfloor\frac{n}{T}\rfloor)\sum_{g|T}gT^2\mu(\frac{T}{g}) \\ =\sum_{T=1}^nT^2\operatorname{sum}^2(\lfloor\frac{n}{T}\rfloor)\sum_{g|T}g\mu(\frac{T}{g}) \\ \]

由狄雷克雷卷积可得：

\[\sum_{d|n}d\mu(\frac{n}{d})=\left(\operatorname{id}*\mu\right)(n)=\phi(n) \]

所以：

\[Ans =\sum_{T=1}^nT^2\operatorname{sum}^2(\lfloor\frac{n}{T}\rfloor)\phi(T)\\ =\sum_{T=1}^nT^2\phi(T)\operatorname{sum}^2(\lfloor\frac{n}{T}\rfloor)\\ \]

\(\operatorname{sum}\) 函数是 \(O(1)\) 的，所以后半部分可以进行数论分块。所以考虑处理前半部分 \(n^2\phi(n)\) 的前缀和，即求:

\[\operatorname{S}(n)=\sum_{i=1}^ni^2\phi(i) \]

因为 \(n \leq 10^{10}\)，所以必须杜教筛。杜教筛式子：

\[\operatorname{g}(1)\operatorname{S}(n)=\sum_{i=1}^n(\operatorname{g}*\operatorname{S})(i)- \sum_{i=2}^n\operatorname{g}(i)\operatorname{S}(\lfloor\frac{n}{i}\rfloor) \]

那么就需要寻找一个函数 \(\operatorname{g}(n)\)，保证我们能快速求出 \(\operatorname{g}(n)\) 和 \((\operatorname{g}*\operatorname{S})(n)\) 的前缀和。
显然可以取 \(\operatorname{g}(n)=\operatorname{id}^2(n)=n^2\)，所以：

\[(\operatorname{g}*\operatorname{S})(n) =\sum_{d|n}\operatorname{g}(d)\operatorname{S}(\frac{n}{d})\\ =\sum_{d|n}d^2(\frac{n}{d})^2\phi(\frac{n}{d})\\ =\sum_{d|n}n^2\phi(\frac{n}{d})\\ =n^2\sum_{d|n}\phi(\frac{n}{d})\\ \]

因为 \(d|n\) 是在枚举 \(n\) 的因数，而 \(\dfrac{n}{d}\) 实际上也是 \(n\) 的因数，因此：

\[(\operatorname{g}*\operatorname{S})(n)=n^2\sum_{d|n}\phi(d)\\ \]

又因为狄利克雷卷积：

\[\phi*\operatorname{I}=\operatorname{id} \]

所以有：

\[(\operatorname{g}*\operatorname{S})(n)=n^3\\ \]

梳理一下，为了杜教筛，现在我们还需要 \(O(1)\) 解决这两个式子的计算：

\[\sum_{i=1}^n\operatorname{g}(n)=\sum_{i=1}^ni^2\\ \sum_{i=1}^n(\operatorname{g}*\operatorname{S})(i)=\sum_{i=1}^ni^3 \]

先考虑平方和：

下列式子均为展开左侧、合并同类项所得。

\[(n+1)^3-n^3=3n^2+3n+1\\ n^3-(n-1)^3=3(n-1)^2+3(n-1)+1\\ (n-1)^3-(n-2)^3=3(n-2)^2+3(n-2)+1\\ (n-2)^3-(n-3)^3=3(n-3)^2+3(n-3)+1\\ ......\\ 3^3-2^3=3*2^2+3*2+1\\ 2^3-1^3=3*1^2+3*1+1 \]

把这 \(n\) 个式子加起来，再把对应次数的项合并。
左侧：

\[(n+1)^3-n^3+n^3-(n-1)^3+(n-1)^3-......-2^3+2^3-1^3\\ =(n+1)^3-1 \]

右侧：

\[3(n^2+(n-1)^2+......+1^2)+ 3(n+(n-1)+......+1)+ (1+1+......+1)\\ =3\sum_{i=1}^ni^2+3\sum_{i=1}^ni+n\\ =3\sum_{i=1}^ni^2+3\dfrac{n(n+1)}{2}+n\\ \]

所以整体是：

\[(n+1)^3-1=3\sum_{i=1}^ni^2+3\dfrac{n(n+1)}{2}+n\\ 3\sum_{i=1}^ni^2=(n+1)^3-1-3\dfrac{n(n+1)}{2}-n\\ 3\sum_{i=1}^ni^2=n^3+3n^2+3n+1-1-\dfrac{3}{2}n^2-\dfrac{3}{2}n-n\\ 3\sum_{i=1}^ni^2=n^3+\dfrac{3}{2}n^2+\dfrac{1}{2}n\\ 3\sum_{i=1}^ni^2=\dfrac{2n^3+3n^2+n}{2}\\ \sum_{i=1}^ni^2=\dfrac{2n^3+3n^2+n}{6}\\ \sum_{i=1}^ni^2=\dfrac{n(2n^2+3n+1)}{6}\\ \sum_{i=1}^ni^2=\dfrac{n(n+1)(2n+1)}{6}\\ \]

再考虑立方和：

下列式子均为展开左侧、合并同类项所得。

\[(n+1)^4-n^4=4n^3+6n^2+4n+1\\ n^4-(n-1)^4=4n^3+6n^2+4n+1\\ ......\\ 2^4-1^4=4*1^3+6*1^2+4*1+1 \]

把这 \(n\) 个式子加起来，再把对应次数的项合并。
左侧：

\[(n+1)^4-n^4+n^4-(n-1)^4+(n-1)^4-......-2^4+2^4-1^4\\ =(n+1)^4-1 \]

右侧：

\[4(n^3+(n-1)^3+......+1^3)+ 6(n^2+(n-1)^2+......+1^2)+ 4(n+(n-1)+......+1)+ (1+1+......+1)\\ =4\sum_{i=1}^ni^3+6\sum_{i=1}^ni^2+4\sum_{i=1}^ni+n\\ =4\sum_{i=1}^ni^3+6\dfrac{n(n+1)(2n+1)}{6}+4\dfrac{n(n+1)}{2}+n\\ =4\sum_{i=1}^ni^3+n(n+1)(2n+1)+2n(n+1)+n\\ \]

所以整体是：

\[(n+1)^4-1=4\sum_{i=1}^ni^3+n(n+1)(2n+1)+2n(n+1)+n\\ 4\sum_{i=1}^ni^3=(n+1)^4-1-n(n+1)(2n+1)-2n(n+1)-n\\ 4\sum_{i=1}^ni^3=n^4+4n^3+6n^2+4n+1-1-2n^3-3n^2-n-2n^2-2n-n\\ 4\sum_{i=1}^ni^3=n^4+2n^3+n^2\\ 4\sum_{i=1}^ni^3=n^2(n^2+2n+1)\\ \sum_{i=1}^ni^3=\frac{n^2(n^2+2n+1)}{4}\\ \sum_{i=1}^ni^3=\frac{n^2(n+1)^2}{4}\\ \sum_{i=1}^ni^3=\left(\frac{n(n+1)}{2}\right)^2\\ \sum_{i=1}^ni^3=\left(\sum_{i=1}^ni^2\right)^2\\ \]

再次梳理一下我们推好的所有式子。

\[Ans= \sum_{T=1}^nT^2\phi(T)\operatorname{sum}^2(\lfloor\frac{n}{T}\rfloor)\\ \sum_{i=1}^n\operatorname{g}(n)=\sum_{i=1}^ni^2=\dfrac{n(n+1)(2n+1)}{6}\\ \sum_{i=1}^n(\operatorname{g}*\operatorname{S})(i)=\sum_{i=1}^ni^3 =\left(\frac{n(n+1)}{2}\right)^2=\operatorname{sum}^2(n) \]

有 \(/2\) 和 \(/6\) 的操作，因为题目保证输入的 \(p\) 为质数，费马小定理求逆元即可。
前缀和取余后可能会出现 \(\operatorname{S}(r) < \operatorname{S}(l-1)\) 的情况，所以记得 \(+p\) 再取余。
因为我们全程假设 \(n\le m\) 因为 \(\gcd\) 是对称的，所以记得在开头进行 swap。
代码：

// Problem: P3768 简单的数学题
// Contest: Luogu
// URL: https://www.luogu.com.cn/problem/P3768
// Memory Limit: 250 MB
// Time Limit: 4000 ms
// 
// Powered by CP Editor (https://cpeditor.org)

#include <bits/stdc++.h>
#include <bits/extc++.h>
#define INF 0x7fffffff
#define MAXN 3981071
#define eps 1e-9
#define foru(a,b,c)	for(int a=b;a<=c;a++)
#define RT return 0;
#define db(x)	cout<<endl<<x<<endl;
#define LL long long
#define int LL
#define LXF int
#define RIN rin()
#define HH printf("\n")
using namespace std;
inline LXF rin(){
	LXF x=0,w=1;
	char ch=0;
	while(ch<'0'||ch>'9'){ 
	if(ch=='-') w=-1;
	ch=getchar();
	}
	while(ch>='0'&&ch<='9'){
	x=x*10+(ch-'0');
	ch=getchar();
	}
	return x*w;
}
int p,n;
int phi[MAXN];
vector<int> prim;
bitset<MAXN> vis;
__gnu_pbds::cc_hash_table<int,int> sphi;
int inv2,inv6;
int ksm(int a,int b){
	int ret=1;
	a%=p;
	while(b){
		if(b&1)	ret=(ret*a)%p;
		b>>=1,a=a*a%p;
	}
	return ret;
}
int sum(int n){
	n%=p;
	return n*(n+1)%p*inv2%p;
}
int sumf(int n){
	n%=p;
	return n*(n+1)%p*(2*n%p+1)%p*inv6%p;
}
void work(){
	vis[1]=1;
	phi[1]=1;
	for(int i=2;i<=MAXN-10;i++){
		if(!vis[i]){
			prim.emplace_back(i);
			phi[i]=(i-1)%p;
		}
		for(int j=0;j<prim.size()&&i*prim[j]<=MAXN-10;j++){
			vis[i*prim[j]]=1;
			if(i%prim[j]==0){
				phi[i*prim[j]]=phi[i]*prim[j]%p;
				break;
			}
			phi[i*prim[j]]=phi[i]*(prim[j]-1)%p;
		}
		phi[i]=(phi[i]*i%p*i%p+phi[i-1])%p;
	}
}
int getphi(int n){
	if(n<=MAXN-50)	return phi[n];
	if(sphi.find(n)!=sphi.end())	return sphi[n];
	int ret=sum(n)*sum(n)%p,d=0;
	for(int l=2,r;l<=n;l=r+1){
		r=min(n,n/(n/l));
		ret-=getphi(n/l)*(sumf(r)-sumf(l-1)+p)%p;
		ret=(ret%p+p)%p;
	}
	return sphi[n]=ret;
}
signed main(){
	p=RIN,n=RIN;
	inv2=ksm(2,p-2);
	inv6=ksm(6,p-2);
	work();
	int ans=0;
	for(int l=1,r;l<=n;l=r+1){
		r=min(n,n/(n/l));
		ans+=sum(n/l)*sum(n/l)%p*(getphi(r)-getphi(l-1)+p)%p;
		ans%=p;
	}
	cout<<ans;
	return 0;
}