关于 GCD

GCD

正常的 gcd

  相信大家都会辗转相除吧 qwq。

更相减损

  首先有一个引理：$\gcd (a,b)=\gcd (a-b,b)$

  证明： 令 $a=k_1\gcd (a,b),b=k_2\gcd (a,b)$ 显然有 $\gcd (k_1,k_2)=1$。然后 $\gcd (a-b,a)=\gcd ((k_1-k_2)\gcd (a,b),k_2\gcd (a,b))=\gcd (a,b)\gcd (k_1-k_2,k_2)$，现在我们就只需要证明 $\gcd (k_1-k_2,k_2)=1$ 就行了。

  所以我们考虑反证法，假设 $k_1-k_2=ax,k_2=bx$ 其中 $x=\gcd (k_1-k_2,k_2)\ne 1$。那么我们发现 $k_1=k_1-k_2+k_2=(a+b)x$，那么 $k_1$ 和 $k_2$ 就有了共同的不等于 $1$ 的因数，所以 $\gcd (k_1,k_2)\ne 1$ 这就和我们的前提矛盾了，所以 $\gcd (k_1-k_2,k_2)=1$。

  证毕。

辗转相除

  有了刚才的那个引理，我们就能很自然的推出下面这个式子：

$$\gcd (a,b)=\gcd (b,a\mathrm{mod}b)$$

  这个就是欧几里得求 $\gcd$ 的核心式子了（正确性不用解释了吧qwq），于是我们就能写出求 $\gcd$ 的非常简便的，只有一行的，正确的，毋庸置疑的，形而上的，代码了：

int gcd(int a, int b) { return b ? gcd(b, a % b) : a; }

  这个代码的复杂度就显然是 $O(\log a)$ 的，因为每次至少除以 $2$。

更快的 gcd

  这里是 $O(值域)$ 预处理 $O(1)$ 查询的 $\gcd$。

几个引理

$\mathrm{A}$

  引理一：可以将值域中的每一个数 $x$ 分解为 $x=a\times b\times c$，满足要么 $a,b,c\le \sqrt{x}$，要么 $a,b,c\in prime$，其中 $prime$ 表示质数集。

  证明：我们考虑找到 $x$ 的最小值因子 $p$，并找到 $\frac{x}{p}$ 的分解，然后把 $p$ 乘到 $\frac{x}{p}$ 的分解中的最小的那个数上面并成为 $x$ 的分解，这样就能满足上述条件（因为 $x$ 是质数的情况非常显然，所以只考虑 $x$ 不是质数时）。

  考虑归纳法证明上述做法的正确性，假设 $\frac{x}{p}$ 的分解是 { a , b , c } \{a, b, c\} {a,b,c} 且 $a\le b\le c$，那么 $a\le \sqrt[3]{\frac{x}{p}}$，那么（这里只需要证明 $ap\le \sqrt{x}$ 就可以了）：

$$ap\le p\sqrt[3]{\frac{x}{p}}$$

  分两种情况：

1. $p\le \sqrt[4]{x}$ 时，有 $ap\le \sqrt[4]{x}\sqrt[3]{\frac{x}{\sqrt[4]{x}}}=\sqrt{x}$ 这时显然成立。
2. $p>\sqrt[4]{x}$ 时，还要分两种情况考虑：
   1）$a=1$，那么显然就是 $ap=p\le \sqrt{x}$
   2）$a\ne 1$，我们令 $\frac{n}{p}$ 的最小质因子为 $q$，那么显然有 $p\le q\le a\le b\le c$，于是 $pabc\ge p^4>\sqrt[4]{x}=x$，这就和前面 $qabc=x$ 的前提矛盾了，所以只能有 $p\le \sqrt[4]{x}$。

$\mathrm{B}$

  第二个引理：对于一个询问 $\gcd (x,y)$，我们已经把 $x$ 拆分成了 $a,b,c$，那么进行以下操作得到的答案与 $\gcd (x,y)$ 相同：

i = gcd(a, y), y /= i;
j = gcd(b, y), y /= j;
k = gcd(c, y), y /= k;
ans = i * j * k;

  这个玩意儿的正确性就非常显然了，所以就不放证明了qwq。

