RSA加解密

本文大部分参考自阮一峰的RSA算法原理(一) 、RSA算法原理(二)

原文最后的推导部分有些地方一笔带过，自己重新加了注解，方便理解。文章末尾增加了快速幂和扩展欧几里得算法求 ax + by = 1的推导。送给有需要的人。

1. 产生背景

以前的加密方法都是对称加密算法：加密和解密的双方都使用同一种规则(密钥)，这决定了加密方必须把加密密钥告诉解密方，否则无法解密。密钥的保存和传递成了一个难题。

1976年，有人提出了"非对称加密算法"：

1. 乙方生成两把密钥（公钥和私钥）。公钥是公开的，任何人都可以获得，私钥则是保密的。
2. 甲方获取乙方的公钥并用它对信息加密。
3. 乙方得到加密后的信息，用私钥解密。

如果公钥加密的信息只有对应的私钥才能解开，那么只要私钥不泄露，通信就是安全的。

1977年，三位数学家Rivest、Shamir 和 Adleman 设计了一种算法，实现了非对称加密。该算法用他们三人的名字命名，叫做RSA算法。直到现在，RSA算法一直是最广为使用的"非对称加密算法"。这种算法非常可靠，密钥越长，它就越难破解。如今4096bit的公钥长度是很安全的。

2.数学推导

2.1 质数和互质

• 质数：又称为素数，是一个大于1的自然数，除了1和它自身外，不能被其他自然数整除
• 合数：比1大但不是素数的数称为合数。 注：1和0既非质数也非合数
• 互质关系：**如果两个正整数，除了1以外，没有其他公因子，我们就称这两个数是互质关系。

注：不是质数也可以构成互质关系**。比如：15和32。

由互质关系不难得出以下结论：

1. 任意两个质数构成互质关系，比如13和31。 ----因为质数互不相同，公因子只有1

2. 一个数是质数，另一个数只要不是前者的倍数，两者就构成互质关系，比如3和10。 ----因为只有后一个数是前一个质数的倍数时，两者才不构成互质关系

3. 如果两个数中较大的数是质数，则两者构成互质关系，比如97和57。----较大的为质数，则不可能是较小数(1除外)的倍数，那么两者构成互质关系

4. 1和任意一个自然数是都是互质关系，比如1和99。 ----由互质关系的定义得出

5. p是大于1的整数，则p和p-1构成互质关系，比如2和1，57和56。----p> 1,最小为2；即从1开始，相邻的两整数构成互质

6. p是大于1的奇数，则p和p-2构成互质关系，比如17和15。 ----p为奇数，p最小为3，即p-2最小为1,即相邻的奇数构成互质

2.2 欧拉函数

对于任意正整数n，求在小于等于n的正整数之中有多少个与n构成互质关系？ 计算这个值的方法就叫做欧拉函数，以φ(n)表示。 比如在1到8之中，与8形成互质关系的正整数有1、3、5、7，故 φ(n) = 4。

φ(n) 的计算方法并不复杂，但是为了得到最后那个公式，需要一步步讨论。

• 情况1: n=1，则 φ(1) = 1 。因为1与任何数（包括自身）都构成互质关系。
• 情况2：n是质数，则 φ(n)=n-1 。因为质数与小于它的每一个数都构成互质关系(由互质关系得出的第3点结论)。 比如5与1、2、3、4都构成互质关系。
• 情况3: n是质数的某一个次方，即 n = $p^k$ (p为质数，k为大于等于1的整数)，则：

比如 φ(8) = φ($2^3$) =$2^3$ - $2^2$ = 8 -4 = 4。

因为从1开始、小于等于$p^k$的范围内： 所有的含p的整数都与$p^k$不构成互质关系：

p、2p、3p... p^k^

提取p，变为:

1、2、3、... p^k-1^

这些数的总数为：p^k-1^

把它们去除，剩下的就是与n互质的数。即与n互质的总数为

$p^k$ - $p^k-1^$

上面的式子还可以写成下面的形式：

可以看出，上面的第二种情况是 k=1 时的特例。

