#莫比乌斯反演，杜教筛#洛谷 3768 简单的数学题

题目

求$$\sum _{i=1}^n\sum _{j=1}^{n}ijgcd(i,j)(\mathrm{m}\mathrm{o}\mathrm{d}p)$$

---

分析

一波推式子
$$=\sum _{i=1}^n\sum _{j=1}^{n}ij\sum _{k|i,k|j}\phi (k)$$
$$=\sum _{k=1}^n\phi (k)k^2\sum _{i=1}^{\lfloor \frac{n}{k}\rfloor }\sum _{j=1}^{\lfloor \frac{n}{k}\rfloor }ij$$
$$=\sum _{k=1}^n\phi (k)k^2(\sum _{i=1}^{\lfloor \frac{n}{k}\rfloor }i{)}^2$$
$$=\sum _{k=1}^n\phi (k)k^2\sum _{i=1}^{\lfloor \frac{n}{k}\rfloor }{i}^3$$
把右边视为$g(\lfloor \frac{n}{k}\rfloor )$，左边视为$f(k)$
那么其实这个式子可以变成$$\sum _{k=1}^{n}f(k)g(\lfloor \frac{n}{k}\rfloor )$$
后面这个东西整除分块，前面套杜教筛，时间复杂度$O(n^{\frac{2}{3}})$
比如说我们要求$f$,用$h=(f\ast g)$去配对，所以要找到一个合适的$g$去卷积,后面的$g$是新的$g$
那么$h(x)={\sum }_{i|x}f(i)\ast g(x/i)$
那$g(x)$可以用$x^2$代替
那么$h(x)={\sum }_{i|x}\phi (i)x^2=x^3$
那就可以按照杜教筛的方式解决

---

代码

#include <cstdio>
#include <map>
#define rr register
using namespace std;
typedef long long ll; map<ll,ll>uk;
const int N=4700000; int v[N|1],cnt,prime[330001]; ll n,mod,inv,phi[N|1],ans;
inline ll mo(ll x,ll y){return x+y>=mod?x+y-mod:x+y;} 
inline ll g(ll n){n%=mod; return n*(n+1)%mod*(n<<1|1)%mod*inv%mod;}
inline ll sqr(ll n){return n*n%mod;}
inline ll cub(ll n){n%=mod; return sqr((n*(n+1)>>1)%mod);}
inline void init(int n){
    phi[1]=1;
    for (rr int i=2;i<=n;++i){
        if (!v[i]) phi[i]=(prime[++cnt]=v[i]=i)-1;
        for (rr int j=1;j<=cnt&&prime[j]*i<=n;++j){
            v[i*prime[j]]=prime[j];
            if (i%prime[j]) phi[i*prime[j]]=phi[i]*phi[prime[j]];
            else {phi[i*prime[j]]=phi[i]*prime[j]; break;}
        }
    }
    for (rr int i=2;i<=n;++i) phi[i]=mo(phi[i-1],phi[i]*sqr(i)%mod);
}
inline ll f(ll n){
    if (n<=N) return phi[n];
    if (uk.find(n)!=uk.end()) return uk[n];
    rr ll ans=cub(n);
    for (rr ll l=2,r;l<=n;l=r+1)
        r=n/(n/l),ans=mo(ans-f(n/l)*(mo(g(r)-g(l-1),mod))%mod,mod);
    return uk[n]=mo(ans,mod);
}
inline ll ksm(ll x,ll y){
    rr ll ans=1;
    for (;y;y>>=1,x=x*x%mod)
        if (y&1) ans=ans*x%mod;
    return ans;
}
signed main(){
    scanf("%lld%lld",&mod,&n);
    init(n<N?n:N); inv=ksm(6,mod-2);
    for (rr ll l=1,r;l<=n;l=r+1)
        r=n/(n/l),ans=mo(ans,cub(n/l)*(mo(f(r)-f(l-1),mod))%mod);
    return !printf("%lld",mo(ans,mod));
}