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

本文链接：https://blog.youkuaiyun.com/sugar_free_mint/article/details/99697554

本文探讨了一个复杂的数学问题，即求解特定形式的双重求和表达式，涉及数论函数如欧拉函数φ(k)。通过数学转换和优化，采用杜教筛算法实现高效计算，最终给出了解决方案的C++代码实现。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

题目

求 $∑i=1n∑j=1nijgcd(i,j)(modp)\sum_{i=1}^n\sum_{j=1}^nijgcd(i,j)\pmod p$

分析

一波推式子
$=∑i=1n∑j=1nij∑k∣i,k∣jφ(k)=\sum_{i=1}^n\sum_{j=1}^nij\sum_{k|i,k|j}\varphi(k)$
$=∑k=1nφ(k)k2∑i=1⌊nk⌋∑j=1⌊nk⌋ij=\sum_{k=1}^n\varphi(k)k^2\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{k}\rfloor}ij$
$=∑k=1nφ(k)k2(∑i=1⌊nk⌋i)2=\sum_{k=1}^n\varphi(k)k^2(\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}i)^2$
$=∑k=1nφ(k)k2∑i=1⌊nk⌋i3=\sum_{k=1}^n\varphi(k)k^2\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}i^3$
把右边视为 $g(⌊nk⌋)g(\lfloor\frac{n}{k}\rfloor)$ ，左边视为 $f (k)$
那么其实这个式子可以变成 $∑k=1nf(k)g(⌊nk⌋)\sum_{k=1}^nf(k)g(\lfloor\frac{n}{k}\rfloor)$
后面这个东西整除分块，前面套杜教筛，时间复杂度 $O(n23)O(n^{\frac{2}{3}})$
比如说我们要求 $f$ ,用 $h = (f * g)$ 去配对，所以要找到一个合适的 $g$ 去卷积,后面的 $g$ 是新的 $g$
那么 $h(x)=∑i∣xf(i)∗g(x/i)h(x)=\sum_{i|x}f(i)*g(x/i)$
那 $g (x)$ 可以用 $x^2$ 代替
那么 $h(x)=∑i∣xφ(i)x2=x3h(x)=\sum_{i|x}\varphi(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));
}