洛谷 - P1829 - Crash的数字表格 - 莫比乌斯反演

求：
\(S(n,m)=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}lcm(i,j)\)

显然：
\(S(n,m)=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}\frac{ij}{gcd(i,j)}\)

枚举g：
\(S(n,m)=\sum\limits_{g=1}^{n}\frac{1}{g}\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}ij[gcd(i,j)==g]\)

除以g：
\(S(n,m)=\sum\limits_{g=1}^{n}g\sum\limits_{i=1}^{\lfloor\frac{n}{g}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{g}\rfloor}ij[gcd(i,j)==1]\)

记：
\(S_1(n,m)=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}ij[gcd(i,j)==1]\)

原式：
\(S(n,m)=\sum\limits_{g=1}^{n}gS_1(\lfloor\frac{n}{g}\rfloor,\lfloor\frac{m}{g}\rfloor)\)

化简\(S_1(n,m)\)，显然：
\(S_1(n,m)=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}ij\sum\limits_{k|gcd(i,j)}\mu(k)\)

枚举k：
\(S_1(n,m)=\sum\limits_{k=1}^{min}\mu(k)\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}ij[k|gcd(i,j)]\)

显然：
\(S_1(n,m)=\sum\limits_{k=1}^{min}\mu(k)\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}ij[k|i][k|j]\)

这种时候可以除以k：
\(S_1(n,m)=\sum\limits_{k=1}^{min}\mu(k)k^2\sum\limits_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{k}\rfloor}ij[1|i][1|j]\)

即：
\(S_1(n,m)=\sum\limits_{k=1}^{min}\mu(k)k^2\sum\limits_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{k}\rfloor}ij\)

记：
\(S_2(n,m)=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}ij\)

原式：
\(S_1(n,m)=\sum\limits_{k=1}^{min}\mu(k)k^2S_2(\lfloor\frac{n}{k}\rfloor,\lfloor\frac{m}{k}\rfloor)\)

显然：
\(S_2(n,m)=\sum\limits_{i=1}^{n}i\sum\limits_{j=1}^{m}j\)

即：
\(S_2(n,m)=\frac{1}{4}n(n+1)m(m+1)\)

时间复杂度：
求\(S_2(n,m)\)是\(O(1)\)，分块求\(S_1(n,m)\)是\(O(n^{\frac{1}{2}})\)（大概），分块求\(S(n,m)\)是\(O(n)\)（大概）。
还需要线性筛出：\(\sum\limits_{k=1}^{min}\mu(k)k^2\)

---

线性筛已经足够了，还好写，不过为什么不试一波杜教筛呢？

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;

const int mod=20101009;
const int MAXN=1e7;

int pri[MAXN+1];
int &pritop=pri[0];
int A[MAXN+1];
int pk[MAXN+1];

void sieve(int n=MAXN) {
    pk[1]=1;
    A[1]=1;
    for(int i=2; i<=n; i++) {
        if(!pri[i]) {
            pri[++pritop]=i;
            pk[i]=i;
            ll tmp=1ll*i*i;
            if(tmp>=mod)
                tmp%=mod;
            tmp=-tmp;
            if(tmp<mod)
                tmp+=mod;
            A[i]=tmp;
        }
        for(int j=1; j<=pritop; j++) {
            int &p=pri[j];
            int t=i*p;
            if(t>n)
                break;
            pri[t]=1;
            if(i%p) {
                pk[t]=p;
                ll tmp=1ll*A[i]*A[p];
                if(tmp>=mod)
                    tmp%=mod;
                A[t]=tmp;
            } else {
                pk[t]=pk[i]*p;
                if(pk[t]==t) {
                    A[t]=0;
                } else {
                    ll tmp=1ll*A[t/pk[t]]*A[pk[t]];
                    if(tmp>=mod)
                        tmp%=mod;
                    A[t]=tmp;
                }
                break;
            }
        }
    }
    for(int i=1; i<=n; i++) {
        A[i]+=A[i-1];
        if(A[i]>=mod)
            A[i]-=mod;
    }
}

inline int qpow(ll x,int n) {
    ll res=1;
    while(n) {
        if(n&1) {
            res*=x;
            if(res>=mod)
                res%=mod;
        }
        x*=x;
        if(x>=mod)
            x%=mod;
        n>>=1;
    }
    return res;
}

const int inv4=qpow(4,mod-2);

inline int S2(int n,int m) {
    ll res=1ll*n*(n+1);
    if(res>=mod)
        res%=mod;
    res*=m;
    if(res>=mod)
        res%=mod;
    res*=(m+1);
    if(res>=mod)
        res%=mod;
    res*=inv4;
    if(res>=mod)
        res%=mod;
    return res;
}

inline int S1(int n,int m) {
    ll res=0;
    int nm=min(n,m);
    for(int l=1,r; l<=nm; l=r+1) {
        int tn=n/l;
        int tm=m/l;
        r=min(n/tn,m/tm);
        ll tmp=A[r]-A[l-1];
        if(tmp<0)
            tmp+=mod;
        tmp*=S2(tn,tm);
        if(tmp>=mod)
            tmp%=mod;
        res+=tmp;
    }
    if(res>=mod)
        res%=mod;
    return res;
}

inline int s1(int l,int r) {
    ll res=(1ll*(l+r)*(r-l+1))>>1;
    if(res>=mod)
        res%=mod;
    return res;
}

