BZOJ 2301 [HAOI2011]Problem b (容斥＋莫比烏斯反演＋分塊優化 詳解)

2301: [HAOI2011]Problem b

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

Description

Input

Output

Sample Input

Sample Output

14

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);   
    }
}