• 情况4: n可分解成两个互质的整数之积：

　　n = p1 × p2

则

　　φ(n) = φ(p1p2) = φ(p1)φ(p2)

即积的欧拉函数等于各个因子的欧拉函数之积。比如，φ(56)=φ(8×7)=φ(8)×φ(7)=4×6=24。

这一条的证明要用到"中国剩余定理"，这里就不展开了，只简单说一下思路：

如果a与p1互质(a<p1)，b与p2互质(b<p2)，c与p1p2互质(c<p1p2)，则c与数对 (a,b) 是一一对应关系。

由于a的值有φ(p1)种可能，b的值有φ(p2)种可能，则数对 (a,b) 有φ(p1)φ(p2)种可能，而c的值有φ(p1p2)种可能，所以φ(p1p2)就等于φ(p1)φ(p2)。

• 情况5:因为任意一个大于1的正整数，都可以写成一系列质数的积。

根据第4条的结论，得到

再根据第3条的结论，得到

也就等于

这就是欧拉函数的通用计算公式。

比如，1323的欧拉函数，计算过程如下：

2.3 欧拉定理

欧拉函数的用处在于欧拉定理。"欧拉定理"指的是：

如果两个正整数a和n互质，则n的欧拉函数 φ(n) 可以让下面的等式成立：

也就是说，a的φ(n)次方被n除的余数为1。或者说，a的φ(n)次方减去1，可以被n整除。

比如，3和7互质，而7的欧拉函数φ(7)等于6，所以3的6次方（729）减去1，可以被7整除（728/7=104）。

欧拉定理的证明比较复杂，这里就省略了。我们只要记住它的结论就行了。

欧拉定理可以大大简化某些运算。比如，7和10互质，根据欧拉定理：

已知 φ(10)等于4 (4个数为：1、3、7、9)，所以马上得到7的4倍数次方的个位数肯定是1。

因此，7的任意次方的个位数（例如7的222次方），心算就可以算出来。---- 7^220^的个位数为1，故7^222^的个位数为9。

欧拉定理有一个特殊情况。

假设正整数a与质数p互质，因为质数p的φ(p)等于p-1，则欧拉定理可以写成

这就是著名的费马小定理。它是欧拉定理的特例。 ----区别： 欧拉定理强调两个数是互质关系即可；而费马小定理要求其中一个为质数。

欧拉定理是RSA算法的核心。理解了这个定理，就可以理解RSA。

2.4 模反元素

如果两个正整数a和n互质，那么一定可以找到整数b，使得 ab-1 被n整除，或者说ab被n除的余数是1。

这时，b就叫做a的"模反元素"。

比如，3和11互质，那么3的模反元素就是4，因为 (3 × 4)-1 可以被11整除。显然，模反元素不止一个， 4加减11的整数倍都是3的模反元素 {...,-18,-7,4,15,26,...}，即如果b是a的模反元素，则 b+kn 都是a的模反元素。

欧拉定理可以用来证明模反元素必然存在。

可以看到: a的 φ(n)-1 次方，就是a的模反元素。

3. 密钥生成的步骤

通过一个例子来理解RSA算法。假设爱丽丝要与鲍勃进行加密通信，她该怎么生成公钥和私钥呢？

3.1. 随机选择两个不相等的质数p和q。

爱丽丝选择了61和53。（实际应用中，这两个质数越大，就越难破解。）

3.2. 计算p和q的乘积n。

n = 61 * 53 = 3233

n的长度就是密钥长度。3233写成二进制是1100 1010 0001，一共有12位，所以这个密钥就是12位。实际应用中，RSA密钥一般是1024位，重要场合则为2048位。

3.3. 计算n的欧拉函数φ(n)。

根据公式：φ(n) = φ(p*q) = φ(p)*φ(q) = (p-1)(q-1)

爱丽丝算出φ(3233)等于60×52，即3120。

3.4. 随机选择一个整数e，条件是1< e < φ(n)，且e与φ(n) 互质。

