1.约定
x%y为x取模y，即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）
不用多说了吧，这是除模运算的性质！ ( X +
Y )%Z
将（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的整数倍……晕，又来了。∴Z*(Z*I*J+I*A+I*B) %Z=0 （X%Z=0）
概据引理，（3）式（ ( X +
Y )%Z ）可化简为：(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次乘法。
plaincopy
- /*计算n^p*/
- unsigned power(unsigned n,unsigned p)
- {
- for(int i=0;i<p;i++) n*=n;
- return n;
- }
该死的乘法，是时候优化一下了！把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
优化后的算法如下：
plaincopy
- unsigned Power(unsigned n,unsigned p)
- {
- unsigned main=n; //用main保存结果
- unsigned odd=1; //odd用来计算“剩下的”乘积
- while (p>1)
- {//一直计算，直到指数小于或等于1
- if((p%2)!=0)
- {// 如果指数p是奇数，则说明计算后会剩一个多余的数，那么在这里把它
- 乘到结果中
- odd*=main; //把“剩下的”乘起来
- }
- main*=main; //主体乘方
- p/=2; //指数除以2
- }
- return main*odd; //最后把主体和“剩下的”乘起来作为结果
- }
够完美了吗？不，还不够！看出来了吗？main是没有必要的，并且我们可以有更快的代码来判断奇数。要知道除法或取模运算的效率很低，所以我们可以利用偶数的一个性质来优化代码，那就是偶数的二进制表示法中的最低位一定为0!
完美版：
plaincopy
- 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; // 最后把主体和“剩下的”乘起来作为结果
- }
========================================================
5."蒙格马利”快速幂模算法
后面我们会用到这样一种运算：(X^Y)%Z。但问题是当X和Y很大时，只有32位的整型变量如何能够有效的计算出结果？
考虑上面那份最终的优化代码和再上面提到过的积模分解公式，我想你也许会猛拍一下脑门，吸口气说：“哦，我懂了！”。
下面的讲解是给尚没有做出这样动作的同学们准备的：
X^Y可以看作Y个X相乘，即然有积模分解公式，那么我们就可以把Y个X相乘再取模的过程分解开来，比如：(17^25)%29则可分解为：( ( 17 * 17 ) % 29 * ( 17 * 17 ) % 29 * ……
如果用上面的代码将这个过程优化，那么我们就得到了著名的“蒙格马利”快速幂模算法：
plaincopy
- unsigned Montgomery(unsigned n, unsigned p, unsigned m)
- { // 快速计算 (n ^ p) % 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^p)%m的值
unsigned k=1;
n%=m;
while(p!=1)
{
if(0!=(p&1)) //如果p是奇数
{ k=(k*n)%m; } // 上面有个推论 注意下面这一步很多都没提到这点
n=(n*n)%m; // 5^8 mod 3 == (5*5)^4 mod 3 == ((5*5)mod 3)^4 mod 3
p>>=1;
}
return (n*k)%m;
}
=====================================================
6.怎么判断一个数是否为素数？
1)笨蛋的作法：
plaincopy
- bool IsPrime(unsigned n)
- {
- if (n<2)
- {
- //小于2的数即不是合数也不是素数
- throw 0;
- }
- for (unsigned i=2;i<n;++i)
- {
- //和比它小的所有的数相除，如果都除不尽，证明素数
- if (n%i==0)
- {
- //除尽了，则是合数
- return false;
- }
- }
- return true;
- }
一个数去除以比它的一半还要大的数，一定除不尽，所以还用判断吗？？
2)下面是小学生的做法：
plaincopy
- bool IsPrime(unsigned n)
- {
- if (n<2)
- {
- //小于2的数即不是合数也不是素数
- throw 0;
- }
- for(unsigned i=2;i<n/2+1;++i)
- {
- // 和比它的一半小数相除，如果都除不尽，证明素数
- if ( 0 == n % i )
- {
- // 除尽了，合数
- return false;
- }
- }
- return true; // 都没除尽，素数
- }
一个合数必然可以由两个或多个质数相乘而得到。那么如果一个数不能被比它的一半小的所有的质数整除，那么比它一半小的所有的合数也一样不可能整除它。建立一个素数表是很有用的。
3)下面是中学生的做法：
plaincopy
- 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进行整除测试，将会大大提高通过率，方法如下:
6)下面是专家的做法：
plaincopy
- 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这三个数，我其实是故意的，我只是想说明并费马测试并不完全可靠。
#include<iostream>
using namespace std;<pre name="code" class="cpp" style="font-size: 14px; line-height: 26px;">//快速计算(n^p)%m的值
typedef unsigned long LL;LL Montgomery(LL n,LL p,LL m)
{
LL k=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 IsPrime3(LL n)
{ if ( n < 2 ) { // 小于2的数即不是合数也不是素数 throw 0; } static LL aPrimeList[] = { 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 }; const int nListNum = sizeof(aPrimeList)
/ sizeof(LL); for (int i=0;i<nListNum;++i) { // 按照素数表中的数对当前素数进行判断//对于两个互质的素数N和P，必有：N^(P-1)%P=1 , 这里的N是素数表里的,P是我们要进行测试的//当然由于 是素数推出 ==1 而==1不能推出是素数 但是我们进行多组测试(50组)的话还是能保证其概率的if(n==aPrimeList[i])return true;//(根据下面的注意1) if (1!=Montgomery(aPrimeList[i],n-1,n))
// 蒙格马利算法 {//<>1说明不互质 那么就是合数 直接返回//当然这边有个漏洞 如果n就是上面素数表里面的如5 5^4%5!=1 return false 所以要在上面就判断是否是素数表里的(注意1) return false; } } return true;}int main(){LL n;while(cin>>n){if(IsPrime3(n))cout<<"It is a prime number."<<endl;
else cout<<"It is not a prime number"<<endl;}return 0;}
等等，什么？好像有点怪，看一下这个数29341，它等于13 * 37 * 61，显然是一个合数，但是竟通过了测试！！哦，抱歉，我忘了在素数表中加入13，37，61这三个数，我其实是故意的，我只是想说明并费马测试并不完全可靠。
现在我们发现了重要的一点，费马定理是素数的必要条件而非充分条件。这种不是素数，但又能通过费马测试的数字还有不少，数学上把它们称为卡尔麦克数，现在数学家们已经找到所有10 ^ 16以内的卡尔麦克数，最大的一个是9585921133193329。我们必须寻找更为有效的测试方法。数学家们通过对费马小定理的研究，并加以扩展，总结出了多种快速有效的素数测试方法，目前最快的算法是拉宾米勒测试算法，下面介绍拉宾米勒测试。
================================================================
7.拉宾米勒测试
拉宾米勒测试是一个不确定的算法，只能从概率意义上判定一个数可能是素数，但并不能确保。算法流程如下:
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;
}
我们发现用随机数种子的话 会比较好 修改如下:
#include<iostream>
#include<ctime>
using namespace std;
typedef unsigned long long LL; //快速计算(n^p)%m的值
LL Montgomery(LL n,LL p,LL m)
{
LL k=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 Miller_Rabin(LL x){
if(x==2) return true;
if(!(x&1)) return false; //除了2的偶数直接返回
int t=10;
LL a,ans;
while(t--){
//对于两个互质的素数a和x，必有：a^(x-1)%x=1 , 这里的a是随机出来的,x是我们要进行测试的
//当然由于 是素数推出 ==1 而==1不能推出是素数 但是我们进行多组测试(50组)的话还是能保证其概率的
a=rand()%(x-2)+2;
ans=Montgomery(a,x-1,x);
if(ans!=1) return false; }
return true;
}
int main(){
LL n;
srand(time(NULL));//置时间种子
while(cin>>n){
if(Miller_Rabin(n))cout<<"It is a prime number."<<endl;
else cout<<"It is not a prime number"<<endl;
}
return 0;
}
另外一种情况 大素数
文字代码：
1：重复MAX次运算
2：在1~n中取得随机数a
3：计算a^(n-1)mod n？=1，在这个计算里，注意到n可能很大，所以a^(n-1)可能越界，就想到用快速幂来边乘，边取模，但是又发现在n很大的时候，a*a都有可能溢出，所以想到了用快速幂的方法，进行快速积取模，边加边取模。这里的两个快速可避免溢出(较上面的方法可解决溢出)
4：在3中可得到的数如果为1，则在循环未结束前继续从2开始操作，否则直接返回0，表示n是合数
5：如果上面的循环能完整做完，说明n已经是强伪素数，我们可以返回1，判定为素数
#include<iostream>
#include<ctime>
using namespace std;
typedef unsigned __int64 ULL;
ULL MMOD(ULL a,ULL b,ULL c)
{
ULL ans;
for(ans=0;b;b>>=1)
{
if(b&1)
{
ans=ans+a;
while(ans>=c)
ans-=c;
}
a<<=1;
while(a>=c) a-=c;
}
return ans;
} ULL QMod(ULL a,ULL b,ULL c)
{
ULL ans;
for(ans=1;b;b>>=1)
{
if(b&1)
ans=MMOD(ans,a,c);
a=MMOD(a,a,c); }
return ans;
}
bool Miller_Rabin(ULL p)
{
if(p<2) return false;
if(p==2) return true;
int t=5;
ULL ans,a;
while(t--)
{
a=rand()%(p-1)+1;
ans=QMod(a,p-1,p);
if(ans!=1)
return false;
}
return true;
}
int main(){
ULL n;
srand(time(NULL));
while(~scanf("%I64u",&n)){
if(Miller_Rabin(n))
printf("It is a prime number.\n");
else
printf("It is not a prime number.\n");
}
return 0;
}
版权声明：本文为博主原创文章，未经博主允许不得转载。