一些gcd计数问题

2022-06-24 13:33:33

数论什么的全都忘光了吧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)的贡献就可以了，存一个μ的前缀和就行了。

 #include <iostream>

 #include <cstdio>

 #include <cstring>

 #include <cmath>

 #include <algorithm>

 #define MAXN 50005

 using namespace std;

 int T, a, b, d;

 int prim[MAXN], pp, isprim[MAXN], miu[MAXN];

 int main(){

     scanf("%d", &T);

     miu[] = ;

     for(int i = ; i <= ; i ++){

         if(!isprim[i]) prim[++ pp] = i, miu[i] = -;

         for(int j = ; j <= pp && i*prim[j] <= ; j ++){

             isprim[i*prim[j]] = ;

             if(i%prim[j] == ) break;

             miu[i*prim[j]] = -miu[i];

         }

     }

     for(int i = ; i <= ; i ++) miu[i] += miu[i-];

     while(T --){

         scanf("%d%d%d", &a, &b, &d); a /= d, b /= d;

         int now = , aa = a, bb = b, na = , nb = , ans = ;

         while(){

             if(na < nb){

                 ans += (miu[na]-miu[now-]) * aa * bb;

                 now = na+; if(now > a) break;

                 aa = a/now;

                 na = a/aa;

             }else{

                 ans += (miu[nb]-miu[now-]) * aa * bb;

                 if(na == nb){

                     now = na+; if(now > a) break;

                     aa = a/now;

                     na = a/aa;

                 }

                 now = nb+; if(now > b) break;

                 bb = b/now;

                 nb = b/bb;

             }

         }

         printf("%d\n", ans);

     }

 //    system("pause");

     return ;

 }

bzoj2301: [HAOI2011]Problem b]

求有多少对数满足 gcd(x,y)=k, a<=x<=b, c<=y<=d

上一题然后随便容斥一下，，没什么好说的

 #include <iostream>

 #include <cstdio>

 #include <cstring>

 #include <cmath>

 #include <algorithm>

 #define MAXN 50005

 using namespace std;

 int T, a, b, c, d, k;

 int prim[MAXN], pp, isprim[MAXN], miu[MAXN];

 int calc(int a, int b){ if(a <=  || b <= ) return ;

     int now = , aa = a, bb = b, na = , nb = , ans = ;

         while(){

             if(na < nb){

                 ans += (miu[na]-miu[now-]) * aa * bb;

                 now = na+; if(now > a) break;

                 aa = a/now;

                 na = a/aa;

             }else{

                 ans += (miu[nb]-miu[now-]) * aa * bb;

                 if(na == nb){

                     now = na+; if(now > a) break;

                     aa = a/now;

                     na = a/aa;

                 }

                 now = nb+; if(now > b) break;

                 bb = b/now;

                 nb = b/bb;

             }

         }

     return ans;

 }

 int main(){

     scanf("%d", &T);

     miu[] = ;

     for(int i = ; i <= ; i ++){

         if(!isprim[i]) prim[++ pp] = i, miu[i] = -;

         for(int j = ; j <= pp && i*prim[j] <= ; j ++){

             isprim[i*prim[j]] = ;

             if(i%prim[j] == ) break;

             miu[i*prim[j]] = -miu[i];

         }

     }

     for(int i = ; i <= ; i ++) miu[i] += miu[i-];

     while(T --){

         scanf("%d%d%d%d%d", &a, &b ,&c, &d, &k); a --; c --; a /= k, b /= k, c /= k, d /= k;

         printf("%d\n", calc(b, d) - calc(a, d) - calc(c, b) + calc(a, c));

     }

     return ;

 }

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)，随便分块做一下就可以啦！

 #include <iostream>

 #include <cstring>

 #include <cmath>

 #include <algorithm>

 #include <cstdio>

 #define MAXN 10000005

 using namespace std;

 int T, n, m;

 int prim[MAXN], pp, miu[MAXN];

 long long G[MAXN];

 bool isprim[MAXN];

 int main(){

     scanf("%d", &T);

     miu[] = ;

     for(int i = ; i <= ; i ++){

         if(!isprim[i]) prim[++ pp] = i, miu[i] = -, G[i] = ;

         for(int j = ; j <= pp && i * prim[j] <= ; j ++){

             isprim[i*prim[j]] = ;

             if(i%prim[j] == ){G[i*prim[j]] = miu[i]; break;}

             miu[i*prim[j]] = -miu[i]; G[i*prim[j]] = miu[i]-G[i];

         }

     }

     for(int i = ; i <= ; i ++) G[i] += G[i-];

     while(T --){

         scanf("%d%d", &n, &m);

         int now = , nn = n, mm = m, nen = n/nn, nem = m/mm; long long ans = ;

         while(){

             if(nen < nem){

                 ans += (G[nen]-G[now-])*nn*mm;

                 now = nen+; if(now > n) break;

                 nn = n/now;

                 nen = n/nn;

             }else{

                 ans += (G[nem]-G[now-])*nn*mm;

                 if(nen == nem){

                     now = nen + ; if(now > n) break;

                     nn = n/now;

                     nen = n/nn;

                 }

                 now = nem + ; if(now > m) break;

                 mm = m/now;

                 nem = m/mm;

             }

         }

         printf("%lld\n", ans);

     }

     return ;

 }

