[HAOI2011][bzoj2301] Problem b [莫比乌斯反演+容斥原理+分块前缀和优化]

Posted 2020-10-25 Orion545

第一周练习 莫比乌斯反演 T2

题面：

传送门

有洛谷就尽量放洛谷链接呗，界面友好一点

思路：

和HDU1695比较像，但是这一回有50000组数据，直接莫比乌斯反演慢慢加的话会T

 

先解决一个前置问题：怎么处理a,c不是1的情况？

很简单，容斥原理搞之

我们设f(x,y)代表gcd(i,j)==e(1<=i<=x,1<=j<=y)的无序数对(i,j)的个数

那么本题答案相当于f(d,b)-f(c-1,b)-f(a-1,d)+f(a-1,c-1)

 

再来看反演超时的问题

我们注意到原反演过程中，f(1)==mu(i)*(d/i)*(b/i)

（对为什么这么做不太清楚的同学可移步上面HDU1695的那个链接）

对于后两项，在i很大的时候其实他们的值是基本不变动的，变化的只有mu[i]

那么我们可以利用这个过程

每一次，我们搜寻当前节点i的下一个“后两项的乘积改变了的”节点j

j的求法是min(d/(d/i).b/(b/i))，就是反过来求变化区间的大小，i越大，变化需要的时间越久

然后我们预处理mu的时候同时把mu的前缀和算出来，在上述情况下把(d/i)*(b/i)的值乘上sum[j]-sum[i-1]

循环结束以后，把i变成j+1，然后开始下一个循环，直到i的值超过min(b,d)

这就是莫比乌斯反演中的分块前缀和优化

顺便说一下，对于欧拉函数也可以利用这个优化，具体可以看这篇博客

 

Code：

 1 #include<iostream>
 2 #include<cstdio>
 3 #include<algorithm>
 4 #include<cstring>
 5 #define ll long long
 6 using namespace std;
 7 bool vis[100010];int mu[100010],pri[100010],cnt=0,sum[100010];
 8 void init(){
 9     int i=1,j,k;
10     mu[i]=sum[i]=1;
11     for(i=2;i<=100000;i++){
12         if(!vis[i]) pri[++cnt]=i,mu[i]=-1;
13         for(j=1;j<=cnt;j++){
14             k=i*pri[j];if(k>100000) break;
15             vis[k]=1;
16             if(i%pri[j]==0){mu[k]=0;break;}
17             else mu[k]-=mu[i];
18         }
19         sum[i]=sum[i-1]+mu[i];
20     }
21 }
22 inline int read(){
23     int re=0,flag=1;char ch=getchar();
24     while(ch>\'9\'||ch<\'0\'){
25         if(ch==\'-\') flag=-1;
26         ch=getchar();
27     }
28     while(ch>=\'0\'&&ch<=\'9\') re=(re<<1)+(re<<3)+ch-\'0\',ch=getchar();
29     return re*flag;
30 }
31 inline void swap(int &l,int &r){l^=r;r^=l;l^=r;}
32 ll f(int l,int r){
33     if(l>r) swap(l,r);ll re=0,i,j=0;
34     for(i=1;i<=l;i=j+1){
35         j=min(l/(l/i),r/(r/i));
36         re+=(ll)(sum[j]-sum[i-1])*(l/i)*(r/i);
37     }
38 //    cout<<"f "<<l<<ends<<r<<ends<<re<<endl;
39     return re;
40 }
41 int main(){
42     int a,b,c,d,e,T=read();
43     init();
44     while(T--){
45         a=read();b=read();c=read();d=read();e=read();
46         if(!e){printf("%d\\n",0);continue;}
47         b/=e;d/=e;a--;c--;a/=e;c/=e;
48         printf("%lld\\n",f(b,d)-f(a,d)-f(c,b)+f(a,c));
49     }
50 }