大整数幂求模问题

      一、 问题描述：
              计算 (a^power) % m , 其中power 是非负的大整数， a, m 为大于1 的整数。

      二、 问题分析： 
        很显然， 由于 power 是大整数，因此，必须考虑到幂计算的溢出问题。怎么避免溢出呢？ 可以通过降低幂次、逐次取模来实现。一个自然的想法是，将power 分成两个整数之和， power = n1 + n2，则 a^power = (a^n1) * (a^n2) .    通常采用二分法， n1 = n2 或 |n1-n2| = 1 。这就涉及到 (a * b) % m 的计算 。
   不难证明：   (a * b) % m = ((a % m) * (b % m)) % m。 
   证明如下：   设 a = pm + r1, b = qm + r2 , 则 r1 = a % m, r2 = b % m , 
                          则  (a * b) % m = (r1 * r2) % m = ((a % m) * (b % m)) % m.   证毕。

   特别地，当 a = b 时， (a^2) % m = ((a % m)^2) % m ;   这是一个简单却又关键性的结论。

三、 算法设计：

        算法1： 分治策略:
（1）若power 为奇数： 令 power = 2k+1, k 为非负整数  , 则
                   a^(2k+1) = (a^k)^2 *a ; a^(2k+1) % m = ((a^k % m)^2 % m * a) % m
  (2)   若power 为偶数:    令 power = 2k, k为非负整数， 则
                   a^(2k) = (a^k)^2 ; a^(2k) % m = (a^k % m)^2 % m
  (3)   若power == 1 : 返回 a % m ； 若 power == 0: 返回 1 % m。

       据以上(1)/(2)/(3) 条， 即可设计出相应的递归求解程序。时间复杂度为T(n) = T(n/2) + C = O(logn) ，空间复杂度为 O(1)， 不足之处在于有一定的递归调用开销。

        算法2： 整数的二进制分解
将大整数 power 按照二进制进行分解：
            power = a[N] * 2^N + a[N-1] * 2^(N-1) + … + a[1] * 2 + a[0]
其中： a[i] 取值为 0 或 1 ( i=0,1,.., N)，则
           a^power = a^(a[N] * 2^N) * a^(a[N-1] * 2^(N-1)) * … * a^(a[1] * 2) * a^a[0]
很显然: 
          (1) 若 a[i] = 0, 则 a[i] * 2^i = 0 , a^(a[i]*2^i) = 1, 该乘积项可以消去；
          (2) 若 a[i] = 1, 则 a[i] * 2^i = 2^i ， a^(2^i) % m = (a^(2^(i-1)) % m)^2 % m. 
                令 temp = a^(2^(i-1)) % m , 则 a^(2^i) % m = (temp * temp) % m。
比如,  power = 26 = 16 + 8 + 2 = (11010)2, 则 a^26 = a^(2^4 + 2^3 + 2^1); 
计算 a^power 实际上是计算 power 的二进制表示中所有位为1对应的乘积项。
                (a^26) % m = ((a^16 %m) * (a^8 %m) * (a^2 % m)) %m
而， a^8 % m = ((a^4 %m) × (a^4%m)) % m  是可以用动态规划法逐次求解的。 简单起见， 将 (a^i) % m 称为 “取模乘积项”。

算法描述：
令 temp = a^(2^i) % m ， i 是 power 的当前二进制位所对应的位置，
                                        temp 表示 power 的当前二进制位所对应的（取模）乘积项
STEP1:   初始化 temp 为 a % m  ,  result = 1;
STEP2:   对 power 进 行二进制分解。 若 power >=1 , 则进行模2运算：否则转至 STEP3
[1]  若余数为1, 则该位置上的二进制为1, 乘积中需要加入此时的 temp 项 :  result = (result * temp) % m；
           下一个二进制位对应的乘积项为 temp = (temp * temp) % m
[2]  若余数为0, 则该位置上的二进制为0,乘积中不需要加入 temp 项， result 保持不变，
           下一个二进制位对应的乘积项为 temp = (temp * temp) % m
STEP3： 所有的二进制位对应的乘积项都已计算，算法结束。

比如， result = (3^26) % 5 的计算过程如下：26 = （11010)2 ;
（1）初始化:  temp = 3^1 % 5 = 3;, result = 1 ; 
  (2)   检测 26 的最低位二进制为0,  则 不加入乘积项，result = 1,   temp =(3^2) % 5 = (temp * temp) % 5 = 4
  (3)   检测 26 的次低位二进制为1, 则 加入乘积项， result = (result * temp) % 3 = 4 , temp = (3^4) % 5 = (4*4) % 5 = 1;
  (4)   检测 26 的次低位二进制为0, 则 不加入乘积项， result = 4, temp = (3^8) % 5 = (1*1) % 5 = 1;
  (5)   检测 26 的次低位二进制为1, 则 加入乘积项， result = (result * temp) % 5 = 4, temp = (3^16) % 5 = 1;
  (6)   检测 26 的次低位二进制为1, 则 加入乘积项， result = (result * temp) % 5 = 4, temp = (3^32) % 5 = 1.
   由于所有位都检测完毕， 故 3^26 % 5 = 4.  由上可知，
   3^26 % 5 = ((3^16) % 5)) * ((3^8) % 5) * ((3^2) % 5) % 5.  其中 3^16 % 5, 3^8 % 5, 3^2 % 5 是通过动态规划法逐渐求得的。

