BZOJ 3518 点组计数 ——莫比乌斯反演

要求$ans=\sum_{i=1}^n \sum_{j=1}^m (n-i)(m-j)(gcd(i,j)-1)$

可以看做枚举矩阵的大小，然后左下右上必须取的方案数。

这是斜率单增的情况

然后大力反演即可。

最后$ans=ans*2+C(n,3)*m+C(m,3)*n$

$\Theta (n \log n)$

#include <cstdio>

#include <cstring>

#include <iostream>

#include <algorithm>

using namespace std;

#define F(i,j,k) for (int i=j;i<=k;++i)

#define D(i,j,k) for (int i=j;i>=k;--i)

#define ll long long

#define md 1000000007

#define inf 0x3f3f3f3f

#define maxn 50005

ll vis[maxn],mu[maxn],pr[maxn],top;

void init1()

{

    mu[1]=1;

    F(i,2,maxn)

    {

        if (!vis[i])

        {

            mu[i]=-1;

            pr[++top]=i;

        }

        F(j,1,top)

        {

            if ((ll)i*pr[j]>=maxn) break;

            vis[i*pr[j]]=1;

            if (i%pr[j]==0) {mu[i*pr[j]]=0;break;}

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

        }

    }

}

ll f1[maxn],f2[maxn],f3[maxn],ans=0;

ll Sum(ll n)

{

    n=(((n+1)*n)>>1)%md;

    return n;

}

void solve(ll n,ll m)

{

    if (n>m) swap(n,m);

    ll ret=0;

    F(d,1,n)

    {

        ll tmp=0;

        F(p,1,n/d)

        {

            tmp+=mu[p]*(n/p/d)*(m/p/d)*m*n; tmp%=md;

            tmp+=mu[p]*d*d*p*p*Sum(n/p/d)*Sum(m/p/d); tmp%=md;

            tmp-=mu[p]*m*d*p*(m/p/d)*Sum(n/p/d); tmp%=md;

            tmp-=mu[p]*n*d*p*(n/p/d)*Sum(m/p/d); tmp%=md;

        }

        ret+=tmp*(d-1);

    }

    ans=(2*ret)%md;

}

ll n,m;

ll C(ll n)

{

    n%=md;

    return (n*(n-1)*(n-2)/6)%md;

}

int main()

{

    init1();//init2();

    scanf("%lld%lld",&n,&m);

    solve(n,m);

    printf("%lld\n",(ans+(n*C(m))%md+(m*C(n))%md)%md);

}

巴特西

BZOJ 3518 点组计数 ——莫比乌斯反演

最新文章

热门文章