P3768 简单的数学题（莫反+杜教筛）

· · 个人记录

题面

一点也不简单。

思路：

\sum^{n}_{i=1}\sum^{n}_{j=1}ijgcd(i,j) $=\sum^{n}_{d=1}d^3\sum^{\frac{n}{d}}_{i=1}\sum^{\frac{n}{d}}_{j=1}ij\sum^{}_{k|gcd(i,j)}\mu(k)$根据$\epsilon=\mu*1$反演。 $=\sum^{n}_{d=1}d^3\sum^{\frac{n}{d}}_{k=1}\mu(k)k^2\sum^{\frac{n}{dk}}_{i=1}\sum^{\frac{n}{dk}}_{j=1}ij$交换求和号$\times$2，考虑每个$\mu(k)$的贡献。 $=\sum^{n}_{d=1}d^3\sum^{\frac{n}{d}}_{k=1}\mu(k)k^2\frac{(\frac{n}{dk}(1+\frac{n}{dk}))^2}{4}

设sum(i)=\frac{(i*(i+1))^2}{4}，则原式为=\sum^{n}_{d=1}d^3\sum^{\frac{n}{d}}_{k=1}\mu(k)k^2sum(\frac{n}{dk})

$=\sum^{n}_{t=1}sum(\frac{n}{t})t^2\sum^{}_{d|t}d\mu(\frac{t}{d})$进一步化简。 $=\sum^{n}_{t=1}sum(\frac{n}{t})t^2\phi(t)$根据$\mu*id=\phi$反演。 考虑如何用杜教筛求积性函数$f(i)=i^2\phi(i)$的前缀和。 令$g(i)=i^2$，套杜教筛的式子，得到$s(n)=\sum^{n}_{i=1}i^3-\sum^{n}_{i=2}i^2s(\frac{n}{i})

时间复杂度O(n^{\frac{2}{3}})

代码：

#include<bits/stdc++.h>
#define int long long
using namespace std;
int n,p,phi[3000005],sum[3000005],prime[3000005],inv2,inv6,tot;
bool vis[3000005];
unordered_map<int,int> q;
int qpow(int x,int y,int mod){
    int ans=1;
    while(y){if(y&1) ans=ans*x%mod;x=x*x%mod,y>>=1;}
    return ans;
}
int sum1(int x){x%=p;return x*(x+1)%p*inv2%p;}
int sum2(int x){x%=p;return x*(x+1)%p*(2*x+1)%p*inv6%p;}
int sum3(int x){x%=p;return sum1(x)*sum1(x)%p;}
void euler(){
    phi[1]=1;
    for(int i=2;i<=3000000;i++){
        if(!vis[i]) vis[i]=1,phi[i]=i-1,prime[++tot]=i;
        for(int j=1;j<=tot&&i*prime[j]<=3000000;j++){
            vis[i*prime[j]]=1;
            if(i%prime[j]==0){phi[i*prime[j]]=phi[i]*prime[j];break;}
            phi[i*prime[j]]=phi[i]*(prime[j]-1);
        }
    }
    for(int i=1;i<=3000000;i++) sum[i]=(sum[i-1]+phi[i]*i%p*i%p)%p;
}
int calc(int x){
    if(x<=3000000) return sum[x];
    if(q.find(x)!=q.end()) return q[x];
    int s=sum3(x);
    for(int l=2,r,t;l<=x;l=r+1){
        t=x/l,r=x/t;
        s=((s-calc(t)*((sum2(r)-sum2(l-1))%p+p)%p)%p+p)%p;
    }
    return q[x]=s;
}
signed main(){
    cin>>p>>n;euler();inv2=qpow(2,p-2,p),inv6=qpow(6,p-2,p);int ans=0;
    for(int l=1,r,t;l<=n;l=r+1) t=n/l,r=n/t,ans=(ans+(sum3(t)*(((calc(r)-calc(l-1))%p+p)%p))%p)%p;
    cout<<ans<<endl;
    return 0;
}