简单的数学题
求\(\sum_{i=1}^n\sum_{j=1}^nijgcd(i,j)\ mod\ p,n\leq 10^{10},5×10^8≤p≤1.1×10^9\)且p为质数。
解
\[ans=\sum_{i=1}^n\sum_{j=1}^nijgcd(i,j)=\sum_{d=1}^nd\sum_{i=1}^n\sum_{j=1}^nij(gcd(i,j)==d)\]
于是设
\[f(d)=\sum_{i=1}^n\sum_{j=1}^nij(gcd(i,j)==d)\]
\[dc(x)=\sum_{i=1}^xi=\frac{(1+x)x}{2}\]
\[F(d)=\sum_{i=1}^n\sum_{j=1}^nij(d|gcd(i,j))=dc(n/d)^2d^2\]
由Mobius反演定理我们有
\[f(d)=\sum_{d|x}dc(n/x)^2x^2\mu(x/d)\]
代入原式，即
\[ans=\sum_{d=1}^nd\sum_{d|x}dc(n/x)^2x^2\mu(x/d)=\]
\[\sum_{x=1}^ndc(n/x)^2x^2\sum_{d|x}d\mu(x/d)\]
显然,后面的函数为积性函数，即\(\mu*id=\phi\),于是
\[\sum_{x=1}^ndc(n/x)^2x^2\varphi(x)\]
注意到\(x^2\)不能整除分块，于是考虑和\(\varphi(x)\)一起维护，所以设\(C(x)=x^2\varphi(x)\),接下来考虑用杜教筛求前缀和，显然只能待定函数法了，于是设\(A=B*C\)，所以
\[\sum_{i=1}^nA(i)=\sum_{i=1}^n\sum_{d|i}B(i/d)d^2\varphi(d)\]
为了抵掉分数形式，自然想到令\(B(x)=x^2\),于是我们有
\[\sum_{i=1}^n\sum_{d|i}B(i/d)d^2\varphi(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;
}