BZOJ 3529 [Sdoi2014]数表 【莫比乌斯反演】
Description
有一张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\n",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 }