欧几里得算法(辗转相除法)及扩展欧几里得算法

· · 个人记录

欧几里得算法

简介

欧几里得,古希腊人,数学家,被誉为“几何之父”。有著作《几何原本》。

欧几里得算法,又称辗转相除法,用于计算两个正整数 a,b 的最大公约数。在中国古代有与之相同的算法,名曰“更相减损术”,该术记录于《九章算术》中。但提出时间晚于欧几里得。

最大公约数(\text{Greatest Common Divisor}\gcd),指两数最大的公共因数(约数)。

因数,又名约数。若 a\div b 余数为零(即 b 能整除 a ,记作 b|a),则称 ba 的一个因数,也称 ba 的一个因子。

算法内容

两个整数的最大公约数等于其中较小的那个数和两数相除所得余数的最大公约数。

即:若 a>b ,则 \gcd(a,b)=\gcd(b,a\mod b)

算法证明

假设 a>ba,b 均为正整数。

以下, x|y 表示 x 能整除 y (或 y 能被 x 整除)

a=rb+s ,其中 r,s 均为正整数且 s<b ,则 s=a\mod b

a,b 的任意一个公因数 d ,即对于任意一个满足 d|a,d|b 的正整数 d

$\therefore d|(a-rb)$ ,即 $d|s$ ; 也就是说,任意一个 $a,b$ 的公因数,都是 $r$ 的因数; 即 $d|\gcd(a,b)\Rightarrow d|\gcd(b,a\mod b)$ 。 假设正整数 $d$ 满足 $d|b,d|s$ ,则 $d|(rb)$ ,则 $d|(rb+s)$ ,即 $r|a$ ; 也就是说,任意一个 $b,r$ 的公因数,都是 $a$ 的因数; 即 $d|\gcd(a,b)\Leftarrow d|\gcd(b,a\mod b)$ 。 综上 $d|\gcd(a,b)\Leftrightarrow d|\gcd(b,a\mod b)$ 。 由于 $d$ 是任意取的,那不妨取 $d=\gcd(a,b)$ ,可以得到 $\gcd(a,b)|\gcd(b,a\mod d)$ , 再取 $d=\gcd(b,a\mod b)$ ,又可以得到 $\gcd(b,a\mod b)|\gcd(a,b)$ 。 所以可以得到 $\gcd(a,b)=\gcd(b,a\mod b)$ ,得证。 ## 算法应用 用于求正整数 $a,b$ 的最大公约数。 因为 $b<a,a\mod b<b$ ,所以使用一次欧几里得算法就会使问题变为求更小的两个数的最大公约数,从而缩小了问题规模。 而当什么时候发现 $a\mod d=0$ ,意味着此时的 $b$ 就是 $a,b$ 的最大公约数。 所以可以写出如下代码,来求 $a,b$ 的最大公约数。 可以使用递归: ```cpp int gcd(int a,int b) { return b?gcd(b,a%b):a; } ``` 或者使用循环: ```cpp int gcd(int a,int b) { int c; for(;b;a%=b) c=a,a=b,b=c; return a; } ``` 或者压行: ```cpp int gcd(int a,int b){int c;for(;b;a%=b,c=a,a=b,b=c);return a;} ``` ## 算法时间复杂度 考虑一下极限状态下,该算法的运行次数。 极限状态其实就是,每次取模都只是相当于做了个减法。 假设现在要计算 $a,b$ 的最大公约数, 设 $a_0=a,a_1=b,a_2=a_0\mod a_1,\cdots,a_n=a_{n-2}\mod a_{n-1}

其中 a_n=0

a_k=a_{k-2}\mod a_{k-1}\leqslant a_{k-2}-a_{k-1}

