1 /*计算n^p*/
2 unsigned power(unsigned n,unsigned p)
3 {
4     for(int i=0;i<p;i++) n*=n;
5     return n;
6 }

unsigned Power(unsigned n,unsigned p)
{unsigned main=n; //用main保存结果unsigned odd=1; //odd用来计算“剩下的”乘积while (p>1) {//一直计算，直到指数小于或等于1if((p%2)!=0){// 如果指数p是奇数，则说明计算后会剩一个多余的数，那么在这里把它乘到结果中odd*=main; //把“剩下的”乘起来
        }main*=main; //主体乘方p/=2; //指数除以2
   }return main*odd; //最后把主体和“剩下的”乘起来作为结果
}

//完美版：
unsigned Power(unsigned n, unsigned p)
{ // 计算n的p次方unsigned odd = 1; //oddk用来计算“剩下的”乘积while (p > 1){ // 一直计算到指数小于或等于1if (( p & 1 )!=0){ // 判断p是否奇数，偶数的最低位必为0odd *= n; // 若odd为奇数，则把“剩下的”乘起来
      }n *= n; // 主体乘方p /= 2; // 指数除以2
     }return n * odd; // 最后把主体和“剩下的”乘起来作为结果
}

 1 unsigned Montgomery(unsigned n, unsigned p, unsigned m)
 2 { // 快速计算 (n ^ e) % m 的值，与power算法极类似
 3     unsigned r = n % m; // 这里的r可不能省
 4     unsigned k = 1;
 5     while (p > 1)
 6     {
 7         if ((p & 1)!=0)
 8         {
 9             k = (k * r) % m; // 直接取模
10         }
11         r = (r * r) % m; // 同上
12         p /= 2;
13     }
14     return (r * k) % m; // 还是同上
15 }

unsigned Montgomery(unsigned n,unsigned p,unsigned m)
{ //快速计算(n^e)%m的值unsignedk=1;n%=m;while(p!=1){if(0!=(p&1))k=(k*n)%m;n=(n*n)%m;p>>=1;}return(n*k)%m;
}

 1 怎么判断一个数是否为素数？
 2
 3 笨蛋的作法：
 4 bool IsPrime(unsigned n)
 5 {
 6     if (n<2)
 7     { //小于2的数即不是合数也不是素数
 8     throw 0;
 9     }
10     for (unsigned i=2;i<n;++i)
11     { //和比它小的所有的数相除，如果都除不尽，证明素数
12         if (n%i==0)
13         {//除尽了，则是合数
14             return false;
15         }
16     }
17     return true;
18 }
19
20 一个数去除以比它的一半还要大的数，一定除不尽，所以还用判断吗？？

 1 下面是小学生的做法：
 2 bool IsPrime(unsigned n)
 3 {
 4     if (n<2)
 5     {//小于2的数即不是合数也不是素数
 6         throw 0;
 7     }
 8     for(unsigned i=2;i<n/2+1;++i)
 9     { // 和比它的一半小数相除，如果都除不尽，证明素数
10         if ( 0 == n % i )
11         { // 除尽了，合数
12             return false;
13         }
14     }
15     return true; // 都没除尽，素数
16 }
17
18 一个合数必然可以由两个或多个质数相乘而得到。那么如果一个数不能被比它的一半小的所有的质数整除，那么比它一半小的所有的合数也一样不可能整除它。建立一个素数表是很有用的。

