简单的数学题

简单的数学题

求\(\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;
}

标签：pt,sum,简单,dc,il,yyb,数学题,ll
来源： https://www.cnblogs.com/a1b3c7d9/p/10801586.html