费马小定理 素数判定 蒙哥马利算法

约定： x%y为x取模y，即x除以y所得的余数，当x 1) {//一直计算，直到指数小于或等于1 if((p%2)!=0) {// 如果指数p是奇数，则说明计算后会剩一个多余的数，那么在这里把它乘到结果中 odd*=main; //把“剩下的”乘起来 } main*=main; //主体乘方 p/=2; //指数除以2 } return main*odd; //最后把主体和“剩下的”乘起来作为结果 } 够完美了吗？不，还不够！看出来了吗？main是没有必要的，并且我们可以有更快的代码来判断奇数。要知道除法或取模运算的效率很低，所以我们可以利用偶数的一个性质来优化代码，那就是偶数的二进制表示法中的最低位一定为0! 完美版： unsigned Power(unsigned n, unsigned p) { // 计算n的p次方 unsigned odd = 1; //oddk用来计算“剩下的”乘积 while (p > 1) { // 一直计算到指数小于或等于1 if (( p & 1 )!=0) { // 判断p是否奇数，偶数的最低位必为0 odd *= n; // 若odd为奇数，则把“剩下的”乘起来 } n *= n; // 主体乘方 p /= 2; // 指数除以2 } return n * odd; // 最后把主体和“剩下的”乘起来作为结果 } ======================================================== 蒙格马利”快速幂模算法 后面我们会用到这样一种运算：(X^Y)%Z 问题是当X和Y很大时，只有32位的整型变量如何能够有效的计算出结果？ 考虑上面那份最终的优化代码和再上面提到过的积模分解公式，我想你也许会猛拍一下脑门，吸口气说：“哦，我懂了！”。 下面的讲解是给尚没有做出这样动作的同学们准备的。X^Y可以看作Y个X相乘，即然有积模分解公式，那么我们就可以把Y个X相乘再取模的过程分解开来，比如：(17^25)%29则可分解为：( ( 17 * 17 ) % 29 * ( 17 * 17 ) % 29 * …… 如果用上面的代码将这个过程优化，那么我们就得到了著名的“蒙格马利”快速幂模算法： unsigned Montgomery(unsigned n, unsigned p, unsigned m) { // 快速计算 (n ^ e) % m 的值，与power算法极类似 unsigned r = n % m; // 这里的r可不能省 unsigned k = 1; while (p > 1) { if ((p & 1)!=0) { k = (k * r) % m; // 直接取模 } r = (r * r) % m; // 同上 p /= 2; } return (r * k) % m; // 还是同上 } 上面的代码还可以优化。下面是蒙格马利极速版： 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; } ===================================================== 怎么判断一个数是否为素数？ 笨蛋的作法： bool IsPrime(unsigned n) { if (n<2) { //小于2的数即不是合数也不是素数 throw 0; } for (unsigned i=2;i 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 >= 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; }