习题：简单的数学题（杜教筛）

题目

传送门

思路

\[\begin{aligned}\sum_{i=1}^{n}\sum_{j=1}^{n}ijgcd(i,j)&=\sum_{i=1}^{n}\sum_{j=1}^{n}ij\sum_{d|gcd(i,j)}\varphi(d)\\&=\sum_{d=1}^{n}\varphi(d)(\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}id)^2\\&=\sum_{d=1}^{n}\varphi(d)d^2(\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}i)^2\end{aligned} \]

此时，整个就可以用一个整除分块来处理，现在主要考虑的为如何快速地求\(\sum_{d=1}^{n}\varphi(d)d^2\)

设\(f(n)=\varphi(n)n^2\)，显然\(f(n)\)是一个积性函数，尝试使用杜教筛

构造\(g(i)=i^2\)，那么\(h(n)=(f*g)(n)=\sum_{d|n}\varphi(d)d^2(\frac{n}{d})^2=n^3\)

设\(F(n)=\sum_{i=1}^{n}f(i),H(n)=\sum_{i=1}^{n}h(i)\)，根据结论，有

\(F(n)=\frac{H(n)-\sum_{d=2}^{n}g(d)F(\lfloor\frac{n}{d}\rfloor)}{g(1)}\)

同时，因为是在整除分块的时候进行杜教筛，所以最后总的时间复杂度为$O(n^{\frac{2}{3}}) $

代码

#include<unordered_map>
#include<iostream>
#include<cstdio>
#include<cmath>
using namespace std;
#define int long long
long long n,mod;
int len=30000000,inv6,inv2;
long long f[30000005];
int lenp,pri[5000005],phi[50000005];
bool vis[5000005];
unordered_map<long long,long long> mf;
long long qkpow(long long a,long long b)
{
    if(b==0)
        return 1;
    if(b==1)
        return a;
    long long t=qkpow(a,b/2);
    t=t*t%mod;
    if(b&1)
        t=t*a%mod;
    return t;
}
long long summ1(long long n)
{
    return n%mod*(n+1)%mod*inv2%mod;
}
long long summ2(long long n)
{
    return 1ll*n%mod*(n+1)%mod*((2*n+1)%mod)%mod*inv6%mod;
}
long long summ3(long long n)
{
    long long temp=(1ll*n%mod*(n+1)%mod*inv2%mod);
    return temp*temp%mod;
}
long long getphi(long long n)
{
    if(n<=len&&f[n])
        return f[n];
    if(mf.count(n))
        return mf[n];
    long long temp=summ3(n);
    for(long long l=2,r;l<=n;l=r+1)
    {
        r=n/(n/l);
        temp=(temp-(summ2(r)-summ2(l-1))%mod*getphi(n/l)%mod)%mod;
    }
    temp=(temp%mod+mod)%mod;
    if(n<=len)
        return f[n]=temp;
    return mf[n]=temp;
}
void prepa(long long n)
{
    len=n;
    f[1]=1;
    for(int i=2;i<=n;i++)
    {
        if(vis[i]==0)
        {
            pri[++lenp]=i;
            phi[i]=i-1;
        }
        for(int j=1;j<=lenp&&pri[j]*i<=n;j++)
        {
            vis[i*pri[j]]=1;
            if(i%pri[j]==0) 
            {
                phi[i*pri[j]]=pri[j]*phi[i];
                break;
            }
            phi[i*pri[j]]=phi[i]*phi[pri[j]];
        }
        f[i]=(f[i-1]+phi[i]*i%mod*i)%mod;
    }
}
signed main()
{
    ios::sync_with_stdio(false);
    cin>>mod>>n;
    prepa(pow(n,2.0/3));
    inv6=qkpow(6,mod-2);
    inv2=qkpow(2,mod-2);
    long long ans=0;
    for(int l=1,r;l<=n;l=r+1)
    {
        r=n/(n/l);
        //cout<<"debug:"<<l<<' '<<r<<'\n';
        ans=(ans+1ll*(getphi(r)-getphi(l-1))*summ1(n/l)%mod*summ1(n/l)%mod)%mod;
    }
    ans=(ans%mod+mod)%mod;
    cout<<ans;
    return 0;
}
/*
998244353 20
*/