Description
给一个整数(n)，让你求
[sumlimits_{i=1}^n sumlimits_{j=1}^n ijgcd(i,j)
]
对一个大质数(p)取模。
保证(n le 10^{10},5 imes 10^{8} le p le 1.1 imes 10^9)，(p)为质数
Solutions
先来推柿子好了，枚举(gcd)的取值，有
[egin{aligned}
Ans&=sumlimits_{k} ksumlimits_{i=1}^nsumlimits_{j=1}^n ij[gcd(i,j)=k] \
&=sumlimits_{k} k^3 sumlimits_{i=1}^{lfloorfrac{n}{k}
floor}sumlimits_{j=1}^{lfloorfrac{n}{k}
floor} ij[gcd(i,j)=1]
end{aligned}
]
考虑求(Sum(n)=sumlimits_{i=1}^nsumlimits_{j=1}^n ij[gcd(i,j)=1])
推柿子
[egin{aligned}
Sum(n)&=sumlimits_{i=1}^nsumlimits_{j=1}^n ijsumlimits_{dmid gcd(i,j)} mu(d) \
&= sumlimits_{d} mu(d) sumlimits_{i=1}^{lfloor n/d
floor}idsumlimits_{j=1}^{lfloor n/d
floor} jd \
&= sumlimits_{d} mu(d)d^2 s(lfloorfrac{n}{d}
floor)^2
end{aligned}
]
其中(s(n)=1+2+cdots+n=frac{n(n+1)}{2})。
所以
[egin{aligned}
Ans&=sumlimits_{k} k^3 Sum(lfloorfrac{n}{k}
floor) \
&=sumlimits_{k} k^3 sumlimits_{d} mu(d)d^2 s(lfloorfrac{lfloor n/k
floor}{d}
floor) \
&=sumlimits_{k} k^3 sumlimits_{d} mu(d)d^2 s(lfloorfrac{n}{kd}
floor)
end{aligned}
]
枚举(T=kd)，有
[egin{aligned}
Ans&=sumlimits_{T} s(lfloorfrac{n}{T}
floor) sumlimits_{dmid T} mu(d)d^2(frac{T}{d})^3 \
&=sumlimits_{T} s(lfloorfrac{n}{T}
floor)T^2sumlimits_{dmid T}mu(d) imes frac{T}{d}
end{aligned}
]
令(f(n)=n^2sumlimits_{dmid n} mu(d) imes frac{n}{d})，如果能快速求出(f)的前缀和的话我们对上面的柿子数论分块就好了。
观察到后面的和式是(mu)与(id)的Dirichlet卷积的形式，假设
[F(n)=sumlimits_{dmid n} mu(d) imes frac{n}{d}
]
根据莫比乌斯反演的结论，必有
[n=sumlimits_{dmid n} F(d)
]
可以得到(F(n)=varphi(n))，所以(f(n)=n^2varphi(n))，我们想快速求出(f)的前缀和，(nle 10^{10})，线筛又死了。
可以考虑杜教筛(djs?)，令(S(n)=sumlimits_{i=1}^n f(i))，我们想找到另一个积性函数(g)，让(f*g)好看一点，我们知道欧拉函数有一个很美妙的性质(sumlimits_{dmid n} varphi(n)=n)，所以为了把(f)中的(n^2)消掉，配(g(n)=n^2)即可，有
[h(n)=(f*g)(n)=sumlimits_{dmid n} f(d)g(frac{n}{d})=sumlimits_{dmid n} n^2varphi(d)
]
即(h(n)=n^2sumlimits_{dmid n} varphi(n)=n^3)，
愉快的套柿子
[S(n)=sumlimits_{i=1}^n h(i)-sumlimits_{i=2}^n g(i)S(lfloorfrac{n}{i}
floor)=sumlimits_{i=1}^n i^3-sumlimits_{i=2}^n i^2S(lfloor frac{n}{i}
floor)
]
用一点小学奥数的知识，我们知道(1^3+2^3+cdots+n^3=(1+2+cdots+n)^2)，(1^2+2^3+cdots+n^2=frac{n(n+1)(2n+1)}{6})，所以
[S(n)=s(n)^2-sumlimits_{i=2}^n i^2S(lfloorfrac{n}{i}
floor)
]
后面的显然可以数论分块，于是处理(f)的前缀和的话杜教筛就好了。
再写一遍答案的柿子
[Ans=sumlimits_{T} s(lfloorfrac{n}{T}
floor)f(T)
]
对(T)数论分块+杜教筛就没了。
#include <bits/stdc++.h>
using namespace std;
#define rep(i,a,b) for (int i=(a);i<(b);++i)
#define per(i,a,b) for (int i=(a)-1;i>=(b);--i)
#define mp make_pair
#define pb push_back
#define fi first
#define se second
#define all(x) (x).begin(),(x).end()
#define SZ(x) ((int)(x).size())
typedef double db;
typedef long long ll;
typedef pair<int,int> PII;
typedef vector<int> VI;
const int maxn=5e6,N=maxn+10;
int mod,i2,i6;
inline int add(int x,int y) {return (x+=y)>=mod?x-mod:x;}
inline int sub(int x,int y) {return (x-=y)<0?x+mod:x;}
inline int normal(ll x) {return x<0?x+mod:x;}
inline int fpow(int x,int y) {
int ret=1; for(;y;y>>=1,x=1ll*x*x%mod)
if(y&1) ret=1ll*ret*x%mod;
return ret;
}
inline int sqr(int x) {return 1ll*x*x%mod;}
inline void initmod(int P) {
mod=P;
i2=fpow(2,mod-2),i6=fpow(6,mod-2);
}
inline int ss1(ll n) {n%=mod;return 1ll*n*(n+1)%mod*i2%mod;}
inline int ss2(ll n) {n%=mod;return 1ll*n*(n+1)%mod*(2*n+1)%mod*i6%mod;}
inline int s1(ll l,ll r) {return sub(ss1(r),ss1(l-1));}
inline int s2(ll l,ll r) {return sub(ss2(r),ss2(l-1));}
int p[N],pn,phi[N];
bool vis[N];
void sieve(int n) {
phi[1]=1;
rep(i,2,n+1) {
if(!vis[i]) {phi[i]=i-1;p[pn++]=i;}
for(int j=0;j<pn&&i*p[j]<=n;j++) {
vis[i*p[j]]=1;
if(i%p[j]==0) {phi[i*p[j]]=phi[i]*p[j];break;}
else phi[i*p[j]]=phi[i]*phi[p[j]];
}
}
rep(i,1,n+1) phi[i]=add(phi[i-1],1ll*i*i%mod*phi[i]%mod);
}
map<ll,int> fsum;
int Sf(ll n) {
if(n<=maxn) return phi[n];
if(fsum.count(n)) return fsum[n];
int ans=sqr(ss1(n));
for(ll l=2,r=0;l<=n;l=r+1) {
r=n/(n/l);
ans=sub(ans,1ll*s2(l,r)*Sf(n/l)%mod);
}
return fsum[n]=ans;
}
int main() {
#ifdef LOCAL
freopen("a.in","r",stdin);
#endif
int p; ll n; scanf("%d%lld",&p,&n);
initmod(p),sieve(maxn);
int ans=0;
for(ll l=1,r=0;l<=n;l=r+1) {
r=n/(n/l);
ans=add(ans,1ll*sub(Sf(r),Sf(l-1))*sqr(ss1(n/l))%mod);
}
printf("%d
",ans);
return 0;
}