求解最大公约数

问题：求两个数a,ba,b的最大公约数。

我们把a,ba,b的最大公约数记为gcd(a,b)gcd(a,b)。

朴素算法

朴素的算法是，求出a,ba,b的素因子分解，取两个数共同的素因子，即为最大公约数：

a=p1a1p2a2...pnana=p_1^{a_1}p_2^{a_2}...p_n^{a_n}

b=p1b1p2b2...pnbnb=p_1^{b_1}p_2^{b_2}...p_n^{b_n}

gcd(a,b)=p1min(a1,b1)p2min(a2,b2)...pnmin(an,bn)gcd(a,b)=p_1^{min(a_1,b_1)}p_2^{min(a_2,b_2)}...p_n^{min(a_n,b_n)}

a1,a2,...,an,b1,b2,...,bn≥0a_1,a_2,...,a_n,b_1,b_2,...,b_n \ge 0（任意数的0次幂等于1，一般我们会省略掉）

例如：588000=25∗3∗53∗72588000=2^5*3*5^3*7^2,540=22∗33∗5∗70540=2^2*3^3*5*7^0

它们的最大公约数为：60=22∗3∗5∗7060=2^2*3*5*7^0

这个算法需要先对a,ba,b做素因子分解，计算量较大，且实现较为复杂，实际中我们不会使用这种方式计算最大公约数。

辗转相除法

有一种更高效的算法可以求最大公约数，叫作辗转相除法（也叫做欧几里得算法）。

在介绍这个算法之前，需要引入一个定理：

(a,b)的公约数与(b,a%b)的公约数完全相同(a%b表示a mod b，即取余操作)(a,b)的公约数与(b,a\%b)的公约数完全相同(a\%b表示a\ mod\ b，即取余操作)

证明：

  1. 如果a<ba<b，则(b,a%b)=(b,a)(b,a\%b)=(b,a)，只是相当于交换了两个数的位置，显然公约数相同
  2. 如果a=ba=b，则(b,a%b)=(b,0)=(a,0)(b,a\%b)=(b,0)=(a,0)，由于任意数都是00的约数，显然公约数也是相同的
  3. 如果a>ba>b，我们需要从两个方向证明
    • 如果dd是(a,b)(a,b)的公约数，则dd也是(b,a%b)(b,a\%b)的公约数：

      a=k∗b+c(0<c<b)c=a%b记a=ka∗d,b=kb∗dc=a−k∗b=ka∗d−k∗kb∗d=(ka−k∗kb)∗d(ka−k∗kb)为整数→d为c=a%b的约数→d为(b,a%b)的公约数\begin{align} a=k*b+c(0<c<b) \\ c=a\%b \\ 记a=k_a*d,b=k_b*d \\ c=a-k*b=k_a*d-k*k_b*d=(k_a-k*k_b)*d \\ (k_a-k*k_b)为整数 \\ \rightarrow d为c=a\%b的约数 \\ \rightarrow d为(b,a\%b)的公约数 \end{align}

    • 如果dd是(b,a%b)(b,a\%b)的公约数，则dd也是(a,b)(a,b)的公约数：

      a=k∗b+c(0<c<b)c=a%b记b=kb∗d,c=kc∗da=k∗b+c=k∗kb∗d+kc∗d=(k∗kb+kc)∗d(k∗kb+kc)为整数→d为a的约数→d为(a,b)的公约数\begin{align} a=k*b+c(0<c<b) \\ c=a\%b \\ 记b=k_b*d,c=k_c*d \\ a=k*b+c=k*k_b*d+k_c*d=(k*k_b+k_c)*d \\ (k*k_b+k_c)为整数 \\ \rightarrow d为a的约数 \\ \rightarrow d为(a,b)的公约数 \end{align}

由于(a,b)(a,b)的公约数与(b,a%b)(b,a\%b)的公约数相同，因此它们的最大公约数也是相同的

我们可以利用上述的定理得到一个递推算法，算法步骤如下：

g=gcd(a,b)=gcd(b,a%b)b⇒a2,a%b⇒b2g=gcd(a2,b2)=gcd(b2,a2%b2)...bi−1⇒ai,ai−1%bi−1⇒big=gcd(ai,bi)\begin{align} g=gcd(a,b)=gcd(b,a\%b) \\ b\rArr a_2,a\%b\rArr b_2 \\ g=gcd(a_2,b_2)=gcd(b_2,a_2\%b_2) \\ ... \\ b_{i-1}\rArr a_i,a_{i-1}\%b_{i-1}\rArr b_i \\ g=gcd(a_i,b_i) \end{align}

由于ai,bia_i,b_i在不断减小，所以递推式必定会在bi=0b_i=0时结束，此时aia_i即为a,ba,b的最大公约数，不存在比aa还要大的约数了。

代码实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
// 递归版本
long long gcd(long long a, long long b) {
  if (b==0) return a;
  return gcd(b, a%b);
}

// 更精简的写法
long long gcd(long long a, long long b) { return b ? gcd(b, a % b) : a; }

// 迭代版本
long long gcd(long long a, long long b) {
  while (b != 0) {
    long long pa = a;
    a = b;
    b = pa % b;
  }
  return a;
}

