如何判断一个数是否为素数
怎么判断一个数是否为素数？
笨蛋的作法：
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;
}
一个数去除以比它的一半还要大的数，一定除不尽，所以还用判断吗？？
下面是小学生的做法：
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; // 都没除尽，素数
}
一个合数必然可以由两个或多个质数相乘而得到。那么如果一个数不能被比它的一半小的所有的质数整除，那么比它一半小的所有的合数也一样不可能整除它。建立一个素数表是很有用的。
下面是中学生的做法：
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;
}
