最大公约数

暴力

时间复杂度$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)))$

Stein算法:辗转相除结合更相减损术

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

a偶b偶,$gcd(a,b) = 2*gcd(a/2, b/2) = 2*gcd(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(std::abs(a - b), std::min(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

http://littledva.cn/article-41/

https://xuanwo.io/2015/03/11/number-theory-gcd/

最小公倍数

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("%d\n", (x % b + b) % b);
    return 0;
}

其他还有

最后修改:2022 年 03 月 03 日
如果觉得我的文章对你有用,请随意赞赏