爱丽丝就在1到3120之间，随机选择了17。（实际应用中，常常选择65537。）

3.5. 计算e对于φ(n)的模反元素d。

所谓"模反元素"就是指有一个整数d，可以使得 e*d 被φ(n) 除的余数为1, 即

ed ≡ 1(mod φ(n))

这个式子等价于

ed = kφ(n) + 1

于是，找到模反元素d，实质上就是对下面这个二元一次方程求解。

ed - kφ(n) = 1

已知 e=17, φ(n)=3120，

17*d - 3120*k = 1

这个方程可以用"扩展欧几里得算法"求解，此处省略具体过程。总之，爱丽丝算出一组整数解为 (x,y)=(2753,15)，即 d=2753。

至此所有计算完成。

3.6. 将n和e封装成公钥，n和d封装成私钥。

在爱丽丝的例子中，n=3233，e=17，d=2753，所以公钥(n,e)就是 (3233,17)，私钥(n, d)就是（3233, 2753）。

实际应用中，公钥和私钥的数据都采用ASN.1格式表达。

3.7. RSA算法的可靠性

回顾上面的密钥生成步骤，一共出现六个数字：

   p
　 q
　 n
　 φ(n)
　 e
　 d

这六个数字之中，公钥用到了两个（n和e），其余四个数字都是不公开的。其中最关键的是d，因为n和d组成了私钥，一旦d泄漏，就等于私钥泄漏。

那么，有无可能在已知n和e的情况下，推导出d？

（1）ed≡1 (mod φ(n))。只有知道e和φ(n)，才能算出d。

（2）φ(n)=(p-1)(q-1)。只有知道p和q，才能算出φ(n)。

（3）n=pq。只有将n因数分解，才能算出p和q。

结论：如果n可以被因数分解，d就可以算出，也就意味着私钥被破解。

可是，大整数的因数分解，是一件非常困难的事情。目前，除了暴力破解，还没有发现别的有效方法。维基百科这样写道：

  "对极大整数做因数分解的难度决定了RSA算法的可靠性。换言之，对一极大整数做因数分解愈困难，RSA算法愈可靠。

　　假如有人找到一种快速因数分解的算法，那么RSA的可靠性就会极度下降。但找到这样的算法的可能性是非常小的。今天只有短的RSA密钥才可能被暴力破解。到2008年为止，世界上还没有任何可靠的攻击RSA算法的方式。

  只要密钥长度足够长，用RSA加密的信息实际上是不能被解破的。"

举例来说，你可以对3233进行因数分解（61×53），但是你没法对下面这个整数进行因数分解。

  12301866845301177551304949
  58384962720772853569595334
  79219732245215172640050726
  36575187452021997864693899
  56474942774063845925192557
  32630345373154826850791702
  61221429134616704292143116
  02221240479274737794080665
  351419597459856902143413

它等于这样两个质数的乘积：

  33478071698956898786044169
  84821269081770479498371376
  85689124313889828837938780
  02287614711652531743087737
  814467999489

　　　　×

  36746043666799590428244633
  79962795263227915816434308
  76426760322838157396665112
  79233373417143396810270092
  798736308917

事实上，这大概是人类已经分解的最大整数（232个十进制位，768个二进制位）。比它更大的因数分解，还没有被报道过，因此目前被破解的最长RSA密钥就是768位。

3.8. 加密和解密

有了公钥和密钥，就能进行加密和解密了。

3.8.1 加密要用公钥 (n,e)

假设鲍勃要向爱丽丝发送加密信息m，他就要用爱丽丝的公钥 (n,e) 对m进行加密。这里需要注意，m必须是整数（字符串可以取ascii值或unicode值），且m必须小于n。

所谓"加密"，就是算出下式的c：

m^e^ ≡ c (mod n)

爱丽丝的公钥是 (3233, 17)，鲍勃的m假设是65，那么可以算出下面的等式：

65^17^ ≡ 2790 (mod 3233)

于是，c等于2790，鲍勃就把2790发给了爱丽丝。

3.8.2 解密要用私钥(n,d)

