BZOJ3529: [Sdoi2014]数表(莫比乌斯反演,离线)
Description
有一张 n×m 的数表,其第 i 行第 j 列(1 <= i <= n, 1 <= j <= m)的数值为
能同时整除 i 和 j 的所有自然数之和。给定 a , 计算数表中不大于 a 的数之和。
Input
输入包含多组数据。
输入的第一行一个整数Q表示测试点内的数据组数
接下来Q行,每行三个整数n,m,a(|a| < =10^9)描述一组数据。
1 < =N.m < =10^5 , 1 < =Q < =2×10^4
Output
对每组数据,输出一行一个整数,表示答案模2^31的值。
Sample Input
2
4 4 3
10 10 5
4 4 3
10 10 5
Sample Output
20
148
148
解题思路:
这道题就是让我们求${\sum_{i=1}^{N}}{\sum_{j=1}^{M}}{\sigma(gcd(i,j))}({\sigma(gcd(i,j))}<=a)$
a比较让人恶心,考虑将${\sigma(gcd(i,j))}$按次序加入答案,直接统计,就是一种离线的做法。
那么就不需要考虑a了,答案就变成了(设n<=m)
${\sum_{i=1}^{n}}{\sum_{j=1}^{m}}{\sigma(gcd(i,j))}$
$={\sum_{d=1}^{n}}{\sigma(d)}{\sum_{d|i}}{n}{\sum_{d|j}^{m}}[gcd(i,j)==d]$
$={\sum_{d=1}^{n}}{\sigma(d)}{\sum_{i=1}^{\left \lfloor {\frac{n}{d}} \right \rfloor}}{\sum_{j=1}^{\left \lfloor {\frac{n}{d}} \right \rfloor}}{[gcd(i,j)==1]}$
$={\sum_{d=1}^{n}}{\sigma(d)}{\sum_{k=1}^{\left \lfloor {\frac{n}{d}} \right \rfloor}}{\mu(k)}{\left \lfloor {\frac{n}{dk}} \right \rfloor}{\left \lfloor {\frac{m}{dk}} \right \rfloor}$
设T=dk
$={\sum_{T=1}^{n}}{\sum_{d|T}}{\sigma(d)}{\mu(\frac{T}{d})}{\left \lfloor {\frac{n}{T}} \right \rfloor}{\left \lfloor \frac{m}{T} \right \rfloor}$
按顺序加入${\sigma(d)}$就好了
由于具有循环性质,在倍数位置加上就好了
取模卡我好长时间,没想到是这种方法。
代码:
1 #include<cstdio> 2 #include<cstring> 3 #include<algorithm> 4 typedef long long lnt; 5 const int N=300010; 6 struct int_2{int a,b;bool friend operator < (int_2 x,int_2 y){if(x.a!=y.a)return x.a<y.a;return x.b<y.b;}}F[N]; 7 int prime[N]; 8 int miu[N]; 9 bool vis[N]; 10 int cnt; 11 int line[N]; 12 struct qust{ 13 int n,m,a; 14 int no; 15 int ans; 16 }q[N]; 17 int lowbit(int x) 18 { 19 return x&(-x); 20 } 21 void update(int pos,int v) 22 { 23 while(pos<N) 24 { 25 line[pos]+=v; 26 pos+=lowbit(pos); 27 } 28 return ; 29 } 30 int query(int pos) 31 { 32 int ans=0; 33 while(pos) 34 { 35 ans+=line[pos]; 36 pos-=lowbit(pos); 37 } 38 return ans; 39 } 40 void gtp(void) 41 { 42 miu[1]=1; 43 for(int i=2;i<N;i++) 44 { 45 if(!vis[i]) 46 { 47 prime[++cnt]=i; 48 miu[i]=-1; 49 } 50 for(int j=1;j<=cnt&&prime[j]*i<N;j++) 51 { 52 int x=i*prime[j]; 53 vis[x]=true; 54 if(i%prime[j]==0) 55 { 56 miu[x]=0; 57 break; 58 } 59 miu[x]=-miu[i]; 60 } 61 } 62 for(int i=1;i<N;i++) 63 { 64 for(int j=i;j<N;j+=i) 65 F[j].a+=i; 66 F[i].b=i; 67 } 68 return ; 69 } 70 bool cmp(qust x,qust y) 71 { 72 return x.a<y.a; 73 } 74 bool cmq(qust x,qust y) 75 { 76 return x.no<y.no; 77 } 78 int main() 79 { 80 gtp(); 81 int T; 82 scanf("%d",&T); 83 for(int i=1;i<=T;i++) 84 { 85 scanf("%d%d%d",&q[i].n,&q[i].m,&q[i].a); 86 q[i].no=i; 87 } 88 std::sort(q+1,q+T+1,cmp); 89 std::sort(F+1,F+N); 90 for(int i=1,j=1;i<=T;i++) 91 { 92 for(;j<N&&F[j].a<=q[i].a;j++) 93 { 94 for(int k=F[j].b;k<N;k+=F[j].b) 95 update(k,F[j].a*miu[k/F[j].b]); 96 } 97 int n=q[i].n,m=q[i].m; 98 if(n>m) 99 std::swap(n,m); 100 for(int u=1,v;u<=n;u=v+1) 101 { 102 v=std::min(n/(n/u),m/(m/u)); 103 q[i].ans+=(query(v)-query(u-1))*(n/u)*(m/u); 104 } 105 } 106 std::sort(q+1,q+T+1,cmq); 107 for(int i=1;i<=T;i++) 108 printf("%d\n",q[i].ans&(0x7fffffff)); 109 return 0; 110 }