51Nod1222：最小公倍数计数 （莫比乌斯反演）

题目传送门：http://www.51nod.com/onlineJudge/questionCode.html#!problemId=1222

---

题目分析：这题又是skywalkert博客里的练习题，我一开始自己手推原式子推出了一个$O(nlog^2(n))$的东西（没有用$\mu$），完全没办法化简，还是要跑去看题解，后来发现别人的思路和我一开始就不是一个方向的……也许我数论还是太弱了吧。

---

题目要我们求lcm(x,y)在[a,b]内的数对(x,y)的数量，其中x<=y。我们先用差分，设F(n)表示lcm(x,y)在[1,n]内的数对(x,y)的数量(x<=y)，然后用F(b)-F(a-1)即可。
然后我们发现x<=y这个限制很难做，于是我们考虑G(n)=lcm(x,y)在[1,n]内的数对(x,y)的数量(不要求x<=y)，则F(n)=(G(n)+n)/2。
现在变成要求：

$$G(n)=\sum_{x=1}^{n}\sum_{y=1}^{n}[\frac{xy}{(x,y)}<=n]$$
接下来我们枚举d=(x,y)，令 $i=\frac{x}{d},j=\frac{y}{d}$得到：
$$G(n)=\sum_{d=1}^{n}\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{d}\rfloor}[ijd<= n]*[(i,j)=1]$$
接下来我们套用莫比乌斯反演的公式 $\sum_{d|n}\mu(d)=[n=1]$，将最右边的那个条件限制(i,j)=1改为枚举：
$$G(n)=\sum_{d=1}^{n}\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{d}\rfloor}[ijd<= n]\sum_{k|i,k|j}\mu(k)$$
然后我们把对k的枚举提到第二层（即d之后），重新令 $x=\frac{i}{k},y=\frac{j}{k}$，得到：
$$G(n)=\sum_{d=1}^{n}\sum_{k=1}^{\lfloor\frac{n}{d}\rfloor}\mu(k)\sum_{x=1}^{\lfloor\frac{n}{dk}\rfloor}\sum_{y=1}^{\lfloor\frac{n}{dk}\rfloor}[dk^2xy<=n]$$
我们可以把最外层的两层枚举看成是 先枚举d，然后再枚举d的倍数k，那么我们可以变成 先枚举倍数k，再枚举有多少个符合条件的约数d，也就是说 $\sum_{d=1}^{n}\sum_{k=1}^{\lfloor\frac{n}{d}\rfloor}\mu(k)$和 $\sum_{k=1}^{n}\mu(k)\sum_{d=1}^{\lfloor\frac{n}{k}\rfloor}$是一样的。
同时，我们把 $[dk^2xy<=n]$变形为 $[dxy<=\lfloor\frac{n}{k^2}\rfloor]$，得到：
$$G(n)=\sum_{k=1}^{n}\mu(k)\sum_{d=1}^{\lfloor\frac{n}{k}\rfloor}\sum_{x=1}^{\lfloor\frac{n}{kd}\rfloor}\sum_{y=1}^{\lfloor\frac{n}{kd}\rfloor}[dxy<=\lfloor\frac{n}{k^2}\rfloor]$$

---

推到这里之后，我们就可以看到这道题的巧妙之处了。关键在于最右边的$[dxy<=\lfloor\frac{n}{k^2}\rfloor]$，这意味着我们的k只需要枚举到$\sqrt{n}$就可以了，再枚举下去也不会对答案产生贡献。同时我们的d只需要枚举到$\lfloor\frac{n}{k^2}\rfloor$，x只需枚举到$\lfloor\frac{n}{k^2d}\rfloor$，y只需枚举到$\lfloor\frac{n}{k^2dx}\rfloor$即可。再枚举下去也不会对答案产生贡献，这样原式就变成了：

$$G(n)=\sum_{k=1}^{\sqrt{n}}\mu(k)H(\lfloor\frac{n}{k^2}\rfloor)$$
其中：
$$H(n)=\sum_{a=1}^{n}\sum_{b=1}^{\lfloor\frac{n}{a}\rfloor}\sum_{c=1}^{\lfloor\frac{n}{ab}\rfloor}1$$
那么我们怎么快速地求H(n)呢？我们不妨令a<=b<=c，然后将a从1枚举到 $n^{\frac{1}{3}}$，b从a枚举到 $\sqrt{\frac{n}{a}}$，直接算出c的个数为num。根据组合数学，当 $a<b<c$的时候，将6*num加进答案；当 $a=b<c$或 $a<b=c$的时候，将3*num加进答案；当a=b=c的时候，将num加进答案。
那么这样的时间复杂度是多少呢？网上神犇说可以用积分证得约为 $O(n^{\frac{2}{3}})$（然而我并不会证，我们学校搞数学竞赛的人好像也不会证……），总之6组数据的总时间为33s（然而出题人用10s，我的慢好多）。

---

CODE：

#include<iostream>
#include<string>
#include<cstring>
#include<cmath>
#include<cstdio>
#include<cstdlib>
#include<stdio.h>
#include<algorithm>
using namespace std;

const int maxn=400100;
typedef long long LL;

int miu[maxn];
bool vis[maxn];
int prime[maxn];
int cur=0;

LL a,b;

void Preparation()
{
    miu[1]=1;
    for (int i=2; i<maxn; i++)
    {
        if (!vis[i]) prime[++cur]=i,miu[i]=-1;
        for (int j=1; j<=cur && i*prime[j]<maxn; j++)
        {
            int k=i*prime[j];
            vis[k]=true;
            if (i%prime[j]) miu[k]=-miu[i];
            else
            {
                miu[k]=0;
                break;
            }
        }
    }
}

LL Get(LL x)
{
    LL temp=0;
    int A=(int)pow(x,1.0/3.0)+1;
    for (int a=1; a<=A; a++)
    {
        LL B=x/(long long)a;
        B=(long long)floor( sqrt( (long double)B )+0.0000001 );
        for (LL b=a+1; b<=B; b++)
        {
            LL c=x/( (long long)a*b );
            c-=b;
            if (c>0) temp+=( c*6LL );
            if ( (long long)a*b*b<=x ) temp+=3LL;
        }
    }
    for (int a=1; a<=A; a++)
    {
        LL c=x/( (long long)a*(long long)a );
        c-=a;
        if (c>0) temp+=( c*3LL );
        if ( (long long)a*(long long)a*(long long)a<=x ) temp++;
    }
    return temp;
}

LL Work(LL n)
{
    int sn=(int)floor( sqrt( (long double)n )+0.0000001 );
    LL temp=0;
    for (int i=1; i<=sn; i++)
    {
        LL x=n;
        x/=(long long)i;
        x/=(long long)i;
        if (!x) break;
        temp+=( (long long)miu[i]*Get(x) );
    }
    temp+=n;
    temp/=2LL;
    return temp;
}

int main()
{
    freopen("c.in","r",stdin);
    freopen("c.out","w",stdout);

    Preparation();
    cin>>a>>b;
    LL ans=Work(b)-Work(a-1);
    cout<<ans<<endl;

    return 0;
}