有些想补充的内容，但不好直接在初始的那一篇里改。因为那里讲得太细致了，是一步步讲得，要想再塞点别的东西进去就杂乱无章或喧宾夺主了。。所以重开一篇，后续有什么问题，都在这里更新。不是说细节我都明白透了（毕竟我也没实际写过RSA，只看过别人写的...看的还不仔细），只是自己懂一点也可以记录一点，顺便和大家分享。
背景就说到这，下面是我自己想补充的内容！
=================================
在第一篇中，从问题的起因，到寻求解决的步骤，是一步步沿着设计者的足迹，得到了算法的实现。每一步都寻求解决一个问题，很自然地到达最终的结果。这样顺水推舟地走，也许比较轻松，但可能走到最后，便忘记了开始，没办法一览顺流而下的整条河的概貌。
这里，我想仅就算法实现，再梳理一下步骤。
前一篇文章中，我粗略描述了蒙哥马利乘的接口（稍微改一下传参，BIGNUM还是不好作为返回值...得作为参数传进去再输出来）：
[1] int bn_mont_Mult (BIGNUM Result, BIGNUM A, BIGNUM B, BIGNUM X);
名称：蒙哥马利乘
功能：输出 Result = A * B [* R^(-1) mod X ]并且在前一篇文章里，实现 A^E mod X 大数模幂运算，写的是如下（半吊子）代码：
（在此仅改了传参规范）
int bn_PowMod_old (BIGNUM Result, BIGNUM A, BIGNUM E, BIGNUM X){ // Result = pow(A^E) mod XBIGNUM A0, Ret;int K, L;int *a = NULL, *e = NULL;K = bn_Expand(A, 2, a);bn_LShift(A0, A, K); // A0 = A*R, *R 实现为移位运算Ret = A0;L = bn_Expand(E, 2, e);for(int i = L; i >= 0; i--){ // 利用了简化幂运算的“平方-乘算法”,从高位开始bn_mont_Mult(Ret, Ret, Ret, X);if(e[i]==1){bn_mont_Mult(Ret, Ret, A0, X);}}bn_RShift(Ret, Ret, K); // Ret = Ret*R', 移回来return Ret;
}这段代码里，有几个问题：
- 我写了一个 bn_LShift 和 bn_RShift 用来表示“蒙哥马利变换”和“蒙哥马利逆变换”，实际上不会那么处理。。大数都是切割成一块块的，移位咋移啊。。
- 可以看到，bn_mont_Mult 里的 R 的长度（或值），实际上是依赖于入参 A 的长度的，也就是说并不是一个写死的值，所以实际上需要传入 A 的长度。
- 其实还有一个问题。。这个之前没有知友看出来。那就是 Ret 没有赋初值 = =。。。好严重的错误有木有....我的锅
那么我们试着改一下，首先不妨（暂时）假设我们相应地修改好了接口，并且方便起见，将bn_mont_Mult 相同两个大数相乘的情况封装为 bn_mont_Square 平方接口，即：
[2.a] int bn_mont_Mult (BIGNUM Result, BIGNUM A, BIGNUM B, uint32_t alen, BIGNUM X);
名称：蒙哥马利乘
功能：输出 Result = A * B [* R^(-1) mod X ], R = 2^(-alen)
备注：在此假设 alen 是 A 的比特长度[3.b] int bn_set_mont(BIGNUM R, uint32_t alen, BIGNUM X)
名称：获取蒙哥马利参数 R
功能：输出大数 R = 10000...(alen 个 0) mod X[2.c] int bn_mont_trans(BIGNUM Result, BIGNUM A, uint32_t alen, BIGNUM X)
名称：蒙哥马利变换
功能：输出 Result = A [* R mod X ][2.d] int bn_mont_inv(BIGNUM Result, BIGNUM A, uint32_t alen, BIGNUM X)
名称：蒙哥马利逆变换
功能：输出 Result = A [* R^(-1) mod X ][2.e] int bn_mont_Square (BIGNUM Result, BIGNUM A, uint32_t alen, BIGNUM X);
名称：蒙哥马利平方
功能：输出 Result = A * A [* R^(-1) mod X ], R = 2^(-alen)
实现：
int bn_mont_Square (BIGNUM Result, BIGNUM A, uint32_t alen, BIGNUM X) {return bn_mont_Mult(Result, A, A, alen, X);
}Ret 的初值应该是什么呢？感觉应该是1...但考虑到我们用的是蒙哥马利乘，也就意味着我们乘任何两个数，进了 bn_mont_Mult 接口都会先除以一个 R；至于普通乘法里乘1，意味着被乘数的值不变，那么作为输入蒙哥马利域的数，怎样的数才满足：任意一个数和它做 bn_mont_Mult 后值不变呢？显然，R是(唯一)一个这样的角色；想起来也很自然，因为蒙哥马利乘相当于我们把做乘法的位数“向右平移”了 R 嘛...那1可不是得“向左平移”个R，来把它抵消掉。。废话不多说了，模幂的代码相应可修改为：
(A0 改成 AR 吧...代表A*R)
int bn_PowMod (BIGNUM Result, BIGNUM A, BIGNUM E, BIGNUM X){ // Result = pow(A^E) mod XBIGNUM AR, Ret;int K, L;int *a = NULL, *e = NULL;K = bn_Expand(A, 2, a);bn_mont_trans(AR, A, K, X); // AR = A*R mod X, 变换至蒙哥马利域bn_set_mont(Ret, K, X); // Ret = R = 10000...(K个0) mod XL = bn_Expand(E, 2, e);for(int i = L; i >= 0; i--){ // 利用了简化幂运算的“平方-乘算法”,从高位开始bn_mont_Square(Ret, Ret, X);if(e[i]==1){bn_MultMod(Ret, Ret, AR, X);}}bn_mont_inv(Ret, Ret, K); // Ret = Ret*R^(-1) mod X, 由蒙哥马利域变回return Ret;
}这样好一些了，至少完整了。。但实际上还可以精简。事实上 bn_mont_trans 这个大数乘法可以也用蒙哥马利乘代替，好比说，我们实现了如下3.几的一组接口，而不是2.几的：
[3.a] int bn_mont_Mult (BIGNUM Result, BIGNUM A, BIGNUM B, uint32_t alen, BIGNUM X);
名称：蒙哥马利乘;双因数
功能：输出 Result = A * B [* R^(-1) mod X ], R = 2^(-alen)
备注：在此假设 alen 是 A 的比特长度[3.b] int bn_set_mont(BIGNUM R, uint32_t alen, BIGNUM X)
名称：获取蒙哥马利参数 R
功能：输出大数 R = 10000...(alen 个 0) mod X[3.c] int bn_mont_inv(BIGNUM Result, BIGNUM A, uint32_t alen, BIGNUM X)
名称：蒙哥马利逆变换;单因数,相当于 B = 1
功能：输出 Result = A [* R^(-1) mod X ]
实现：
int bn_mont_inv(BIGNUM Result, BIGNUM A, uint32_t alen, BIGNUM X) {BIGNUM one = 1; // 反正就是个赋值...领会就好return bn_mont_Mult(Result, A, one, alen, X);
}[3.d] int bn_mont_Square (BIGNUM Result, BIGNUM A, uint32_t alen, BIGNUM X);
名称：蒙哥马利平方
功能：输出 Result = A * A [* R^(-1) mod X ], R = 2^(-alen)
实现：
int bn_mont_Square (BIGNUM Result, BIGNUM A, uint32_t alen, BIGNUM X) {return bn_mont_Mult(Result, A, A, alen, X);
}正如前文所述，若调用蒙哥马利乘来实现使 A 赋值变到蒙哥马利域上，也就是使 AR = A*R mod X，则除 A 之外的另一个因数需为 R^2，即：
bn_mont_Mult (AR, A, R^2, alen, X)
才能使得：
AR = A*R^2 [*R^(-1) mod X]
于是代码修改为：
int bn_PowMod (BIGNUM Result, BIGNUM A, BIGNUM E, BIGNUM X){ // Result = pow(A^E) mod XBIGNUM AR, Ret, R2;int K, L;int *a = NULL, *e = NULL;K = bn_Expand(A, 2, a);bn_set_mont(R2, 2*K, X); // R2 = 10000...(2*K个0) mod Xbn_MultMod(AR, A, R2, alen, X); // A0 = A*R^2 [*R^(-1) mod X] = A*R mod Xbn_mont_inv(Ret, R2, K); // Ret = R = R2*R^(-1) mod XL = bn_Expand(E, 2, e);for(int i = L; i >= 0; i--){ // 利用了简化幂运算的“平方-乘算法”,从高位开始bn_mont_Square(Ret, Ret, X);if(e[i]==1){bn_MultMod(Ret, Ret, AR, X);}}bn_mont_inv(Ret, Ret, K); // Ret = Ret*R^(-1) mod X, 由蒙哥马利域变回return Ret;
}这样应该就比较接近真实的 RSA 蒙哥马利模幂的接口实现了。。至少看那些开源代码应该能明白哪一步大体在干啥了。
至于剩下的细节，比如说，如何保证模乘的因数比模数小，具体在哪一步取模，等这些细节问题日后有机会再谈吧（其实目前我也搞不清...主要没有需求让我写 RSA ...国际上那么多优秀的人优化了那么多年，怎么也轮不着自己实现（哈哈哈...）不过确实还是挺有意思的算法，不是吗？
)