[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，，，，（好无聊

 #include <iostream>

 #include <cstdio>

 #include <cmath>

 #include <algorithm>

 #include <cstring>

 #define ll long long

 using namespace std;

 int m, N, bin[], cbin[], cnt;

 ll ans;

 void dfs(int x, int s, int phin){

     if(x == cnt+){

         ans += phin*(m/s); return;

     }

     dfs(x+, s, phin);

     s*=bin[x],phin*=(bin[x]-);

     for(int i = ; i <= cbin[x]; i ++){

         dfs(x+, s, phin);

         s *= bin[x]; phin *= bin[x];

     }

 }

 int main(){

     cin >> N; m = N;

     for(int i = ; i*i <= N; i ++)if(N%i==){

         bin[++cnt] = i;

         while(N%i==)cbin[cnt] ++, N/=i;

     }

     if(N)bin[++cnt]=N, cbin[cnt] = ;

     dfs(, , );

     printf("%lld\n", ans);

     return ;

 }

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

 #include <iostream>

 #include <cstdio>

 #include <cmath>

 #include <cstring>

 #include <algorithm>

 #define MAXN 100005

 #define lowbit(x) (x&(-x))

 using namespace std;

 int T, N, prim[MAXN], pp, isprim[MAXN], miu[MAXN], cm[MAXN], fs[MAXN], A[MAXN], ans[MAXN];

 inline void Add(int x, int c){

     for(; x <= N; x += lowbit(x)) A[x] += c;

 }

 inline int Sum(int x){

     int ret = ;

     for(; x; x -= lowbit(x)) ret += A[x];

     return ret;

 }

 struct FF{int f, x;}F[MAXN];

 struct point{int n, m, a, num;}p[];

 bool cmp(point x, point y){return x.a < y.a;}

 bool cmpf(FF a, FF b){return a.f < b.f;}

 inline void Insert(int x, int c){

     for(int i = ; i*x <= N; i ++) Add(i*x, c*miu[i]);

 }

 int calc(int a, int b){

     int now = , aa = a, bb = b, na = , nb = , ans = , pre = , la;

         while(){

             if(na < nb){

                 la = Sum(na);

                 ans += (la-pre) * aa * bb;

                 pre = la;

                 now = na+; if(now > a) break;

                 aa = a/now;

                 na = a/aa;

             }else{

                 la = Sum(nb);

                 ans += (la-pre) * aa * bb;

                 pre = la;

                 if(na == nb){

                     now = na+; if(now > a) break;

                     aa = a/now;

                     na = a/aa;

                 }

                 now = nb+; if(now > b) break;

                 bb = b/now;

                 nb = b/bb;

             }

         }

     return ans;

 }

 int main(){

     scanf("%d", &T);

     miu[] = F[].f = F[].x = ;

     for(int i = ; i <= ; i ++){

         if(!isprim[i]) prim[++ pp] = i, miu[i] = -, F[i].f = i+, cm[i] = i, fs[i] = ;

         for(int j = ; j <= pp && i*prim[j] <= ; j ++){

             isprim[i*prim[j]] = ;

             if(i%prim[j] == ){

                 miu[i*prim[j]] = ;

                 F[i*prim[j]].f = F[i].f + fs[i]*cm[i]*prim[j];

                 cm[i*prim[j]] = cm[i]*prim[j];

                 fs[i*prim[j]] = fs[i];

                 break;

             }

             miu[i*prim[j]] = -miu[i];

             F[i*prim[j]].f = F[i].f*(+prim[j]);

             cm[i*prim[j]] = prim[j], fs[i*prim[j]] = F[i].f;

         }

         F[i].x = i;

     }

     for(int i = ; 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));

     sort(p + , p + T + , cmp);

     sort(F + , F + +, cmpf);

     int jj = ;

     for(int i = ; i <= T; i ++){

         while(jj <=  && F[jj].f <= p[i].a){Insert(F[jj].x, F[jj].f); jj ++;}

         ans[p[i].num] = calc(p[i].n, p[i].m)&((<<)-);

     }

     for(int i = ; i <= T; i ++) printf("%d\n", ans[i]);

     return ;

 }

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)*d²}*p }

