BZOJ 2301 [HAOI2011]Problem b（莫比乌斯反演）

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

咱也不知道莫比乌斯函数，莫比乌斯反演怎么来的，直接用就得了...

①莫比乌斯函数：

else表示质因子分解后有质因子出现一次以上

莫比乌斯函数跟着欧拉筛跑一遍就可以求出来

②莫比乌斯反演的两种写法：

   

其中f(i)是题目询问的东西，你发现直接求求不出来。

但是你发现你可以轻松算出i的所有因子的f函数的和（写法1），

或者轻松算出i所有倍数的f函数的和（写法2），就可以利用一堆F（），和莫比乌斯函数来反推出这个f（i）

③这道题的思路：

对于a<=x<=b,c<=y<=d,询问gcd（x，y） =  k的对数，

我们可以把它转化为：

求1<=x<=m,1<=y<=n,询问gcd（x，y）= k的对数，假设用函数写作ff(m,n,k)

我们的答案就是 ff(b,d,k) - ff(a-1,d,k) - ff(b,c-1) + ff(a-1,c-1)

因此我们的问题就是怎么求    f(i)  ——  1<=x<=m,1<=y<=n,gcd（x，y）= i 的个数

我们看莫比乌斯反演两种F哪种可以简单的求出呢？

写法1：还是那个范围，求所有gcd(x,y) = i的因数  的 对数 和

写法2：还是那个范围，求所有gcd(x,y) = i的倍数  的 对数 和

我们发现写法2求的F(i)可以直接求出，因为只要x，y都是i的倍数，求出来的gcd(x,y)都是i或者i的倍数

因此x，y的对数为

于是对于每次询问的k,n,m

于是每次询问可以通过枚举k的倍数 O（N）的解决

但是还是太慢了

我们发现对于每次枚举的k的倍数，会存在某一段连续枚举的区间[x*k,y*k],被n整除的值都相同

所以用到下面讲的这个玩意儿

④整除分块：

举个例子，对于500这个数

除数251~500，整除结果都等于1

除数167~250，整除结果都等于2

......

我们把整除结果相同的数归入一个块

对于一个数i,:

n/(n/i)可以得到这个块的右界

因此你可以对于当前的i求得右界，这段区间整除的结果都是一样的直接整块处理。然后下一次直接让i变成上一个块右界+1即可

这里有个结论：块最多块

因此枚举块复杂度是O（sqrt(N)）级别的

⑤回到这题，怎么在这题中实现？

这题的问题在于两个是两个整除分块乘积，由于n,m大小不同，整除分块的右界不同，其中一个过了当前块的右界那个整除块的数就变了，所以我们枚举当前i的右界时，取m，n两个的右界里较小的，这样保证这一段乘积都一样。

这样枚举的复杂度为,因为最多这么多个右界

此时我们再看前面的莫比乌斯函数，当前i，我们求得右界last=min(n/(n/i),m/(m/i))

(i,last)每个k的倍数乘积都是,这段区间的莫比乌斯函数范围是(i/k,last/k)，我们莫比乌斯函数提前预处理前缀和

那么(i,last)这段区间共有（presum[last/k] - presum[i/k-1]）个

枚举完区间就OK了。

⑥后来，不过已经不重要了

发现题目转换成

求1<=x<=n/k,1<=y<=m/k,询问gcd(x,y) = 1貌似可以降个几倍的复杂度，因为之前总要枚举完区间min(n,m)，现在枚举到min(n/k,m/k)就行了

不过也就差个几倍，因为到后面块越来越大跳的很快的。

代码：

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

using namespace std;

#define ll long long
#define for1(i,a,b) for (int i=a;i<=b;i++)
#define for0(i,a,b) for (int i=a;i<b;i++)
#define rof1(i,a,b) for (int i=a;i>=b;i--)
#define rof0(i,a,b) for (int i=a;i>b;i--)
#define pb push_back
#define fi first
#define se second
#define debug(x) printf("----Line %s----\n",#x)
#define pt(x,y) printf("%s = %d\n",#x,y)
#define INF 0x3f3f3f3f
#define df(x) ll x;scanf("%I64d",&x)
#define df2(x,y) ll x,y;scanf("%I64d %I64d",&x,&y)
#define mod 1000000007
#define duozu(T) int T;scanf("%d",&T);while (T--)

const int N = 5e4+5;

int miu[N];
bool vis[N];
int prime[N>>2];

void euler()///欧拉筛求莫比乌斯函数
{
    int tot = 0;
    miu[1] = 1;
    for (int i=2;i<=N-5;i++){
        if (!vis[i]){
            prime[tot++] = i;
            miu[i] = -1;
        }
        for (int j=0;j<tot && prime[j]*i<=N-5;j++){
            vis[i*prime[j]] = true;
            if (i%prime[j]==0){
                miu[i*prime[j]] = 0;
                break;
            }
            else {
                miu[i*prime[j]] = -miu[i];
            }
        }
    }
}

ll query(int n,int m)
{
    ll ans = 0;
    int last;
    if (n>m) swap(n,m);
    for (int i=1;i<=n;i=last+1){
        last = min(n/(n/i),m/(m/i));
        ans += 1LL*(n/i)*(m/i)*(miu[last]-miu[i-1]);
    }
    return ans;
}

ll query(int n,int m,int k)
{
    ll ans = 0;
    int last;
    if (n>m) swap(n,m);
    for (int i=1;i<=n;i=last+1){
        last = min(n/(n/i),m/(m/i));
        ans += 1LL*(n/i)*(m/i)*(miu[last/k]-miu[(i-1)/k]);
    }
    return ans;
}

int main()
{
    //freopen("C:/Users/DELL/Desktop/input.txt", "r", stdin);
    //freopen("C:/Users/DELL/Desktop/output.txt", "w", stdout);
    euler();
    for1(i,2,N-5) miu[i] += miu[i-1];


    duozu(T){
        int a,b,c,d,k;
        scanf("%d %d %d %d %d",&a,&b,&c,&d,&k);
        ll ans = 0;

        ans = query(b/k,d/k) - query(b/k,(c-1)/k) - query((a-1)/k,d/k) + query((a-1)/k,(c-1)/k);//这样快

        //ans = query(b,d,k) - query(b,c-1,k) - query(a-1,d,k) + query(a-1,c-1,k);//这样慢

        printf("%lld\n",ans);
    }

    return 0;
}