数论什么的全都忘光了吧QAQ 做了几道简单的题练习一下。
bzoj1101: [POI2007]Zap
求有多少对数满足 gcd(x,y)=d, 1<=x<=a, 1<=y<=b
首先那个d肯定是要除下去的， 变成 gcd(x,y) = 1, 1<=x<=a, 1<=y<=b这样一个问题。
这样就变成了最经典的莫比乌斯反演， 设f(d)为 gcd(x,y) = d, 1<=x<=a, 1<=y<=b时满足要求的x,y的个数，F(d) = Σ(d|n)f(n) = [a/d] * [b/d]， 那么f(1) = Σ μ(d) * [a/d] * [b/d]
复杂度是O(n)的， 但是多组询问会T，这时候想到了bzoj1257 CQOI2007]余数之和sum的思想：对于一个N， [n/d]的取值是sqrt(n)级别的（显然当d <= sqrt(n)，只有sqrt(n)个d所以有sqrt(n)中取值， 当d > sqrt(n), n/d < sqrt(n)所以只有sqrt(n)中取值）， 所以只要在 a/d或是b/d中有任何一个变化了以后再计算对f(1)的贡献就可以了，存一个μ的前缀和就行了。
1 #include <iostream> 2 #include <cstdio> 3 #include <cstring> 4 #include <cmath> 5 #include <algorithm> 6 #define MAXN 50005 7 using namespace std; 8 int T, a, b, d; 9 int prim[MAXN], pp, isprim[MAXN], miu[MAXN]; 10 int main(){ 11 scanf("%d", &T); 12 miu[1] = 1; 13 for(int i = 2; i <= 50000; i ++){ 14 if(!isprim[i]) prim[++ pp] = i, miu[i] = -1; 15 for(int j = 1; j <= pp && i*prim[j] <= 50000; j ++){ 16 isprim[i*prim[j]] = 1; 17 if(i%prim[j] == 0) break; 18 miu[i*prim[j]] = -miu[i]; 19 } 20 } 21 for(int i = 2; i <= 50000; i ++) miu[i] += miu[i-1]; 22 while(T --){ 23 scanf("%d%d%d", &a, &b, &d); a /= d, b /= d; 24 int now = 1, aa = a, bb = b, na = 1, nb = 1, ans = 0; 25 while(1){ 26 if(na < nb){ 27 ans += (miu[na]-miu[now-1]) * aa * bb; 28 now = na+1; if(now > a) break; 29 aa = a/now; 30 na = a/aa; 31 }else{ 32 ans += (miu[nb]-miu[now-1]) * aa * bb; 33 if(na == nb){ 34 now = na+1; if(now > a) break; 35 aa = a/now; 36 na = a/aa; 37 } 38 now = nb+1; if(now > b) break; 39 bb = b/now; 40 nb = b/bb; 41 } 42 } 43 printf("%d ", ans); 44 } 45 // system("pause"); 46 return 0; 47 }
bzoj2301: [HAOI2011]Problem b]
求有多少对数满足 gcd(x,y)=k, a<=x<=b, c<=y<=d
上一题然后随便容斥一下，，没什么好说的
1 #include <iostream> 2 #include <cstdio> 3 #include <cstring> 4 #include <cmath> 5 #include <algorithm> 6 #define MAXN 50005 7 using namespace std; 8 int T, a, b, c, d, k; 9 int prim[MAXN], pp, isprim[MAXN], miu[MAXN]; 10 int calc(int a, int b){ if(a <= 0 || b <= 0) return 0; 11 int now = 1, aa = a, bb = b, na = 1, nb = 1, ans = 0; 12 while(1){ 13 if(na < nb){ 14 ans += (miu[na]-miu[now-1]) * aa * bb; 15 now = na+1; if(now > a) break; 16 aa = a/now; 17 na = a/aa; 18 }else{ 19 ans += (miu[nb]-miu[now-1]) * aa * bb; 20 if(na == nb){ 21 now = na+1; if(now > a) break; 22 aa = a/now; 23 na = a/aa; 24 } 25 now = nb+1; if(now > b) break; 26 bb = b/now; 27 nb = b/bb; 28 } 29 } 30 return ans; 31 } 32 int main(){ 33 scanf("%d", &T); 34 miu[1] = 1; 35 for(int i = 2; i <= 50000; i ++){ 36 if(!isprim[i]) prim[++ pp] = i, miu[i] = -1; 37 for(int j = 1; j <= pp && i*prim[j] <= 50000; j ++){ 38 isprim[i*prim[j]] = 1; 39 if(i%prim[j] == 0) break; 40 miu[i*prim[j]] = -miu[i]; 41 } 42 } 43 for(int i = 2; i <= 50000; i ++) miu[i] += miu[i-1]; 44 while(T --){ 45 scanf("%d%d%d%d%d", &a, &b ,&c, &d, &k); a --; c --; a /= k, b /= k, c /= k, d /= k; 46 printf("%d ", calc(b, d) - calc(a, d) - calc(c, b) + calc(a, c)); 47 } 48 return 0; 49 }
bzoj2820: YY的GCD
求1<=x<=N, 1<=y<=M且gcd(x,y)为质数有多少对
天啊窝一定是世界上最后一个知道 [n/(x*y)] = [[n/x]/y] 的人！！ QAQ
很显然可以直接枚举一个质数然后然后想原来的题那样做，但是T*600000*n的复杂度肯定无法接受。
ans = ∑(p为质数)∑(d) μ(d)*[n/pd]*[m/pd]
容易想到像之前那样把n和m分块处理一下， 可以做到T*600000*sqrt(n)， 但显然还是无法接受QAQ
觉得p为质数这种东西放在外层循环肯定就必须枚举了， 但是如果放在内层可能有奇怪的事情发生？
ans = ∑(d)∑(p为质数) μ(d)*[n/pd]*[m/pd]
感觉很靠谱的样子！但是这里n除了除以一个p以外还除以了一个d， 让内层无法分块了， 所以考虑设一个数s = pd， 那么
ans = ∑(s)[n/s][m/s]∑(p为质数)μ(s/p)
随便算一算发现G(s) = ∑(p为质数)μ(s/p) 的这个函数是可以在线性筛的过程中算出来的！
所以问题转换为 ans = ∑(s)[n/s][m/s]G(s)， 随便分块做一下就可以啦！
1 #include <iostream> 2 #include <cstring> 3 #include <cmath> 4 #include <algorithm> 5 #include <cstdio> 6 #define MAXN 10000005 7 using namespace std; 8 int T, n, m; 9 int prim[MAXN], pp, miu[MAXN]; 10 long long G[MAXN]; 11 bool isprim[MAXN]; 12 int main(){ 13 scanf("%d", &T); 14 miu[1] = 1; 15 for(int i = 2; i <= 10000000; i ++){ 16 if(!isprim[i]) prim[++ pp] = i, miu[i] = -1, G[i] = 1; 17 for(int j = 1; j <= pp && i * prim[j] <= 10000000; j ++){ 18 isprim[i*prim[j]] = 1; 19 if(i%prim[j] == 0){G[i*prim[j]] = miu[i]; break;} 20 miu[i*prim[j]] = -miu[i]; G[i*prim[j]] = miu[i]-G[i]; 21 } 22 } 23 for(int i = 2; i <= 10000000; i ++) G[i] += G[i-1]; 24 while(T --){ 25 scanf("%d%d", &n, &m); 26 int now = 1, nn = n, mm = m, nen = n/nn, nem = m/mm; long long ans = 0; 27 while(1){ 28 if(nen < nem){ 29 ans += (G[nen]-G[now-1])*nn*mm; 30 now = nen+1; if(now > n) break; 31 nn = n/now; 32 nen = n/nn; 33 }else{ 34 ans += (G[nem]-G[now-1])*nn*mm; 35 if(nen == nem){ 36 now = nen + 1; if(now > n) break; 37 nn = n/now; 38 nen = n/nn; 39 } 40 now = nem + 1; if(now > m) break; 41 mm = m/now; 42 nem = m/mm; 43 } 44 } 45 printf("%lld ", ans); 46 } 47 return 0; 48 }
[bzoj2705: [SDOI2012]Longge的问题]
求 sigma(gcd(i,n), 1<=i<=n<2^32)
很简单的题。ans = ∑(p|n) phi(n/p)*p
显然无法线性筛预处理phi, 那么就用 phi(n) = n * ∏((pi-1)/pi)这个式子， dfs一遍就好啦。
对了这道题虽然它说 n <= 1<<32 但实际上没有n>=1<<31的情况诶！除了ans以外全部开int有机会从8ms降到4ms，，，，（好无聊
1 #include <iostream> 2 #include <cstdio> 3 #include <cmath> 4 #include <algorithm> 5 #include <cstring> 6 #define ll long long 7 using namespace std; 8 int m, N, bin[33], cbin[33], cnt; 9 ll ans; 10 void dfs(int x, int s, int phin){ 11 if(x == cnt+1){ 12 ans += phin*(m/s); return; 13 } 14 dfs(x+1, s, phin); 15 s*=bin[x],phin*=(bin[x]-1); 16 for(int i = 1; i <= cbin[x]; i ++){ 17 dfs(x+1, s, phin); 18 s *= bin[x]; phin *= bin[x]; 19 } 20 } 21 int main(){ 22 cin >> N; m = N; 23 for(int i = 2; i*i <= N; i ++)if(N%i==0){ 24 bin[++cnt] = i; 25 while(N%i==0)cbin[cnt] ++, N/=i; 26 } 27 if(N)bin[++cnt]=N, cbin[cnt] = 1; 28 dfs(1, 1, 1); 29 printf("%lld ", ans); 30 return 0; 31 }
3529: [Sdoi2014]数表
有一张N×m的数表，其第i行第j列（1 < =i < = n，1 < =j < =m）的数值为能同时整除i和j的所有自然数之和。给定a，计算数表中不大于a的数之和。
感觉这题挺神的。
先不考虑a， 如何计算数表中所有数之和呢？
设F(x) 为能整除x的所有自然数之和， 这个可以线性筛时胡乱搞搞算出来。
ans = sigema{F(gcd(i,j)) }
= sigema(p) { sigema(d){miu(d)*[n/dp]*[m/dp]} * F(p)}
然后像bzoj2820: YY的GCD 那样， 试图把可以sqrt(n)计算的[n/dp]什么的挪到外层以降低复杂度
[n/dp] 和[m/dp] 这两项以外包含了/p的， 只要乘上p以后就可以挪到外边了
ans = sigema(p){ sigema(p|d){miu(d/p)*[n/d]*[m/d]} * F(p)}
= sigema(d){ [n/d]*[m/d] * sigema(p|d){miu(d/p) * F(p)} }
好啦！现在只需要能够预处理出来 G(x) = sigema(p|x){miu(x/p)*F(p)} 就可以分块前缀和啦！
因为约数的期望是O(logn)的所以只要暴力一下就可以做到O(nlogn)的复杂度啦！
现在考虑原问题。。。。。。。。。感觉a这个条件都要被遗忘了。
如果F(x) > a的话， 直接让F(x) = 0 就可以了！ 所以只要把询问离线， 按a排一个序， 再把F(x)排一个序， 一个一个加进去， 加的时候用树状数组维护前缀和就可以了！好机智啊！
复杂度 O(nlog^2n + q*sqrt(n)*logn)
话说我知道现在还会在取模意义下做减法的时候不再加上一个mod ，，，，，，被自己蠢哭啦QAQ
1 #include <iostream> 2 #include <cstdio> 3 #include <cmath> 4 #include <cstring> 5 #include <algorithm> 6 #define MAXN 100005 7 #define lowbit(x) (x&(-x)) 8 using namespace std; 9 int T, N, prim[MAXN], pp, isprim[MAXN], miu[MAXN], cm[MAXN], fs[MAXN], A[MAXN], ans[MAXN]; 10 inline void Add(int x, int c){ 11 for(; x <= N; x += lowbit(x)) A[x] += c; 12 } 13 inline int Sum(int x){ 14 int ret = 0; 15 for(; x; x -= lowbit(x)) ret += A[x]; 16 return ret; 17 } 18 struct FF{int f, x;}F[MAXN]; 19 struct point{int n, m, a, num;}p[20005]; 20 bool cmp(point x, point y){return x.a < y.a;} 21 bool cmpf(FF a, FF b){return a.f < b.f;} 22 inline void Insert(int x, int c){ 23 for(int i = 1; i*x <= N; i ++) Add(i*x, c*miu[i]); 24 } 25 int calc(int a, int b){ 26 int now = 1, aa = a, bb = b, na = 1, nb = 1, ans = 0, pre = 0, la; 27 while(1){ 28 if(na < nb){ 29 la = Sum(na); 30 ans += (la-pre) * aa * bb; 31 pre = la; 32 now = na+1; if(now > a) break; 33 aa = a/now; 34 na = a/aa; 35 }else{ 36 la = Sum(nb); 37 ans += (la-pre) * aa * bb; 38 pre = la; 39 if(na == nb){ 40 now = na+1; if(now > a) break; 41 aa = a/now; 42 na = a/aa; 43 } 44 now = nb+1; if(now > b) break; 45 bb = b/now; 46 nb = b/bb; 47 } 48 } 49 return ans; 50 } 51 int main(){ 52 scanf("%d", &T); 53 miu[1] = F[1].f = F[1].x = 1; 54 for(int i = 2; i <= 100000; i ++){ 55 if(!isprim[i]) prim[++ pp] = i, miu[i] = -1, F[i].f = i+1, cm[i] = i, fs[i] = 1; 56 for(int j = 1; j <= pp && i*prim[j] <= 100000; j ++){ 57 isprim[i*prim[j]] = 1; 58 if(i%prim[j] == 0){ 59 miu[i*prim[j]] = 0; 60 F[i*prim[j]].f = F[i].f + fs[i]*cm[i]*prim[j]; 61 cm[i*prim[j]] = cm[i]*prim[j]; 62 fs[i*prim[j]] = fs[i]; 63 break; 64 } 65 miu[i*prim[j]] = -miu[i]; 66 F[i*prim[j]].f = F[i].f*(1+prim[j]); 67 cm[i*prim[j]] = prim[j], fs[i*prim[j]] = F[i].f; 68 } 69 F[i].x = i; 70 } 71 for(int i = 1; i <= T; i ++) scanf("%d%d%d", &p[i].n, &p[i].m, &p[i].a), p[i].num = i, N = max(N, max(p[i].n,p[i].m)); 72 sort(p + 1, p + T + 1, cmp); 73 sort(F + 1, F + 100000+1, cmpf); 74 int jj = 1; 75 for(int i = 1; i <= T; i ++){ 76 while(jj <= 100000 && F[jj].f <= p[i].a){Insert(F[jj].x, F[jj].f); jj ++;} 77 ans[p[i].num] = calc(p[i].n, p[i].m)&((1<<31)-1); 78 } 79 for(int i = 1; i <= T; i ++) printf("%d ", ans[i]); 80 return 0; 81 }
BZOJ 2154 Crash的数字表格
求Σ[1<=i<=n]Σ[1<=j<=m]LCM(i,j) mod 20101009
这题没什么技巧，直接倒一下就可以了。
ans = sigema(p) {i*j/p ((i,j)=p) }
设S(x) = x*(x-1)/2
ans = sigema(p) { sigema(d){miu(d)*S(n/pd)*S(m/pd)*d2}*p }
然后分块套分块做一遍就好了 O(sqrt(n)*sqrt(n)) = O(n)
1 #include <iostream> 2 #include <cstdio> 3 #include <cstring> 4 #include <cmath> 5 #include <algorithm> 6 #define ll long long 7 #define mod 20101009 8 #define MAXN 10000005 9 using namespace std; 10 int N, n, m, prim[MAXN], isprim[MAXN], pp; 11 ll sum[MAXN], miu[MAXN]; 12 ll SS(ll x, ll y){ 13 return (((x*(x+1)/2)%mod) * ((y*(y+1)/2)%mod) )%mod; 14 } 15 ll F(int a, int b){ 16 int now = 1, aa = a, bb = b, na = 1, nb = 1; 17 ll ans = 0; 18 while(1){ 19 if(na < nb){ 20 (ans += SS(a/now, b/now) * ((miu[na]-miu[now-1]+mod)%mod) %mod )%=mod; 21 now = na+1; if(now > a) break; 22 aa = a/now; 23 na = a/aa; 24 }else{ 25 (ans += SS(a/now, b/now) * ((miu[nb]-miu[now-1]+mod)%mod) %mod ) %= mod; 26 if(na == nb){ 27 now = na+1; if(now > a) break; 28 aa = a/now; 29 na = a/aa; 30 } 31 now = nb+1; if(now > b) break; 32 bb = b/now; 33 nb = b/bb; 34 } 35 } 36 return ans; 37 } 38 ll calc(int a, int b){ 39 int now = 1, aa = a, bb = b, na = 1, nb = 1; 40 ll ans = 0; 41 while(1){ 42 if(na < nb){ 43 (ans += F(a/now, b/now) * ((sum[na]-sum[now-1]+mod)%mod) %mod ) %=mod; 44 now = na+1; if(now > a) break; 45 aa = a/now; 46 na = a/aa; 47 }else{ 48 (ans += F(a/now, b/now) * ((sum[nb]-sum[now-1]+mod)%mod) %mod )%=mod; 49 if(na == nb){ 50 now = na+1; if(now > a) break; 51 aa = a/now; 52 na = a/aa; 53 } 54 now = nb+1; if(now > b) break; 55 bb = b/now; 56 nb = b/bb; 57 } 58 } 59 return ans; 60 } 61 int main(){ 62 scanf("%d%d", &n, &m); N = min(n, m); 63 miu[1] = 1; 64 for(int i = 2; i <= N; i ++){ 65 if(!isprim[i]) prim[++ pp] = i, miu[i] = -1; 66 for(int j = 1; j <= pp && i*prim[j]<=N; j ++){ 67 isprim[i*prim[j]] = 1; 68 if(i%prim[j] == 0) break; 69 miu[i*prim[j]] = -miu[i]; 70 } 71 } 72 for(int i = 1; i <= N; i ++) sum[i] = (sum[i-1]+i)%mod, miu[i] = (miu[i-1]+miu[i]*i*i)%mod; 73 cout << calc(n, m)<<endl; 74 return 0; 75 }
2693: jzptab
上一题的加强版， 1000组数据。
又是这样！窝打了很多东西后电脑崩了都没存下来QAQ
思路不难，考虑把上一题的那个式子继续化简， 像以前的方法一样通过把里面的那个多项式都乘上一个p或是用换元把pd换成S（其实这两种方法本质上是一样的）， 然后就可以成功地把[n/p]这个非常棒的可以sqrt(n)搞出来的东西放到外层多项式了，现在只要考虑如何O(n)地预处理出内层多项式。
G(d) = sigema(i|d) {d/i * i2 * μ(i)} = d* sigema(i|d) {d*μ(d)}
用莫比乌斯反演推一下容易得出积性函数的因数和也是积性函数， 所以G(d)是积性函数， 线性筛一遍随便搞一搞就可以了
1 #include <iostream> 2 #include <cstdio> 3 #include <cstring> 4 #include <cmath> 5 #include <algorithm> 6 #define mod 100000009 7 #define MAXN 10000005 8 #define ll long long 9 using namespace std; 10 int T, N, M, prim[MAXN], pp, miu[MAXN], isprim[MAXN]; 11 ll G[MAXN]; 12 ll SS(ll a, ll b){return (((a*(a+1)/2)%mod)*((b*(b+1)/2)%mod))%mod;} 13 int calc(int a, int b){ 14 int now = 1, aa = a, bb = b, na = 1, nb = 1, ans = 0; 15 while(1){ 16 if(na < nb){ 17 (ans += ((G[na]-G[now-1]) * SS(aa, bb))%mod)%=mod; 18 now = na+1; if(now > a) break; 19 aa = a/now; 20 na = a/aa; 21 }else{ 22 (ans += ((G[nb]-G[now-1]) * SS(aa, bb))%mod)%=mod; 23 if(na == nb){ 24 now = na+1; if(now > a) break; 25 aa = a/now; 26 na = a/aa; 27 } 28 now = nb+1; if(now > b) break; 29 bb = b/now; 30 nb = b/bb; 31 } 32 } 33 return ans; 34 } 35 int main(){ 36 miu[1] = 1; G[1] = 1; 37 for(int i = 2; i <= 10000000; i ++){ 38 if(!isprim[i]) prim[++ pp] = i, miu[i] = -1, G[i] = (1-i+mod)%mod; 39 for(int j = 1; j <= pp && i*prim[j] <= 10000000; j ++){ 40 isprim[i*prim[j]] = 1; 41 if(i%prim[j] == 0){ 42 G[i*prim[j]] = G[i]; break; 43 } 44 miu[i*prim[j]] = -miu[i]; 45 G[i*prim[j]] = (G[i]*G[prim[j]])%mod; 46 } 47 } 48 for(int i = 1; i <= 10000000; i ++) G[i] = (G[i-1] + G[i]*i)%mod; 49 scanf("%d", &T); 50 while(T--){ 51 int a, b; scanf("%d%d", &a, &b); 52 printf("%d ", (calc(a, b)+mod)%mod); 53 } 54 return 0; 55 }
总结
上面的这些题其实都是一个问题换一下形式而已，可以起名叫做“gcd计数问题”？无非是看到要求一个gcd相关的东西， 然后不好做， 所以要莫比乌斯反演一下把(x,y)=d转化为(d|x 且 d|y)， 这就是一个直接用除法取整可以解决的问题了，而除法取整有一个非常好的性质就是不同的 f(x) = x/i 的数量是 sqrt(x)级别的， 如果可以把他提出来且把它里面的东西的前缀和预处理出来就可以O(sqrt(n))地解决这个问题。有的时候 [x/ij] 被套在了内层于是只能O(sqrt(n))地算出内层部分， 不能接受， 于是就要想办法把[x/ij]换到外层， 内层整体乘上j就可以了，然后再考虑如何O(n)地计算被换到内层的部分的前缀和就可以了。