RSA的大数运算
余年寄山水
659次浏览
2021年01月19日 21:34
最佳经验
本文由作者推荐
dnf黄金粉末-说尺子
俺曾经查阅了网上找得到的各种用于实现
RSA
的大数运算库，然而最终还是决定自 己
动手写一个。因为凡是效率高速度快的代码（
crypto++
、
mira cl
、
freelip
、
rsaref
等）
，要么使
用的数据结构过于复杂，
要么编码风格杂乱无章，
俺的水平和耐心都实在是有限，
以 至于无
法读懂这些东西。
而俺读得懂的一些代码，
其实现方式却又过于幼稚，
效率极低速度一塌糊
涂。
俺觉得像俺这样的人不在少数，
于是决心写一个清晰易懂，< br>效率也过得去的东西奉献给
大家。
这个函数库刚做好的时候，生成
1024
位的随机密钥耗时大约
5
分钟，俺认为是可以接
受的。但后来找到一个叫
tE!
的老外用
m iracl
库写的
RsaTools
，发现其生成
1024
位的密钥
耗时不超过三秒钟！
于是俺针对俺的代码开始了艰苦的优化工作，
希望能达到甚至超过 这一
水平
?
一周之后
1024
位密钥的平均生成时间已经降至
5
秒左右，但是单单依靠优化代码来
进一步提高速度也非常困难了。于是俺开始借助金山词 霸来查阅能够通过
找到的一
切与
RSA
算法相关的论 文，但是网上关于
RSA
算法的论述绝大多数都是用于硬件实现的，
将其算法流程用软 件设计语言来实现极其繁琐
?
而且俺发现这样做下去俺只会离自己的初衷
越来越远:
俺的代码将不再清晰易懂
?
所以俺一度准备放弃
?
准备放弃之后，
心态平静了许多，
再回头去看那些原来不太能够理解的
RSA
算法原理，
却发现其实也不是那么高深莫测，不急不躁地慢慢看，慢慢想，突然就一下子全明白 了。一
番改进之后，
现在这个版本的函数库同样具有非常简单而清晰的结构，
速度也不 算慢，
生成
1024
位的密钥在俺
PIII
900
的笔记本上平均耗时不超过两秒。程序使用
C++
编写，可在
VC6.0
下直接编译通过，希望大家喜欢。如果发现
Bug
或者有好的修改建议，俺将非常感
谢您能够给俺一个
。
最后，感谢看雪论坛，感谢雪兄多次热心相助，俺在此学到了很多知识，当然
还要乘机拍拍马屁，感谢俺家甜甜的支持！
afanty@
RSA
原理
:
选取两个不同的大素数
p
、
q
，并计算
N=p*q
选取小素数
d
，并计算
e
，使
d*e % (p-1)(q-1)=1
对于任意
A
若
B=A**d % N
则
A=B**e % N
可 见
d
、
e
形成了非对称秘钥关系，加密者用公钥
d
加密，解 密者可用私钥
e
解密，第三者
即使拦截了密文
B
、公钥
d< br>和
N
，在不知道
p
、
q
的前提下，无法推算出
e
，从而无法获得
明文
A
。当
N
取非常大的值时，将其因 式分解成
p
、
q
是非常困难的，例如当
N
为
102 4 bit
时，据分析，需动用价值数千万美金的大型计算机系统并耗费一年的时间
?
RSA
密钥的选取和加解密过程都非常简洁，在算法上主要要实现四个问题：
1 ?
如何处理大数运算
2
、如何求解同余方程
XY % M = 1
3 ?
如何快速进行模幂运算
4 ?
如何获取大素数
实际上，
在实现
RSA
算 法的过程中大家会发现后三个问题不是各自独立的，
它们互有关联，
环环相套，相信届时你会意 识到：
RSA
算法是一种“优美”的算法！
RSA
依赖大数运算，目前主流
RSA
算法都建立在
1024
位的大数运算之上。而大多数
的编译器只能支持到
64
位的整数运算，即我们在运算中所使用的整数必须小于等于
64
位，
即：0xffffffffffffffff
，也就是
18446744
，这远远达不 到
RSA
的需要，于是需要
专门建立大数运算库来解决这一问题。
最简单的办法是将大数当作数组进行处理，数组 的各元素也就是大数每一位上的数字，
通常采用最容易理解的十进制数字
0
—
9
。然后对“数字数组”编写加减乘除函数。但是这
样做效率很低，
因为二进制为1024
位的大数在十进制下也有三百多位，
对于任何一种运算，
都需要在两个有 数百个元素的数组空间上多次重循环，
还需要许多额外的空间存放计算的进
退位标志及中间结果 。
另外，
对于某些特殊的运算而言，
采用二进制会使计算过程大大简化，
而这 种大数表示方法转化成二进制显然非常麻烦，
所以在某些实例中则干脆采用了二进制数
组的方法 来记录大数，当然这样效率就更低了
?
一个有效的改进方法是将大数表示为一个
n
进制数组，
对于目前的
32
位系统而言
n
可
以取值为
2
的
32
次方，即
0x1 00000000
，假如将一个二进制为
1024
位的大数转化成
0x100 00000
进制，
就变成了
32
位，
而每一位的取值范围不再是二进 制的
0
—
1
或十进制的
0
—
9
，
而是
0-0xffffffff
，
我们正好可以用一个
32
位的DWORD
（如：
无符号长整数，
unsigned
long
）
类型来表示该值。所以
1024
位的大数就变 成一个含有
32
个元素的
DWORD
数组，
而针对
DWORD
数组进行各种运算所 需的循环规模至多
32
次而已。而且
0x100000000
进
制与二进制，对于计算机来说，几乎是一回事，转换非常容易
?
例如大数
18446744
，等于
0xffffffff
ffffffff
，就相当于十进制的
99
：有两
位，每位都是
0 xffffffff
。而
18446744
等于
0x0000 00000000
，就
相当于十进制的
100
：有三位，第一位是
1
，其它两位都是
0
，如此等等。在实际应用中，
“数字数组”的排列顺序采 用低位在前高位在后的方式，这样，大数
A
就可以方便地用数
学表达式来表示其值：
A=Sum[i=0 to n] (A[i]*r**i)
，
r=0x100000000
，
0<=A[i]< r
其中
Sum
表示求和，
A[i]
表示用以记录
A
的数组的第
i
个元素，
**
表示乘方。
任何整数运算最终都能分解成数字与数字之间的运算，在
0x100000000
进 制下其“数
字”最大达到
0xffffffff
，其数字与数字之间的运算，结果也必 然超出了目前
32
位系统的字
长。在
VC++
中，存在一个
__int64
类型可以处理
64
位的整数，所以不用担心这一问题，而
在 其它编译系统中如果不存在
64
位整形，就需要采用更小的进制方式来存储大数，例如
16
位的
WORD
类型可以用来表示
0x10000
进制。
但效率更高的办法还是采用
32
位的
DWORD
类型，只不过将
0x100000000
进制改成
0x 40000000
进制，这样两个数字进行四则运算的最
大结果为
0x3fffffff * 0x3fffffff
，小于
0xffffffffff ffff
，可以用一个双精度浮点类型（
double
，
52
位有效 数字）来储存这一中间结果，只是不能简单地用高位低位来将中间结果拆分成两
个“数字”
。< br>
设有大数
A
、
B
和
C
，其中< br>A>=B
：
A=Sum[i=0 to p](A[i]*r**i)
B=Sum[i=0 to q](B[i]*r**i)
C=Sum[i=0 to n](C[i]*r**i)
r=0x100000000
，
0<=A[i],B[ i],C[i]
p>=q
则当
C=A+B
、C=A-B
、
C=A*B
时，我们都可以通过计算出
C[i]
来 获得
C
：
C=A+B
，
显然
C[i]
不总是等于
A[i]+B[i]
，
因为
A[i]+B[i]
可能
>0xffffffff
，
而
C[i]
必须
<=0xf fffffff
，
这时就需要进位，当然在计算
C[i-1]
时也可能产生了 进位，所以计算
C[i]
时还要加上上次的
进位值。如果用一个
64
位变量
result
来记录和，另一个
32
位变量
carry
来记录进位，则有：
carry = 0
FOR i FROM 0 TO p
result=A[i]+B[i]+carry
C[i]=result%0x100000000
carry=result/0x100000000
IF carry=0 n=p
Else: n = p + 1
C=A-B
，同 理
C[i]
不总是等于
A[i]-B[i]
，因为
A[i]-B[i ]
可能
<0
，而
C[i]
必须
>=0
，这时就需< br>要借位，同样在计算
C[i-1]
时也可能产生了借位，所以计算
C[i]时还要减去上次的借位值
:
carry = 0
FOR i FROM 0 TO p
IF A[i]-B[i]-carry>=0
C[i]=A[i]-B[i]-carry
carry = 0
Else
C[i]=0x100000000+A[i]-B[i]-carry
carry = 1
n = p
WHILE C[n]=0 n=n-1
C=A*B
，首先我们需要观察日常做乘法所用的“竖式计算”过程：
A3
A2
A1
A0
*
B2
B1
B0
-------------------------------
=
A3B0 A2B0 A1B0 A0B0
+
A3B1 A2B1 A1B1 A0B1
+ A3B2 A2B2 A1B2 A0B2
-------------------------------
=
C5
C4
C3
C2
C1
C0
可以归纳出：
C[i]=Sum[j=0 to q](A[i-j]*B[j])
，其中
i-j
必须
>=0
且
<=p
。当然这一结论 没有考
虑进位，
虽然计算
A[i-j]*B[j]
和
Sum
的时候都可能发生进位，
但显然这两种原因产生的进位
可以累加成一个进位值
?
最终可用如下算法完成乘法
:
n = p + q - 1
carry = 0
For i FROM 0 TO n
Result = carry
For j FROM 0 TO q
IF 0<=i-j<=p result=result+A[i-j]*B[j]
C[i]=result%0x100000000
carry=result/0x100000000
IF carry!=0
n = n + 1
c [n] = carry
对于
C=A/B
，由于无法将
B
对
A
“试商”< br>，我们只能转换成
B[q]
对
A[p]
的试商来得到一个
近似 值，所以无法直接通过计算
C[i]
来获得
C
，只能一步步逼近
C< br>。由于：
B*(A[p]/B[q]-1)*0x100000000**(p-q)
令：
X=0
，重复
A=A-X*B
，
X=X+(A[p]/B[q] -1)*0x100000000**(p-q)
，直到
A则：
X=A/B
，且此时的
A=A%B
注意对于任意大数
A*0x100000000**k
，都等价于将
A
的数组中的各元素左移
k
位，不必
计算；同样，
A/0x1000 00000**k
则等价于右移。
在
RSA
算法中，往往要在已知
A
、
N
的情况下，求
B
，使得
(A*B)%N=1
。即相当于求解
B
、
M
都是未知数的二元一次不定方程
A*B-N*M=1
的最小整数解。
而针对不定方程
ax-by=c
的最小整数解，
古今中外都进行过详尽的研 究，
西方有著名的
欧几里德算法，即辗转相除法，中国有秦九韶的“大衍求一术”
。事 实上二元一次不定方程
有整数解的充要条件是
c
为
a
、
b< br>的最大公约数。即当
c=1
时，
a
、
b
必须互质。而 在
RSA
算
法里由于公钥
d
为素数，素数与任何正整数互质，所以可 以通过欧几里德方程来求解私钥
e?
欧几里德算法是一种递归算法，比较容易理解：
例如：
11x-49y=1
，求
x
（
a
）
11 x - 49 y = 1
49%11=5 ->
（
b
）
11 x -
5 y = 1
11%5 =1 ->
（
c
）
x -
5 y = 1
令
y=0
代入（
c
）得
x=1
令
x=1
代入（
b
）得
y=2
令
y=2
代入（
a
）得
x=9
同理可使用递归算法求得任意
ax-by=1
（
a
、
b
互质）
的解。
实 际上通过分析归纳将递归
算法转换成非递归算法就变成了大衍求一术
:
x=0,y=1
While a! = 0
i = y
y = X - y * a / b
X = i
i = a
a=b%a
b = i
IF x<0 x=x+b
RETURN x
模幂运算是
RSA
的核心算法，
最直接地决定了
RSA
算法的性能。
针对快速模幂运算
这一课题，
西方现代数学家提出了大量的解决方案，< br>通常都是先将幂模运算转化为乘模运算
?
例如求
D=C**15 % N
，由于：
a*b % n = (a % n)*(b % n) % n
，所以：
C1 =C*C
% N
=C**2
% N
C2 =C1*C
% N
=C**3
% N
C3 =C2*C2 % N
=C**6
% N
C4 =C3*C
% N
=C**7
% N
C5 =C4*C4 % N
=C**14 % N
C6 =C5*C
% N
=C**15 % N
即：
对于
E=15
的幂模运算可分解为
6
个乘 模运算，
归纳分析以上方法可以发现对于任
意
E
，都可采用以下算法计算D=C**E % N
：
d = 1
While e >= 0
IF E%2=0
C=C*C % N
e = e / 2
Else
D=D*C % N
e = e - 1
RETURN D
继续分析会发现，要知道
E
何时能整除
2
，并不需要反复进行减一或除二的操作，只
需验证
E
的二进制各位是
0
还是
1
就可以了，
从左至右或从右至 左验证都可以，
从左至右
会更简洁，设
E=Sum[i=0 to n](E[i]*2**i)
，
0<=E[i]<=1
，则：