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

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

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

原式

$=\sum_{i=1}^N\sum_{j=1}^M\frac{i*j}{gcd(i.j)}$

$=\sum_{d=1}^{N}\sum_{i=1}^{\lfloor\frac{N}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{M}{d}\rfloor}dij[gcd(i,j)==1]$

$=\sum_{d=1}^{N}\sum_{i=1}^{\lfloor\frac{N}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{M}{d}\rfloor}dij\sum_{k|gcd(i,j)}\mu(k)$

$=\sum_{d=1}^{N}\sum_{i=1}^{\lfloor\frac{N}{dk}\rfloor}\sum_{j=1}^{\lfloor\frac{M}{dk}\rfloor}dijk^2\sum_{k=1}^{\lfloor\frac{N}{d}\rfloor} \mu(k)$

$=\sum_{d=1}^{N}d\sum_{k=1}^{\lfloor\frac{N}{d}\rfloor} k^2 \mu(k)\sum_{i=1}^{\lfloor\frac{N}{dk}\rfloor}i\sum_{j=1}^{\lfloor\frac{M}{dk}\rfloor}j$

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

$\sum_{d=1}^{N}d\sum_{k=1}^{\lfloor\frac{N}{d}\rfloor} k^2\mu(k)\sum_{i=1}^{\lfloor\frac{N}{dk}\rfloor}i\sum_{j=1}^{\lfloor\frac{M}{dk}\rfloor}j$中的$\lfloor\frac{N}{d}\rfloor$对于不同的$d$也是可以整除分块的;然后对于$k^2\mu(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[<<],*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=,fix=; register char ch; while(!isdigit(ch=getchar())) fix=ch=='-'?-:fix;

        do ret=ret*+(ch^); while(isdigit(ch=getchar())); return ret*fix;

    }

}using Fread::g;

const int M=,N=;

int mu[N],pri[N/],sum[N],n,m,cnt;

bool v[N]; ll Inv;

inline void MU(int n) { mu[]=;

    for(R i=;i<=n;++i) {

        if(!v[i]) pri[++cnt]=i,mu[i]=-;

        for(R j=;j<=cnt&&i*pri[j]<=n;++j) {

            v[i*pri[j]]=true;

            if(i%pri[j]==) break;

            mu[i*pri[j]]=-mu[i];

        }

    } for(R i=;i<=n;++i) sum[i]=(ll)(sum[i-]+(ll)(mu[i]+M)*i%M*i)%M;

}

inline int S(int x,int y) {return (ll)x*(x+)%M*Inv%M*y%M*(y+)%M*Inv%M;}

inline int F(int n,int m) {register ll ret=; n>m?swap(n,m):void();

    for(R l=,r;l<=n;l=r+) {

        r=min(n/(n/l),m/(m/l));

        ret=(ret+(ll)(sum[r]-sum[l-]+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(); MU(n); register ll ans=;Inv=(M+)/;

    for(R l=,r;l<=n;l=r+) {

        r=min(n/(n/l),m/(m/l));

        (ans+=(ll)(r-l+)*(l+r)%M*Inv%M*F(n/l,m/l)%M)%=M;

    } printf("%lld\n",ans);

}

2019.06.09

巴特西

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

最新文章

热门文章