「学习笔记」基于值域预处理的快速 GCD 算法

Description

给定 \(m\) 个数 \(a_1,a_2,\cdots,a_m\)，再给定 \(n\) 个数 \(b_1,b_2,\cdots,b_m\)。

试求出所有 \(\gcd(a_i,b_j)\) 的值，其中 \(1\le i,j\le m\)。

（已知 \(\forall 1\le i\le n,a_i,b_i\in[1,n]\)）

Analysis

设 \(x\in[1,n]\)，若可以将 \(x\) 分为三个数的乘积 \(x=a×b×c(a\le b\le c)\)，并且 \(a,b,c\) 或为质数，或不超过 \(\sqrt{n}\)。

则称有序三元组 \((a,b,c)\) 是 \(x\) 的一个分解。

如何求任意一个数 \(x\) 的分解呢？

考虑线筛，若 \(x\) 为质数，显然 \((1,1,x)\) 是它的一个分解；

若 \(x\) 为合数，设 \(p\) 是 \(x\) 的最小质因子， \((a_0,b_0,c_0)\) 是 \(x/p\) 的一个分解；

下证 \(a_0p,b_0,c_0\) 排序后是 \(x\) 的一个分解，

• 若 \(p\le\sqrt[4]{x}\)，结合 \(a_0\le\sqrt[3]{x/p}\)，可得 \(a_0p\le\sqrt{x}\)。

• 若 \(p>\sqrt[4]{x}\)，当 \(a_0=1\) 时，\(a_0p=p\)，它是一个质数，结论显然成立；当 \(a_0\ne 1\) 时，取 \(x/p\) 的最小质因子 \(q\)，因为 \(q\) 也是 \(x\) 的质因子，所以 \(p\le q\)，又由于 \(a_0,b_0,c_0\)是 \(x/p\) 不为1的因数，从而 \(q\le a_0\le b_0\le c_0\)，\(p×a_0×b_0×c_0>x\)，矛盾，不存在此情况。

于是 \(a_0p,b_0,c_0\) 排序后是 \(x\) 的一个分解。

容易在 \(O(n)\) 的时间内预处理出 \(\sqrt{n}\) 内的数两两的最大公因数，下面考虑求 \(\gcd(x,y)\)：

引理：若 \(r|a,r|b\)，则 \(\gcd(a,b)=r\gcd(a/r,b/r)\)。

这是显然的。

设 \((a,b,c)\) 是 \(x\) 的一个分解，令 \(p_1=y,r_1=\gcd(a,p_1)\)；

则 \(r_1|a,r_1|p_1\)，从而\(\gcd(x,y)=r_1\gcd(x/r_1,p_1/r_1)\)；

再令 \(p_2=p_1/r_1,r_2=\gcd(b,p_2)\)，同上有 \(\gcd(x,y)=r_1r_2\gcd(x/(r_1r_2),p_2/r_2)\)；

再令 \(p_3=p_2/r_2,r_3=\gcd(c,p_3)\)，同上有 \(\gcd(x,y)=r_1r_2r_3\gcd(x/(r_1r_2r_3),p_3/r_3)\)；

易知此时 \(\gcd(x/(r_1r_2r_3),p_3/r_3)=1\)，故 \(\gcd(x,y)=r_1r_2r_3\)。

由于我们预处理出了 \(\sqrt{n}\) 内所有数两两的 GCD，故上述算法过程可以做到 \(O(1)\)。

Code

Talking is cheap, show you the code.

#include <cstdio>
#include <algorithm>
#include <iostream>

using namespace std;
const int N = 5000;
const int M = 1e6;
const int T = 1e3;
const int mod = 998244353;
int n, tot, a[N + 1], b[N + 1], prime[M + 1], fac[M + 1][3], gcd[T + 1][T + 1];
int A[N + 1];
bool check[M + 1];

int read()
{
	int ans = 0; char ch;
	for(; !isdigit(ch = getchar()); );
	for(; isdigit(ch); ans = (ans << 1) + (ans << 3) + (ch ^ 48), ch = getchar());
	return ans;
}
void sieve()
{
	fac[1][0] = fac[1][1] = fac[1][2] = 1;
	for(int i = 2; i <= M; i++)
	{
		if(!check[i])
		{
			prime[++tot] = i;
			fac[i][0] = fac[i][1] = 1, fac[i][2] = i;
		}
		for(int j = 1; j <= tot && i * prime[j] <= M; j++)
		{
			int k = i * prime[j];
			check[i * prime[j]] = 1;
			fac[k][0] = fac[i][0] * prime[j];
			fac[k][1] = fac[i][1], fac[k][2] = fac[i][2];
			//sort
			if(fac[k][0] > fac[k][1]) swap(fac[k][0], fac[k][1]);
			if(fac[k][1] > fac[k][2]) swap(fac[k][1], fac[k][2]);
			//
			if(i % prime[j] == 0)
				break;
		}
	}
}
void init()
{
	for(int i = 1; i <= T; i++)
	{
		gcd[i][0] = gcd[0][i] = i;
		for(int j = 1; j <= i; j++)
			gcd[i][j] = gcd[j][i] = gcd[j][i % j];
	}
}
int GCD(int a, int b)
{
	int ret = 1;
	for(int i = 0, r; i < 3; i++)
	{
		if(fac[a][i] > T)
		{
			if(b % fac[a][i]) r = 1;
				else r = fac[a][i];
		}
		else
			r = gcd[fac[a][i]][b % fac[a][i]];
		b /= r;
		ret = ret * r;
	}
	return ret;
}
int main()
{
	n = read();
	for(int i = 1; i <= n; i++)
		a[i] = read();
	for(int i = 1; i <= n; i++)
		b[i] = read();
	sieve();
	init();
	for(int i = 1; i <= n; i++)
	{
		for(int j = 1, pw = i; j <= n; j++, pw = 1ll * pw * 1ll * i % mod)
			A[i] = (1ll * A[i] + 1ll * pw * 1ll * GCD(a[i], b[j])) % mod;
		printf("%d\n", A[i]);
	}
	return 0;
}