传送门

题意：给定$p,N$,求

$$\sum _{i=1}^N\sum _{j=1}^{N}ijgcd(i,j)modp$$

$p$为质数，在$1e9$左右

$N\le 1e10$

神仙题

前置芝士：杜教筛

懒得重新写，所以顺便讲一下

杜教筛可以在非线性时间求积性函数的前缀和

设所求函数为$f(n)$,$S(n)为前缀和$

$$S(n)=\sum _{i=1}^{n}f(i)$$

构造一个可以快速计算前缀和的函数$g(n)$，且和$f(n)$的狄利克雷卷积$(f\ast g)(n)$的前缀和可以快速求出

考虑

$$\sum _{i=1}^n(f\ast g)(i)$$

暴力拆开

$$\sum _{i=1}^n\sum _{d\mid i}f(d)g(\frac{i}{d})$$

枚举$d$

$$\sum _{d=1}^{n}g(d)\sum _{d\mid i}f(\frac{i}{d})$$

$$\sum _{d=1}^{n}g(d)S(\lfloor \frac{n}{d}\rfloor )$$

我们要求的是$S(n)$,所以可以求出$g(1)S(n)$

即

$$\sum _{d=1}^{n}g(d)S(\lfloor \frac{n}{d}\rfloor )-\sum _{d=2}^{n}g(d)S(\lfloor \frac{n}{d}\rfloor )$$

换成上面

$$\sum _{i=1}^n(f\ast g)(i)-\sum _{d=2}^{n}g(d)S(\lfloor \frac{n}{d}\rfloor )$$

左边直接求，右边整除分块套递归

复杂度是$O(N^{\frac{3}{4}})$

如果线性筛出$N^{\frac{2}{3}}$以内的并用$map$记忆化，可以达到$O(N^{\frac{2}{3}})$

正文

$$\sum _{i=1}^N\sum _{j=1}^{N}ijgcd(i,j)$$

套路

$$\sum _{d=1}^N\sum _{i=1}^N\sum _{j=1}^N[gcd(i,j)=d]ijd$$

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

右边可以反演

设

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

$$F(d)=\sum _{d\mid n}f(n)=\sum _{i=1}^N\sum _{j=1}^N[d\mid i][d\mid j]ij$$

记$sum(n)=\frac{n(n+1)}{2}$

$$F(d)=d^2(sum\lfloor \frac{N}{d}\rfloor {)}^2$$

$$f(d)=\sum _{d\mid n}F(n)\mu (\frac{n}{d})=\sum _{d\mid n}{n}^2(sum\lfloor \frac{N}{n}\rfloor {)}^2\mu (\frac{n}{d})$$

原来求的是

$$\sum _{d=1}^{N}df(d)$$

$$\sum _{d=1}^{N}d\sum _{d\mid n}{n}^2(sum\lfloor \frac{N}{n}\rfloor {)}^2\mu (\frac{n}{d})$$

枚举$n$

$$\sum _{n=1}^N{n}^2(sum\lfloor \frac{N}{n}\rfloor {)}^2\sum _{d\mid n}d\mu (\frac{n}{d})$$

右边就是$\phi$

$$\sum _{n=1}^N(sum\lfloor \frac{N}{n}\rfloor {)}^2{n}^2\phi (n)$$

左边是个整除分块，所以求出右边就可以了

#undef f

设$f(n)=n^2\phi (n)$

显然这是个积性函数，考虑杜教筛

$$(f\ast g)(n)=\sum _{d\mid n}g(d)f(\frac{n}{d})=\sum _{d\mid n}g(d)\frac{n^2}{d^2}\phi (\frac{n}{d})$$

令$g(n)=n^2$就可以消掉

则

$$(f\ast g)(n)=\sum _{d\mid n}{n}^2\phi (\frac{n}{d})=n^2\sum _{d\mid n}\phi (d)=n^3$$

众所周知

$$1^2+2^2+3^2+...+n^2=\frac{(n)(n+1)(2n+1)}{6}$$

$$1^3+2^3+3^3+...+n^3=(1+2+3+...+n{)}^2$$

两个前缀和都可以$O(1)$算出

这个问题就解决了

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cctype>
#include <map>
#define MAXN 10000005
using namespace std;
typedef long long ll;
int MOD,inv2,inv6;
inline int add(const int& x,const int& y){return x+y>=MOD? x+y-MOD:x+y;}
inline int dec(const int& x,const int& y){return x-y<0? x-y+MOD:x-y;}
inline int mul(const int& x,const int& y){return (ll)x*y%MOD;}
inline int qpow(int a,int p)
{int ans=1;while (p){if (p&1) ans=mul(ans,a);a=mul(a,a);p>>=1;}return ans;
}
const int N=10000000;
int np[MAXN],pl[1000005],cnt;
int phi[MAXN],f_sum[MAXN];
void init()
{np[1]=phi[1]=1;for (int i=2;i<=N;++i){if (!np[i]) phi[pl[++cnt]=i]=i-1;int x;for (int j=1;(x=i*pl[j])<=N;++j){np[x]=1;if (i%pl[j]==0){phi[x]=mul(phi[i],pl[j]);break;}phi[x]=mul(phi[i],pl[j]-1);}}for (int i=1;i<=N;i++) f_sum[i]=add(f_sum[i-1],mul(phi[i],mul(i,i)));
}
inline int calc(const int& n){return mul(mul(n,n+1),inv2);}
inline int calc2(const int& n){return mul(mul(n,n+1),mul(2*n+1,inv6));}
inline int calc3(const int& n){int t=calc(n);return mul(t,t);}
map<int,int> m;
int getsum(ll n)
{if (n<=N) return f_sum[n];if (m[n]) return m[n];int ans=calc3(n%MOD);for (ll l=2,r;l<=n;l=r+1){r=n/(n/l);int tmp=getsum(n/l);ans=dec(ans,mul(dec(calc2(r%MOD),calc2((l-1)%MOD)),tmp));}return m[n]=ans;
}
//int gcd(int a,int b)
//{
//	return b? gcd(b,a%b):a;
//}
int main()
{ll n;scanf("%d%lld",&MOD,&n);inv2=(MOD+1)>>1;inv6=qpow(6,MOD-2);init();	int ans=0;
//	for (int i=1;i<=n;i++)
//		for (int j=1;j<=n;j++)
//			ans=(ans+(ll)i*j%MOD*gcd(i,j))%MOD;
//	cerr<<ans<<'\n';
//	ans=0;for (ll l=1,r;l<=n;l=r+1){r=n/(n/l);ans=add(ans,mul(calc3((n/l)%MOD),dec(getsum(r),getsum(l-1))));}printf("%d\n",ans);return 0;
}