“扩展欧几里得”及“线性同余方程”

“扩展欧几里得”及“线性同余方程”

五月 20, 2020

“扩展欧几里得”及“线性同余方程”

欧几里得算法

有两个数a , b,我们要求gcd(a,b),怎么做?枚举因子显然过于笨重,那该怎么做?

欧几里得有个十分有用的定理:gcd(a,b)=gcd(b,a%b)gcd(a,b)=gcd(b,a\%b) 这里省略证明

这样我们可以十分迅速的求解

1
2
3
4
5
long long gcd(long long a,long long b)
{
if(b == 0)return a;
else return gcd(b,a%b);
}

那啥玩意是扩展欧几里得呢???


扩展欧几里得

现在我们知道a , b的最大公约数是 gcd,那么利用裴蜀定理,我们可以得出:一定可以找到 x , y使得:ax+by=gcda*x + b*y = gcd,现在我们只要找到一组特解x0,y0,就可以用其表示出整个方程的通解:

x=x0+(b/gcd)ty=y0(a/gcd)ttZx=x_0+(b/gcd)*t \\y=y_0-(a/gcd)*t\\ t\in Z

那问题转化到,如何求特解x0和y0了。

观察欧几里得算法,易得出其停止条件为a=gcd , b=0a=gcd\space ,\space b = 0,此时为最终状态,这时我们会有a1+b0=gcda*1+b*0=gcd

我们是不是可以从最终状态反推到最初状态呢?


推导

假设初始状态:ax+by=gcda*x+b*y=gcd,下一组解为bx1+(a%b)y1=gcdb*x_1+(a\%b)*y_1=gcd

由欧几里得定理:gcd(a,b)=gcd(b,a%b)gcd(a,b)=gcd(b,a\%b)

易得:ax+by=bx1+(a%b)y1a*x+b*y=b*x_1+(a\%b)*y_1

我们可以推出a%b=a(a/b)ba\%b = a-(a/b)*b(“/”为整除向下取整)

所以可以得出:

bx1+(a%b)y1=bx1+(a(a/b)b)y1=bx1+ay1(a/b)by1=ay1+b(x1(a/b)y1)b*x_1+(a\%b)*y_1=b*x_1+(a-(a/b)*b)*y_1\\ =b*x_1+a*y_1-(a/b)*b*y_1\\ =a*y_1+b*(x_1-(a/b)*y_1)

即:ax+by=ay1+b(x1(a/b)y1)a*x+b*y=a*y_1+b*(x_1-(a/b)*y_1)

即:

x=y1y=x1(a/b)y1x = y_1\\ y = x_1 - (a/b)*y_1

递归即可求出gcd,和通解x、y。


代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
int exgcd(int a,int b,int &x,int &y)
{
if(b == 0)
{
x = 1;
y = 0;
return a;
}
int gcdd = exgcd(b,a%b,x,y);
int t = x;
x = y;
y = t - (a/b)*y;
return gcdd;
}

我们可以看出,y = t - (a/b)*y 中,t是上次的y,y是上次的x,所以可以简化一下

1
2
int gcdd = exgcd(b,a%b,y,x);
y-=(a/b)*x;

扩展欧几里得的用途

  • 求解形如 ax +by = c 的通解,但是一般没有谁会无聊到让你写出一串通解出来,都是让你在通解中选出一些特殊的解,比如一个数对于另一个数的乘法逆元

乘法逆元(在维基百科中也叫倒数,当然是 mod p后的,其实就是倒数不是吗?):

如果(a*x)%b≡1,且gcd(a,b)=1(a与b互质),则称a关于模b的乘法逆元为x。

因为:ax%b=1a*x\%b=1

根据:a%b=a(a/b)ba\%b = a-(a/b)*b

可以得出:ax%b=axb(ax/b)=ax+b(ax/b)a*x\%b=a*x-b*(a*x/b)=a*x+b*(-a*x/b)

ax/b=y-a*x/b=y,则有ax+py=1a*x+p*y=1

