题目

题目链接：https://www.luogu.com.cn/problem/P3768
输入一个整数 \(n\) 和一个整数 \(p\)，你需要求出：

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

思路

大力化式子

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

\[=\sum^{n}_{d=1}d\sum^{n}_{i=1}\sum^{n}_{j=1}ij[\gcd(i,j)=d] \]

\[=\sum^{n}_{d=1}d^3\sum^{\lfloor\frac{n}{d}\rfloor}_{i=1}\sum^{\lfloor\frac{n}{d}\rfloor}_{j=1}ij[\gcd(i,j)=1] \]

\[=\sum^{n}_{d=1}d^3\sum^{\lfloor\frac{n}{d}\rfloor}_{i=1}\sum^{\lfloor\frac{n}{d}\rfloor}_{j=1}ij\sum^{}_{k|\gcd(i,j)}\mu(k) \]

显然把 \(\mu\) 扔到前面去

\[=\sum^{n}_{d=1}d^3\sum^{n}_{k=1}\mu(k)k^2\sum^{\lfloor\frac{n}{dk}\rfloor}_{i=1}\sum^{\lfloor\frac{n}{dk}\rfloor}_{j=1}ij \]

设 \(g(x)=\sum^{n}_{i=1}\sum^{n}_{j=1}ij\)，则 \(g(x)=(\frac{n(n+1)}{2})^2\)。
那么

\[=\sum^{n}_{d=1}d^3\sum^{n}_{k=1}\mu(k)k^2\times g(\lfloor\frac{n}{dk}\rfloor) \]

设 \(i=dk\)，那么 \(d^3k^2=i^2d\)，

\[=\sum^{n}_{i=1}g(\lfloor\frac{n}{i}\rfloor)i^2\sum^{}_{d|i}\mu(\frac{i}{d})d \]

惊奇的发现最后是一个狄利克雷卷积的形式，\(\sum^{}_{d|i}\mu(\frac{i}{d})d=\varphi(i)\)。

\[=\sum^{n}_{i=1}g(\lfloor\frac{n}{i}\rfloor)i^2\varphi(i) \]

前面的 \(g\) 可以整除分块，这样我们就需要求出 \(i^2\varphi(i)\) 的前缀和。
发现这是一个积性函数，设 \(f(n)=n^2\varphi(n)\)，我们需要找到一个合适的积性函数 \(g\)，然后杜教筛求 \((f*g)\) 的前缀和。
可以发现，取 \(g(n)=n^2\)，那么 \((f*g)(n)=\sum_{d|n}f(d)g(\frac{n}{d})=\sum_{d|n}d^2\varphi(d)\times \frac{n^2}{d^2}=n^2\sum_{d|n}\varphi(d)=n^3\)。
通过 OEIS 发现 \(\sum^{n}_{i=1}i^3=(\sum^{n}_{i=1}i)^2\)，是可以 \(O(1)\) 计算的。
而通过 OEIS 可以再次知道 \(g\) 的前缀和 \(\sum^{n}_{i=1}i^2=\frac{n(n+1)(2n+1)}{6}\)。也是可以 \(O(1)\) 计算的。
那么我们就可以杜教筛求出 \(f\) 的前缀和了。
听 my_dog 大爷说，时间复杂度的证明涉及到微积分。所以复杂度就是 \(O(\text{能过})\) /fad。

代码 
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;

const int N=2e7;
ll n,m,inv,MOD,ans,phi[N],prm[N];
bool v[N];
map<ll,int> sum;

void findprm(int n)
{
	phi[1]=1;
	for (int i=2;i<=n;i++)
	{
		if (!v[i]) prm[++m]=i,phi[i]=i-1;
		for (int j=1;j<=m;j++)
		{
			if (i>n/prm[j]) break;
			v[i*prm[j]]=1; phi[i*prm[j]]=phi[i]*(prm[j]-1);
			if (i%prm[j]==0)
			{
				phi[i*prm[j]]=phi[i]*prm[j];
				break;
			}
		}
	}
}

ll fpow(ll x,ll k)
{
	ll ans=1;
	for (;k;k>>=1,x=x*x%MOD)
		if (k&1) ans=ans*x%MOD;
	return ans;
}

ll calc2(ll x)
{
	ll p=x%MOD;
	return p*(p+1)%MOD*(2*p+1)%MOD*inv%MOD;
}

ll calc3(ll x)
{
	ll p=x%MOD;
	return ((1+p)*p/2%MOD)*((1+p)*p/2%MOD)%MOD;
}

ll getsum(ll n)
{
	if (n<N) return phi[n];
	if (sum[n]) return sum[n];
	ll res=calc3(n%MOD);
	for (ll l=2,r;l<=n;l=r+1)
	{
		r=n/(n/l);
		res=(res-getsum(n/l)*(calc2(r)-calc2(l-1)))%MOD;
	}
	return sum[n]=res;
}

int main()
{
	scanf("%lld%lld",&MOD,&n);
	findprm(N-1);
	for (int i=1;i<N;i++)
		phi[i]=(phi[i]*i%MOD*i+phi[i-1])%MOD;
	inv=fpow(6,MOD-2);
	for (ll l=1,r;l<=n;l=r+1)
	{
		r=n/(n/l);
		ans=(ans+calc3(n/l)*(getsum(r)-getsum(l-1))%MOD)%MOD;
	}
	printf("%lld",(ans%MOD+MOD)%MOD);
	return 0;
}