根据这个不等式的形式,很容易联想到斐波那契数列。 所以问题基本上相当于求 ${\large\frac{a}{\gcd(a,b)}}$ 在斐波那契数列中哪两项之间。 而我们显然可以知道,若设 $F_1=F_2=1,F_n=F_{n-1}+F_{n-2}$ ,则 $F_n\geqslant 2F_{n-2}$ ,所以当 $n$ 稍大些时,便有 $n\leqslant2\log_2F_n$ 。 所以算法的时间复杂度在 $O(log_2n)$ 附近,但实际上远远达不到。 + 下面给出斐波那契数列的前 $50$ 项(每行 $5$ 个): $1,1,2,3,5, 8,13,21,34,55, 89,144,233,377,610, 987,1597,2584,4181,6765, 10946,17711,28657,46368,75025, 121393,196418,317811,514229,832040, 1346269,2178309,3524578,5702887,9227465, 14930352,24157817,39088169,63245986,102334155 165580141,267914296,433494437,701408733,1134903170, 而 $2log_2(12586269025)\approx67>50$ 。 可见斐波那契数列增长极快。 # 扩展欧几里得算法 ## 简介 扩展欧几里德算法($exgcd$),用于求不定方程 $ax+by=m$ 的解 $x,y$ 。基于欧几里得算法,在辗转相除的过程中求解不定方程。 ## 前置知识 ### 1. 裴蜀定理 线性不定方程 $ax+by=m$ 有解,当且仅当 $\gcd(a,b)|m$ ,且有解时一定有无穷多解。 证明: 设 $m_0$ 为 $ax+by$ 的最小正值,为了方便表示,设 $r={\large\lfloor\frac{a}{m_0}\rfloor}$ , 则 $a\mod m_0=a-rm_0=a-r(ax+by)=a(1-rx)+b(-ry)$ , 我们发现 $a\mod m_0$ 也可以写成 $ax'+by'$ 的形式, 考虑到 $0\leqslant a\mod m_0<m_0$ ,而 $m_0$ 是 $ax+by$ 的最小正值, 所以 $a\mod m_0=0$ ,即 $m_0|a$ ; 同理 $m_0|b$ ; 于是 $m_0|\gcd(a,b)$ 。 而 $\gcd(a,b)|a,\gcd(a,b)|b$ ,所以 $\gcd(a,b)|(ax+by)$ , 也就是说 $\gcd(a,b)|m_0$ 。 所以 $m_0=\gcd(a,b)$ 。 即 **$ax+by$的最小正整数值为 $\gcd(a,b)$** 。 若 $m_0|m$ ,则 $ax+by=m$ 一定有解。 若 $m_0\not|\ \ m$ ,则 $ax+by=m$ 有解 $\Leftrightarrow$ $ax+by=m\mod m_0$ 有解,而后者显然无解,所以 $ax+by=m$ 无解。 综上,$ax+by=m$ 有解 $\Leftrightarrow \gcd(a,b)|m$ ,得证。 不妨假设不定方程的最小正整数解为 $\begin{cases}x=x_0\\y=y_0\end{cases}$ , 易知 $\begin{cases}x=x_0+{\large\frac{b}{\gcd(a,b)}}t\\y=y_0-{\large\frac{a}{\gcd(a,b)}}t\end{cases}(\ t$ 为参数,且为整数$)$ 都是不定方程的解。 而且这个参数方程包含了所有的解,证明略。 ### 2. 欧几里得算法 上面已经讲过了。 ## 算法思路 若 $ax+by=m$ 有解。 考虑在递归求最大公因数时,最大公因数一直不变。$ax+by=m$ 有解,则 $(bx+(a\mod b)y=m)$ 也有解。 **那么这些解之间有没有关系呢?** 当然是有的。来看看关系是什么。 考虑递归边界,$a'=\gcd(a,b),b'=0$ , 此时不定方程 $a'x+b'y=m$ 显然有一解 $\begin{cases}x_0={\large\frac m{\gcd(a,b)}}\\y_0=0\end{cases}$ , 我们考虑上一层递归的答案,不妨就假设 $a'=b,b'=a\mod b$ 好了, 已知为 $a'x_0+b'y_0=m$ ,即$bx_0+(a\mod b)y_0=m$ ,现在要把这个式子写成 $ax_1+by_1=m$ 这种形式,来研究 $x_1,y_1$ 与 $x_0,y_0$ 的关系,从而由 $x_0,y_0$ 推算出 $x_1,y_1$ 的值。 $\ \ \ \ \ bx_0+(a\mod b)y_0=m \Leftrightarrow bx_0+(a-{\large\lfloor\frac ab\rfloor}b)y_0=m 所以 $\begin{cases}x_1=y_0\\y_1=x_0-{\large\lfloor\frac ab\rfloor}y_0\end{cases}$ , 于是我们得到了递归时,相邻两层之间答案的关系, 而 $a,b$ 在递归时是已知的,下一层答案的计算又总早于上一层答案的计算。 所以在自下而上回溯时,就可以逐层将答案计算出来啦。 因为每层的 $a,b$ 在回溯时还有用,所以常使用递归来实现扩展欧几里得算法。 ## 算法应用 + 求解不定方程 $ax+by=m$ 。 + 求同余方程 $ax\equiv m\pmod b$ 的最小正整数解 $x$ 。 方法是:将这个同余方程写成不定方程 $ax+by=m$ 的形式,解出一个 $x$ ,则 $((x\mod b+b)\mod b)$ 一定是该同余方程的最小正整数解。至于 $y$ 最后算出来是多少,我们并不关心。 + 求 $a$ 在模 $b$ 意义下的逆元。 方法是:将题意转化为求同余方程 $ax\equiv1\pmod b$ 的最小正整数解 $x$ 。再转化为不定方程的形式。 + 作为各种各样数论题的基础。比如 $\text{noip 2017}$ 某凯的疑惑。 求解不定方程 $ax+by=\gcd(a,b)$ 的代码如下: ```cpp void exgcd(int a,int b,int &x,int &y) { if(!b) {x=1;y=0;return ;} exgcd(b,a%b,y,x); y-=a/b*x; } ``` # 例题 实在没有模板题。。。就用这个来代替吧: [P1082 同余方程](https://www.luogu.org/problem/P1082)。 $AC$ 代码: ```cpp #include<iostream> #include<cstdio> #include<cmath> using namespace std; int a,b,x,y; void exgcd(int p,int q,int &x,int &y) { if(q==0) {x=1;return ;} exgcd(q,p%q,y,x); y-=(p/q*x); } int main() { scanf("%d%d",&a,&b); exgcd(a,b,x,y); printf("%d",(x%b+b)%b); return 0; } ```