数论——扩展欧几里得算法与线性同余方程

数论——扩展欧几里得算法与线性同余方程

一、扩展欧几里得算法

裴蜀(Bézout)定理

对任何整数a、b和它们的最大公约数论——扩展欧几里得算法与线性同余方程数d=gcd(a,b),一定存在整数x,y,使ax+by=d成立。 
  下面证明裴蜀定理:

证明:

  在欧几里得算法的最后一步,即b=0时,显然有一对整数x = 1,y = 0,使得a*1+0*0=gcd(a,0)。

  当b>0时,

  ∵ by + (a mod b) x = d,即 by + (a - a/b * b) x = d(注:a/b下取整)。

  ∴ 整理得 ax + b(y - a/b * x)=d。

  令x' = x,y' = y,于是 ax' + by' = d。所以x'和y'就是满足条件的一组解。

  由欧几里得算法递归过程及数学归纳法,可知裴蜀定理成立。

  证毕!

  由上述证明过程可得整数x和y的计算方法,这种计算方法被称为扩展欧几里得算法

代码实现

  模板题链接:扩展欧几里得算法

  代码如下:

int exgcd(int a,int b,int &x,int &y)
{
    if(!b)
    {
        x=1;y=0;
        return a;
    }
    int d=exgcd(b,a%b,y,x);
    y-=a/b*x;
    return d;
}

ax+by=d的通解

x = x0 + k*(b/d)

y = y0 - k*(b/d)

其中k可以是任意一个整数,x0,y0是ax+by=d的一组特殊解。

  下面证明上式的正确性:

  ax + by = a(x0 + k(b/d)) + b(y0 - k(a/d)) = ax0 +by0 + k(ab/d) - k(ab/d) =d。

  下面证明ax+by=d的通解均是上述形式:

证明:

  ax0 + by0 = d

  ax' + by' =d

  作差可得:a(x0 - x') + b(y0 - y') = 0

  等式两边同除以d得,a/d(x0 - x') + b/d(y0 - y') = 0

  移项得,a/d(x0 - x') = - b/d(y0 - y')

  因为d是a和b的最大公约数,所以(a/d)与(b/d)互质。

  又因为 (b/d) | a/b(x0 - x'),所以 (b/d) | (x0 - x')。

  所以x0 - x' = k * (b/d),即x' = x0 - k * (b/d)。

  由于k可以是任意的整数,所以上式等价于x' = x0 + k * (b/d)。

  证毕!

二、线性同余方程

解法原理

给定整数a,b,m,或者给出无解。

  a * x ≡ b(mod m),等价于 ax + my = b(b一定是gcd(a,m)的倍数)。

  下面给出证明:

证明:

  设ax mod m = b mod m = r

  mk1 + r =ax

  mk2 + r =b

  其中k1,k2都是整数,

  作差得:m(k1 - k2) = ax - b

  所以 ax = mk +b

  证毕!

  根据扩展欧几里得算法,可以得到 ax' + my' = d(d = gcd(a,m))的一组解x'和y'。

  那么在其等式两边同乘上(b/d)即可得:a((b/d)x') + m((b/d)y') = d * (b/d) = b。

  于是我们便得到了线性同余方程的一组解即为x = (b/d)x'和y = (b/d)y'。

  因为只要求x,所以只需返回x即可。

代码实现

  模板题链接:线性同余方程

  根据上述原理,我们只需要用欧几里得算法得出一组解之后,再将x乘上(b/d)即可。

  代码如下:

#include <iostream>
#include <algorithm>
#include <cmath>
#include <cstdio>
using namespace std;
int exgcd(int a,int b,int &x,int &y)
{
    if(!b)
    {
        x=1;y=0;
        return a;
    }
    int d=exgcd(b,a%b,y,x);
    y-=a/b*x;
    return d;
}
int main()
{
    int n;scanf("%d",&n);
    for(int i=1;i<=n;i++)
    {
        int a,b,m;scanf("%d%d%d",&a,&b,&m);
        int x,y;
        int d=exgcd(a,m,x,y);
        if(b%d)
        {
            puts("impossible");
            continue;
        }
        x=(long long)x * (b/d) % m;
        printf("%d
",x);
    }
    return 0;
}