爱丽丝拿到鲍勃发来的2790以后，就用自己的私钥(3233, 2753) 进行解密。可以证明，下面的等式一定成立：

c^d^ ≡ m (mod n)

也就是说，c的d次方除以n的余数为m。现在，c等于2790，私钥(n,d)是(3233, 2753)，那么，爱丽丝算出

2790^2753^ ≡ 65 (mod 3233)

因此，爱丽丝知道了鲍勃加密前的原文就是65。

至此，"加密--解密"的整个过程全部完成。

我们可以看到，如果不知道d，就没有办法从c求出m。而前面已经说过，要知道d就必须分解n，这是极难做到的，所以RSA算法保证了通信安全。

你可能会问，公钥(n,e) 只能加密小于n的整数m，那么如果要加密大于n的整数，该怎么办？有两种解决方法：

• 一种是把长信息分割成若干段短消息，每段分别加密；
• 另一种是先选择一种"对称性加密算法"（比如DES），用这种算法的密钥加密信息，再用RSA公钥加密DES密钥。

3.9 私钥解密的证明

最后，我们来证明，为什么用私钥解密，一定可以正确地得到m。也就是证明下面这个式子：

$c^d$ ≡ m (mod n)

因为，根据加密规则:

$ｍ^e$ ≡ c (mod n)

于是，c可以写成下面的形式：

c = $m^e$ - kn

将c代入到我们要证明的那个解密规则：

$(m^e - kn)^d$ ≡ m (mod n)

可以表示为：

$(m^e - kn)^d$ -hn = m

左边展开多项式：

$(m^e - kn)^d$

=$(m^e - kn)(m^e - kn)...(m^e - kn)$ // d个$(m^e - kn)$相乘

= m^ed^ + knm^e(d-1)^ + $(kn)^2$ * m^e(d-2)^ + ...+ $(kn)^d$ // 从第二项开始，都含有n，可以被n整除。

故：

$(m^e - kn)^d$ - hn

= $m^ed$ + knm^e(d-1)^ + $(kn)^2$ * m^e(d-2)^ + ...+ $(kn)^d$ - hn

= m^ed^ + gn //g为从第二项开始、提取n后的各项之和

所以求证：

$(m^e - kn)^d$ ≡ m (mod n)

等同求证：

m^ed^ ≡ m (mod n)

也可从取模公式上得到以上结论。数学上对取余(取模) 有如下等式：

(a + b) % p = (a % p + b % p) % p ---公式1
(a - b) % p = (a % p - b % p) % p
(a * b) % p = (a % p * b % p) % p
a ^ b % p = ((a % p)^b) % p

参考公式1，所以

(m^ed^ + kn * m^e(d-1)^ + $(kn)^2$ * m^e(d-2)^ + ...+ $(kn)^d$) (mod n)

= (m^ed^ (mod n) + kn * m^e(d-1)^ (mod n) + $(kn)^2$ * m^e(d-2)^ (mod n) + ... + $(kn)^d$ * (mod n)) (mod n) = $m^ed$ (mod n)

故求证：

$(m^e - kn)^d$ ≡ m (mod n)

等同求证：

m^ed^ ≡ m (mod n)

由于

ed ≡ 1(mod φ(n))

所以

ed = h φ(n) + 1

将ed代入, 即求证：

m^hφ(n)+1^ ≡ m (mod n)

接下来，分成两种情况证明上面这个式子。

（1）m与n互质。

根据欧拉定理，此时

m^φ(n)^ ≡ 1 (mod n)

所以

$m^φ(n)^ = kn + 1

所以

m^hφ(n)+1^
= m^hφ(n)^* m
=(kn + 1)^h^*m
=(kn + 1)(kn + 1)...(kn + 1)*m //h个(kn + 1)相乘
=（kn)^h^*m+(kn)^h-1^*m+ ... + (kn)*m + m

根据公式：

(a + b) % p = (a % p + b % p) % p

得到

