【51nod 1222】 最小公倍数计数

**题目来源： Project Euler
基准时间限制：6 秒 空间限制：131072 KB 分值: 640 难度：8级算法题**
定义F(n)表示最小公倍数为n的二元组的数量。
即：如果存在两个数（二元组）X，Y(X <= Y)，它们的最小公倍数为N，则F(n)的计数加1。
例如：F(6) = 5，因为[2,3] [1,6] [2,6] [3,6] [6,6]的最小公倍数等于6。

给出一个区间[a,b]，求最小公倍数在这个区间的不同二元组的数量。
例如：a = 4，b = 6。符合条件的二元组包括：
[1,4] [2,4] [4,4] [1,5] [5,5] [2,3] [1,6] [2,6] [3,6] [6,6]，共10组不同的组合。
Input
输入数据包括2个数：a, b，中间用空格分隔（1 <= a <= b <= 10^11)。
Output
输出最小公倍数在这个区间的不同二元组的数量。
Input示例
4 6
Output示例
10
可以推导,先考虑大小没有限制的二元组情况：

$$F’(n)=\sum_{i=1}^n\sum_{j=1}^n[\frac{ij}{gcd(i,j)}==n]\\= \sum_{d|n}\sum_{d|i}\sum_{d|j}[\frac{ij}{d}==n]\\ =\sum_{d|n}\sum_{i}^{\frac{n}{d}}\sum_{j}^{\frac{n}{d}}[ijd==n]\\ =\sum_{d|n}\sum_{i}^{\frac{n}{d}}\sum_{j}^{\frac{n}{d}}[ij==\frac{n}{d}]\\ =\sum_{d|n}\sum_{i=1}^d\sum_{j=1}^d[ij==d]\\ =\sum_{d|n}\sum_{k|d}1$$
但是这种是对于二元组没有大小关系限制时的答案。
设有$d(n)=\sum_{d|n}$即因子个数函数。

现在：

$$\sum_{i=1}^nF(i)=\sum_{i=1}^n\sum_{j|i}d(j)=\sum_{j=1}^nd(j)\lfloor\frac{n}{j}\rfloor$$

我们来考虑一下

$$G(n)=\sum_{i=1}^nd(i)$$
那么; $$G(n)=\sum_{i=1}d(i)=\sum_{i=1}\lfloor\frac{n}{i}\rfloor$$
然后分块处理：

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<map>
#include<cmath>
#define maxx 5000000
#define mod 1000000007
#define ll long long
using namespace std;

map<ll,ll>g;
ll G(ll n)
{
    if(g[n])return g[n];
    ll ans=0;
    //while((ans+1)*(ans+1)<=n)ans++;
    for(ll i=1,last;i<=n;i=last+1)
    {
        last=n/(n/i);
        ans+=(last-i+1)*(n/i);
    }
    g[n]=ans;
    return ans;
}
map<ll,ll>f;
ll F(ll n)
{
    if(f[n])return f[n];
    ll ans=0;
    for(ll i=1,last;i<=n;i=last+1)
    {
        last=n/(n/i);
        ans+=(G(last)-G(i-1))*(n/i);
    }
    f[n]=ans;
    return ans;
}
int main()
{
    ll a,b;
    cin>>a>>b;
    cout<<(F(b)-F(a-1)+b-a+1)/2<<endl;
    return 0;
}

但是这种推导下来，我能预计会超时，但是居然还wa了，实在不知道为啥，希望有大佬可以指出原因。

换另一种推法：

$$ans(n)=\sum_{i=1}^n\sum_{j=1}^n[\frac{ij}{gcd(i,j)}<=n?1:0]$$
枚举gcd: $$ans(n)=\sum_{d=1}^n\sum_{d|i}^n\sum_{d|j}^n[gcd(i,j)==d][\frac{ij}{d}<=n]\\= \sum_{d=1}^n\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{d}\rfloor}[gcd(i,j)==1][ijd<=n]$$
这里可以引入欧拉函数或者莫比乌斯函数，但是欧拉函数似乎没什么用啊。
用莫比乌斯搞一下： $$=\sum_{d=1}^n\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{k|gcd(i,j)}\mu(k)[ijd<=n]\\= \sum_{d=1}^n\sum_{k=1}^{\lfloor\frac{n}{d}\rfloor}\mu(k)\sum_{k|i}^{\lfloor\frac{n}{d}\rfloor}\sum_{k|j}^{\lfloor\frac{n}{d}\rfloor}[ijd<=n]\\= \sum_{d=1}^n\sum_{d|k}^{n}\mu(k)\sum_{i=1}^{\lfloor\frac{n}{dk}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{dk}\rfloor}[ijk^2d<=n]\\= \sum_{k=1}^n\mu(k)\sum_{d|k}\sum_{i=1}^{\lfloor\frac{n}{dk}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{dk}\rfloor}[ijk^2d<=n]\\= \sum_{k=1}^n\mu(k)\sum_{d=1}^{\lfloor\frac{n}{k}\rfloor}\sum_{i=1}^{\lfloor\frac{n}{dk}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{dk}\rfloor}[ijd<=\frac{n}{k^2}]$$

根据现在的条件约束，我们是可以适当的缩写范围的。

$$ans=\sum_{k=1}^{\sqrt n}\mu(k)\sum_{d=1}^{\lfloor\frac{n}{k^2}\rfloor}\sum_{i=1}^{\lfloor\frac{n}{dk^2}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{dk^2}\rfloor}[ijd<=\frac{n}{k^2}]$$

这里就只要去枚举（i,j,d）,因为这三个的地位是相等的，我们假设i<=j=d,那么i就只需要枚举$x^{\frac{1}{3}}$次
代码：

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<map>
#include<cmath>
#define maxx 316238
#define mod 1000000007
#define ll long long
#define LL long long
using namespace std;

bool isP[maxx];
int prime[maxx];
int cnt;
int mu[maxx];
void init()
{
    mu[1]=1;
    for(int i=2;i<maxx;i++)
    {
        if(!isP[i])
        {
            prime[cnt++]=i;
            mu[i]=-1;
        }
        for(int j=0;j<cnt&&(ll)i*prime[j]<maxx;j++)
        {
            isP[i*prime[j]]=true;
            if(i%prime[j])
                mu[i*prime[j]]=-mu[i];
            else
            {
                mu[i*prime[j]]=0;
                break;
            }

        }
    }
}

ll F(ll x)
{
    if(!x)return 0;
    ll ans=0;
    for(ll k=1;k*k<=x;k++)
    {
        if(mu[k]!=0)
        {
            ll m=x/(k*k),s=0;
            for(ll i=1;i*i*i<=m;i++)
            {
                for(int j=i+1;i*j*j<=m;j++)
                    s+=(m/(i*j)-j)*6+3;
                s+=(m/(i*i)-i)*3+1;
            }
            ans+=mu[k]*s;
        }
    }
    return (ans+x)/2;
}
int main()
{
    init();
    ll a,b;
    cin>>a>>b;
    cout<<F(b)-F(a-1)<<endl;
    return 0;
}