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
",(ans+(n*C(m))%md+(m*C(n))%md)%md);
}