((kn)^h^*m+(kn)^h-1^*m+ ... + (kn)*m + m) (mod n)
= ((kn)^h^*m (mod n) + (kn)^h-1^*m (mod n) + ... + (kn)*m (mod n) + m (mod n)) (mod n)
= m (mod n)

得证。

（2）m与n不是互质关系。

知识点：

两个不同质数的乘积有且只有四个因数：这两个质数、1和它本身

由于n等于质数p和q的乘积，可知 n的因数只有1、n、p、q。 m与n不互质，则双方含有非1公因数。

因 m < n, 所以m不能为n的倍数，故m = kp，k<q 或 m = kq, k<p。

以 m = kp为例，此时k<q , 可知k与q互质 ---- 由质数的推论3 "如果两个数中较大的数是质数，则两者构成互质关系"得出

又p和q互质，故 kp 与q互质。

注：

p,q互质, 且p,k互质, 那么kp 与q 为什么互质呢？

a、b互质的另一种定义它们的最大公约数gcd(a,b)=1,等价于存在整数u、v使得 ua+vb=1。
同时，a、c互质, gcd(a,c)=1, 等价于存在整数s、t使得sa+tc=1。
所以bc=(1-ua)/v * (1-sa)/t。
vtbc=(1-ua)(1-sa)=1-ua-sa+usa^2^。
a*(u+s-usa)+bc*vt=1。
所以gcd(a,bc)=1。

根据欧拉定理：

(kp)^φ(q)^ ≡ 1 (mod q)

即

(kp)^q-1^≡ 1 (mod q), 可得：(kp)^q-1^ = tq + 1， 即 m^q-1^ = tq + 1

又因 ed = hφ(n) + 1, 所以：

m^ed^
=m^hφ(n)+1^
=m^hφ(pq)+1^
=m^hφ(p)*φ(q)+1^
=m^h(p-1)(q-1)+1^
=(m^(q-1)^)^h(p-1)^ * m
=(tq + 1)^h(p-1)^ * m
=(tq + 1)(tq + 1)...(tq + 1)*m //h(p-1)个(tq + 1)相乘
=(tq)^h(p-1)^*m + (tq)^h(p-1)-1^*m + ... + (tq)*m + m
=(tq)^h(p-1)^*kp + (tq)^h(p-1)-1^*kp + ... + (tq)*kp + kp
=(tq)^h(p-1)-1^*tk(pq) + (tq)^h(p-1)-2^*tk(pq) + ...tk(pq) + kp
=(tq)^h(p-1)-1^*tkn + (tq)^h(p-1)-2^*tkn + ...tkn + kp
≡ m (mod n)

得证。

4. e为什么要默认取65537

公钥的组成是(n, e), 其中n的二进制位数就是公钥的长度。 假如我们的公钥采用目前最长的长度4096位。 早期的Openssl为了加密和签名的效率，默认的e取值3。

举个例子，前方战线请示总部是否开战，总部收到消息后电文速回No。如果e取3，会出现什么问题呢？

首先总部查看ASCII表，得出"No"的 ASCII 为 0x4E 0x6F， 所以得到明文 M = 0x4E6F.

接着总部进行机密, 密文 C = (0x4E6F)^e mod n , 计算 (0x4E6F)^3 = 0x75CCE07084F，即 C = 0x75CCE07084F mod n

n是一个4096位的大数，肯定大于0x75CCE07084F的，所以我们可以直接得到密文:

C = 0x75CCE07084F

这时敌军监听拦截到这个传输的密文0x75CCE07084F , 觉得很数字很短，有机会进行破解，于是：

敌军将0x75CCE07084F 转成10进制，得到8095174953039

再尝试着进行开根操作，当进行到开3次根的时候：

$\sqrt[3]{8095174953039}$ = 20079

接着将它转成十六进制 hex(20079)= 0x4E4F，

敌军赶紧查一下ASCII表，得到了No的明文，遂奔向作战指挥部

从上面的操作我们可以看到，由于n特别大，而e特别小，导致整个加解密的过程都可以不用到n，即使公钥的n很大，也有可能容易被破解。

