BZOJ 2301 [HAOI2011]Problem b (容斥＋莫比乌斯反演＋分块优化详解)

最新推荐文章于 2019-09-09 08:19:50 发布

原创最新推荐文章于 2019-09-09 08:19:50 发布 · 3k 阅读

1 ·

CC 4.0 BY-SA版权

文章标签：

#BZOJ #莫比乌斯反演 #容斥

ACM 同时被 3 个专栏收录

635 篇文章

订阅专栏

数论

78 篇文章

订阅专栏

组合数学

56 篇文章

订阅专栏

本文针对HAOI2011 Problemb题目进行了解析，介绍了一种利用容斥原理和分块求和的方法来计算在给定区间内满足特定最大公约数条件的数对数量。

2301: [HAOI2011]Problem b

Time Limit: 50 Sec Memory Limit: 256 MB
Submit: 2096 Solved: 909
[Submit][Status][Discuss]

Description

对于给出的n个询问，每次求有多少个数对(x,y)，满足a≤x≤b，c≤y≤d，且gcd(x,y) = k，gcd(x,y)函数为x和y的最大公约数。

Input

第一行一个整数n，接下来n行每行五个整数，分别表示a、b、c、d、k

Output

共n行，每行一个整数表示满足要求的数对(x,y)的个数

Sample Input

2

2 5 1 5 1

1 5 1 5 2

Sample Output

HINT

100%的数据满足：1≤n≤50000，1≤a≤b≤50000，1≤c≤d≤50000，1≤k≤50000

题目链接：http://www.lydsy.com/JudgeOnline/problem.php?id=2301

题目分析：首先n就是五万，因此每次即使是O(n)的计算总的下来也是O(n^2)了，因此对于每次操作我们要让时间复杂度小于O(n)，题意很清楚

先将原式变形：

a / k <= x / k <= b / k，c / k <= y / k <= d / k，gcd(x / k, y / k) = 1，令cal(b / k，d / k)为1 <= x / k <= b / k，1 <= y / k <= d / k时取出的满足条件的个数，则根据容斥原理有

ans = cal(b / k，d / k) - cal((a - 1) / k，d / k) - cal((c - 1) / k，b / k) + cal((a - 1) / k，(c - 1) / k)，因为1～a-1和1～c-1都不在我们所求的范围内，又减的时候这段区间减了两次，因此要再加上一个，接下来看cal函数，这里要用到分块求和，如果不做优化，就是直接枚举公约数ans += mob[i] * (l / i) * (r / i)但这样会超时，考虑到不能整除的特性，在很大一段区间内(l / i)和(r / i)的值是相同的，举个简单的例子，l = 10，r = 11那么可以看出 i从6到10，l / i和r / i的值都是1，因此考虑分块求和，从i开始最长的相等区间长度为min(l / (l / i)，r / (r / i))，这里如果不能理解的话，比如l / (l / i)，设l / i ＝ p，p表示整除时的值，那么l / p就是从i开始整除值为p的个数了

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
int const MAX = 50005;
int p[MAX], mob[MAX], sum[MAX];
bool prime[MAX];
int a, b, c, d, k;

void Mobius()
{
    int pnum = 0;
    memset(prime, true, sizeof(prime));
    memset(sum, 0, sizeof(sum));
    mob[1] = 1;
    sum[1] = 1;
    for(int i = 2; i < MAX; i++)
    {
        if(prime[i])
        {
            p[pnum ++] = i;
            mob[i] = -1;
        }
        for(int j = 0; j < pnum && i * p[j] < MAX; j++)
        {
            prime[i * p[j]] = false;
            if(i % p[j] == 0)
            {
                mob[i * p[j]] = 0;
                break;
            }
            mob[i * p[j]] = -mob[i];
        }
        sum[i] = sum[i - 1] + mob[i];
    }
}

int cal(int l, int r)
{
    if(l > r)
        swap(l, r);
    int ans = 0;
    for(int i = 1, last = 0; i <= l; i = last + 1)
    {
        last = min(l / (l / i), r / (r / i));
        ans += (l / i) * (r / i) * (sum[last] - sum[i - 1]);
    }
    return ans;
}

int main()
{
    Mobius();
    int n;
    scanf("%d", &n);
    while(n --)
    {
        scanf("%d %d %d %d %d", &a, &b, &c, &d, &k);
        int ans = 0;
        ans += cal(b / k, d / k);
        ans -= cal((a - 1) / k, d / k);
        ans -= cal((c - 1) / k, b / k);
        ans += cal((a - 1) / k, (c - 1) / k);
        printf("%d\n", ans);   
    }
}