「学习笔记」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. 易得$a_0\leq\sqrt[3]{\frac{x}{p}}$
    解$a_0\times p\leq\sqrt{x}$得$p\leq \sqrt[4]{x}$,考虑$p>\sqrt[4]{x}$的情况,由于$x$不是素数,$\frac{x}{p}$的最小质因子$q$即为$x$的第二小质因子,一定$\geq p$,而$a_0,b_0,c_0$都为$\frac{x}{p}$的因子,有$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;
}

本博客采用知识共享署名-相同方式共享 4.0 国际许可协议进行许可
本文链接:https://blog.seniorious.cc/2020/quick-gcd/