从欧几里得讲起……

欧几里得的辗转相除法计算的是两个自然数a和b的最大公约数g,意思是能够同时整除a和b的自然数中最大的一个。
两个数的最大公约数通常写成GCD(a, b),或者简写成(a, b)。
(Greatest Common Divisor)

计算方法……辗转相除!

作为数论中一个很基础的内容,它的形式很漂亮,代码也很干净。

gcd(a,b)=gcd(b,agcd(a,b) = gcd(b,a % b)

证明过程

r=ar = a % b则有a=k×b+ra = k \times b + r故可知r=ak×br = a - k \times b
现假设ddaa,bb的一个公约数
则有dad|a , dbd|b

又由上式可得d(ak×b)d|(a - k \times b)dr    d(ad|r \iff d|(a%b)
由此证明dd也是bb(a(a%b)的公约数
gcd(a,b)=gcd(b,a%b)gcd(a,b) = gcd(b,a\%b)

注:dad|a表示dd整除aa

代码实现

long long gcd(long long a,long long b)
{   return b == 0 ? a : gcd(b,a%b); }

贝祖等式……问题的开端!

贝祖等式(又称裴蜀等式)是一个形如ax+by=max + by = m的方程。
而有整数解(x,y)    m=k×gcd(a,b),kZ(x,y) \iff m = k \times gcd(a,b) , k \in \mathbb Z
裴蜀等式有解时必然有无穷多个解。

补充:特别的,ax+by=1    gcd(a,b)=1ax + by = 1 \iff gcd(a,b)=1

而利用辗转相除法可以得到这一方程的整数特解。
而且不难发现,大多数情况下xxyy总是一正一负。

一段巧妙的证明

ax+by=gcd(a,b)ax + by = gcd(a,b)开始……
由于长得漂亮对称的结构,我们不妨设a>ba \gt b

  • b=0b = 0时,gcd(a,b)=agcd(a,b) = a这个时候我们的等式就变成了ax=aax = a的形态不难看出这时候的一组解(x,y)=(1,0)(x,y) = (1,0)
  • b0b \not= 0时候,我们可以得到
    ax1+by1=gcd(a,b)ax_1 + by_1 = gcd(a,b)
    bx2+(a%b)y2=gcd(b,a%b)bx_2 + (a\%b)y_2 = gcd(b,a\%b)
    又由于gcd(a,b)=gcd(b,a%b)gcd(a,b) = gcd(b,a\%b)
    所以有ax1+by1=bx2+(a%b)y2ax_1 + by_1 = bx_2 + (a\%b)y_2
    定义可得a%b=aab×ba\%b = a - \lfloor {\cfrac{a}{b}} \rfloor \times b将其带入上式

    ax1+by1=bx2+(aab×b)y2ax_1 + by_1 = bx_2 + (a - \lfloor {\cfrac{a}{b}} \rfloor \times b)y_2

    ax1+by1=bx2+ay2(ab×b)y2ax_1 + by_1 = bx_2 + ay_2 - (\lfloor {\cfrac{a}{b}} \rfloor \times b)y_2

    aabb的对应系数相等,所以可推出
    x1=y2x_1 = y_2
    y1=x2aby2y_1 = x_2 - \lfloor {\cfrac{a}{b}} \rfloor y_2

综上,我们推导出了求该等式特解的边界条件和递推过程。

代码实现

注意 gcd , x , y 变量均使用引用传参,这是为了每一步都能改变原始数据。

void exgcd(long long a, long long b, long long &gcd, long long &x, long long &y)
{
    if (b == 0)
    {
        gcd = a;
        x = 1;
        y = 0;
        return;
    }
    exgcd(b, a % b, gcd, x, y);
    long long xt = x, yt = y;
    x = yt;
    y = xt - a / b * yt;
    return;
}

应用举例

求解同余方程的最小正整数解ax1(modb)ax \equiv 1 (\mod b),给出aabb

分析

ax1(modb)    ax+by=1ax \equiv 1 (\mod b) \iff ax + by = 1
其中yy是我们引入的一个整数而且十有八九应该是负数。
但是有一个问题,ax+by=1ax + by = 1又不是ax+by=gcd(a,b)ax + by = gcd(a,b)
这是怎么牵扯上关系的呢?还不是你乱点鸳鸯谱

由最大公约数的定义,总有agcd(a,b)a|gcd(a,b)而且bgcd(a,b)b|gcd(a,b)
如果有x,yZx,y \in \mathbb Z则必定(ax+by)gcd(a,b)(ax + by)|gcd(a,b)
即若有ax+by=max + by = mmgcd(a,b)m|gcd(a,b)
m=1m=1时,就有1gcd(a,b)1|gcd(a,b)
故而gcd(a,b)gcd(a,b)就只能等于11

但是我们求出的xx十有八九是负数也有可能太大,怎么得到最小正整数解呢?
ax+by=1ax + by = 1
ax+by+k×abk×ab=1ax + by + k \times ab - k \times ab= 1
a(x+kb)+b(yka)=1a(x + kb) + b(y - ka)= 1
故而我们使用(x%b+b)%b(x' \% b + b)\%b的方式计算最终结果

这里一开始我不太明白
后来想想,我们要求的是xx'在加上kkbb之后最小的整数结果
前提:x<0,b>0x' \lt 0, b \gt 0
x%bx'\%b可以将其结果限制在[b,0][-b,0]
b<x%b<0-b \lt x'\%b \lt 0时候加bb便是最终结果了
xx'%b = 0时候加bb再对bb取模便是0不会出现错误

这里涉及到负数的取模运算,但这个计算的结果是百家争鸣的,我们先以C++中的结果为准
在C++中,有7%3=1-7 \% 3 = -1
a=3,b=10a = 3, b = 10为例解得x=3,y=1x = -3,y = 1
计算得到(3%10+10)%10=7(-3 \% 10 + 10) \% 10 = 7

代码实现

void exgcd(long long a, long long b, long long &gcd, long long &x, long long &y)
{
    if (b == 0)
    {
        gcd = a;
        x = 1;
        y = 0;
        return;
    }
    exgcd(b, a % b, gcd, x, y);
    long long xt = x, yt = y;
    x = yt;
    y = xt - a / b * yt;
    return;
}
int main()
{
    long long a,b;
    long long g,x,y;
    cin >> a >> b;
    exgcd(a,b,g,x,y);
    x = (x % b + b)%b;
    cout << x << endl;
    return 0;
}