51nod 1222 最小公倍数计数 莫比乌斯反演+杜教筛

题意

定义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组不同的组合。
1 <= a <= b <= 10^11

分析

答案就、可以用$\sum\limits_{i=1}^bF(i)-\sum\limits_{i=1}^{a-1}F(i)$来表示，现在就看如何求$\sum\limits_{i=1}^nF(i)$。
先不考虑二元组的顺序问题，只要最后把答案加上n再/2就好了。
我们要求的是

$$Ans=\sum_{i=1}^n\sum_{j=1}^n[lcm(i,j)<=n]$$
化简一下后可以得到 $$Ans=\sum_{d=1}^n\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{d}\rfloor}[gcd(i,j)==1][ij<=\lfloor\frac{n}{d}\rfloor]$$
大力反演一波后不难得到 $$Ans=\sum_{d=1}^{\sqrt n}\mu(d)*\sum_i\sum_j\sum_r[ijr<=\lfloor\frac{n}{d^2}\rfloor]$$
显然前面那部分可以用杜教筛，后面那部分的话，考虑求有多少 $a<=b<=c$满足 $abc<=n$。
考虑枚举 $a$到$n^{\frac{1}{3}}$，设 $t=n^{\frac{1}{3}}$，那么 $b$只要从$a$枚举到 $\sqrt t$就可以了。
我们再分别考虑三个数相等，两个数相等和都不相等的情况，再乘上对应组合数即可。

跑的飞快。

代码

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<map>
using namespace std;

typedef long long LL;

const int N=10000005;

int prime[N],tot,mu[N];
bool not_prime[N];
map<LL,LL> w;

void get_prime(int n)
{
    mu[1]=1;
    for (int i=2;i<=n;i++)
    {
        if (!not_prime[i]) prime[++tot]=i,mu[i]=-1;
        for (int j=1;j<=tot&&i*prime[j]<=n;j++)
        {
            not_prime[i*prime[j]]=1;
            if (i%prime[j]==0) break;
            mu[i*prime[j]]=-mu[i];
        }
    }
    for (int i=2;i<=n;i++) mu[i]+=mu[i-1];
}

LL get_mu(LL n)
{
    if (n<=10000000) return mu[n];
    if (w[n]) return w[n];
    LL ans=1;
    for (LL i=2,last;i<=n;i=last+1)
    {
        LL t=n/i;
        last=n/t;
        ans-=(last-i+1)*get_mu(t);
    }
    return w[n]=ans;
}

LL get_s(LL n)
{
    LL ans=0;
    for (LL i=1;i*i*i<=n;i++)
    {
        LL t=n/i;
        ans+=(t/i-i)*3+1;
        for (LL j=i+1;j*j<=t;j++) ans+=(t/j-j)*6+3;
    }
    return ans;
}

LL solve(LL n)
{
    LL ans=0;
    for (LL i=1,last;i*i<=n;i=last+1)
    {
        LL t=n/i/i;
        last=sqrt(n/t);
        ans+=(get_mu(last)-get_mu(i-1))*get_s(t);
    }
    return (ans+n)/2;
}

int main()
{
    get_prime(10000000);
    LL a,b;
    cin>>a>>b;
    printf("%lld",solve(b)-solve(a-1));
    return 0;
}