Binary GCD

https://www.cnblogs.com/cmy-blog/p/binary-gcd.html)

Binary GCD

Luogu P5435 基于值域预处理的快速 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\)，就可以得到：

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 快速统计。

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 代码

//#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;
}