然后分块套分块做一遍就好了 O(sqrt(n)*sqrt(n)) = O(n)

 #include <iostream>

 #include <cstdio>

 #include <cstring>

 #include <cmath>

 #include <algorithm>

 #define ll long long

 #define mod 20101009

 #define MAXN 10000005

 using namespace std;

 int N, n, m, prim[MAXN], isprim[MAXN], pp;

 ll sum[MAXN], miu[MAXN];

 ll SS(ll x, ll y){

     return (((x*(x+)/)%mod) * ((y*(y+)/)%mod) )%mod;

 }

 ll F(int a, int b){

     int now = , aa = a, bb = b, na = , nb = ;

     ll ans = ;

     while(){

         if(na < nb){

             (ans += SS(a/now, b/now) * ((miu[na]-miu[now-]+mod)%mod) %mod )%=mod;

             now = na+; if(now > a) break;

             aa = a/now;

             na = a/aa;

         }else{

             (ans += SS(a/now, b/now) * ((miu[nb]-miu[now-]+mod)%mod) %mod ) %= mod;

             if(na == nb){

                 now = na+; if(now > a) break;

                 aa = a/now;

                 na = a/aa;

             }

             now = nb+; if(now > b) break;

             bb = b/now;

             nb = b/bb;

         }

     }

     return ans;

 }

 ll calc(int a, int b){

     int now = , aa = a, bb = b, na = , nb = ;

     ll ans = ;

     while(){

         if(na < nb){

             (ans += F(a/now, b/now) * ((sum[na]-sum[now-]+mod)%mod) %mod ) %=mod;

             now = na+; if(now > a) break;

             aa = a/now;

             na = a/aa;

         }else{

             (ans += F(a/now, b/now) * ((sum[nb]-sum[now-]+mod)%mod) %mod )%=mod;

             if(na == nb){

                 now = na+; if(now > a) break;

                 aa = a/now;

                 na = a/aa;

             }

             now = nb+; if(now > b) break;

             bb = b/now;

             nb = b/bb;

         }

     }

     return ans;

 }

 int main(){

     scanf("%d%d", &n, &m); N = min(n, m);

     miu[] = ;

     for(int i = ; i <= N; i ++){

         if(!isprim[i]) prim[++ pp] = i, miu[i] = -;

         for(int j = ; j <= pp && i*prim[j]<=N; j ++){

             isprim[i*prim[j]] = ;

             if(i%prim[j] == ) break;

             miu[i*prim[j]] = -miu[i];

         }

     }

     for(int i = ; i <= N; i ++) sum[i] = (sum[i-]+i)%mod, miu[i] = (miu[i-]+miu[i]*i*i)%mod;

     cout << calc(n, m)<<endl;

     return ;

 }

2693: jzptab

上一题的加强版， 1000组数据。

又是这样！窝打了很多东西后电脑崩了都没存下来QAQ

思路不难，考虑把上一题的那个式子继续化简，像以前的方法一样通过把里面的那个多项式都乘上一个p或是用换元把pd换成S（其实这两种方法本质上是一样的），然后就可以成功地把[n/p]这个非常棒的可以sqrt(n)搞出来的东西放到外层多项式了，现在只要考虑如何O(n)地预处理出内层多项式。

G(d) = sigema(i|d) {d/i * i² * μ(i)} = d* sigema(i|d) {d*μ(d)}

用莫比乌斯反演推一下容易得出积性函数的因数和也是积性函数，所以G(d)是积性函数，线性筛一遍随便搞一搞就可以了

 #include <iostream>

 #include <cstdio>

 #include <cstring>

 #include <cmath>

 #include <algorithm>

 #define mod 100000009

 #define MAXN 10000005

 #define ll long long

 using namespace std;

 int T, N, M, prim[MAXN], pp, miu[MAXN], isprim[MAXN];

 ll G[MAXN];

 ll SS(ll a, ll b){return (((a*(a+)/)%mod)*((b*(b+)/)%mod))%mod;}

 int calc(int a, int b){

     int now = , aa = a, bb = b, na = , nb = , ans = ;

         while(){

             if(na < nb){

                 (ans += ((G[na]-G[now-]) * SS(aa, bb))%mod)%=mod;

                 now = na+; if(now > a) break;

                 aa = a/now;

                 na = a/aa;

             }else{

                 (ans += ((G[nb]-G[now-]) * SS(aa, bb))%mod)%=mod;

                 if(na == nb){

                     now = na+; if(now > a) break;

                     aa = a/now;

                     na = a/aa;

                 }

                 now = nb+; if(now > b) break;

                 bb = b/now;

                 nb = b/bb;

             }

         }

     return ans;

 }

 int main(){

     miu[] = ; G[] = ;

     for(int i = ; i <= ; i ++){

         if(!isprim[i]) prim[++ pp] = i, miu[i] = -, G[i] = (-i+mod)%mod;

         for(int j = ; j <= pp && i*prim[j] <= ; j ++){

             isprim[i*prim[j]] = ;

             if(i%prim[j] == ){

                 G[i*prim[j]] = G[i]; break;

             }

             miu[i*prim[j]] = -miu[i];

             G[i*prim[j]] = (G[i]*G[prim[j]])%mod;

         }

     }

     for(int i = ; i <= ; i ++) G[i] = (G[i-] + G[i]*i)%mod;

     scanf("%d", &T);

     while(T--){

         int a, b; scanf("%d%d", &a, &b);

         printf("%d\n", (calc(a, b)+mod)%mod);

     }

     return ;

 }

总结

上面的这些题其实都是一个问题换一下形式而已，可以起名叫做“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)地计算被换到内层的部分的前缀和就可以了。

码农公寓