P1829 [国家集训队]Crash的数字表格 / JZPTAB 莫比乌斯反演

又一道。。。分数和取模次数成正比$qwq$

---

求:$sum_{i=1}^Nsum_{j=1}^Mlcm(i,j)$

原式

$=sum_{i=1}^Nsum_{j=1}^Mfrac{i*j}{gcd(i.j)}$

$=sum_{d=1}^{N}sum_{i=1}^{lfloorfrac{N}{d} floor}sum_{j=1}^{lfloorfrac{M}{d} floor}dij[gcd(i,j)==1]$

$=sum_{d=1}^{N}sum_{i=1}^{lfloorfrac{N}{d} floor}sum_{j=1}^{lfloorfrac{M}{d} floor}dijsum_{k|gcd(i,j)}mu(k)$

$=sum_{d=1}^{N}sum_{i=1}^{lfloorfrac{N}{dk} floor}sum_{j=1}^{lfloorfrac{M}{dk} floor}dijk^2sum_{k=1}^{lfloorfrac{N}{d} floor} mu(k)$

$=sum_{d=1}^{N}dsum_{k=1}^{lfloorfrac{N}{d} floor} k^2 mu(k)sum_{i=1}^{lfloorfrac{N}{dk} floor}isum_{j=1}^{lfloorfrac{M}{dk} floor}j$

其中，$sum_{i=1}^{lfloorfrac{N}{dk} floor}isum_{j=1}^{lfloorfrac{M}{dk} floor}j$为等差数列，可以$O(1)$求,并且对于不同的$k$是可以整除分块的；

$sum_{d=1}^{N}dsum_{k=1}^{lfloorfrac{N}{d} floor} k^2mu(k)sum_{i=1}^{lfloorfrac{N}{dk} floor}isum_{j=1}^{lfloorfrac{M}{dk} floor}j$中的$lfloorfrac{N}{d} floor$对于不同的$d$也是可以整除分块的;然后对于$k^2mu(k)$先线性筛出来，再做个前缀和。

所以时间复杂度是$O(N)$的。

#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cstring>
#include<cmath>
#include<cctype>
#include<cstdlib>
#include<vector>
#include<queue>
#include<map>
#include<set>
#define ll long long
#define R register int
using namespace std;
namespace Fread {
    static char B[1<<15],*S=B,*D=B;
    #define getchar() (S==D&&(D=(S=B)+fread(B,1,1<<15,stdin),S==D)?EOF:*S++)
    inline int g() {
        R ret=0,fix=1; register char ch; while(!isdigit(ch=getchar())) fix=ch=='-'?-1:fix;
        do ret=ret*10+(ch^48); while(isdigit(ch=getchar())); return ret*fix;
    }
}using Fread::g;
const int M=20101009,N=10000010;
int mu[N],pri[N/10],sum[N],n,m,cnt;
bool v[N]; ll Inv;
inline void MU(int n) { mu[1]=1;
    for(R i=2;i<=n;++i) { 
        if(!v[i]) pri[++cnt]=i,mu[i]=-1;
        for(R j=1;j<=cnt&&i*pri[j]<=n;++j) {
            v[i*pri[j]]=true;
            if(i%pri[j]==0) break;
            mu[i*pri[j]]=-mu[i];
        }
    } for(R i=1;i<=n;++i) sum[i]=(ll)(sum[i-1]+(ll)(mu[i]+M)*i%M*i)%M;
}
inline int S(int x,int y) {return (ll)x*(x+1)%M*Inv%M*y%M*(y+1)%M*Inv%M;}
inline int F(int n,int m) {register ll ret=0; n>m?swap(n,m):void(0);
    for(R l=1,r;l<=n;l=r+1) {
        r=min(n/(n/l),m/(m/l));
        ret=(ret+(ll)(sum[r]-sum[l-1]+M)*S(n/l,m/l)%M)%M;
    } return ret;
}
signed main() {
#ifdef JACK
    freopen("NOIPAK++.in","r",stdin);
#endif
    n=g(),m=g(); n>m?swap(n,m):void(0); MU(n); register ll ans=0;Inv=(M+1)/2;
    for(R l=1,r;l<=n;l=r+1) {
        r=min(n/(n/l),m/(m/l));
        (ans+=(ll)(r-l+1)*(l+r)%M*Inv%M*F(n/l,m/l)%M)%=M;    
    } printf("%lld
",ans);    
}

---

2019.06.09