浅谈扩展欧几里得算法

什么是拓展欧几里得?简单的说,就是求关于x,y的方程 ax + by = gcd(a,b) 的所有整数解


现在我们来解决四个问题

  • 什么是裴属定理,如何证明裴属定理?

  • 怎么用扩展欧几里得来求ax + by = gcd(a,b) 的特解?

  • 怎么求由特解推出其他的所有解?

  • 扩展欧几里得的三大应用


一、裴蜀定理内容即证明

内容

设a, b是不全为零的整数,则存在整数x, y 使得 ax + by = gcd(a, b).

证明:

  1. a, b中假如有一个为0,则上述定理完全吻合,如a = 0, 则gcd(a, b) = b , 显然成立

  2. 若a, b 都不等于0

    gcd(a, b) = d,则d != 0

    等式两边同时除以d,得:

    a1x + b1y = 1

    再利用辗转相除法:gcd(a, b) --> gcd(b, a % b) --> ……

    设模出来的余数叫做x, 则有:

    gcd(a1, b1) = gcd(b1, r1) = gcd(r1, r2) = …… = gcd(rn-1, rn) = 1

    把辗转相除法中的运算展开,做成带余数的除法,得

    • a1 = x1b + r1

    • b1 = x2r1 + r2

    • r1 = x3r2 + r3

      ……

    • rn-3 = xn-1rn-2 + rn-1

    • rn-2 = xnrn-1 + rn

    • rn-1 = xn+1rn

    不妨令辗转相除法在除到互质的时候退出则 rn = 1 , 由倒数第二行得

    • rn-2 = xnrn-1 + 1 —> 1 = rn-2 - xnrn-1

    由倒数第三行得:**rn-1 = rn-3 - xn-1rn-2**带入上式

    得:1 = (1 + xnxn-1)rn-2 - xnrn-3

    采取同样的方法,逐个消去rn-2,……, r1

    最终会得到1 = f(x)a1 + g(x)b1

    其中f(x)和g(x) 是关于x的函数

    设f(x) = x, g(x) = y,则可以得到1 = a1x + b1y, 得证!

    正常人都不喜欢看上面写的奇奇怪怪的推导,那就看下面我这个具体的推导叭,这里是假设到 r4 = 1,然后往前推导滴

    在这里插入图片描述

二、求特解

我们都知道,欧几里得公式可以由这个式子表达

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

通过这个式子,我们可以不断递推到b = 0, 此时a即为a和b的最大公约数

现在我们对这个式子进行展开,看看有什么好玩的嘛

gcd(a, b) = a * x1 + b * x2

gcd(b, a % b) = b * x2 + (a % b) * y2

其中a % b = a - (a / b) * b

故由第一行的欧几里得公式得:

a * x1 + b * y1 = b * x2 + {a - (a / b) * b}y2

化简得a * x1 + b * y1 = a * y2 + b * {x2 - (a / b) * y2}

根据待定系数法得到

  • x1 = y2

  • y1 = x2 - (a / b) * y2

也就是说,如果有了gcd(b, a % b)的解x2, y2就可以推出gcd(a, b)的解x1, y1

那么这就和求gcd的过程一样,一直推到最后x = 1, y = 0,就可以回溯回去,得到gcd(a, b)的特解x1, y1

特解其实是最小的一个解

二、怎么利用特解推出其他所有整数解

先说结论:

x = x 0 + k ∗ b g c d ( a , b ) x = x0 + k * \frac{b}{gcd(a,b)} x=x0+kgcd(a,b)b

y = y 0 − k ∗ a g c d ( a , b ) y = y0 - k * \frac{a}{gcd(a,b)} y=y0kgcd(a,b)a

任意的x + y = x0 + y0

x0, y0就是方程的特解

对于a * x + b * y = g来说,让x增加b/gcd,让y增加a / gcd

这样就相当于加a * b / gcd,减去a * b / gcd,一加一减,最终不变

三、扩展欧几里得三大应用

  • 用扩展欧几里得算法解不定方程ax+by=c

如果 c % gcd(a, b) == 0,即c 是 gcd(a, b)的整数倍,则方程有解,且解就在上面的基础上乘以 c g c d ( a , b ) \frac{c}{gcd(a, b)} gcd(a,b)c即可

void e_gcd(int a, int b, int &gcd, int &x, int &y)
{
    if (b == 0)
    {
        x = 1;
        y = 0;
        gcd = a;
    }
    else
    {
        e_gcd(b, a % b, gcd, y, x);
        y -= x * (a / b);
    }
}

int main(){
    int a, b, x, y, gcd;
  	cin>>a>>b;
    e_gcd(a, b, gcd, x, y);
    cout<<x<<' '<<y<<endl;
    cout<<gcd<<endl;
    return 0;
}
  • 用扩展欧几里得算法求解模线性方程ax≡b (mod n)

    对于模线性方程ax≡b (mod n)可以化简为ax + ny = b,设d = gcd(a, n) 当且仅当b % d == 0时有解,且有d个解

    b % d == 0时,设ax + ny = d,的特解为x,y

    等式两边同时乘以 b d \frac{b}{d} db,则方程变成 a x b d + n y b d = b ax\frac{b}{d} + ny\frac{b}{d} = b axdb+nydb=b

    则原方程ax + ny = b的解为** x ′ = x ∗ b d , y ′ = y ∗ b d x' = x*\frac{b}{d},y' = y*\frac{b}{d} x=xdb,y=ydb**

    故**ax≡b (mod n)**的特解为 x0 = x * (b / d) % n

    d个解为xi= x0 + i*(n/d) mod n,{i = 0……n-1}

    解的间隔为n/d

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

bool modular_linear_equation(int a, int b, int n){
    int x, y, x0, gcd;
    exgcd(a, n, gcd, x, y);
    int d = gcd;//d为gcd(a, n)
    if(b % d)return false;//无解
    x0 = x * (b/d) % n;//特解
    for(int i = 0; i < d; ++i)cout<<x0 + i * (n / d)<<endl;//所有的解
    return true;
}
  • 用欧几里德算法求乘法逆元:

什么叫乘法逆元?

ax ≡ 1 (mod p)

这里我们称x为a关于p的逆元

上述式子可以还原为:ax + py = 1,根据上面讲的扩展欧几里得算法,我们知道gcd(a, p) = 1时才有解,通过扩展欧几里得算法得到的解x0,利用x0%p就可以得到最小解

为什么呢?

因为通解为x = x0 + p * k

那么,也就是说, a 关于 p 的逆元是一个关于 p 同余的,那么根据最小整数原理,一定存在一个最小的正整数,它是 a 关于p 的逆元,而最小的肯定是在(0 , p)之间的,而且只有一个,这就好解释了。

但是,由于问题的特殊性,有时候我们得到的特解 x0 是一个负数,还有的时候我们的 p 也是一个负数这怎么办?

当 p 是负数的时候,我们取 p 的绝对值就行了,当 x0 是负数的时候,x0% p 的结果仍然是负数(在计算机计算的结果上是这样的,虽然定义的时候不是这样的),这时候,我们仍然让 x0 对abs§ 取模,然后结果再加上abs§ 就行了,于是,我们不难写出下面的代码求解一个数 a 对于另一个数 p 的乘法逆元:

int a, b, x, y, gcd, ans, p;

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

inline int cal(int a, int p){
    exgcd(a, p, gcd, x, y);
    if(1 % gcd != 0)return -1;
    x = x * 1 / gcd;
    p = abs(p);
    ans = x % p;
    if(ans <= 0)ans += p;
    return ans;
}
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值