下面是中学生的做法：
bool IsPrime2(unsigned n)
{if ( n < 2 ){ // 小于2的数即不是合数也不是素数throw 0;}static unsigned aPrimeList[] = { // 素数表1, 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41,43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 113, 193, 241, 257, 337, 353, 401, 433, 449, 577, 593, 641, 673, 769, 881, 929, 977, 1009, 1153, 1201, 1217, 1249, 1297,1361, 1409, 1489, 1553, 1601, 1697, 1777, 1873, 1889, 2017, 2081, 2113, 2129, 2161, 2273, 2417, 2593, 2609, 2657, 2689, 2753, 2801, 2833, 2897, 3041, 3089, 3121, 3137, 3169, 3217, 3313, 3329, 3361, 3457, 3617, 3697, 3761, 3793, 3889, 4001, 4049, 4129, 4177, 4241, 4273, 4289, 4337, 4481, 4513, 4561, 4657, 4673, 4721, 4801, 4817, 4993, 5009, 5153, 5233, 5281, 5297, 5393, 5441, 5521, 5569, 5857, 5953, 6113, 6257, 6337, 6353, 6449, 6481, 6529, 6577, 6673, 6689, 6737, 6833, 6961, 6977, 7057, 7121, 7297, 7393, 7457, 7489, 7537, 7649, 7681, 7793, 7841, 7873, 7937, 8017, 8081, 8161, 8209, 8273, 8353, 8369, 8513, 8609, 8641, 8689, 8737, 8753, 8849, 8929, 9041, 9137, 9281, 9377, 9473, 9521, 9601, 9649, 9697, 9857 };const int nListNum = sizeof(aPrimeList)/sizeof(unsigned);//计算素数表里元素的个数for (unsigned i=2;i<nListNum;++i ){ if(n/2+1<aPrimeList[i]){return true;}if(0==n%aPrimeList[i]){return false;}}/*由于素数表中元素个数是有限的，那么对于用素数表判断不到的数，就只有用笨蛋办法了*/for (unsigned i=aPrimeList[nListNum-1];i<n/2+1;i++ ){ if (0==n%i){ // 除尽了，合数 return false;}}return true;
} 

 还是太糟了，我们现在要做的对于大型素数的判断，那个素数表倒顶个P用！当然，我们可以利用动态的素数表来进行优化，这就是大学生的做法了。但是动态生成素数表的策略又复杂又没有效率，所以我们还是直接跳跃到专家的做法吧：根据上面讲到的费马小定理，对于两个互质的素数N和P，必有：N^(P-1)%P=1 那么我们通过这个性质来判断素数吧，当然，你会担心当P很大的时候乘方会很麻烦。不用担心！我们上面不是有个快速的幂模算法么？好好的利用蒙格马利这位大数学家为我们带来的快乐吧！算法思路是这样的： 对于N，从素数表中取出任意的素数对其进行费马测试，如果取了很多个素数，N仍未测试失败，那么则认为N是素数。当然，测试次数越多越准确，但一般来讲50次就足够了。另外，预先用“小学生”的算法构造一个包括500个素数的数组，先对Q进行整除测试，将会大大提高通过率，方法如下:
bool IsPrime3(unsigned n)
{if ( n < 2 ){ // 小于2的数即不是合数也不是素数throw 0;}static unsigned aPrimeList[] = {2, 3, 5, 7, 11, 17, 19, 23, 29, 31, 41,43, 47, 53, 59, 67, 71, 73, 79, 83, 89, 97};const int nListNum = sizeof(aPrimeList) / sizeof(unsigned);for (int i=0;i<nListNum;++i){ // 按照素数表中的数对当前素数进行判断if (1!=Montgomery(aPrimeList[i],n-1,n)) // 蒙格马利算法
        {return false;}}return true;
}OK，这就专家的作法了。 等等，什么？好像有点怪，看一下这个数29341，它等于13 * 37 * 61，显然是一个合数，但是竟通过了测试！！哦，抱歉，我忘了在素数表中加入13，37，61这三个数，我其实是故意的，我只是想说明并费马测试并不完全可靠。现在我们发现了重要的一点，费马定理是素数的必要条件而非充分条件。这种不是素数，但又能通过费马测试的数字还有不少，数学上把它们称为卡尔麦克数，现在数学家们已经找到所有10 ^ 16以内的卡尔麦克数，最大的一个是9585921133193329。我们必须寻找更为有效的测试方法。数学家们通过对费马小定理的研究，并加以扩展，总结出了多种快速有效的素数测试方法，目前最快的算法是拉宾米勒测试算法，下面介绍拉宾米勒测试。

================================================================
拉宾米勒测试拉宾米勒测试是一个不确定的算法，只能从概率意义上判定一个数可能是素数，但并不能确保。算法流程如下:1.选择T个随机数A，并且有A<N成立。2.找到R和M，使得N=2*R*M+1成立。快速得到R和M的方式：N用二进制数B来表示，令C=B-1。因为N为奇数（素数都是奇数），所以C的最低位为0，从C的最低位的0开始向高位统计，一直到遇到第一个1。这时0的个数即为R，M为B右移R位的值。3.如果A^M%N=1，则通过A对于N的测试，然后进行下一个A的测试4.如果A^M%N!=1，那么令i由0迭代至R，进行下面的测试5.如果A^((2^i)*M)%N=N-1则通过A对于N的测试，否则进行下一个i的测试 6.如果i=r，且尚未通过测试，则此A对于N的测试失败，说明N为合数。7.进行下一个A对N的测试，直到测试完指定个数的A通过验证得知，当T为素数，并且A是平均分布的随机数，那么测试有效率为1 / ( 4 ^ T )。如果T > 8那么测试失误的机率就会小于10^(-5)，这对于一般的应用是足够了。如果需要求的素数极大，或着要求更高的保障度，可以适当调高T的值。下面是代码：bool RabbinMillerTest( unsigned n )
{if (n<2){ // 小于2的数即不是合数也不是素数throw 0;}const unsigned nPrimeListSize=sizeof(g_aPrimeList)/sizeof(unsigned);//求素数表元素个数for(int i=0;i<nPrimeListSize;++i){// 按照素数表中的数对当前素数进行判断if (n/2+1<=g_aPrimeList[i]){// 如果已经小于当前素数表的数，则一定是素数return true;}if (0==n%g_aPrimeList[i]){// 余数为0则说明一定不是素数return false;}}// 找到r和m，使得n = 2^r * m + 1;int r = 0, m = n - 1; // ( n - 1 ) 一定是合数while ( 0 == ( m & 1 ) ){m >>= 1; // 右移一位r++; // 统计右移的次数
    }const unsigned nTestCnt = 8; // 表示进行测试的次数for ( unsigned i = 0; i < nTestCnt; ++i ){ // 利用随机数进行测试，int a = g_aPrimeList[ rand() % nPrimeListSize ];if ( 1 != Montgomery( a, m, n ) ){int j = 0;int e = 1;for ( ; j < r; ++j ){if ( n - 1 == Montgomery( a, m * e, n ) ) {break;}e <<= 1;}if (j == r){return false;}}}return true;
}