洛谷 3768簡單的數學題（莫比烏斯反演+杜教篩）

題目連接

\(Description\)
求\[\sum_{i-1}^n\sum_{j=1}^nijgcd(i,j)\mod p\]

\(n<=10^{10}\)

---

\(Solution\)
\[\sum_{i-1}^n\sum_{j=1}^nijgcd(i,j)\]
\[=\sum_{d=1}^nd^3\sum_{i=1}^{n/d}\sum_{j=1}^{n/d}ij[gcd(i,j)==1]\]
發現跟這道題只差了一個\(d^3\)。最後化簡獲得
\[ans=\sum_{T=1}^n(\frac{\frac{n}{T}\times(\frac{n}{T}+1)}{2})^2\phi(T)T^2\]
令\(f(i)=\phi(i)i^2\)。
已知\(\phi*1=Id\),因此能夠令\(g(i)=i^2\),這樣\(f*g(x)=\sum_{d|x}\phi(d)d^2\frac{x^2}{d^2}=x^3\),\(f*g\)的前綴和就能夠\(O(1)計算了\)。

#include<complex>
#include<cstdio>
#include<map>
#define LL long long
using namespace std;
const int N=7e6+7;
int tot,ans,mod,nn,div6;
LL n;
int prime[N],phi[N];
bool check[N];
map<LL,int>mp;
void Init()
{
    check[1]=phi[1]=1;
    nn=min(n,(LL)N-1);
    for(int i=2;i<=nn;i++)
    {
        if(!check[i])prime[++tot]=i,phi[i]=i-1;
        for(int j=1;j<=tot && i*prime[j]<=nn;j++)
        {
            check[i*prime[j]]=1;
            if(i%prime[j])phi[i*prime[j]]=phi[i]*phi[prime[j]];
            else
            {
                phi[i*prime[j]]=phi[i]*prime[j];
                break;
            }
        }
    }
    for(int i=1;i<=nn;i++)
        phi[i]=(phi[i-1]+1ll*i*i%mod*phi[i]%mod)%mod;
}
int calc1(LL x)
{
    x%=mod;
    return x*(x+1)/2%mod;
}
int calc2(LL x)
{
    x%=mod;
    return x*(x+1)%mod*(x+x+1)%mod*div6%mod;
}
LL Sum(LL x)
{
    if(x<=nn)return phi[x];
    if(mp[x])return mp[x];
    LL res=calc1(x);res=res*res%mod;
    for(LL l=2,r;l<=x;l=r+1)
    {
        r=x/(x/l);
        res=(res-(calc2(r)-calc2(l-1))*Sum(x/l)%mod)%mod;
    }
    return mp[x]=(res+mod)%mod;
}
int Fpow(LL b,int p)
{
    LL res=1;
    for(;p;p>>=1,b=b*b%mod)
        if(p&1)res=res*b%mod;
    return res;
}
int main()
{
    scanf("%d%lld",&mod,&n);
    div6=Fpow(6,mod-2);
    Init();
    for(LL l=1,r;l<=n;l=r+1)
    {
        r=n/(n/l);
        LL tmp=calc1(n/l);
        ans=(ans+tmp*tmp%mod*(Sum(r)-Sum(l-1))%mod)%mod;
    }
    printf("%d\n",(ans+mod)%mod);
    return 0;
}