高阶幂的求余的方法

通常会有如下问法：

 

    有两个数，A和B，A的范围较小，B的范围较大。问A的B次幂的最后n位是多少？ n一般都小于5.

 

或

    有两个数，A和B，他们的范围都很大。问A的B次幂的最后n位是多少？

 

 

该题目主要思路是A的B次幂最后求余即可。但是A和B要不一个很大，要不都很大，是很难求到的。

 

都用到了数学里的同余数的概念。

 

同余数概念如下：

 

采用同余的习惯写法,对三个整数a,b,p, p>0
  a=b (mod p) (其实等式应该是三横,我没有找到HTML的表示方法)
表示a%p=b%p. 同余有个很好的基本性质: 支持加法和乘法. 即:
若 a=b (mod p), c=d (mod p), 则有
  a+c=b+d (mod p), a*c=b*d (mod p)
因此, 若a=b(mod p), 则对任意自然数n, 有
  a^n=b^n (mod p)
上式可以帮助解决楼主的问题(求a^n除以p的余数): 若a是个很大的数,我们先求出它除以p的余数b,然后计算b^n除以p的余数就行了.
但是,如果n也很大, 计算b^n的工作量还是不小. 受启于同余式中的Euler定理, 我们可有如下作法:
考虑b,b^2,b^3,...b^n,...这些数除以p的余数得到的数列S,S(有无穷项)中每一项都是0到p-1之间的某一个数.因此,一定有不相等的两个数m,n使得
  S[m]=S[n]
注意只要S[m]=S[n], 就一定有S[m+1]=S[n+1], S[m+2]=S[n+2], ... s[m+k]=S[n+k], ...
因此, 数列S实际上是循环的,并且循环的周期小于等于p, 因此用程序可很快算出其周期.

 

----------------------------------------------------------------------

 

 

取模的方法有一点a^(2^n)%c=(a^2%c)^(2^(n-1))%c

#include <iostream> using namespace std; long PowerMod(long a, long b, int k) { long tmp = a, ret = 1; while (b) { if (b & 1) ret = ret * tmp % k; tmp = tmp * tmp % k; b >>= 1; } return ret; } int main() { long a; long b; int i,N; cin>>N; for(i=0;i<N;i++) { cin>>a>>b; cout<<PowerMod(a,b,9907)<<endl; } return 0; }

 

 

 

求最后几位，如最后四位，只要模10000就行了。

要算一共有多少位，可以求log10(a^b)=blog10（a）（小数进位），例如5^10

有10*log10(5)=6.98,进位得7。

 

 

 

【解释】

 

zz: http://blog.csdn.net/shuqin1984/article/details/6543677

 

整数的二进制分解

 将大整数 power 按照二进制进行分解：

            power = b[N] * 2^N + b[N-1] * 2^(N-1) + … + b[1] * 2 + b[0]

其中： b[i] 取值为 0 或 1 ( i=0,1,.., N)，则

           a^power = a^(b[N] * 2^N) * a^(b[N-1] * 2^(N-1)) * … * a^(b[1] * 2) * a^b[0]
 很显然:  

          (1) 若 b[i] = 0, 则 b[i] * 2^i = 0 , a^(b[i]*2^i) = 1, 该乘积项可以消去；

          (2) 若 a[i] = 1, 则 b[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 是通过动态规划法逐渐求得的。

 

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