https://www.cnblogs.com/cmy-blog/p/binary-gcd.html)
Binary GCD
给定 $a_1,a_2,\cdots,a_n,b_1,b_2,\cdots,b_n$，对于 $i\in[1,n]$，求 $\displaystyle A_i=\sum_{j=1}^{n}i^j\gcd\left(a_i,b_j\right)$ 对 $998244353$ 取模。
$1\leq n\leq5000$，$1\leq a_i,b_i\leq10^6$。时间限制 1.2s。
以下讨论对于 $\gcd(a,b)$，均认为 $a,b\geq0$。
考虑 $\mathcal O\left(n^2\log V\right)$ 好像不太能过，又不会预处理，因此乱搞。辗转相除虽然是 $\mathcal O(\log V)$ 的，但是除法太慢了，复杂度优化不了就优化常数。
考虑优化辗转相减。
以下记 $\mid$ 表示位运算 |，$\gg$ 表示位运算 >>，$\ll$ 表示位运算 <<。其运算优先级均等同于 C++ 中的运算优先级。
对于 $\gcd(a,b)$：
-
$a=b$，$\gcd(a,b)=a$。
-
$a=0$ 或 $b=0$，则 $\gcd(a,b)=a\mid b$。
-
$a,b$ 均为偶数，则 $\gcd(a,b)=\gcd(a\gg1,b\gg1)\ll1$。
-
$a$ 为偶数，$b$ 为奇数，则 $\gcd(a,b)=\gcd(a\gg1,b)$。
$a$ 为奇数，$b$ 为偶数，则 $\gcd(a,b)=\gcd(a,b\gg1)$。
-
$a,b$ 均为奇数。
考虑 $\gcd(a,b)=\gcd(b,a-b)$，可以发现 $a-b$ 为偶数，因此 $\gcd(a-b\gg1,b)$。
注意需要钦定 $a>b$，否则会出现意想不到的错误（位运算会炸掉）。
记录一下 $\ll1$ 的次数 $p$，就可以得到：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
int gcd(int a,int b){
int p=0;
while(true){
if(a==b){
return a<<p;
}
if(!a||!b){
return (a|b)<<p;
}
if(~a&1){
if(~b&1){
a>>=1;b>>=1;
p++;
}else{
a>>=1;
}
}else if(~b&1){
b>>=1;
}else{
if(a<b){
swap(a,b);
}
a=a-b>>1;
}
}
}
但是这还是不够快，还需要更快。
每次一个一个去除后缀 $0$ 显然是缓慢的，可以利用 GCC 内建函数 __builtin_ctz 快速统计。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
inline int gcd(int a, int b){
if(!a||!b){
return a|b;
}
int az=__builtin_ctz(a),bz=__builtin_ctz(b);
int z=min(az,bz);
b>>=bz;
while(a){
a>>=az;
int diff=abs(b-a);
az=__builtin_ctz(diff);
b=min(a,b);
a=diff;
}
return b<<z;
}
如果调用 $\gcd(a,b)$ 时，保证 $a,b\neq 0$，可以不要 if 判断 $a=0$ 或 $b=0$。
AC 代码
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
//#include<bits/stdc++.h>
#include<algorithm>
#include<iostream>
#include<cstring>
#include<iomanip>
#include<cstdio>
#include<string>
#include<vector>
#include<cmath>
#include<ctime>
#include<deque>
#include<queue>
#include<stack>
#include<list>
using namespace std;
constexpr const int N=5000,P=998244353;
int n,a[N+1],b[N+1];
inline int gcd(int a, int b){
int az=__builtin_ctz(a),bz=__builtin_ctz(b);
int z=min(az,bz);
b>>=bz;
while(a){
a>>=az;
int diff=abs(b-a);
az=__builtin_ctz(diff);
b=min(a,b);
a=diff;
}
return b<<z;
}
int main(){
/*freopen("test.in","r",stdin);
freopen("test.out","w",stdout);*/
ios::sync_with_stdio(false);
cin.tie(0);cout.tie(0);
cin>>n;
for(int i=1;i<=n;i++){
cin>>a[i];
}
for(int i=1;i<=n;i++){
cin>>b[i];
}
for(int i=1;i<=n;i++){
int ans=0,p=1;
for(int j=1;j<=n;j++){
p=1ll*p*i%P;
ans=(ans+1ll*p*gcd(a[i],b[j])%P)%P;
}
cout<<ans<<'\n';
}
cout.flush();
/*fclose(stdin);
fclose(stdout);*/
return 0;
}
