「学习笔记」O(1)求最大公约数

一种\(O(\)值域\()\)时间预处理\(O(1)\)时间求最大公约数(\(\gcd\))的算法

一些约定

1. \(N\)为询问的值域
2. \(Prime\)为全体素数集合
3. 集合\(\{a_1,a_2,\cdots,a_m\}\)是\(n\)的分解，当且仅当\(a_1\times a_2\times\cdots\times a_m=n\)

原理

定理一

内容

可以将值域中的每个\(x\)分解成\(\{a,b,c\}\)，满足\(a,b,c\leq\sqrt{x}\)或\(\in Prime\)（定义这种分解为合法分解） #### 证明 不妨设\(a\leq b\leq c\)若\(c\notin Prime\)且\(c>\sqrt{x}\),则\(c\)可分解为\(\{d,e\}\)且\(d\leq e\)有\(d\leq\sqrt{x}\)，而\(a\times b=\frac{x}{c}<\sqrt{x}\)则有\(n\)的分解\(\{d,a\times b,e\}\)，若\(e>\sqrt{x}\)则可按该规律一直分解直到\(e\in Prime\)或\(\leq\sqrt{x}\) ### 定理二 #### 内容 对于询问\(\gcd(x,y)\),分别考虑\(a,b,c\)对答案的贡献,\(a\)对答案的贡献为\(\gcd(a,y)\)，再将\(y\)除以\(\gcd(a,y)\)（这个因子已经被算过，不能重复计算），再对\(b,c\)干同样的事，三个贡献相乘即为\(\gcd(x,y)\)

证明

易得引理若\(r\mid p, r\mid q\)则\(\gcd(p,q)=r\times\gcd(\frac{p}{r},\frac{q}{r})\)

分别代入$ \[\begin{cases} &p_1=a\times b\times c,q_1=y,r_1=\gcd(a,q_1) \\ &p_2=b\times c,q_2=\frac{q_1}{r_1},r_2=\gcd(b,q_2) \\ &p_3=c,q_3=\frac{q_2}{r_2},r_3=\gcd(c,q_3) \end{cases}\]

$即可

实现

我们发现实现的难点在于如何在\(O(N)\)时间内进行分解，询问部分较为容易

分解

对于\(x\geq2\)，找到\(x\)的最小质因子\(p\)以及\(\frac{x}{p}\)的合法分解\(\{a_0,b_0,c_0\}\)且\(a_0\leq b_0\leq c_0\)，\(x\)的一种合法分解即为\(\{a_0\times p,b_0,c_0\}\)的升序排序

考虑证明： 1. \(x\in Prime\)时显然成立,分解为\(\{1,1,x\}\) 2. 当\(p\le\sqrt[4]{x}\)时将\(a_0\leq\sqrt[3]{\frac{x}{p}}\)带入有\(a_0\times p\le\sqrt{x}\) 3. 考虑\(p>\sqrt[4]{x}\)的情况，\(\left(1.\right)\) \(a_0=1\)，显然有\(a_0\times p=p\le\sqrt{x}\)；\(\left(2.\right)\) \(a\neq1\)，由于\(x\)不是素数，\(\frac{x}{p}\)的最小质因子\(q\)即为\(x\)的第二小质因子，一定\(\geq p\)，而\(a_0,b_0,c_0\)都为\(\frac{x}{p}\)的非\(1\)因子，有\(p\leq q\leq a_0\leq b_0\leq c_0\),\(p\times a_0\times b_0\times c_0>(\sqrt[4]{x})^4=x\)与\(p\times a_0\times b_0\times c_0=x\)相矛盾，故不存在此情况

所以只用跑一次线性筛，用最小质因子更新即可，然后预处理出\(\sqrt{n}\times\sqrt{n}\)的\(\gcd\)数组

代码如下

// fac为合法分解，isp表示是否非质数，pri为质数数组，tot为pri的大小，pre为预处理的gcd数组，M为值域，T为sqrt(M)
void work() {
  fac[1][0] = fac[1][1] = fac[1][2] = 1;
  for (int i = 2; i <= M; ++i) {
    if (!isp[i]) {
      fac[i][0] = fac[i][1] = 1;
      fac[i][2] = i;
      pri[++tot] = i;
    }
    for (int j = 1; pri[j] * i <= M; ++j) {
      int tmp = pri[j] * i;
      isp[tmp] = true;
      fac[tmp][0] = fac[i][0] * pri[j];
      fac[tmp][1] = fac[i][1];
      fac[tmp][2] = fac[i][2];
      if (fac[tmp][0] > fac[tmp][1]) {
        fac[tmp][0] ^= fac[tmp][1] ^= fac[tmp][0] ^= fac[tmp][1];
      }
      if (fac[tmp][1] > fac[tmp][2]) {
        fac[tmp][1] ^= fac[tmp][2] ^= fac[tmp][1] ^= fac[tmp][2];
      }
// 对于整数运算a ^= b ^= a ^= b等价于swap(a, b)这里就是手动进行length = 3的排序
      if (i % pri[j] == 0) {
        break;
      }
    }
  }
  for (int i = 0; i <= T; ++i) {
    pre[0][i] = pre[i][0] = i;
  }
  for (int i = 1; i <= T; ++i) {
    for (int j = 1; j <= i; ++j) {
      pre[i][j] = pre[j][i] = pre[j][i % j];
    }
  }
}

查询

若当前枚举的为\(a\)，只需注意\(a>\sqrt{N}\)时分\(a\mid y\)和\(a\nmid y\)两种情况即可

以下为代码

int gcd(int a, int b) {
// 不想写if-else所以用三目运算符代替并缩进了一下
  int ans = 1;
  for (int i = 0; i < 3; ++i) {
    int tmp = (fac[a][i] > T) ?
                (b % fac[a][i] ?
                   1
                 : fac[a][i]
                )
              : pre[fac[a][i]][b % fac[a][i]];
    b /= tmp;
    ans *= tmp;
  }
  return ans;
}