inline int S(int n,int m) {
    ll res=0;
    int nm=min(n,m);
    for(int l=1,r; l<=nm; l=r+1) {
        int tn=n/l;
        int tm=m/l;
        r=min(n/tn,m/tm);
        ll tmp=s1(l,r);
        tmp*=S1(tn,tm);
        if(tmp>=mod)
            tmp%=mod;
        res+=tmp;
    }
    if(res>=mod)
        res%=mod;
    return res;
}

int main() {
#ifdef Yinku
    freopen("Yinku.in","r",stdin);
#endif // Yinku
    int n,m;
    scanf("%d%d",&n,&m);
    sieve(max(n,m));
    printf("%d\n",S(n,m));
    return 0;
}

---

杜教筛还是非常快的。

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;

const int mod=20101009;
//杜教筛
const int MAXN=pow(1e7,2.0/3.0)+1;
int pri[MAXN+1];
int &pritop=pri[0];
int A[MAXN+1];
int pk[MAXN+1];

void sieve(int n=MAXN) {
    pk[1]=1;
    A[1]=1;
    for(int i=2; i<=n; i++) {
        if(!pri[i]) {
            pri[++pritop]=i;
            pk[i]=i;
            ll tmp=1ll*i*i;
            if(tmp>=mod)
                tmp%=mod;
            tmp=-tmp;
            if(tmp<mod)
                tmp+=mod;
            A[i]=tmp;
        }
        for(int j=1; j<=pritop; j++) {
            int &p=pri[j];
            int t=i*p;
            if(t>n)
                break;
            pri[t]=1;
            if(i%p) {
                pk[t]=p;
                ll tmp=1ll*A[i]*A[p];
                if(tmp>=mod)
                    tmp%=mod;
                A[t]=tmp;
            } else {
                pk[t]=pk[i]*p;
                if(pk[t]==t) {
                    A[t]=0;
                } else {
                    ll tmp=1ll*A[t/pk[t]]*A[pk[t]];
                    if(tmp>=mod)
                        tmp%=mod;
                    A[t]=tmp;
                }
                break;
            }
        }
    }
    for(int i=1; i<=n; i++) {
        A[i]+=A[i-1];
        if(A[i]>=mod)
            A[i]-=mod;
    }
}

inline int qpow(ll x,int n) {
    ll res=1;
    while(n) {
        if(n&1) {
            res*=x;
            if(res>=mod)
                res%=mod;
        }
        x*=x;
        if(x>=mod)
            x%=mod;
        n>>=1;
    }
    return res;
}

const int inv4=qpow(4,mod-2);

inline int S2(int n,int m) {
    ll res=1ll*n*(n+1);
    if(res>=mod)
        res%=mod;
    res*=m;
    if(res>=mod)
        res%=mod;
    res*=(m+1);
    if(res>=mod)
        res%=mod;
    res*=inv4;
    if(res>=mod)
        res%=mod;
    return res;
}

const int inv6=qpow(6,mod-2);

inline int s2(int n) {
    ll res=1ll*n*(n+1);
    if(res>=mod)
        res%=mod;
    res*=(2ll*n+1);
    if(res>=mod)
        res%=mod;
    res*=inv6;
    if(res>=mod)
        res%=mod;
    return res;
}

unordered_map<int,int> GA;
inline int Get_A(int n){
    if(n<=MAXN)
        return A[n];
    if(GA.count(n))
        return GA[n];
    ll res=1;
    for(int l=2,r;l<=n;l=r+1){
        int t=n/l;
        r=n/t;
        ll tmp=s2(r)-s2(l-1);
        if(tmp<0)
            tmp+=mod;
        tmp*=Get_A(t);
        if(tmp>=mod)
            tmp%=mod;
        res-=tmp;
    }
    res%=mod;
    if(res<0)
        res+=mod;
    return GA[n]=res;
}

inline int S1(int n,int m) {
    ll res=0;
    int nm=min(n,m);
    for(int l=1,r; l<=nm; l=r+1) {
        int tn=n/l;
        int tm=m/l;
        r=min(n/tn,m/tm);
        ll tmp=Get_A(r)-Get_A(l-1);
        if(tmp<0)
            tmp+=mod;
        tmp*=S2(tn,tm);
        if(tmp>=mod)
            tmp%=mod;
        res+=tmp;
    }
    if(res>=mod)
        res%=mod;
    return res;
}

inline int s1(int l,int r) {
    ll res=(1ll*(l+r)*(r-l+1))>>1;
    if(res>=mod)
        res%=mod;
    return res;
}

inline int S(int n,int m) {
    ll res=0;
    int nm=min(n,m);
    for(int l=1,r; l<=nm; l=r+1) {
        int tn=n/l;
        int tm=m/l;
        r=min(n/tn,m/tm);
        ll tmp=s1(l,r);
        tmp*=S1(tn,tm);
        if(tmp>=mod)
            tmp%=mod;
        res+=tmp;
    }
    if(res>=mod)
        res%=mod;
    return res;
}

int main() {
#ifdef Yinku
    freopen("Yinku.in","r",stdin);
#endif // Yinku
    int n,m;
    scanf("%d%d",&n,&m);
    sieve();
    printf("%d\n",S(n,m));
    return 0;
}