利用扩展欧几里得可以求出x,且只有gcd(a,b)=1时,才存在逆元x。

  • 用来求线性同余方程

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
#include<iostream>
using namespace std;

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

int move_reverse(int a,int b)
{
int d, x, y;
d = exgcd(a, b, x, y);
if(d != 1)
return -1;
else
return (x%b+b)%b;
}

int main()
{
int a, b;
cin >> a >> b;
cout << move_reverse(a, b) << endl;
return 0;
}

洛谷例题


原理

  • 知道初始条件
  • 当前状态的答案可以从下一状态转移过来
  • 知道最后的状态

线性同余方程

定义

对于形如

axc (mod b):(ax)%b=ca*x≡c \space (mod\space b)\\ 即:(a*x)\%b=c

的方程,称为线性同余方程(Congruence Equation)


性质

  • 方程ax+by=ca*x+b*y=c与方程axc  (mod b)a*x≡c\space\space(mod\space b)是等价的,有整数解的充要条件是gcd(a,b)cgcd(a,b)|c

  • 若方程有整数解,则共有gcd(a,b)gcd(a,b)个解

  • gcd(a,b)=dgcd(a,b)=d,则线性同余方程在[0,b/d1]\left[0,b/d-1\right]上有唯一解。

  • gcd(a,b)=1gcd(a,b)=1,且x0x_0y0y_0为方程ax+by=ca*x+b*y=c的一组解

    则方程的通解为

    x=x0+tby=y0tatZx=x_0+t*b\\ y=y_0-t*a\\ t\in Z


求解方法

根据第一条性质,方程ax+by=ca*x+b*y=c,我们先用扩展欧几里得求出一组x0x_0y0y_0

也就是ax0+by0=gcd(a,b)a*x_0+b*y_0=gcd(a,b)

然后两边同时除以gcd(a,b)gcd(a,b),再乘cc

就得到了方程:acx0/gcd(a,b)+bcy0/gcd(a,b)=ca*c*x_0/gcd(a,b)+b*c*y_0/gcd(a,b)=c

这时我们就找到了方程的一个解。

根据第四条性质,可以求出方程的所有解,但在实际问题中,我们往往被要求求出一个最小整数解,也就是一个特解xx

t=b/gcd(a,b)x=(x%t+t)%tt=b/gcd(a,b)\\ x=(x\%t+t)\%t

此时的yy可以用xx推出:

y=(c/d(a/d)x)/(b/d)y=(c/d-(a/d)*x)/(b/d)

这个x=(x%t+t)%t是为了让x最小


代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
#include<iostream>
#include<cstdio>
#include<cmath>
using namespace std;

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

void CongruenceEquation(int a, int b, int c)
{
int x;
int y;
int d = exgcd(a, b, x, y);
if(c % d != 0)
return;
//求出x推y
x = c * x / d;//乘c除以d得到方程的一个解
int t = b / d;
int min_x = (x % t + t) % t;
int min_y = (c / d - (a / d) * min_x) / (b / d);
cout << min_x << " " << abs(min_y) << endl;
//求出y推x
y = c * y / d;
t = a / d;
min_y = (y % t + t) % t;
min_x = (c / d - (b / d) * min_y) / (a / d);
cout << abs(min_x) << " " << min_y << endl;
}

int main()
{
//(a*x)%b=n
//a*x=n (mod b)
//a*x+b*y=c
int a,b,c;
cin >> a >> b >> c;
CongruenceEquation(a, b, c);
return 0;
}

线性同余方程组

{xa1 (mod m1)xa2 (mod m2)xan (mod mn)\begin{cases} x ≡ a_1 \space (mod\space m_1)\\ x ≡ a_2 \space (mod \space m_2)\\ \vdots\\ x ≡ a_n \space (mod \space m_n)\\ \end{cases}\\

其中m1,m2m_1,m_2\dots\dots不保证互质。

求最小非负整数解XX

即若干个线性同余方程的组合。

求解方式可以参考我们普通的方程组,先解式1,将式1的通解带入式2后化简,再解式2,如此重复即可得满足所有式子的解。

例如