bzoj2301 [HAOI2011]Problem b

2301: [HAOI2011]Problem b

Time Limit: 50 Sec  Memory Limit: 256 MB
Submit: 6192  Solved: 2830
[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

14
3

HINT

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

分析：莫比乌斯反演套路题.n组询问，每次区间长度有50000，那么也就是说每一次单独询问的复杂度至少是O(sqrt(n))的.要求gcd(x,y)=k的数对(x,y)的数量，一个非常常见的套路就是同时除以k,就变成了找gcd(x/k,y/k)=1的数对数，一个比较好的处理方法是欧拉函数，只不过复杂度是O(n)的，不能接受.事实上问题已经被转化为求[1,p]的gcd = 1的数对数，通过容斥原理，可以得到最终的答案为f(b,d) - f(a-1,d) - f(b,c-1) + f(a-1,c-1).关键就是如何快速求出f的值.可以尝试莫比乌斯反演.一般能用莫比乌斯反演做的题都跟gcd有关.

          设f(i)为gcd(x,y)=1的(x,y)的数量，显然我们想求得的是f(1).那么F(i)就表示 i | gcd(x,y)的(x,y)的数量.通过莫比乌斯反演的第二种形式，可以表示为f(i) = Σμ(d/i) * F(d),d是i的倍数.F(i)非常好求，就是(n/i) * (m/i),其中n,m是x,y的取值上限.枚举i的倍数d就可以得到f(i).可是i=1,时间复杂度还是O(n),这要怎么优化呢？换一下枚举的方式，我们不能一个个地去枚举d，而是每次都能跳着枚举.因为F(d) = (n/d) * (m/d), (n/d)有很多结果是相同的，我们可以枚举d,d每次跳到min{(n/d - 1),(m/d - 1)}的d值上去.这样跳的这一段的F值是固定的，只需要算一下μ(d)的前缀和即可,结果最多有根号种，所以复杂度为O(sqrt(n)),就可以通过本题了.

          关于跳数有一点技巧，每次n/d,m/d都是下取整，跳的时候要跳到最大的能取到的值.

          莫比乌斯反演一个非常常用的套路就是枚举除法的结果，从n降到sqrt(n).一些关于gcd的常见变形要记住：gcd(x,y)=k，右边肯定是要变成1的,根据这个可以变形.求i | gcd(x,y)的对数，也有很常用的解法：(n/i) * (m/i).最核心的思路就是从已经会做的简单题中提炼出方法，将原问题不断的变形成几个小问题的结合，从而求解.

#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>

using namespace std;

typedef long long ll;

ll a, b, c, d, k, n, mo[50010], vis[50010], prime[50010], tot, sum[50010];

void init()
{
    mo[1] = 1;
    for (ll i = 2; i <= 50000; i++)
    {
        if (!vis[i])
        {
            prime[++tot] = i;
            mo[i] = -1;
        }
        for (ll j = 1; j <= tot; j++)
        {
            ll t = prime[j] * i;
            if (t > 50000)
                break;
            vis[t] = 1;
            if (i % prime[j] == 0)
            {
                mo[t] = 0;
                break;
            }
            mo[t] = -mo[i];
        }
    }
    for (ll i = 1; i <= 50000; i++)
        sum[i] = sum[i - 1] + mo[i];
}

ll solve(ll l, ll r)
{
    ll res = 0, last = 0;
    for (ll i = 1; i <= min(l, r); i = last + 1)
    {
        last = min(l / (l / i), r / (r / i)); //最大的使得l/last = l/i 并且 r/last = r/i的值
        res += (l / i) * (r / i) * (sum[last] - sum[i - 1]);
    }
    return res;
}

int main()
{
    init();
    scanf("%lld", &n);
    while (n--)
    {
        scanf("%lld%lld%lld%lld%lld", &a, &b, &c, &d, &k);
        printf("%lld\n", solve(b / k, d / k) - solve((a - 1) / k, d / k) - solve(b / k, (c - 1) / k) + solve((a - 1) / k, (c - 1) / k));
    }

    return 0;
}