BZOJ 1101 [POI2007]Zap 【莫比乌斯反演】

Description

FGD正在破解一段密码，他需要回答很多类似的问题：对于给定的整数a,b和d，有多少正整数对x,y，满足x<=a，y<=b，并且gcd(x,y)=d。作为FGD的同学，FGD希望得到你的帮助。

数据范围：多组询问，1<=n<= 50000，每一组 1<=d<=a,b<=50000

Solution

题意：求1<=x<=a，1<=y<=b，满足gcd( x，y ) = d 的（x，y）对数

为了学莫比乌斯反演从黄学长那里找的模板题

最直接的想法当然是

for i 1 to a

for j 1 to b

if( gcd==d ) count++;

我们把它写成求和的形式

$\sum_{i=1}^a\sum_{j=1}^b[\text{gcd}(i,j)==d]$

令 $a'=\left \lfloor a/d \right \rfloor,b'=\left \lfloor b/d \right \rfloor$

$=\sum_{i=1}^{a'}\sum_{j=1}^{b'}[\text{gcd}(i,j)==1]$

莫比乌斯反演中有一个公式

这个公式可以用二项式定理随便推一下

观察以上两个式子

发现二式恰好完美契合了当且仅当 gcd==d 时数量加一的限制条件

于是我们考虑将gcd( i，j )用二式替换掉，并且将原式中 n 变为 gcd ( i，j )

此时得到一个新的式子

$=\sum_{i=1}^{a'}\sum_{j=1}^{b'}\sum_{d|\text{gcd}(i,j)}\mu(d)$

$d|\text{gcd}(i,j) \Leftrightarrow d|i \,\text{and} \,d|j$

于是我们改变枚举方式

在最外层枚举 gcd 可能的因数 d，然后再从 a' ，b' 中任取两个d的倍数

这样的方案数是有⌊a' / d ⌋⌊ b' / d ⌋ 的，并且转化枚举方式后本质与上式相同

于是得到

$=\sum_{d=1}^{a'}\mu(d)\left \lfloor a'/d \right \rfloor \left \lfloor b'/d \right \rfloor$

然后μ也有线性筛的办法，可以学一下

问题已经接近解决，复杂度大约 n*a ，n是数据组数

然而看一下数据50000*50000是过不了的

再仔细想一想，对于⌊a' / d ⌋，一定是多段连续的定值，也就是说对于一段区间上的 d，⌊a' / d ⌋ 的取值不会变化，⌊ b' / d ⌋ 同样

于是可以记录一下μ的前缀和，然后分块来求

但是 ⌊a' / d ⌋ 理论上有 √a' 个取值，⌊ b' / d ⌋ 同样，⌊a' / d ⌋⌊ b' / d ⌋ 取值不是 √a'*b' 个么，复杂度似乎并没有的到优化

实际上由于二者的除数d在同一时刻是相同的，所以在 d 一定的时候 ⌊a' / d ⌋ 仅对应一个 ⌊ b' / d ⌋ 取值仅有 √a'+√b' 个

于是复杂度变成了 n* √a'

问题解决

  1 #include<map>
  2 #include<cmath>
  3 #include<ctime>
  4 #include<queue>
  5 #include<stack>
  6 #include<cstdio>
  7 #include<climits>
  8 #include<iomanip>
  9 #include<cstring>
 10 #include<cstdlib>
 11 #include<iostream>
 12 #include<algorithm>
 13 
 14 #define maxp 50000
 15 #define maxn 50000+5
 16 #define set(a,b) memset(a,(b),sizeof(a))
 17 #define fr(i,a,b) for(ll i=(a),_end_=(b);i<=_end_;i++)
 18 #define rf(i,b,a) for(ll i=(a),_end_=(b);i>=_end_;i--)
 19 #define fe(i,a,b) for(int i=first[(b)],_end_=(a);i!=_end_;i=s[i].next)
 20 #define fec(i,a,b) for(int &i=cur[(b)],_end_=(a);i!=_end_;i=s[i].next)
 21 
 22 using namespace std;
 23 
 24 typedef long long ll;
 25 
 26 int prime[maxn],pri[maxn],tot=0;
 27 int sum[maxn],miu[maxn];
 28 int ans;
 29 int n,m;
 30 
 31 void read()
 32 {
 33 #ifndef ONLINE_JUDGE
 34   freopen("1101.in","r",stdin);
 35   freopen("1101.out","w",stdout);
 36 #endif
 37   //cin >> n ;
 38   scanf("%d",&n);
 39 }
 40 
 41 void write()
 42 {}
 43 
 44 void print()
 45 {
 46   //cout << ans << endl ;
 47   printf("%d\n",ans);
 48 }
 49 
 50 void mobius()
 51 {
 52   miu[1]=1;
 53   fr(i,2,maxp){
 54     if( !prime[i] ) pri[++tot]=i,miu[i]=-1;
 55     int j=1;
 56     while( j<=tot && pri[j]*i<=maxp ){
 57       prime[pri[j]*i]=1;
 58       if( i%pri[j]==0 ){
 59     　　miu[pri[j]*i]=0;
 60     　　break;
 61       }
 62       miu[pri[j]*i]=-miu[i];
 63       j++;
 64     }
 65   }
 66   fr(i,1,maxp)
 67     sum[i]=sum[i-1]+miu[i];
 68 }
 69 
 70 int cal(int x,int y)
 71 {
 72   int res=0,pos,i=1;
 73   if(x>y) swap(x,y);
 74   while(i<=x){
 75     pos=min(x/(x/i),y/(y/i));
 76     res+=(sum[pos]-sum[i-1])*(x/i)*(y/i);
 77     i=pos+1;
 78   }
 79   return res;
 80 }
 81 
 82 void work()
 83 {
 84   mobius();
 85   fr(i,1,n){
 86     int a,b,d;
 87     //cin >> a >> b >> d ;
 88     scanf("%d%d%d",&a,&b,&d);
 89     ans=cal(a/d,b/d);
 90     print();
 91   }
 92 }
 93 
 94 int main()
 95 {
 96   read();
 97   work();
 98   write();
 99   return 0;
100 }

posted @ 2015-07-01 13:29 ST_Saint 阅读(164) 评论(0) 编辑收藏举报

刷新页面返回顶部

ST_Saint

BZOJ 1101 [POI2007]Zap 【莫比乌斯反演】

Description

Solution

公告