《基于值域预处理的快速 GCD》

在 \(\mathrm O(v)\) 预处理后 \(\mathrm O(1)\) 回答两个数的 gcd，其中 \(v\) 是值域。奇技淫巧了属于是。

注意到，若 \(n\) 的所有质因子都不超过 \(\sqrt n\)，则一定能找到三个数 \(a,b,c\) 的乘积等于 \(n\)，且三个数都不超过 \(\sqrt n\)。

证明很简单。考虑将 \(n\) 的所有质因子随意排成一列（允许重复），设为 \(p_{1\sim m}\)，\(P(l,r)=\prod\limits_{i=l}^rp_i\)。找到从左往右最后一个 \(P(1,x)\leq\sqrt n\) 的 \(x\)，显然有：

1. \(P(1,x)\leq\sqrt n\)，但 \(P(x+1,m)\geq\sqrt n\)；
2. \(P(1,x+1)>\sqrt n\)，但 \(P(x+2,m)<\sqrt n\)。

这样只需要把 \(p_{x+1}\) 单独分出来，令 \(a=P(1,x),b=p_{x+1},c=P(x+2,m)\)，即可保证 \(a,b,c\leq\sqrt n\)。

这样容易知道，对任意 \(n\)，都能找到 \(a\leq b\leq c\) 且 \(abc=n\)（此时显然 \(a,b\leq\sqrt n\)），满足 \(c\leq\sqrt n\) 或 \(c\in\mathbb P\)。如果对 \(1\sim v\) 以内每个数 \(x\) 都找到这样的 \(a_x,b_x,c_x\)，那么若求 \(\gcd(x,y)\)，则可对 \(d=(a/b/c)_x\) 依次进行 ans *= gcd(d, y), y /= gcd(d, y)。若 \(d\leq\sqrt n\)，有 \(\gcd(d,y)=\gcd(d,y\bmod d)\)，可以预处理 \(\sqrt v\times\sqrt v\) 的 gcd 表（复杂度 \(\mathrm O(v)\)）并 \(\mathrm O(1)\) 回答；对 \(d\in\mathbb P\)，只需要看是否有 \(d\mid y\) 即可。回答 gcd 共需进行 3 次取模、2 次除法，为大常数 \(\mathrm O(1)\)。

至于如何求出所有 \(a_x,b_x,c_x\)？一个简洁的算法是，设 \(x\) 的最小质因子为 \(p\)，令 \(y=\dfrac xp\)，则将 \(\{a_yp,b_y,c_y\}\) 排序则可得到 \(a_x,b_x,c_x\)，直接线性筛即可，复杂度 \(\mathrm O(v)\)。正确性证明是 dirty works，此处从略。

放一份模板题的代码：

code

constexpr int N = 5010, V = 1e6 + 10;

struct gcd_solver {
	static constexpr int B = 1000;
	int d[B + 10][B + 10];
	struct _tuple {
		int a, b, c;
		_tuple(int a = 1, int b = 1, int c = 1) : a(a), b(b), c(c) {}
		void operator*=(int x) {
			a *= x;
			if(a > b) {
				swap(a, b);
				if(b > c) swap(b, c);
			}
		}
	} tp[V];
	gcd_solver() {
		REP(x, 1, B) REP(y, 0, x) d[x][y] = d[y][x] = y ? d[y][x % y] : x;
		static int prm[V], cnt; static bool vis[V];
		REP(i, 2, 1e6) {
			if(!vis[i]) prm[++cnt] = i, tp[i] = _tuple(1, 1, i);
			REP(j, 1, cnt) {
				int x = prm[j];
				if(i * x > V) break;
				vis[i * x] = true;
				auto &t = tp[i * x] = tp[i]; t *= x;
				if(i % x == 0) break;
			}
		}
	}
	int operator()(int x, int y) {
		auto &t = tp[y];
		int g = d[x % t.a][t.a], ans = g; x /= g;
		g = d[x % t.b][t.b], ans *= g; x /= g;
		if(t.c <= B) ans *= d[x % t.c][t.c];
		else if(x % t.c == 0) ans *= t.c;
		return ans;
	}
} gcd;

int n;
int a[N], b[N];

void mian() {
	n = read();
	REP(i, 1, n) a[i] = read();
	REP(i, 1, n) b[i] = read();
	REP(i, 1, n) {
		int ij = 1, ans = 0;
		REP(j, 1, n) addto(ans, (ll)(ij = (ll)ij * i % P) * gcd(a[i], b[j]) % P);
		prt(ans), pc('\n');
	}
}