大整數冪求模問題

   
      一、 問題描述：

              計算 (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 # 運行

 

/* * 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);}

#include <stdio.h>#include <time.h>#define MAX_SCALE 10long 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. 可以首選分治法和二進位分解法。分治法將所要求解的數值規模減半，而二進位分解則從數值的一個特別角度來尋求問題的解答。