51nod 1222 最小公倍数计数

最新推荐文章于 2021-11-13 15:03:19 发布

原创最新推荐文章于 2021-11-13 15:03:19 发布 · 500 阅读

1 ·

CC 4.0 BY-SA版权

文章标签：

#信息学 #莫比乌斯反演

51nod 同时被 2 个专栏收录

21 篇文章

订阅专栏

莫比乌斯反演

16 篇文章

订阅专栏

本文介绍了一种求解特定区间内满足条件的二元组[x,y]数量的方法，其中条件为a≤lcm(x,y)≤b且x≤y。通过数学推导和编程实现，提供了一个高效的解决方案。

题目大意：

求有多少个二元组[x,y]满足a<=lcm(x, y)<=b且x<=y。
1<=a<=b<= $10^{11}$

题解：

先不考虑x<=y的情况，求出答案以后加上(b - a + 1)再除以2即可。
要满足a<=lcm(x,y)<=b,可以先求出1 ~ a - 1, 1 ~ b的，再相减。
由于我很弱，所以我一步一步来：

$Ans = \sum_{i = 1}^n \sum_{j = 1}^n [lcm(i, j) <= n]$
$= \sum_{i = 1}^n \sum_{j = 1}^n [i * j / gcd(i , j) <= n]$
$= \sum_{d = 1}^n \sum_{i = 1}^n \sum_{j = 1}^n [i * j * d <= n] * [gcd(i, j) = 1]$

设 $f(x) = \sum_{i = 1}^n \sum_{j = 1}^n [i * j * d <= n] * [gcd(i, j) = x]$
$g(x) = \sum_{i = 1}^{n / x} f(i * x)$
这里是 $f(1)$
$f(1) = \sum_{d' = 1}^n μ(d') * g(d')$
$= \sum_{d' = 1}^n μ(d') *\sum_{i' = 1}^n \sum_{j' = 1}^n [d * d'^2 * i' * j' <= n]$
上面这个转换有点难理解。在这里d’相当于原来i,j的最大公约数,所以是 $d'^2$

$Ans = \sum_{d = 1}^n \sum_{d' = 1}^n μ(d') *\sum_{i' = 1}^n \sum_{j' = 1}^n [d * d'^2 * i' * j' <= n]$
$= \sum_{d' = 1}^n μ(d') * \sum_{d = 1}^n \sum_{i = 1}^n \sum_{j = 1}^n[d'^2 * d * i * j <= n]$
缩小范围：
$=\sum_{d' = 1}^{\sqrt{n}} μ(d') * \sum_{d = 1}^n \sum_{i = 1}^n \sum_{j = 1}^n[d'^2 * d * i * j <= n]$

我们可以强行使d<=i<=j，这样的复杂度用微积分证明是 $O(n^{2/3})$
那么求出三个都相同的个数，有两个相同的个数，三个都不同的个数，乘上对应的组合数，就是没有顺序的方案数了。

感谢老曹对中间反演过程指点。

Code：

#include<cmath>
#include<cstdio>
#define fo(i, x, y) for(int i = x; i <= y; i ++)
#define ll long long

const int Maxn = 320000;

bool bz[Maxn + 1]; int mu[Maxn + 1], p[Maxn];
ll a, b;

ll calc(ll n) {
    ll s = 0;
    for(ll a = 1; a * a * a <= n; a ++) {
        for(ll b = a; a * b * b <= n; b ++) {
            s += (n / a / b - b + 1) * 6;
            if(a == b)
                s -= (n / a / b - b) * 3; else
                s -= 3;
        }
        s -= 5;
    }
    return s;
}

ll get_ans(ll n) {
    if(!n) return 0;
    ll s = 0;
    fo(i, 1, sqrt(n))
        if(mu[i])
            s += mu[i] * calc(n / i / i);
    return (s + n) / 2;
}

int main() {
    mu[1] = 1;
    fo(i, 2, Maxn) {
        if(!bz[i]) p[++ p[0]] = i, mu[i] = -1;
        fo(j, 1, p[0]) {
            int k = i * p[j];
            if(k > Maxn) break;
            bz[k] = 1;
            if(i % p[j] == 0) {
                mu[k] = 0; break;
            }
            mu[k] = -mu[i];
        }
    }
    scanf("%lld %lld", &a, &b);
    printf("%lld\n", get_ans(b) - get_ans(a - 1));
}