算法与数据结构-数论之蒙哥马利模乘

        对于乘模运算 A*B%N，如果A、B都是1024位的大数，先计算A*B，再% N，就会产生2048位的中间结果，如果不采用动态内存分配技术就必须将大数定义中的数组空间增加一倍，这样会造成大量的浪费，因为在绝大多数情况下不会用到那额外的一倍空间，而采用动态内存分配技术会使大数存储失去连续性而使运算过程中的循环操作变得非常繁琐。所以模乘运算的首要原则就是要避免直接计算A*B。

设A=Sum[i=0 to k](A*r**i)，r=0x10000000，0<=A<r，则：C = A*B = Sum[i=0 to n](A*B*r**i) %N
可以用一个循环来处理：
C=0;
FOR i=n to 0
C=C*r
C=C+A*B
RETURN C

这样将一个多位乘法转换成了一系列单位乘法和加法，由于：
a*b %n = (a%n * b%n) %n
a+b %n = (a%n + b%n) %n

所以，对于求C=A*B %N，我们完全可以在求C的过程中完成：
C=0;
FOR i=n to 0
C=C*r %N
C=C+A*B %N
RETURN C

这样产生的最大中间结果是A*B 或C*r，都不超过1056位，空间代价会小得多，但是时间代价却加大了，因为求模的过程由一次变成了多次。对于孤立的乘模运算而言这种时间换空间的交易还是值得的，但是对于反复循环的乘模运算，这种代价就无法承受，必须另寻出路。

不难发现，模乘过程中复杂度最高的环节是求模运算，因为一次除法实际上包含了多次加法、减法和乘法，如果在算法中能够尽量减少除法甚至避免除法，则算法的效率会大大提高。
设A=Sum[i=0 to k](A*2**i)，0<=A<=1，则：
C= A*B = Sum[i=0 to k](A*B*2**i)

可用循环处理为：
C=0
FOR i FROM k TO 0
C=C*2
C=C+A*B
RETURN C

若令 C'= A*B*2**(-k)，则：
C'= Sum[i=0 to k](A*B*2**(i-k))

用循环处理即：
C'=0
FOR i FROM 0 TO k
C'=C'+A*B
C'=C'/2
RETURN C'

通过这一算法求A*B*2**(-k)是不精确的，因为在循环中每次除以2都可能有余数被舍弃了，但是可以通过这一算法求A*B*2**(-k) %N的精确值，方法是在对C'除2之前，让C'加上C'[0]*N。由于在RSA中N是两个素数的积，总是奇数，所以当C'是奇数时，C'[0]=1，C'+C'[0]*N 就是偶数，而当C'为偶数时C'[0]=0，C'+C'[0]*N还是偶数，这样C'/2 就不会有余数被舍弃。又因为C'+N %N = C' %N，所以在计算过程中加若干次N，并不会影响结果的正确性。

可以将算法整理如下：
C'=0
FOR i FROM 0 TO k
C'=C'+A*B
C'=C'+C'[0]*N
C'=C'/2
IF C'>=N C'=C'-N
RETURN C'

由于在RSA中A、B总是小于N，又0<=A,C'[0]<=1，所以：
C' = (C'+A*B+C'[0]*N)/2
C' < (C'+2N)/2
2C' < C'+2N
C' < 2N

既然C'总是小于2N，所以求C' %N 就可以很简单地在结束循环后用一次减法来完成，即在求A*B*2**(-k) %N的过程中不用反复求模，达到了我们避免做除法的目
的。当然，这一算法求得的是A*B*2**(-k) %N，而不是我们最初需要的A*B %N。但是利用A*B*2**(-k)我们同样可以求得A**E %N。

设R=2**k %N，R'=2**(-k) %N，E=Sum[i=0 to n](E*2**i)：
A'=A*R %N
X=A'
FOR i FROM n TO 0
X=X*X*R' %N
IF E=1 X=X*A'*R' %N
X=X*1*R' %N
RETURN X

最初：
X = A*R %N，

开始循环时：
X = X*X*R' %N
= A*R*A*R*R' %N
= A**2*R %N

反复循环之后：
X = A**E*R %N

最后：
X = X*1*R' %N
= A**E*R*R' %N
= A**E %N

如此，我们最终实现了不含除法的模幂算法，这就是著名的蒙哥马利算法，而X*Y*R' %N 则被称为“蒙哥马利模乘”。以上讨论的是蒙哥马利模乘最简单，最容易理解的二进制形式。蒙哥马利算法的核心思想在于将求A*B %N转化为不需要反复取模的A*B*R' %N，但是利用二进制算法求1024位的A*B*R' %N，需要循环1024次之多，我么必然希望找到更有效的计算A*B*R' %N的算法。

考虑将A表示为任意的r进制：
A = Sum[i=0 to k](A*r**i) 0<=A<=r

我们需要得到的蒙哥马利乘积为：
C'= A*B*R' %N R'=r**(-k)

则以下算法只能得到C'的近似值
C'=0
FOR i FROM 0 TO k
C'=C'+A*B
C'=C'/r
IF C'>=N C'=C'-N
RETURN C'

因为在循环中每次C'=C'/r 时，都可能有余数被舍弃。假如我们能够找到一个系数 q，使得(C' + A*B + q*N) %r =0，并将算法修改为：
C'=0
FOR i FROM 0 TO k
C'=C'+A*B+q*N
C'=C'/r
IF C'>=N C'=C'-N
RETURN C'

则C'的最终返回值就是A*B*R' %N的精确值，所以关键在于求q。由于：
(C' + A*B + q*N) %r =0
==> (C' %r + A*B %r + q*N %r) %r =0
==> (C'[0] + A*B[0] + q*N[0]) %r =0

若令N[0]*N[0]' %r =1，q=(C'[0]+A*B[0])*(r-N[0]') %r，则：
(C'[0] + A*B[0] + q*N[0]) %r
= (C'[0]+A*B[0] - (C'[0]+A*B[0])*N[0]'*N[0]) %r) %r
= 0

于是我们可以得出r为任何值的蒙哥马利算法：
m=r-N[0]'
C'=0
FOR i FROM 0 TO k
q=(C'[0]+A*B[0])*m %r
C'=(C'+A*B+q*N)/r
IF C'>=N C'=C'-N
RETURN C'

如果令 r=0x100000000，则 %r 和 /r 运算都会变得非常容易，在1024位的运算中，循环次数k 不大于32，整个运算过程中最大的中间变量C'=(C'+A*B+q*N)< 2*r*N < 1057位，算法效率就相当高了。唯一的额外负担是需要计算 N[0]'，使N[0]*N[0]' %r =1，而这一问题前面已经用欧几里德算法解决过了，而且在模幂运算转化成反复模乘运算时，N是固定值，所以N[0]'只需要计算一次，负担并不大。