题意:求$(sum_{i=1}^{n}sum_{j=1}^{n}ijgcd(i,j))mod p$(p为质数,n<=1e10)
很显然,推式子。
$sum_{i=1}^{n}sum_{j=1}^{n}ijgcd(i,j)$
=$sum_{d=1}^{n}sum_{i=1}^{n}sum_{j=1}^{n}ijd[gcd(i,j)==d]$
=$sum_{d=1}^{n}d^3sum_{i=1}^{lfloor frac{n}{d} floor}sum_{j=1}^{lfloor frac{n}{d} floor}ij[gcd(i,j)==d]$
=$sum_{d=1}^{n}d^3sum_{i=1}^{lfloor frac{n}{d} floor}mu(i)i^2S({lfloor frac{n}{id} floor})^2,S(n)=(n+1)*n/2$
=$sum_{T=1}^{n}S({lfloor frac{n}{T} floor})^2sum_{d|T}d^3(frac{T}{d})^2mu(frac{T}{d})$
=$sum_{T=1}^{n}S({lfloor frac{n}{T} floor})^2T^2sum_{d|T}dmu(frac{T}{d})$
由$mu*id=varphi $可得$sum_{T=1}^{n}S({lfloor frac{n}{T} floor})^2T^2varphi (T)$
前面整除分块,只需要预处理$T^2varphi(T)$ 前缀和即可。
由于n有1e10那么大,就需要用到非线性的求前缀和的方法,这里用到杜教筛,见代码。
#include <bits/stdc++.h> #define ll long long using namespace std; const int N=1e6+5; int pri[N],tot,phi[N],sum[N]; bool p[N]; ll n,MD,ans,inv6; unordered_map<ll,int> w; void init() { phi[1]=1; for(int i=2;i<N;i++) { if(!p[i]) phi[i]=i-1,pri[tot++]=i; for(int j=0;j<tot&&pri[j]*i<N;j++) { p[i*pri[j]]=true; if(i%pri[j]==0) { phi[i*pri[j]]=phi[i]*pri[j]; break; } else phi[i*pri[j]]=phi[i]*phi[pri[j]]; } } for(int i=1;i<N;i++) sum[i]=(sum[i-1]+1LL*i*i%MD*phi[i]%MD)%MD; } ll pre_3(ll x) { x%=MD; ll t=x*(x+1)/2%MD; return t*t%MD; } ll pre_2(ll x) { x%=MD; return x*(x+1)%MD*(2*x+1)%MD*inv6%MD; } int quick_pow(int x,int y) { int ans=1; while(y) { if(y&1) ans=1LL*ans*x%MD; y>>=1; x=1LL*x*x%MD; } return ans; } int cal(ll x) { if(x<N) return sum[x]; if(w[x]) return w[x]; ll ans=pre_3(x); for(ll l=2,r;l<=x;l=r+1) { r=x/(x/l); ans=(ans-(pre_2(r)-pre_2(l-1)+MD)%MD*cal(x/l)%MD+MD)%MD; } return w[x]=ans; } int main() { scanf("%lld%lld",&MD,&n); inv6=quick_pow(6,MD-2),init(); for(ll l=1,r;l<=n;l=r+1) { r=n/(n/l); ans=(ans+pre_3(n/l)*(cal(r)-cal(l-1)+MD)%MD)%MD; } printf("%lld ",ans); return 0; }