【BZOJ2154】Crash的数字表格，数论练习之二维LCM（莫比乌斯反演）

Time:2016.06.18
Author:xiaoyimi
转载注明出处谢谢
看不见的分割线
传送门
思路：
设n≤m
$\displaystyle \sum^n_{i=1}\sum^m_{j=1}lcm(i,j)$
$\displaystyle =\sum^n_{i=1}\sum^m_{j=1}\frac{ij}{gcd(i,j)}$
$\displaystyle =\sum^n_{i=1}\sum^m_{j=1}\sum^n_{d=1}[d|i][d|j][gcd(i,j)=d]\frac{ij}{d}$
令i=di,j=dj
$\displaystyle =\sum^n_{d=1}\sum^\frac nd_{i=1}\sum^\frac md_{j=1}[gcd(i,j)=1]dij$
$\displaystyle =\sum^n_{d=1}d\sum^\frac nd_{i=1}\sum^\frac md_{j=1}\sum^\frac nd_{p=1}[p|i][p|j]μ(p)ij$
$\displaystyle =\sum^n_{d=1}d\sum^\frac nd_{p=1}μ(p)*p^2*\sum^\frac n{dp}_{i=1}\sum^\frac m{dp}_{j=1}ij$
$\displaystyle =\sum^n_{d=1}d\sum^\frac nd_{p=1}μ(p)*p^2*\frac{[\frac n{dp}][\frac m{dp}]([\frac n{dp}]+1)([\frac m{dp}]+1)}4$
两次分块求解，复杂度为$O((\sqrt n)^2=n)$
具体做法就是先通过n/d,m/d的取值对d分块，即在每一个块内对于不同的d来说，n/d,m/d的取值是分别相等的，然后再通过一个函数S(n,m)求得S(n/d,m/d)，其中S函数也是可以分块求得的，此处不再赘述
$\displaystyle S(n,m)$
$\displaystyle =\sum_{i=1}^n\sum_{j=1}^m[gcd(i,j)=1]ij$
$\displaystyle =\sum_{d=1}^nμ(d)*d^2*\frac{[\frac n{d}][\frac m{d}]([\frac n{d}]+1)([\frac m{d}]+1)}4$
注意：
$μ(d)*d^2$可以处理前缀和；
有很多地方会炸int甚至long long
代码：

#include<cstdio>
#include<iostream>
#include<cstring>
#define M 10000005
#define mo 20101009
#define LL long long
using namespace std;
int n,m,ans;
int prime[M>>3],mu[M];
bool vis[M];
int cal(int n,int m)
{
    int up,ans=0;
    for (int i=1;i<=n;i=up+1)
    {
        up=min(n/(n/i),m/(m/i));
        ans=(ans+((LL)(n/i)*(n/i+1)/2%mo)*((LL)(m/i)*(m/i+1)/2%mo)%mo*(mu[up]-mu[i-1])%mo)%mo;
    }
    return ans;
}
main()
{
    scanf("%d%d",&n,&m);
    if (n>m) swap(n,m);
    mu[1]=1;
    for (int i=2;i<=m;i++)
    {
        if (!vis[i]) prime[++prime[0]]=i,mu[i]=-1;
        for (int j=1;j<=prime[0];j++)
            if (i*prime[j]<=m)
            {
                vis[i*prime[j]]=1;
                if (i%prime[j]) mu[i*prime[j]]=-mu[i];
                else {mu[i*prime[j]]=0;break;}
            }
            else break;
    }
    for (int i=2;i<=m;i++) mu[i]=(mu[i-1]+(LL)mu[i]*i*i%mo)%mo;
    int up;
    for (int i=1;i<=n;i=up+1)
    {
        up=min(n/(n/i),m/(m/i));
        ans=(ans+(LL)(up-i+1)*(up+i)/2%mo*cal(n/i,m/i)%mo)%mo;
    }
    printf("%d",(ans+mo)%mo);
}