最大公约数

暴力

时间复杂度$O(min (a, b))$

辗转相除

辗转相除法，又名欧几里得除法，目的是求出两个正整数的最大公约数。

两个正整数a和b（a>b），它们的最大公约数等于a除以b的余数c和b之间的最大公约数。比如10和25，25除以10商2余5,那么10和25的最大公约数，等同于10和5的最大公约数。

int gcd(int a, int b) {
  if (!b || a == b) return a;
  if (a > b) return gcd(b, a % b);
  else return gcd(a, b % a);
}

时间复杂度近似$O(log min(a, b))$

更相减损术

辗转相除法涉及取模运算，当整数较大时，取模运算的效率比较差。这时可以看更相减损术。

更相减损术， 出自于中国古代的《九章算术》，也是一种求最大公约数的算法。两个正整数a和b（a>b），它们的最大公约数等于a-b的差值c和较小数b的最大公约数。比如10和25，25减去10的差是15,那么10和25的最大公约数，等同于10和15的最大公约数。

int gcd(int a, int b) {
  if (a == b) return a;
  if (a > b) return gcd(b, a - b);
  else return gcd(a, b - a);
}

时间复杂度$O(max (a, b)))$

辗转相除结合更相减损术

更相减损术不稳定，运算次数更大，尤其是当两个数比较悬殊的时候。考虑两者结合：

a偶b偶，$gcd(a,b) = 2gcd(a/2, b/2) = 2gcd(a>>1, b>>1)$

a偶b奇，$gcd(a,b) = gcd(a/2, b) = gcd(a>>1, b)$

a奇b偶，$gcd(a,b) = gcd(a, b/2) = gcd(a, b>>1)$

a奇b奇，利用更相减损术运算一次，$gcd(a,b) = gcd(b, a-b)$， 此时a-b必然是偶数，又可以继续进行移位运算。

int gcd(int a, int b) {
  if (a == b) return a;
  if (a & 1) {
    if (b && 1) return gcd(a, std::abs(a - b));
    else return gcd(a, b >> 1);
  }
 	else {
    if (b & 1) return gcd(a >> 1, b);
    else return 2 * gcd(a >> 1, b >> 1);
  }
}

时间复杂度$O(log max (a, b)))$

参考

https://houbb.github.io/2017/08/23/math-03-common-gcd-03

https://zhuanlan.zhihu.com/p/31824895

最小公倍数

a * b / gcd(a, b)

扩展欧几里得

原理

可用于求 $ a x + b y = gcd(a, b) $ 的一组解$(x_1, y_1)$。
$$
a x_1 + b y_1 = gcd(a, b) tag{1}
$$

下面是简要求解过程。

用b与$a % b$来替代a和b，得到
$$
bx_2 + (a % b) y_2 = gcd (b, a % b) tag{2}
$$

根据欧几里得算法以及$a % b$的算数表达形式
$$
begin{align}
gcd(a, b) &= gcd (b, a % b) tag{3}
a % b &= a - lfloor frac{a}{b} rfloor b tag{4}
end{align}
$$
联立式子$(2)$、$(3)$和$(4)$
$$
a y_2 + b left [ x_2 - lfloor frac{a}{b} rfloor y_2 right ] = gcd (a, b) tag{5}
$$
联立式子$(1)$和$(5)$
$$
begin{align}
x_1 &= y_2
y_1 &= x_2 - lfloor frac{a}{b} rfloor y_2
end{align}
$$
由此写出以下程序。

int exgcd(int a,int b,int &x,int &y)
{
    if (b==0) {
        x = 1;
        y = 0;
        return a;
    }
    int r = exgcd(b, a%b, x, y);
    int t = y;
    y = x - (a / b) * y;     //2的情况
    x = t;
    return r;
}

通解

如果 $(x_0, y_0)$ 是 $a x + b y = c$ 一组解，那么通解可以表示为
$$
begin{align}
x &= x_0 + frac{b}{c} k
y &= y_0 - frac{a}{c} k
end{align}
$$
事实上
$$
a x + b y = a left( x_0 + frac{b}{c} k right) + b left( y_0 - frac{a}{c} k right) = a x_0 + b y_0 = c
$$

参考

https://zhuanlan.zhihu.com/p/100567253

https://zhuanlan.zhihu.com/p/42707457

https://www.desgard.com/algo/docs/part2/ch02/3-ext-euclidean/

https://ksmeow.moe/euclid/

https://www.cnblogs.com/wkfvawl/p/9350867.html

例子

以Acwing 203 同余方程为例，$ax equiv 1 (mod m)$ 可以写成 $a x - 1 = m y$，即 $ax - my = 1$。形式上等同于 $a x + b y = c, b = -m, c = 1$。

根据扩展欧几里得算法可以解出一个$x_0$，那么根据通项解可知，就是在$x_0$上加减若干个 $frac{b}{c}m = m$ 。由于$x_0$可能是负数，或者是大于$m$的正数，通过$(x % b + b) % b$使其成为$(0, m)$之间的正整数。

#include <cstdio>
using namespace std;

int exgcd(int a, int b, int &x, int &y) {
    if (!b) {
        x = 1, y = 0;
        return a;
    }
    int gcd = exgcd(b, a % b, x, y), temp = y;
    y = x - a / b * y;
    x = temp;
    return gcd;
}
int main() {
    int a, b, x, y;
    scanf("%d%d", &a, &b);
    exgcd(a, b, x, y);
    printf("%dn", (x % b + b) % b);
    return 0;
}

其他还有

• Leetcode365 水壶问题