扩展欧几里得

扩展欧几里德算法

欧几里德算法

有两个数 a b，现在，我们要求 a b 的最大公约数，怎么求？枚举他们的因子？不现实，当 a b 很大的时候，枚举显得那么的naïve ，那怎么做？ 

欧几里德不愧是一个令我们钦佩的数学家提出了许多十分又有用的定理，其中就包括欧几里得算法： \(gcd(a, b) = gcd(b , a\%b)\) ，这样，我们可以在几乎是 \(log\) 的时间内解出 \(a\) 和 \(b\) 的最大公约数了。

代码:

int gcd(int x, int y)
{
  if(y == 0)
    return x;
  else
    return(y, x % y);
}

扩展欧几里德

首先我们存在两个数\(a, b\)那么一定存在\(x,y\) ，使得: \(ax + by = gcd(a,b)\) 现在让我们求解出\(x,y\)使得： $ax + by = gcd(a,b) $的通解 \(x 和 y\),此时用到的算法就是扩展欧几里得算法。
讲实现和原理之前，我们一定要深刻理解我们要求什么，及\(a,b,x,y\)这四个量的意义，首先\(x\),\(y\)它的值是由\(a\),\(b\)确定的，所以\(x,y\)是因变量(注意是变化的)，而根据原有的\(a,b\)直接求\(x,y\)是很难的，所以我们只得采用扩展欧几里得。

原理及证明

这是一个不定方程，有多解是一定的，但是只要我们找到一组特殊的解 \(x_{end}\) 和 \(y_{end}\) 那么，我们就可以用 \(x_0\) 和 \(y_0\) 表示出通解：
       \(x_{通} = x_{end}+(\frac{b}{gcd})*t\)
       \(~ y_{通} ~= ~y_{end} ~ – ~ (\frac{a}{gcd})*t\)
    我们令 \(B = \frac{b}{gcd}\) ， \(A = \frac{a}{gcd}\) ， 那么，\(A\)和\(B\)一定是互素的，所以最小的系数一定就是\(A\)和\(B\)了。

实现

现在，我们知道了有\(x\)和\(y\)使得 ：\(ax + by = gcd(a,b)\) ， 那么，求出这个特解\(x\)和\(y\)的方法就是使用扩展欧几里得算法。首先我们因为由上文可知\(gcd(a,b) = gcd(b, a \% b)\)，所以我们就可以从\(ax+by=gcd(a,b)\)继续向下推，可得\(bx_{now} + (a\%b)y_{now} = gcd(b, (a\%b)) = gcd(a,b) = ax_0 + by_0\)，此时我们要用下标将\(x_0，y_0\)和\(x_{now}，y_{now}\)分开,因为他们所在式子中的\(a,b\)是不一样的，所以他们的值也是不一样的(这点非常重要)
    我们观察到：欧几里德算法停止的状态是： \(a= gcd\) ， \(x_{now}=1,y_{now}=0\) ，那么这时，我们就会有： \(a*1 + b * 0 = gcd(a,b)\)
    可这是最终状态，我们仅仅只是求出了最终状态下的解\(x\)和\(y\),但此时的\(a\),\(b\)已经物是人非了,所以我们必须从最终状态回溯到最初的状态。
    假设当前我们已知\(a\)和\(b\)及其\(gcd\)，想要求出\(x_0\)和\(y_0\)使得\(ax_{0}+ by_{0} = gcd(a,b)\),而我们已经求出了下一组\(x_{now}\)和\(y_{now}\) 使得： \(bx_{now} + (a\%b)y_{now} = gcd(b,a\%b)\)， 那么是否可以寻找到它们之间的关系。
    我们知道： \(a\%b = a - (\frac{a}{b})*b\)(\(\frac{a}{b}\)指的是整除)。
那么，我们还可以继续向下推：
  \(gcd(b,(a\%b))\)
  $=bx_{now} +(a-(\frac{a}{b})*b)y_{now} $

​ \(= bx_{now} + ay_{now}– (\frac{a}{b})*by_{now}\)

\(= ay_{now} + b(x_{now} – \frac{a}{b}y_{now})\)

​ \(=ax_0+by_0\)
你是否找到什么关系？
        \(x_0 = y_{now}\)

\(y_0 = x_{now} – \frac{a}{b}y_{now}\)

以上就是扩展欧几里德算法的原理和实现。
代码:

void exgcd(ll a,ll b,ll &x,ll &y)
{
    if(!b) x=1,y=0;
    else
    {
        exgcd(b,a%b,y,x);
        y-=(a/b)*x;
    }
}