出现上面这种情况的原因是明文进行指数e运算的时候，得到的值小于n，从而变得加密过程用不到n，所以解密过程也可以用不到n。所以需要把e变大，但又不能无限大，影响运算效率。比如现在大多数RSA算法的默认e大小65537(0x10001)，这样指数运算后的值肯定大于4096位的大小， 这样n在加解密过程中就必须参与到。

5. 快速幂

对于a的n次方 f(n) = $a^n$ ，a为底数，n为指数。随着n的递增，f(n)呈现指数增长。

快速幂的思想是：

• 底数a做平方，可使总体的幂运算次数降为原来的1/2。即 $a^n$ = (a*a)^n/2^。

• 利用公式 (a * b) % p = (a % p * b % p) % p，将幂运算后取余，转化为各个乘数取余、再乘积，最后再取余。

比如直接计算 2^100^ 的后三位。

一般思路： 2^100^ % 1000

/**
*
* @param base 底数
* @param power 指数
* @return 幂运算结果的后三位， 例如12345，则输出345。
*/
private long normalPower(long base, long power) {
    int result = 1;
    for (int i=1; i<= power; i++) {
        result *= base;
    }
    return result % 1000;
}

@Test
public void normalPowerTest() {
    long base = 2, power = 100;
    long result = normalPower(base, power);
    System.out.println(String.format("normalPower(%d, %d) = %d", base, power, result));
}

输出：

normalPower(2, 100) = 0

分析：幂运算结果很大，long型也无法存这么大的数据, 发生数据越界，输出0。

• 改进1: 利用公式 (a * b) % p = (a % p * b % p) % p，转化为各个乘数取余、再乘积，最后再取余。

private long quickPower(long base, long power) {
    int result = 1;
    for (int i=1; i<= power; i++) {
        result *= base;
        result %= 1000;
    }
    return result % 1000;
}

@Test
public void normalPowerTest() {
    long base = 2, power = 100;
    long result = normalPower(base, power);
    System.out.println(String.format("normalPower(%d, %d) = %d", base, power, result));

    long start = System.currentTimeMillis();
    result = quickPower(base, power);
    long time = System.currentTimeMillis() - start;
    System.out.println(String.format("quickPower(%d, %d) = %d, time: %d", base, power, result, time));

    power = 1000000000;
    start = System.currentTimeMillis();
    result = quickPower(base, power);
    time = System.currentTimeMillis() - start;
    System.out.println(String.format("quickPower(%d, %d) = %d, time: %d", base, power, result, time));
}

输出：

normalPower(2, 100) = 0
quickPower(2, 100) = 376, time: 0
quickPower(2, 1000000000) = 376, time: 4916

可见：转化为各个乘数取余、再乘积，最后再取余，可以得到正确结果。但如果指数很大时，运算效率较低，比较耗时。因为时间复杂度O(N), N为指数。

5.1 快速幂优化

快速幂的核心思想就是每一步都把指数分为两半，相应的底数做平方运算。这样能把很大的指数不断变小，所需循环次数也变小。

举例。

3^10 = 3 * 3 * 3 * 3 * 3 * 3 * 3 * 3 * 3 * 3;
3^10 = (3 * 3) * (3 * 3) * (3 * 3) * (3 * 3) * (3 * 3);
3^10 = 9^5;
2^10000 = 4^5000; //底数做一次平方，指数就从10000变成了5000， 减少了5000次循环操作。
9^5 = 9^(4+1) = 9^4 * (9^1); //指数为奇数时，先拆分成偶数+1的形式。偶数部分可继续执行缩半操作
= (9*9)^2 * (9^1)
= (81^2) * (9^1)
= 6561 * 9 //幂运算结果为当所有指数都为奇数时底数的乘积。及 6561^1 * 9 ^1 = 6561*9

改进2: 利用快速幂优化