完整的程序如下： 
程序在 Ubuntu10.10 gcc4.4.5 环境下编译运行通过。
  $ gcc -g -Wall bintmode.c runtime.c -o bintmode # 编译连接
  $ bintmode or gdb bintmode # 运行

[cpp] view plaincopyprint?
/*
* bintmode.c :  此程序计算 (a^power) % mod. power 是大整数
* 基本公式： (a*b) mod m = ((a mod m) * (b mod m)) mod m
*/ 
 
#include <stdio.h> 
#include <assert.h> 
 
int modMRec(int, int, int); 
int modMDyn(int, int ,int); 
void testModMRec(int); 
void testModMDyn(int); 
 
void testValid(int (*fun)(int a, int power, int mode)); 
void testInvalid(int (*fun)(int a, int power, int mode)); 
 
int main() 
{ 
   printf(" ******** 使用递归技术求解： ******** \n"); 
   testValid(modMRec); 
   /* testInvalid(modMRec);  */  
   defRuntime(testModMRec); 
 
   printf(" ******** 使用二进制分解技术求解： ******** \n"); 
   testValid(modMDyn); 
   /* testInvalid(modMDyn);  */  
   defRuntime(testModMDyn); 
 
   return 0; 
} 
 
/*
* modMRec: 递归求解 (a^power) % mod , power 是个大整数
*/ 
int modMRec(int a, int power, int mod) 
{ 
   assert(a>=1 && power >=0 && mod >=1); 
   if (power == 0) { 
       return 1 % mod; 
   } 
   if (power == 1) { 
       return a % mod; 
   } 
   if (power % 2 != 0) { 
       int temp = modMRec(a, power/2, mod); 
       return  (temp * temp * a) % mod; 
   }  
   else { 
       int temp = modMRec(a, power/2, mod); 
       return  (temp * temp) % mod; 
   } 
} 
 
/*
* modMDyn: 使用二进制分解技术求解 (a^power) % mod , power 是个大整数
*/ 
int modMDyn(int a, int power, int mod) 
{ 
   assert(a>=1 && power >=0 && mod >=1); 
   int result = 1; 
   int temp = a % mod; 
   while (power >= 1) { 
      if (power % 2 != 0) { 
         result = (result * temp) % mod;    
      } 
      temp = (temp * temp) % mod; 
      power /= 2;    
   } 
   return result; 
} 
 
void testValid(int (*fun)(int a, int power, int mode)) 
{ 
   int base[5] = {2,3,5,7,11}; 
   int i,k; 
 
   for (i=0; i < 5; i++) { 
     for (k = 0; k < 20; k++) { 
        printf("%d ^ %d mod %d = %d .\n", base[i], k, 10, (*fun)(base[i], k, 10)); 
     }     
   } 
} 
 
void testInvalid(int (*fun)(int a, int power, int mode)) 
{  
   printf("0 ^ 3 mod 5"); 
   (*fun)(0, 3, 5);   
   printf("3 ^ -1 mod 5"); 
   (*fun)(3, -1, 5); 
   printf("3 ^ 5 mod 0"); 
   (*fun)(3, 5, 0); 
} 
 
void testModMRec(int n) 
{ 
   modMRec(2, n, 10); 
}   
 
void testModMDyn(int n) 
{ 
   modMDyn(2, n , 10); 
} 

[cpp] view plaincopyprint?
#include <stdio.h> 
#include <time.h> 
 
#define MAX_SCALE 10 
 
long defScale[MAX_SCALE] = {1, 10, 100, 1000, 10000, 100000, 1000000, 10000000, 100000000, 1000000000}; 
 
/* runtime: 测量 fun 指向的函数的运行时间
* fun: 指向测试函数的指针 ；
* scale:  问题规模数组 
*/ 
void runtime(void (*fun)(long scale), long* scale, int len); 
 
void defRuntime(void (*fun)(long scale)); 
 
void runtime(void (*fun)(long scale), long *scale, int len) 
{ 
   int i; 
   clock_t start, end; 
   for (i = 0; i < len; i++) { 
      start = clock(); 
      (*fun)(scale[i]); 
      end = clock(); 
      printf("scale: %12ld\t run time: %8.4f (ms). \n", scale[i], ((double)(end-start)) * 1000 / CLOCKS_PER_SEC); 
   } 
} 
 
void defRuntime(void (*fun)(long scale)) 
{ 
   runtime(fun, defScale, MAX_SCALE); 
} 


四、总结： 关于科学计算的算法
1. 往往要先查阅相关资料，了解相应的数学性质及结论，并对问题进行化简；例如本例中使用 (a*b)%m = ((a%m)*(b%m))%m 的模性质，避免了乘法溢出的问题。
2. 可以首选分治法和二进制分解法。分治法将所要求解的数值规模减半，而二进制分解则从数值的一个特别角度来寻求问题的解答。

评论

2 楼 刘彦明 2013-03-29  

其实这种做法我还不太熟悉、你一下就看懂了。能给我讲下不？

1 楼 三断笛 2013-03-28  

好东西啊，检测素数会要幂模，有了博主的这些就爽了