高速幂取模

快速幂取模

问题定义：

数论中经常出现的一个问题是对一个数的幂取模，也称为模取幂，即求a^b mod n。如果计算量较小，可以直接计算出a^b的值，再作模n运算。但是如果a和b的值都非常大，a^b的值用计算机难以表示，或者即使可以用大数运算的方式用计算机表示，也会因为耗时过长难以应用。基于模运算的基本性质，可以设计出一种算法，快速求解这一问题。这种方法为“快速幂取模”，也称为“反复平方法“。

算法描述：[1]

算法原理：

算法基础在于模运算的基本性质：  (a*b)%n = ( (a%n) * (b%n) ) %n

算法实现的原理是：

b = b(k) * 2^k + b(k-1)*2^(k-1) + ... + b(1) * 2^1 +b(0)*2^0   (其中< b(k),b(k-1)...b(1),b(0) > 是b 的二进制位)

a ^ b = a ^ (b(k)*2^k + b(k-1)*2^(k-1) + ... + b(1)*2 + b(0) )

         = [a ^ (b(k)*2^k)]  * [a^(b(k-1)*2^(k-1))] *...*[a^(b(1)*2) ]* a

         = p(k) * p(k-1) * ... * p(1) * p(0)

注意b(i)是二进制位，所以有 b(i) = 1 或 b(i) =0，即p(i) = a ^ (2^i) 或者 p(i) = 0

p(i) = a^(2^i) 即是通过对 a不断取平方得到， 这由上述算法描述中 步骤6完成，即 (a^(2^(i-1)))^2 = (a^(2^i))

若p(i) != 0，p(i)有一个上述的生成过程， 算法中d保持着现有的p(i)， 步骤9引入a，开始生成 p(j)，生成的过程是在每次的步骤6中完成由a->a^2->(a^2)^2 = a^(2^2) = ... = a^(2^i)的转化。

算法实现1：

#define MAX 0x40000000
int expMod ( int x, int y, int k )
{
    int mask = MAX;
    int modRes = 1;
    int times = 31;

    while (!(mask&y)) {
	mask = (mask>>1);
	--times;
    }
    x %= k;

    while (times--) {
	modRes = (modRes * modRes)%k;
	if (mask&y) {
	    modRes = (modRes * x)%k;
	}
	mask = (mask>>1);
    }
    return modRes;
}

算法实现2：

算法导论上在给出上述算法的同时提出了问题，即对b的二进制位从右向左扫描的算法如何实现，下面是我的一个实现，比实现1要简单

#define MASK 0xffffffff

int expMod2 ( int x, int y, int k )
{
    int tx = x;
    int modRes = 1;
    tx %= k;

    while (y&MASK) {
        if (y&1) 
            modRes = modRes * tx % k;
        y = (y>>1);
        tx = tx * tx % k;
    }
    return modRes;
}

算法应用：

素数测试[2]

RSA公钥加密算法[3]

参考资料：

[1] 算法导论31.6节

[2] 拉宾米勒素数测试：http://www.nocow.cn/index.php/%E6%8B%89%E5%AE%BE%E7%B1%B3%E5%8B%92%E7%B4%A0%E6%95%B0%E6%B5%8B%E8%AF%95

[3] RSA(algorithm) Wikipedia