private long quickPower2(long base, long power) {
    long result = 1;
    while (power > 0) {
        if (power % 2 == 0) {
            power /= 2; //指数缩小一半
            base = base * base % 1000; //底数做平方运算
        }else {
            power -= 1; //指数为奇数时先减1，转化为偶数
            result = result * base % 1000; //收集指数为1的结果
            power /= 2;
            base = base * base % 1000;
        }
    }
    return result;
}

或简写为：

private long quickPower3(long base, long power) {
    long result = 1;
    while (power > 0) {
        if (power % 2 == 1) {
            result = result * base % 1000; //收集指数为1的结果
        }
        power /= 2; //指数缩小一半, power为奇数和偶数都一样
        base = base * base % 1000; //底数做平方运算
    }
    return result;
}

输出：

normalPower(2, 100) = 0
quickPower(2, 100) = 376, time: 0
quickPower(2, 1000000000) = 376, time: 4911
quickPower2(2, 1000000000) = 376, time: 0

6. 扩展欧几里得算法求 ax + by = 1

6.1 欧几里得算法

欧几里得算法又称辗转相除法，非负整数a，b, 计算a，b的最大公约数gcd(a,b)时，有计算公式：

gcd(a,b) = gcd(b,a mod b)。

证明：

设d为a和b的公约数，则a、b可表示为: a = md, b = nd。
a mod b 可表示为 a - kb， a和kb都为d的倍数，kb = knd。
则 a - kb = (m - nk)d。
而 b = nd, 所以 d也为b 和 a mod b的公约数。最大公约数也一定相等。
所以有: gcd(a, b) = gcd(b, a mod b)。

6.1.1 欧几里得算法java算法实现

//b不等于0时，继续gcd运算
//b等于0时，此时a的值就是a和b的最大公约数
public static int gcd(int a, int b) {
    return b != 0 ? gcd(b, a%b) : a;
}

6.2 扩展欧几里得算法

扩展欧几里得算法：

对于正整数a和b，必存在有整数x和y，使得ax + by = gcd(a, b)。

我们知道，欧几里得算法总是把gcd(a,b)转化为求解gcd(b,a%b)，而当b变为0时返回a，此时的a就等于gcd。也就是说，欧几里得算法结束时变量a中存放的是gcd，变量b中存放的是0，因此此时显然有 a × 1 + b × 0 = gcd成立，此时有x=1、y=0成立。

不妨利用上面的欧几里得算法的过程来计算x和y。目前已知的是递归边界成立时为x=1、y=0，需要想办法反推出最初始的x和y。

分析：

ax + by = gcd(a, b); 设 r = a % b = a - (a / b) * b; // a / b 为a除以b的整数部分
因 gcd(a, b) = gcd(b, a % b)
则 gcd(a, b) = gcd(b, r), 而gcd(b, r)可表示为 bx' + ry';

故 bx' + ry' = bx' + (a - (a / b) b) y'
= bx' + ay' - (a / b)by'
= ay' + b(x' - (a / b) y')
则 ax + by = ay' + b(x' - (a / b) y')

故有：
x = y';
y = x' - (a / b) y';

结果表示 求解使得ax + by = gcd(a, b)的x和y可以通过辗转相除法转化求解x'和y'。

继续, 因为有：

gcd(b, r) = gcd(r, b % r)， 设r' = b % r, 则 gcd(b, r) = gcd(r, r') = x''r + y''r';
用r' = b % r = b - (b / r) r 带入：
x''r + y''r'
= x''r + y''(b - (b / r) r)
= by'' + r(x'' - (b / r)y'')
而 gcd(b, r) = bx' + ry'

故有：
x' = y'';
y' = x'' - (b / r)y'';

结果表示 求解x'和y'可以通过辗转相除法转化求解x''和y'' 。

一直递归。

中止条件：根据欧几里得，这个过程会有一个尽头：b = 0

第k次递归时，b变为了0, a$x^k$ + b$y^k$ = 𝑔𝑐𝑑(𝑎,𝑏),

即 a$x^k$ = 𝑔𝑐𝑑(𝑎,𝑏)

为使等式成立，可以令$x^k$=1, $y^k$=0(当然也可以为其他值)。此时gcd(a, b) = a。

