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

一、扩展欧几里得算法

裴蜀（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 = x₀ + k*(b/d)

y = y₀ - k*(b/d)

其中k可以是任意一个整数，x_0,y₀是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的通解均是上述形式：

证明：

　　ax₀ + by₀ = d

　　ax' + by' =d

　　作差可得：a(x₀ - x') + b(y₀ - y') = 0

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

　　移项得，a/d(x₀ - x') = - b/d(y₀ - y')

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

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

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

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

　　证毕！

二、线性同余方程

解法原理

给定整数a,b,m，求一个整数x满足a∗x≡b(mod m)，或者给出无解。

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

　　下面给出证明：

证明：

　　设ax mod m = b mod m = r

　　mk₁ + r =ax

　　mk₂ + r =b

　　其中k₁,k₂都是整数，

　　作差得：m(k₁ - k₂) = 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\n",x);
    }
    return 0;
}