简单的数学题

I.简单的数学题

在做这题之前，我们先来见一位老朋友：

$\sum\limits_{i=1}^n\sum\limits_{j=1}^n\gcd(i,j)$

我们在1.v.[NOI2010]能量采集中就已经接触到了这道题。当时我们运用了$\sum\limits_{d|n}\mu(d)=[n=1]$的公式。现在，我们要运用另一个公式$\sum\limits_{d|n}\varphi(d)=n$。

$\begin{aligned}&\sum_{i=1}^n\sum_{j=1}^n\gcd(i,j)\\=&\sum_{i=1}^n\sum_{j=1}^n\sum_{x|\gcd(i,j)}\varphi(x)\\=&\sum_{i=1}^n\sum_{j=1}^n\sum_{x|i,x|j}\varphi(x)\\=&\sum_{x=1}^n\varphi(x)\left\lfloor\dfrac{n}{x}\right\rfloor^2\end{aligned}$

前面的东西可以线性筛前缀和，后面的东西可以整除分块，这样我们也得到一种和之前运用的$\mu$复杂度一致的运算。

而现在这道题呢？

我们要求$\sum\limits_{i=1}^n\sum\limits_{j=1}^nij\gcd(i,j)$。

再来！

$\begin{aligned}&\sum_{i=1}^n\sum_{j=1}^nij\gcd(i,j)\\=&\sum_{i=1}^n\sum_{j=1}^nij\sum_{x|\gcd(i,j)}\varphi(x)\\=&\sum_{i=1}^n\sum_{j=1}^nij\sum_{x|i,x|j}\varphi(x)\\=&\sum_{x=1}^n\varphi(x)x^2\sum_{i=1}^{n/x}\sum_{j=1}^{n/x}ij\\=&\sum_{x=1}^n\varphi(x)x^2(\sum_{i=1}^{n/x}i)^2\end{aligned}$

后面的$(\sum_{i=1}^{n/x}i)^2$可以直接套等差数列求和公式并整除分块。前面的我们考虑杜教筛。

设$f(x)=\varphi(x)x^2$，我们要求$S(n)=\sum_{i=1}^n\varphi(i)i^2$。

我们想找一个合适的$g$。

考虑$(f*g)(n)=\sum\limits_{d|n}\varphi(d)d^2g(\dfrac{n}{d})$。

我们想把这个恶心的$d^2$消掉。这样的话，我们不如设$g=id^2$。这样的话，$(f*g)(n)=\sum\limits_{d|n}\varphi(d)d^2(\dfrac{n}{d})^2=\sum\limits_{d|n}\varphi(d)n^2=n^3$！！！

我们好像找对了。

再来看式子：

$g(1)S(n)=\sum_{d=1}^n(f*g)(d)-\sum_{d=2}^ng(d)S(\dfrac{n}{d})$

代入，得$S(n)=\sum_{i=1}^ni^3-\sum_{d=2}^ni^2S(\dfrac{n}{d})$

我们有公式：

$\sum\limits_{i=1}^ni^2=\dfrac{n(n+1)(2n+1)}{6}$

$\sum\limits_{i=1}^ni^3=\dfrac{n^2(n+1)^2}{4}$

这样我们就可以杜教筛了。

我们回到一开始的式子：$\sum_{x=1}^n\varphi(x)x^2(\sum_{i=1}^{n/x}i)^2$。后半部分是$O(\sqrt{n})$的整除分块，前半部分要用杜教筛，这样，看起来复杂度为$O(n^{\frac{2}{3}}*\sqrt{n})=O(n^{\frac{7}{6}})$。

Emm?这东西一看就不像能过的样子呀？

首先，因为杜教筛是基于记忆化搜索实现的，多次使用复杂度明显是低于$O(n^{\frac{2}{3}})$的。至于低多少，就要看出题人是否是用脚造数据了。

同时，因为整除分块的过程中，杜教筛要筛的$n$不会总是$n$级别的。总复杂度应该低于$O(\sum_{i=1}^n(\dfrac{n}{i})^{\frac{2}{3}})$

并且，因为杜教筛是基于分块思想的，我们一开始预处理的部分也可以大于$n^{\frac{2}{3}}$，这样单次杜教筛复杂度就会低于$O(n^{\frac{2}{3}})$。这样一番优化下来，这个算法复杂度应该是亚线性，甚至直接是$O(n^{\frac{2}{3}})$的（其实是我不会证）

或者换一种说法，它的复杂度是$O(\text{能过})$。

代码：

#pragma GCC optimize(3)
#include<bits/stdc++.h>
using namespace std;
#define int long long
const int N=4000000;
int p,phi[N+10],inv6,pri[N+10],n;
inline int ksm(int x,int y){
	register int res=1;
	for(;y;x=(1ll*x*x)%p,y>>=1)if(y&1)res=(1ll*res*x)%p;
	return res;
}
inline void init(){
	inv6=ksm(6,p-2);
	phi[1]=1;
	for(register int i=2;i<=N;i++){
		if(!pri[i])pri[++pri[0]]=i,phi[i]=i-1;
		for(register int j=1;j<=pri[0]&&i*pri[j]<=N;j++){
			pri[i*pri[j]]=true;
			if(!(i%pri[j])){phi[i*pri[j]]=(1ll*phi[i]*pri[j])%p;break;}
			phi[i*pri[j]]=phi[i]*phi[pri[j]];
		}
	}
	for(register int i=1;i<=N;i++)phi[i]=(phi[i-1]+1ll*phi[i]*i%p*i%p)%p;
}
inline int sqrsum(int x){
	x%=p;
	return 1ll*x*(x+1)%p*(2*x+1)%p*inv6%p;
}
inline int cubsum(int x){
	x%=p;
	return (1ll*((1ll*x*(x+1)/2)%p)*((1ll*x*(x+1)/2)%p))%p;
}
unordered_map<int,int>mp;
inline int djs(int x){
	if(x<=N)return phi[x];
	if(mp[x])return mp[x];
	register int res=cubsum(x);
	for(register int l=2,r;l<=x;l=r+1)r=x/(x/l),res=(res-(1ll*(sqrsum(r)-sqrsum(l-1)+p)%p*djs(x/l)%p)+p)%p;
	return mp[x]=res;
}
inline int solve(){
	register int res=0;
	for(register int l=1,r;l<=n;l=r+1)r=n/(n/l),res=(1ll*(djs(r)-djs(l-1)+p)%p*cubsum(n/l)%p+res)%p;
	return res;
}
signed main(){
	scanf("%lld%lld",&p,&n),init();
	printf("%lld\n",solve());
	return 0;
}