这就找到了一组可行解，再一层层倒退回去，就能得到原始方程的一组整数解(x, y)。

6.2.1 扩展欧几里得算法Java实现

int x, y;
private int exgcd(int a, int b) {
    if (b == 0) {
        x = 1;
        y = 0;
        return a;
    }
    int d = exgcd(b, a % b); //递归计算exGcd（b,a%b）
    int temp = x;
    x = y;
    y = temp - (a / b)*y; //更新y = x(old) - a/b*y(old)
    return d; //d是gcd
}

最后得到的x,y 即为要求的一组整数解。

在得到这样一组解之后，就可以通过下面的式子得到全部解（其中K为任意整数）：

$x' = x + K * \frac{b}{gcd}$

$y' = y + K * \frac{a}{gcd}$

简单证明一下：

假设新的解为 x + s1、 y − s2，即有 a ∗ ( x + s1 ) + b ∗ ( y − s2 ) = gcd成立，

通过代入 ax + by = gcd 可以得到 as1 = bs2，

于是 s1/s2 = b/a 成立。

为了让s1和s2尽可能小，可以让分子和分母同时除以一个尽可能大的数，同时保证它们仍然是整数。显然，由于 b/gcd 与 a/gcd互质，因此gcd是允许作为除数的最大数，于是

s1/s2 = b/a = b/gcd / (a/gcd),

得s1和s2的最小取值就是 b/gcd 与 a/gcd。

证毕

也就是说，x和y的所有解分别以 b/gcd ​与 a/gcd​为周期。

6.2.2 x的最小非负整数解

那么x的最小非负整数解是什么呢？直观来说就是 x%(b/gcd)。但是由于通过exGcd函数计算出来的x、y可以为负数，因此实际上 x%(b/gcd)会得到一个负数，例如 (−15)。考虑到即便x是负数， x % (b/gcd)的范围也是在( − b/gcd, 0)，因此对于任意整数来说，

( x % (b/gcd) + b/gcd ) % (b/gcd) 

才是对应的最小非负整数解。

如果gcd=1，即ax+by=1时，全部解的公式简化为下式，且x的最小非负整数解也可以简化为

( x % b + b ) % b

故x为最小非负整数解为：

int exgcd = exgcd(a, b);
x = ( x % b + b ) % b; // x的最小非负整数解
y = (1 - 11*x)/20;

6.2.3 ax + by = c 是否一定存在整数解

二元一次方程：ax + by = c 是否一定存在整数解呢？ 不一定。

• 1）当a == 0 或 b == 0的时候，方程转为一元一次方程，可直接求得。

• 2）当c不是 gcd(a,b)的倍数的时候，方程无解。

因为a和b都是gcd(a, b)的倍数，所以ax+by(a, b的线性组合)也是gcd(a, b)的倍数。如果c不是gcd(a, b)的倍数，那么两边不相等，无解。这条结论叫做"裴蜀定理"。

• 3）当c是gcd(a, b)的倍数，设g = gcd(a, b)，则用ax + by = g求解出的一组解（x0，y0）,有ax + by = c 求出的一组解为(x0c/g, y0/g)

6.2.4 例子

求7d + 12k = 1的整数解。

再用代码测试一下：

@Test
public void exgcdTest() {
    int a = 7;
    int b = 12;
    //7x + 12y = 1
    exgcd = exgcd(a, b);
    x = ( x % b + b ) % b; // x的最小非负整数解
    y = (1 - 7*x)/12;
    System.out.println(String.format("exgcd = %d, 使%dx + %dy = 1的x的最小非负整数解为: x=%d, y=%d", exgcd, a, b, x, y));
}

int x, y;
private int exgcd(int a, int b) {
    if (b == 0) {
        x = 1;
        y = 0;
        return a;
    }

    int d = exgcd(b, a % b);
    int temp = x;
    x = y;
    y = temp - (a / b)*y;
    return d;
}

输出：

exgcd = 1, 使7x + 12y = 1的x的最小非负整数解为: x=7, y=-4