求(sum_{i=1}^nsum_{j=1}^nijgcd(i,j) mod p,nleq 10^{10},5×10^8≤p≤1.1×10^9)且p为质数。
解
[ans=sum_{i=1}^nsum_{j=1}^nijgcd(i,j)=sum_{d=1}^ndsum_{i=1}^nsum_{j=1}^nij(gcd(i,j)==d)
]
于是设
[f(d)=sum_{i=1}^nsum_{j=1}^nij(gcd(i,j)==d)
]
[dc(x)=sum_{i=1}^xi=frac{(1+x)x}{2}
]
[F(d)=sum_{i=1}^nsum_{j=1}^nij(d|gcd(i,j))=dc(n/d)^2d^2
]
由Mobius反演定理我们有
[f(d)=sum_{d|x}dc(n/x)^2x^2mu(x/d)
]
代入原式，即
[ans=sum_{d=1}^ndsum_{d|x}dc(n/x)^2x^2mu(x/d)=
]
[sum_{x=1}^ndc(n/x)^2x^2sum_{d|x}dmu(x/d)
]
显然,后面的函数为积性函数，即(mu*id=phi),于是
[sum_{x=1}^ndc(n/x)^2x^2varphi(x)
]
注意到(x^2)不能整除分块，于是考虑和(varphi(x))一起维护，所以设(C(x)=x^2varphi(x)),接下来考虑用杜教筛求前缀和，显然只能待定函数法了，于是设(A=B*C)，所以
[sum_{i=1}^nA(i)=sum_{i=1}^nsum_{d|i}B(i/d)d^2varphi(d)
]
为了抵掉分数形式，自然想到令(B(x)=x^2),于是我们有
[sum_{i=1}^nsum_{d|i}B(i/d)d^2varphi(d)=sum_{i=1}^ni^3=(frac{(1+n)n}{2})^2
]
于是套杜教筛基本式，我们不难有,设(S(n)=sum_{i=1}^nC(i))
[S(n)=(frac{(1+n)n}{2})^2-sum_{d=2}^nS(n/d)d^2
]
(d^2)前缀和很好处理，进行杜教筛基本操作即可。
参考代码：
#include <iostream>
#include <cstdio>
#define il inline
#define ri register
#define ll long long
#define control 6000000
using namespace std;
template<class f1,class f2>
struct chain{
chain *next;f1 x;f2 y;
};
template<class f1,class f2,class f3,const f3 key>
struct unordered_map{
chain<f1,f2>*head[key],*pt;
il void insert(f1 x,f2 y){
pt=new chain<f1,f2>,pt->x=x,pt->y=y;
x%=key,pt->next=head[x],head[x]=pt;
}
il chain<f1,f2>* find(ll x){
for(pt=head[x%key];pt!=NULL;pt=pt->next)
if(x==pt->x)return pt;
return NULL;
}
};
bool check[control+1];
ll prime[500000],pt,phi[control+1],
yyb,inv6;
il void prepare();
il ll pow2(ll),djs(ll),dc(ll),
pow(ll,ll),p2sum(ll),p3sum(ll);
int main(){
ll n,i,j,ans(0);
scanf("%lld%lld",&yyb,&n),inv6=pow(6,yyb-2),
prepare();
for(i=1;i<=n;i=j+1)
j=n/(n/i),(ans+=pow2(dc(n/i))*(djs(j)-djs(i-1))%yyb)%=yyb;
printf("%lld",(ans+yyb)%yyb);
return 0;
}
il ll pow(ll x,ll y){
ll ans(1);while(y){
if(y&1)ans=ans*x%yyb;
x=x*x%yyb,y>>=1;
}return ans;
}
il ll p2sum(ll n){
return n%=yyb,n*(n+1)%yyb*(n*2+1)%yyb*inv6%yyb;
}
il ll p3sum(ll n){
return n%=yyb,pow2(dc(n));
}
il ll dc(ll n){
return n%=yyb,(n*(n+1)>>1)%yyb;
}
unordered_map<ll,ll,ll,999931>sxr;
il ll djs(ll n){
if(n<=control)return phi[n];
if(sxr.find(n)!=NULL)return sxr.find(n)->y;
ll i,j,ans(p3sum(n));
for(i=2;i<=n;i=j+1)
j=n/(n/i),(ans-=djs(n/j)*(p2sum(j)-p2sum(i-1))%yyb)%=yyb;
return sxr.insert(n,ans),ans;
}
il ll pow2(ll x){
return x*x%yyb;
}
il void prepare(){
ll i,j;check[1]|=phi[1]|=true;
for(i=2;i<=control;++i){
if(!check[i])prime[++pt]=i,phi[i]=i-1;
for(j=1;j<=pt&&prime[j]*i<=control;++j){
check[i*prime[j]]|=true;
if(!(i%prime[j])){phi[i*prime[j]]=phi[i]*prime[j];break;}
phi[i*prime[j]]=phi[i]*phi[prime[j]];
}
}for(i=1;i<=control;++i)
(((phi[i]*=pow2(i))%=yyb)+=phi[i-1])%=yyb;
}