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

2019-04-13 15:45发布

站内文章 / 模拟电子
16224 0 1246
转载于http://blog.csdn.net/arvonzhang/article/details/8564836


1.约定

x%y为x取模y，即x除以y所得的余数，当x象都为整数。
x^y表示x的y次方。乘方运算的优先级高于乘除和取模，加减的优先级最低。
见到x^y/z这样，就先算乘方，再算除法。
A/B，称为A除以B，也称为B除A。
若A%B=0，即称为A可以被B整除，也称B可以整除A。
A*B表示A乘以B或称A乘B，B乘A，B乘以A……都一样。 复习一下小学数学
公因数：两个不同的自然数A和B，若有自然数C可以整除A也可以整除B，那么C就是A和B的公因数。 公倍数：两个不同的自然数A和B，若有自然数C可以被A整除也可以被B整除，那么C就是A和B的公倍数。 互质数：两个不同的自然数，它们只有一个公因数1，则称它们互质。 费马是法国数学家，又译“费尔马”，此人巨牛，他的简介请看下面。不看不知道，一看吓一跳。 费马人物简介：http://baike.baidu.com/view/6303430.htm?fromId=5194&redirected=seachword　　

2.费马小定理：

有N为任意正整数，P为素数，且N不能被P整除（显然N和P互质），则有：N^P%P=N(即：N的P次方除以P的余数是N)。 但是我查了很多资料见到的公式都是这个样子： (N^(P-1))%P=1后来分析了一下，两个式子其实是一样的，可以互相变形得到。 原式可化为：(N^P-N)%P=0(即：N的P次方减N可以被P整除，因为由费马小定理知道N的P次方除以P的余数是N)把N提出来一个，N^P就成了你N*(N^(P-1))，那么(N^P-N)%P=0可化为： (N*(N^(P-1)-1))%P=0
请注意上式，含义是：N*(N^(P-1)-1)可以被P整除 又因为N*(N^(P-1)-1)必能整除N（这不费话么!）
所以，N*(N^(P-1)-1)是N和P的公倍数，小学知识了^_^ 又因为前提是N与P互质，而互质数的最小公倍数为它们的乘积，所以一定存在 正整数M使得等式成立：N*(N^(P-1)-1)=M*N*P
两边约去N，化简之：N^(P-1)-1=M*P
因为M是整数，显然：N^(P-1)-1)%P=0即：N^(P-1)%P=1
============================================

3.积模分解公式

先有一个引理，如果有：X%Z=0，即X能被Z整除，则有：(X+Y)%Z=Y%Z
设有X、Y和Z三个正整数，则必有：(X*Y)%Z=((X%Z)*(Y%Z))%Z 想了很长时间才证出来，要分情况讨论才行： 1.当X和Y都比Z大时，必有整数A和B使下面的等式成立：
X=Z*I+A（1）
Y=Z*J+B（2）
不用多说了吧，这是除模运算的性质！
将（1）和（2）代入(X*Y)modZ得：((Z*I+A)(Z*J+B))%Z乘开，再把前三项的Z提一个出来，变形为：(Z*(Z*I*J+I*A+I*B)+A*B)%Z（3） 因为Z*(Z*I*J+I*A+I*B)是Z的整数倍……晕，又来了。
概据引理，（3）式可化简为：(A*B)%Z又因为：A=X%Z，B=Y%Z，代入上面的式子，就成了原式了。 2.当X比Z大而Y比Z小时，一样的转化：
X=Z*I+A
代入(X*Y)%Z得：
(Z*I*Y+A*Y)%Z
根据引理，转化得：(A*Y)%Z
因为A=X%Z，又因为Y=Y%Z，代入上式，即得到原式。
同理，当X比Z小而Y比Z大时，原式也成立。 3.当X比Z小，且Y也比Z小时，X=X%Z，Y=Y%Z，所以原式成立。
=====================================================

4.快速计算乘方的算法