扩展欧几里得算法

如果gg为(a,b)(a,b)的最大公约数，则必定存在一对整数x,yx,y，使得g=a∗x+b∗yg=a*x+b*y，我们可以用扩展欧几里得算法来计算x,yx,y的值

方法步骤如下：在使用辗转相除法计算(a,b)(a,b)的最大公约数时，假设最终迭代了nn次，最终得到(an,bn=0)(a_n,b_n=0)的最大公约数为g=ang=a_n，我们可以得到(xn=1,yn=0)(x_n=1,y_n=0)：

g=an∗1+bn∗0g=a_n*1+b_n*0

根据数学归纳法，当我们知道满足以下等式的(xi,yi)(x_i,y_i)的值时

g=ai∗xi+bi∗yig=a_i*x_i+b_i*y_i

我们可以通过辗转相除法的递推关系倒推得到(xi−1,yi−1)(x_{i-1},y_{i-1})的值，具体做法如下：

ai=bi−1a_i=b_{i-1}

bi=ai−1%bi−1=ai−1−⌊ai−1bi−1⌋∗bi−1b_i=a_{i-1}\%b_{i-1}=a_{i-1}- \lfloor \frac{a_{i-1}}{b_{i-1}} \rfloor*b_{i-1}

⇒g=ai∗xi+bi∗yi=bi−1∗xi+(ai−1−⌊ai−1bi−1⌋∗bi−1)∗yi=ai−1∗yi+bi−1∗(xi−⌊ai−1bi−1⌋∗yi)\rArr g=a_i*x_i+b_i*y_i=b_{i-1}*x_i+(a_{i-1}- \lfloor \frac{a_{i-1}}{b_{i-1}} \rfloor*b_{i-1})*y_i\\=a_{i-1}*y_i+b_{i-1}*(x_i-\lfloor \frac{a_{i-1}}{b_{i-1}} \rfloor*y_i)

⇒xi−1=yi,yi−1=xi−⌊ai−1bi−1⌋∗yi\rArr x_{i-1}=y_i,y_{i-1}=x_i-\lfloor \frac{a_{i-1}}{b_{i-1}} \rfloor*y_i

由此递推式，如果我们已知xn,ynx_n,y_n，便可以计算得到xn−1,yn−1x_{n-1},y_{n-1}，不断重复此过程，最终可以得到一对x,yx,y，使得g=a∗x+b∗yg=a*x+b*y

代码实现

1
2
3
4
5
6
7
8
9
10
// 递归版本
long long exgcd(long long a, long long b, long long& x, long long& y) {
  if (b==0) {
    x = 1; y = 0;
    return a;
  }
  long long g = exgcd(b, a%b, y, x);
  y = y - a/b * x;
  return g;
}

迭代版本

上述的递推式需要先知道xi,yix_i,y_i，才能够计算xi−1,yi−1x_{i-1},y_{i-1}，看起来无法使用非递归方式迭代实现，然而，我们通过观察可以得到：

{xi−1=yiyi−1=xi−⌊ai−1bi−1⌋∗yi\left \{ \begin{aligned} x_{i-1} & = y_i \\ y_{i-1} & = x_i-\lfloor \frac{a_{i-1}}{b_{i-1}} \rfloor*y_i \end{aligned} \right.

这可以写成矩阵的形式：

[xi−1yi−1]=[011−⌊ai−1bi−1⌋]∗[xiyi]\begin{bmatrix} x_{i-1} \\ y_{i-1} \end{bmatrix} = \begin{bmatrix} 0 & 1 \\ 1 & -\lfloor \frac{a_{i-1}}{b_{i-1}} \rfloor \end{bmatrix} * \begin{bmatrix} x_{i} \\ y_{i} \end{bmatrix}

将迭代过程中的所有项展开后：

[xy]=[011−⌊ab⌋]∗[011−⌊a1b1⌋]∗...∗[011−⌊an−1bn−1⌋]∗[xn=1yn=0]\begin{bmatrix} x \\ y \end{bmatrix} = \begin{bmatrix} 0 & 1 \\ 1 & -\lfloor \frac{a}{b} \rfloor \end{bmatrix} * \begin{bmatrix} 0 & 1 \\ 1 & -\lfloor \frac{a_{1}}{b_{1}} \rfloor \end{bmatrix} *...* \begin{bmatrix} 0 & 1 \\ 1 & -\lfloor \frac{a_{n-1}}{b_{n-1}} \rfloor \end{bmatrix} * \begin{bmatrix} x_{n} = 1 \\ y_{n} = 0 \end{bmatrix}

代码实现
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
// 迭代版本
long long exgcd(long long a, long long b, long long& x, long long& y) {
  long long w[2][2];
  w[0][1] = w[1][0] = 0;
  w[0][0] = w[1][1] = 1;
  while (b != 0) {
    x = w[0][0];
    w[0][0] = w[0][1];
    w[0][1] = x - a/b * w[0][1];
    y = w[1][0];
    w[1][0] = w[1][1];
    w[1][1] = y - a/b * w[1][1];


    long long pa = a;
    a = b;
    b = pa % b;
  }
  x = w[0][0];
  y = w[1][0];
  return a;
}