算法的基本流程

  对于一个询问 $\gcd (a,b)$，我们就考虑用引理 $\mathrm{B}$ 来进行计算，那么我们就可以先考虑 $\gcd (x,a)$，答案分三种情况讨论：

1. $x\le \sqrt{n}$，考虑手动递归一层（辗转相除）：$\gcd (x,a)=\gcd (x,a\mathrm{mod}x)$，这里就有 $x\le \sqrt{n}$ 且 $a\mathrm{mod}x\le \sqrt{n}$，所以我们可以考虑打一个 $\sqrt{n}\times \sqrt{n}$ 的表，记录 $\gcd (i,j)$ 的值，然后就可以 $O(1)$ 查询了。
2. $x>\sqrt{n}\wedge x\mid a$，这个时候答案显然为 $x$
3. $x>\sqrt{n}\wedge x\nmid a$，这个时候，根据引理 $\mathrm{A}$，我们知道 $x$ 此时必定是一个质数，所以必定有 $\gcd (x,a)=1$。

  现在我们发现，如果把 $\sqrt{n}\times \sqrt{n}$ 的 $\gcd$ 表打出来并且解决数的拆分，然后就可以 $O(1)$ 计算 $\gcd$ 了。

  拆分的话就按照我们的 $\mathrm{A}$ 引理的证明中提到的方案直接构造就好了，时间复杂度时 $O(n)$ 的（线性筛）。

代码

#include<bits/stdc++.h>
using namespace std;
#define in read()
#define MAXN 5050
#define MAXM 1001000
#define MAXK 1010
#define mod 998244353

inline int read(){
	int x = 0; char c = getchar();
	while(c < '0' or c > '9') c = getchar();
	while('0' <= c and c <= '9'){
		x = x * 10 + c - '0'; c = getchar();
	}
	return x;
}

int n = 0; int tot = 0;
int a[MAXN] = { 0 };
int b[MAXN] = { 0 };
int fac[MAXM][3] = { 0 };                          // 分解 
int pre[MAXK][MAXK] = { 0 };                       // 预处理的 gcd 数组 
int vis[MAXM] = { 0 };
int pri[MAXM] = { 0 };

void prework(int N, int K){
	fac[1][0] = fac[1][1] = fac[1][2] = 1;
	for(int i = 2; i <= N; i++){
		if(!vis[i]) pri[++tot] = i, fac[i][0] = fac[i][1] = 1, fac[i][2] = i;
		for(int j = 1; j <= tot and i * pri[j] <= N; j++){
			int temp = i * pri[j];
			vis[temp] = 1;
			fac[temp][0] = fac[i][0] * pri[j];
			fac[temp][1] = fac[i][1];
			fac[temp][2] = fac[i][2];
			if (fac[temp][0] > fac[temp][1]) swap(fac[temp][0], fac[temp][1]);
	      	if (fac[temp][1] > fac[temp][2]) swap(fac[temp][1], fac[temp][2]);
			if(i % pri[j] == 0) break;
		}
	}
	for(int i = 0; i <= K; i++) pre[i][0] = pre[0][i] = i;
	for(int i = 1; i <= K; i++)
		for(int j = 1; j <= i; j++)
			pre[i][j] = pre[j][i] = pre[j][i % j];
}

int gcd(int x, int y){
	int ans = 1;
	for(int i = 0; i < 3; i++){
		int d = 0;
		if(fac[x][i] <= 1e3) d = pre[fac[x][i]][y % fac[x][i]];
		else if(y % fac[x][i] != 0) d = 1;
		else d = fac[x][i];
		y /= d, ans *= d;
	}
	return ans;
}

int main(){
	prework(1e6, 1e3);
	n = in;
	for(int i = 1; i <= n; i++) a[i] = in;
	for(int i = 1; i <= n; i++) b[i] = in;
	for(int i = 1; i <= n; i++){
		int ans = 0;
		for(int j = 1, now = i; j <= n; j++, now = 1ll * now * i % mod)
			ans = (ans + (1ll * now * gcd(a[i], b[j]) % mod) % mod) %mod;
		cout << ans << '\n';
	}
	return 0;
}