如计算2^13，则传统做法需要进行12次乘法。 [cpp] view plain copy
  1. /*计算n^p*/  
  2. unsigned power(unsigned n,unsigned p)  
  3. {  
  4.     for(int i=0;i
  5.     return n;  
  6. }  
该死的乘法，是时候优化一下了！把2*2的结果保存起来看看，是不是成了： 4*4*4*4*4*4*2
再把4*4的结果保存起来：16*16*16*2 
一共5次运算，分别是2*2、4*4和16*16*16*2 这样分析，我们算法因该是只需要计算一半都不到的乘法了。
为了讲清这个算法，再举一个例子2^7：2*2*2*2*2*2*2 
两两分开：(2*2)*(2*2)*(2*2)*2 
如果用2*2来计算，那么指数就可以除以2了，不过剩了一个，稍后再单独乘上它。 再次两两分开，指数除以2： ((2*2)*(2*2))*(2*2)*2
实际上最后一个括号里的2 * 2是这回又剩下的，那么，稍后再单独乘上它 现在指数已经为1了，可以计算最终结果了：16*4*2=128 优化后的算法如下： [cpp] view plain copy
  1. unsigned Power(unsigned n,unsigned p)   
  2. {  
  3.    unsigned main=n; //用main保存结果  
  4.    unsigned odd=1; //odd用来计算“剩下的”乘积  
  5.    while (p>1)   
  6.    {//一直计算，直到指数小于或等于1  
  7.         if((p%2)!=0)  
  8.         {// 如果指数p是奇数，则说明计算后会剩一个多余的数，那么在这里把它  
  9. 乘到结果中  
  10.             odd*=main; //把“剩下的”乘起来  
  11.         }  
  12.         main*=main; //主体乘方  
  13.         p/=2; //指数除以2  
  14.    }  
  15.    return main*odd; //最后把主体和“剩下的”乘起来作为结果  
  16. }  
够完美了吗？不，还不够！看出来了吗？main是没有必要的，并且我们可以有更快的代码来判断奇数。要知道除法或取模运算的效率很低，所以我们可以利用偶数的一个性质来优化代码，那就是偶数的二进制表示法中的最低位一定为0! 完美版： [cpp] view plain copy
  1. unsigned Power(unsigned n, unsigned p)   
  2. { // 计算n的p次方  
  3.     unsigned odd = 1; //oddk用来计算“剩下的”乘积  
  4.     while (p > 1)  
  5.     { // 一直计算到指数小于或等于1  
  6.        if (( p & 1 )!=0)  
  7.       { // 判断p是否奇数，偶数的最低位必为0  
  8.              odd *= n; // 若odd为奇数，则把“剩下的”乘起来  
  9.       }  
  10.       n *= n; // 主体乘方  
  11.       p /= 2; // 指数除以2  
  12.      }  
  13.     return n * odd; // 最后把主体和“剩下的”乘起来作为结果  
  14. }  

========================================================

5."蒙格马利”快速幂模算法

后面我们会用到这样一种运算：(X^Y)%Z。但问题是当X和Y很大时，只有32位的整型变量如何能够有效的计算出结果？ 考虑上面那份最终的优化代码和再上面提到过的积模分解公式，我想你也许会猛拍一下脑门，吸口气说：“哦，我懂了！”。 下面的讲解是给尚没有做出这样动作的同学们准备的： X^Y可以看作Y个X相乘，即然有积模分解公式，那么我们就可以把Y个X相乘再取模的过程分解开来，比如：(17^25)%29则可分解为：( ( 17 * 17 ) % 29 * ( 17 * 17 ) % 29 * …… 如果用上面的代码将这个过程优化，那么我们就得到了著名的“蒙格马利”快速幂模算法： [cpp] view plain copy
  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. }  

上面的代码还可以优化。下面是蒙格马利极速版： [cpp] view plain copy
  1. unsigned Montgomery(unsigned n,unsigned p,unsigned m)  
  2. { //快速计算(n^p)%m的值  
  3.       unsignedk=1;  
  4.       n%=m;  
  5.      while(p!=1)  
  6.      {  
  7.          if(0!=(p&1))k=(k*n)%m;  
  8.          n=(n*n)%m;  
  9.          p>>=1;  
  10.     }  
  11.     return(n*k)%m;  
  12. }  

=====================================================

6.怎么判断一个数是否为素数？

1)笨蛋的作法：

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

2)下面是小学生的做法：

[cpp] view plain copy
  1. bool IsPrime(unsigned n)  

热门文章