BZOJ 2154 Crash的数字表格 ——莫比乌斯反演

求$sum_{i=1}^nsum_{j=1}^n lcm(i,j)$

枚举因数

$ans=sum_{d<=n} F(d) * d$

$F(d)$表示给定范围内两两$sum_{gcd(i,j)=d} i*j $

令$f(p)=Sum(lfloor n/p floor) Sum(lfloor m/p floor) * p^2$

那么 $f(i)=sum_{i mid n}F(n)$

反演得到$F(i)=sum_{i mid n} mu(n/i) f(n)$

那么我们代入就得到了

$ans=sum_{d<=n}d*sum_{i<=lfloor n/d floor} i^2 *mu(i)* Sum(lfloor frac {n}{i*d} floor,lfloor frac {m}{i*d} floor)$

然后外面分块一次，里面分块一次

时间复杂度$Theta (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 inv 10050505LL
#define maxn 10000005
#define md 20101009LL
  
int mu[maxn],pr[maxn],top=0;
ll ps[maxn];
bool vis[maxn];
  
int n,m;
 
void init()
{
    memset(vis,false,sizeof vis);
    mu[1]=1;ps[1]=1;
    F(i,2,n)
    {
        if (!vis[i]) mu[i]=-1,pr[++top]=i;
        F(j,1,top)
        {
            if (pr[j]*i>n) break;
            vis[pr[j]*i]=true;
            if (i%pr[j]==0) {mu[i*pr[j]]=0;break;}
            mu[i*pr[j]]=-mu[i];
        }
        ps[i]=(ps[i-1]+((ll)mu[i]*i*i))%md;
    }
}
  
ll sum(int n,int m)
{
    n=((ll)n*(n+1)/2)%md;
    m=((ll)m*(m+1)/2)%md;
    return ((ll)n*m)%md;
}
  
ll Function(int n,int m)
{
    if (n>m) swap(n,m);
    ll ret=0;
    for (int i=1,last=0;i<=n;i=last+1)
    {
        last=min(n/(n/i),m/(m/i));
        ret=(ret+((sum(n/i,m/i))*(ps[last]-ps[i-1]+md)%md)%md)%md;
    }
    return ret;
}
  
ll S(int n)
{
    return ((1LL+n)*n/2)%md;
}
  
ll solve(int n,int m)
{
    if (n>m) swap(n,m);
    ll ret=0;
    for (int i=1,last=0;i<=n;i=last+1)
    {
        last=min(n/(n/i),m/(m/i));
        ret=(ret+((ll)Function(n/i,m/i))*(S(last)-S(i-1))%md+md)%md;
    }
    return ret;
}
  
  
int main()
{
    scanf("%d%d",&n,&m);
    if (n<m) swap(n,m);
    init();
    printf("%lld
",solve(n,m));
}