蒙哥马利对自然是的快速求幂的原理在北邮人论坛上面找到了，通俗易懂。
"
快速幂运算的好处是减少了运算量，极大地提高了速度。例如A^65535（A的65535次幂），原始算法要做65535－1＝65534次乘法，而快速幂运算只需要做(16－1)×2＝30次乘法。
原理其实很好懂，假设要计算A^B，即底数是A，指数是B。把B写成二进制形式，拿4位来举例：B＝b4b3b2b1（二进制）。
先用B＝1111（二进制）来做解释。显然A^B ＝ A × A^2 × A^4 × A^8。
又显然，
A^2 ＝ A × A，
A^4 ＝ A^2 × A^2，
A^8 ＝ A^4 × A^4。
假如B的二进制位数更多，则依此类推。
上面这段看懂了吗？如果看懂的，就应该能够写出B＝1111（二进制）情况下，A^B的快速幂运算程序。
然后再来看B＝b4b3b2b1（二进制）的情况。这时候A^B ＝ A1 × A2 × A3 × A4。
其中，
A1 ＝ A（如果b1＝1）或者1（如果b1＝0），
A2 ＝ A^2（如果b2＝1）或者1（如果b2＝0），
A3 ＝ A^4（如果b3＝1）或者1（如果b3＝0），
A4 ＝ A^8（如果b4＝1）或者1（如果b4＝0）。
假如B的二进制位数更多，则依此类推。
这一段看懂了吗？如果看懂的，我不用继续解释，能不能把上面B＝1111（二进制）的程序改一改，加入对b4、b3、b2、b1等于0或者1的条件分支判断，让它变成A^B的快速幂运算程序？
"
根据上面的解释我实现了算法：
#include
/* a^b蒙哥马利快速幂运算算法 */
long long power(long long a, long long b)
{
if (b == 0)
return 1;
long long total = 1;
while (b != 0)
{
if ((b & 1) != 0)
{
total *= a;
}
a *= a;
b >>= 1;
}
return total;
}
int main()
{
long long n, t, k;
while (scanf("%lld%lld", &n, &t) != 2)
{
k = power(n, t);
printf("%lld\n", k);
}
return 0;
}
解释如下（没看懂，数学功底不够）：
“
如果直接计算C=A*B mod N，必须做很费时的mod N运算。从我前面讲的快速幂运算算法可知，乘法运算要进行多次，所以mod N运算也需要多次，比较慢。
假如我们做个转换，A-->A'，B-->B'，然后计算C'=A'*B' mod N，这时mod N运算可以快得多（当然这个转换是需要技巧的）。等到所有乘法运算都做完之后，最后结果再做一次转换，C'-->C，就转回到原本需要的结果了。
这个跟快速乘法运算有点类似：先用快速傅立叶变换，A-->A'，B-->B'，然后乘法变加法，C'=A'+B'，最后来个傅立叶逆变换，C'-->C，那么得到的C=A*B。加法当然比乘法快得多。
然后就是蒙哥马利算法之中，那个神奇的转换了：
A'=A*2^n mod N，
B'=B*2^n mod N，
显然，乘以2^n就是左移n位。至于mod N嘛……我相信n足够大，并且取某些特定值的时候，mod N可以变成不是除法的运算。例如，变成少量乘法、加减法。
由于N是奇数（通常还是大素数或者两个素数的乘积），所以2^n与N的最大公约数是1。2^n有乘法逆元（mod N域内），表示成2^-n（mod N），实际上逆元也是一个正整数，2^n乘这个整数会得1（mod N）。求逆元的算法不讨论了，欧几里得算法是最常规的一种。
显然，上述逆转换就是：
C=C'*2^n-1 mod N
我是现炒现卖的，原理大致是这些。具体的实现细节肯定还有很多问题，例如n值怎样选取，我就没有查到。貌似这里就衍生出不同的技巧性算法。
”
通过查论坛我得出一个结论：北邮确实比南邮好不少。
阅读(1059) | 评论(0) | 转发(0) |
