【HAOI2011】Problem b

· · 个人记录

学高考一段时间后我的OI实力居然越来越强了
原题：
对于给出的n个询问，每次求有多少个数对(x,y)，满足a≤x≤b，c≤y≤d，且gcd(x,y) = k，gcd(x,y)函数为x和y的最大公约数。
100%的数据满足：1≤n≤50000，1≤a≤b≤50000，1≤c≤d≤50000，1≤k≤50000

（我的推导似乎和大部分人的不太一样
这个是距离高考还有10天写的……
因为快高考了所以放飞自我，自习课不想复习拿反演玩一玩，发现这个挺简单的啊，以前怎么学不会呢
大概是学高考增强了我的推公式的能力吧

啊不扯了，直接开始题解 首先必须要想到的思路是四项容斥简化问题
问题转化为求有对少个数对 (x,y) 满足 1\leq x \leq n, 1\leq y \leq m 且 gcd(x,y)=1
看过其他反演题的话可能容易想到， \sum \sum [gcd(i,j)=1] 是比较好算的
那么可以把n和m都除等于k，问题就转化为求 \displaystyle\sum_{i=1}^{n} \sum_{j=1}^{m} [gcd(i,j)=1]
因为考虑到对于原来的满足条件的一个数对 (i,j) ，因为i和j的gcd是k，所以可以写成这样 (i'\cdot k,j'\cdot k) ，并且i'和j'的gcd是1
显然若 i'\cdot k\leq n ，则 i'\leq \lfloor \frac{n}{k} \rfloor，j'同理
现在的问题已经相对简单，剩下的就是推推推了~

ans(n,m)=\sum_{i=1}^{n}\sum_{j=1}^{m}[gcd(i,j)=1] =\sum_{i=1}^{n}\sum_{j=1}^{m}\sum_{d|gcd(i,j)}\mu(d)

这里应用到了

\sum_{d|n}\mu(d)=\begin{cases}0& n=1\\ 1& n\neq1 \end{cases}

的性质
因为 d|gcd(i,j) 这个表现形式不好，既不利于继续化简，亦不利于代码实现，所以现在换一种表示方式

ans(n,m)=\sum_{d}^{min(n,m)}\sum_{i=1}^{\lfloor \frac{n}{d} \rfloor}\sum_{j=1}^{\lfloor \frac{m}{d}\rfloor}\mu(d)

先枚举所有可能的gcd，再枚举所有因子有这个gcd的i和j，不难想象和先枚举ij再枚举d所枚举到的东西是一样的，这个似乎是莫比乌斯反演里常用的转化方法

ans(n,m)=\sum_{d}^{min(n,m)}\mu(d)\sum_{i=1}^{\lfloor \frac{n}{d} \rfloor}\sum_{j=1}^{\lfloor \frac{m}{d}\rfloor}1

因为d和ij的枚举情况无关所以可以提到前面来
后面的东西就可以秒算了

ans(n,m)=\sum_{d}^{min(n,m)}\mu(d)\cdot \lfloor \frac{n}{d} \rfloor \cdot \lfloor \frac{m}{d}\rfloor

至此公式推导完毕
我的推导和popoqqq的不太一样，他是构造了f和F并使用反演的公式，我是根据反演的性质来转化问题，现在看来两种方法本质是一样的
我的思路受《crash的表格》这道题的影响比较深
直接 O(n) 计算会T（因为有n组询问，需要用莫比乌斯反演常用的 O(\sqrt{n}) 求答案的计算方式，这个在po姐的ppt里解释得很清楚，不再赘述
代码： 

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
using namespace std;
int rd(){int z=0,mk=1;  char ch=getchar();
    while(ch<'0'||ch>'9'){if(ch=='-')mk=-1;  ch=getchar();}
    while(ch>='0'&&ch<='9'){z=(z<<3)+(z<<1)+ch-'0';  ch=getchar();}
    return z*mk;
}
int n,m,n1,n2,m1,m2;
int mu[51000],mus[51000];  //mu's!!!!!!
int prm[51000],prt=0;  bool prg[51000];
void gtmu(){
    memset(prg,0,sizeof(prg));
    mu[1]=mus[1]=1;
    for(int i=2;i<=50000;++i){
        if(!prg[i])  prm[++prt]=i,mu[i]=-1;
        for(int j=1;prm[j]*i<=50000 && j<=prt;++j){
            prg[prm[j]*i]=true;
            if(!(i%prm[j])){  mu[i*prm[j]]=0;  break;}
            mu[prm[j]*i]=-mu[i];
        }
        mus[i]=mus[i-1]+mu[i];
    }
}
int sm(int x,int y){  return (x+1)*x*(y+1)*y/4;}
int cclt(int x,int y){
    if(x>y)  swap(x,y);
    int tmp,bwl=0;
    for(int i=1;i<=x;i=tmp+1){
        tmp=min(x/(x/i),y/(y/i));
        //bwl+=sm(x/i,y/i)*(mus[tmp]-mus[i-1]);
        bwl+=(x/i)*(y/i)*(mus[tmp]-mus[i-1]);
    }
    return bwl;
}
int main(){freopen("ddd.in","r",stdin);
    gtmu();
    cin>>n;
    while(n --> 0){  //趋向于
        //cin>>n1>>n2>>m1>>m2>>m;
        n1=rd(),n2=rd(),m1=rd(),m2=rd(),m=rd();
        printf("%d\n",cclt(n2/m,m2/m)+cclt((n1-1)/m,(m1-1)/m)-cclt((n1-1)/m,m2/m)-cclt(n2/m,(m1-1)/m));
        //注意-1
    }
    return 0;
}