BZOJ 3529 [Sdoi2014]数表【莫比乌斯反演】

有一张N×m的数表，其第i行第j列（1 < =i < =n，1 < =j < =m）的数值为
能同时整除i和j的所有自然数之和。给定a，计算数表中不大于a的数之和。

HINT

1 < =N．m < =10^5 ， 1 < =Q < =2×10^4

Solution

先假设a的限制不存在

正面把答案强行写出来是这样的

其中F(i)为i的约数和，可以线性筛处理

然后慢慢化简其中的各个部分

令g(i)为1<=x<=n,1<=y<=m,gcd(x,y)=i的数对(x,y)的个数

则有

用莫比乌斯反演代换

此时只需要一个的前缀和就可以继续愉快地分块了

然而此时还有a的限制，怎么继续做前缀和？

如果我们可以每次求前缀和的时候都只选择F( i )比a小的就好了

此时想到其实还有一种前缀和的方法树状数组，并且可以以加入的顺序控制前缀和

于是自然地离线询问，按a值排序，按F( i )排序后依次加入树状数组

取模用了PoPoQQQ的自然溢出

  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 100000
 15 #define maxn 100000+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 struct sigma{
 27   int number,F;
 28 }s[maxn];
 29 
 30 struct query{
 31   int n,m,a;
 32   int number;
 33 }q[maxn];
 34 
 35 int bit[maxn];
 36 int prime[maxn],pri[maxn],miu[maxn];
 37 int ans[maxn],tot=0,st=1;
 38 int T;
 39 
 40 void read()
 41 {
 42 #ifndef ONLINE_JUDGE
 43   freopen("3529.in","r",stdin);
 44   freopen("3529.out","w",stdout);
 45 #endif
 46   //cin >> T ;
 47   scanf("%d",&T);    
 48 }
 49 
 50 void write()
 51 {
 52   fr(i,1,T)
 53     //cout << (ans[i]&0x7fffffff) << endl ;
 54     printf("%d
",ans[i]&0x7fffffff);
 55 }
 56 
 57 void get()
 58 {
 59   miu[1]=1;
 60   fr(i,2,maxp){
 61     if( !prime[i] ) pri[++tot]=i,miu[i]=-1;
 62     int j=1;
 63     while( j<=tot && pri[j]*i<=maxp ){
 64       prime[i*pri[j]]=1;
 65       if( i%pri[j]==0 ){
 66     miu[i*pri[j]]=0;
 67     break;
 68       }
 69       miu[i*pri[j]]=-miu[i];
 70       j++;
 71     }
 72   }
 73   fr(i,1,maxp){
 74     fr(j,1,maxp/i)
 75       s[i*j].F+=i;
 76     s[i].number=i;
 77   }
 78 }
 79 
 80 bool comp(sigma a,sigma b)
 81 {
 82   return a.F<b.F;
 83 }
 84 
 85 bool cmp(query a,query b)
 86 {
 87   return a.a<b.a;
 88 }
 89 
 90 void add(int x,int w)
 91 {
 92   while( x<=maxp ){
 93     bit[x]+=w;
 94     x+=x&-x;
 95   }
 96 }
 97 
 98 int sum(int x)
 99 {
100   int res=0;
101   while( x>0 ){
102     res+=bit[x];
103     x-=x&-x;
104   }
105   return res;
106 }
107 
108 int calc(int x,int y)
109 {
110   if( x>y ) swap(x,y);
111   int res=0,pos;
112   int i=1;
113   while( i<=x ){
114     pos=min((x/(x/i)),(y/(y/i)));
115     res+=(sum(pos)-sum(i-1))*(x/i)*(y/i);
116     i=pos+1;
117   }
118   return res;
119 }
120 
121 void work()
122 {
123   fr(i,1,T){
124     //cin >> q[i].n >> q[i].m >> q[i].a ;
125     scanf("%d%d%d",&q[i].n,&q[i].m,&q[i].a);
126     q[i].number=i;
127   }
128   get();
129   sort(q+1,q+T+1,cmp);
130   sort(s+1,s+maxp+1,comp);
131   fr(i,1,T){
132     while( st<=maxp && s[st].F<=q[i].a ){
133       fr(j,1,maxp/s[st].number)
134     add(s[st].number*j,s[st].F*miu[j]);
135       st++;
136     }
137     ans[q[i].number]=calc(q[i].n,q[i].m);
138   }
139 }
140 
141 int main()
142 {
143   read();
144   work();
145   write();
146   return 0;
147 }

BZOJ 3529 [Sdoi2014]数表 【莫比乌斯反演】

Description

HINT

Solution

BZOJ 3529 [Sdoi2014]数表【莫比乌斯反演】