题目描述: 在RSA密码体系中,欧几里得算法是加密或解密运算的重要组成部分。它的基本运算过程就是解 (x*a) % n = 1 这种方程。 其中,x,a,n皆为正整数。现在给你a和n的值(1 < a,n < 140000000),请你求出最小的满足方程的正整数解x(保证有解). 如:a = 1001, n = 3837,则输出23。
示例: 输入: a = 1001 n = 3837
输出:23
分析: 此题使用扩展欧几里得算法求解。
欧几里得算法: 欧几里德算法又称辗转相除法,用于计算两个整数a,b的最大公约数 gcd(a,b)。基本算法:设 a = qb + r,其中a,b,q,r都是整数,则
gcd(a, b) = gcd(b, r),即 gcd(a, b) = gcd(b, a%b)。
证明: a = qb + r
推理1: 如果 r = 0,那么 a 是 b 的倍数,此时显然 b 是 a 和 b 的最大公约数。
推理2: 如果 r ≠ 0,任何整除 a 和 b 的数必定整除 a - qb = r,而且任何同时整除 b 和 r 的数必定整除 qb + r = a,所以 a 和 b 的公约数集合与 b 和r 的公约数集合是相同的。特别的,a 和 b 的最大公约数是相同的。
例如: a = i * k, b = j * k
那么 r = i * k - j * q * k = (i - j*q) * k
a, b, r 都可以整除 k 。
python 参考代码如下:
1
2
3
4
5
6
7
8
def gcd(a, b):
return a if b == 0 else gcd(b, a%b)
def gcd(a, b):
if b == 0:
return a
else:
return gcd(b, a%b)
扩展欧几里得算法: 在求 a , b 的最大公约数 m = g c d ( a , b ) 的同时,求出贝祖等式 ax + by = m的一个解 ( x , y ) 。
对贝祖等式有兴趣的可以百度了解一下。
ax + by = gcd(a, b) ——- (1)
bx + a%by = gcd(b, a%b) —–(2)
a%b = a - a//b * b —–(3)
gcd(a, b) == gcd(b, a%b) —-(4)
所以就有: \(ax + by = bx + a\%b y = bx + (a-a//b*b)y = bx + ay - a//b*by = ay + b(x-a//b*y)\) 于是就有:
x = y
y = x - a // b*y
递归版参考代码:可以和上面的 朴素欧几里得算法对比学习。
1
2
3
4
5
6
7
8
9
10
11
12
# arr[0] 相当于 x, arr[1] 相当于 y
def exgcd(a, b, arr):
# 推理1 终止条件
if b == 0:
arr[0] = arr[1] = 1, 0
return a
r = exgcd(b, a%b, arr[0], arr[1])
t = arr[1]
arr[1] = arr[0] - (a//b) * arr[1]
arr[0] = t
return r
真正难的是 非递归扩展欧几里得算法,这个才是真的有点点难度。需要消耗一点点的思考吧。
先给出算法步骤:
- [余数为0?] : 如果 r = 0 ,算法终止,此时有
am + bn = d
, 正如所求.
假定 m = 1769 , n = 551, 我们相继得到:
a’ | a | b’ | b | c | d | q | r |
---|---|---|---|---|---|---|---|
1 | 0 | 0 | 1 | 1769 | 551 | 3 | 116 |
0 | 1 | 1 | -3 | 551 | 116 | 4 | 87 |
1 | -4 | -3 | 13 | 116 | 87 | 1 | 29 |
-4 | 5 | 13 | -16 | 87 | 29 | 3 | 0 |
我们可以算出 gcd(1769, 551) = 29
.
现在我们已经了解了算法的步骤,现在再来分析一下,此算法为什么是对的?
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
执行第 1 步后: a'= b = 1, a = b' = 0, c = m, d = n;
执行 第 4 步以前:
a'm + b'n = c -------- 等式(1)
am + bn = d ---------等式(2)
a'm + b'n = qd+r ----- c = qd + r
a'm + b'n - q(am+bn) = r ----- d = am + bn
(a'-qa)m + (b'-qb)n = r --------等式(3)
执行第 4 步: c=d, d=r, t=a', a'=a, a=t-qa, t=b', b'=b, b=t-qb
等式(2) 变成 a'm + b'n = c
等式(3) 变成 am + bn = d
(PS: 这一步好好的看看,里面各个变量转换的顺序不能错。)
所以最后我们又变成了等式(1)和等式(2),然后就又可以了新一步的循环了。
为了得到 等式(1) 和 等式 (2) 我们需要当 a' = b = 1, a = b' = 0, c = m, d = n 这样才符合要求。
为什么要得到等式(1) 和 等式(2) 呢?
因为 gcd(m, 0) = m, gcd(0, n) = n。
也就是说:
根据gcd(m, 0) = m, gcd(0, n) = n 这个原理,
我们初始化 a' = b = 1, a = b' = 0, c = m, d = n
得到 等式(1) 和 等式(2), 再根据 c = qd+r, d = am+bn 公式,得到等式 (3)。
为了进行下一组循环,使用第4步,又得到等式(1)和等式(2)。
非递归扩展欧几里得算法的代码可以看下面这个代码:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
def extgcd(m, n):
# 用 ta, tb 表示 a' 和 b'
ta = b = 1
a = tb = 0
c, d = m, n
q, r = c//d, c%d
while r:
taa = ta - q * a
tbb = tb - q * b
ta, tb = a, b
a, b = taa, tbb
c, d = d, r
q, r = c//d, c%d
# print(ta, tb, a, b, c, d, q, r)
return d, a, b
现在我们可以得到 ax + by = m 的一组解了。 但是对于我们的题目 (x * a) % n = 1 而言还有一点距离。
又 因为 a%n = a - a//n * n 令 a // n = y 则 a % n = a - y*n, 所以推导出 (x * a) % n = x * a - y * n = 1
又 k 可以为负数或者是正数,上面的公式又可以转变为 x*a + y * n = 1 , 这个时候,就和我们的 ax+by=m, 有一点点像了。
所以对于我们这个题目,先通过扩展欧几里得算法得到一组解。
假设现在得到的一组解为 (x0, b0)
因为要求最小正整数解,然而我们求出来的那一组解,可能不是最小的正整数解。于是就需要变化一下。
那么一直一组解如何求另外一组解呢?
公式:
x = x0 - b*t
y = y0 + a*t
其中 t 为整数。
我们的目的是得到下一组的解,现在我们来简单验证一下,这个公式为啥可以?
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
ax0 + by0 = 1
带入公式:
a(x0-n*t) + b(y0+n*t) = 1
变形:
ax0 - a*b*t + by0 + b*a*t = 1
合并公因式:
ax0 + by0 + (b*a*t - a*b*t) = 1
结果:
ax0 + by0 = 1
这时 x0 和 y0 的值,随着 t 的变化而变化,t可以是负整数,也可以是正整数。
代码:
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
def extgcd(m, n):
# 用 ta, tb 表示 a' 和 b'
ta = b = 1
a = tb = 0
c, d = m, n
q, r = c//d, c%d
while r:
taa = ta - q * a
tbb = tb - q * b
ta, tb = a, b
a, b = taa, tbb
c, d = d, r
q, r = c//d, c%d
# print(ta, tb, a, b, c, d, q, r)
return d, a, b
_, x0, y0 = extgcd(a, n)
t = 1
while x0 <= 0:
x0 = x0 + n*t
y0 = y0 - a*t
t += 1
t=1
while True:
if (x0 - n * t) <= 0:
break
x0 = x0 - n * t
y0 = y0 + a * t
t += 1
print(x0)