大素数判断_fermat素性测试+Miller-Rabin素性测试

一、朴素的判断一个数是否为素数：

原理：若一个数为合数，那么必然存在这样的两个数：2<=a<=sqrt(n) <=b<n，使得n=a*b。

解法：从 2 到 sqrt(n) 枚举，若存在数字 a 为数 n 的因子，那么数字 n 即为合数。若不存在，则数字 n 为偶数。

代码：

isprime(i)==1  表示数字 i 为质数

bool isprime(int n)
{
  if(n<2)return 0;
  if(n==2)return 1;
  if(n%2==0)return 0;
  /*
    去掉n为偶数的情况，那么n的因子就只能是奇
    数，下方枚举时，就可以只枚举奇数因子。 
  */
  int i,j,k;
  k=(int)sqrt(n*1.0);
  /*
    这里需要特别注意，sqrt里面的参数为实数，如果参数为整数，
    可能会报错。记得我有一次在hdoj上做了一题，用c++可以编译
    通过，用G++就不行，找了好久，才发现是这个参数问题。
  */ 
  for(i=3;i<=k;i+=2)
    if(n%i==0)return 0;
  return 1;  
}

以下的讲解为简化版本，更加详细的请参考：http://www.matrix67.com/blog/archives/234 

其次，fermat素性测试与Miller-Rabin素性测试的结果均为不确定的，即不保证完全正确。保证完全正确的算法有：AKS算法 （发表的相应论文名：ＰＲＩＭＥ　ｉｓ　ｉｎ　Ｐ）。

二、fermat素性测试

原理：费马小定理：

           假如p是质数，且Gcd(a,p)=1，那么 a(p-1) ≡1（mod p）。即：假如a是整数，p是质数，且a,p互质(即两者只有一个公约数1)，那么a的(p-1)次方除以p的余数恒等于1。

思路：那么，根据费马小定理的逆定理，如果们找到一个 a ，使得 a ^(n-1)%n !=1，是否就可以确定数字 n 不是质数呢？很遗憾，这是不行的，费马小定理的逆定理并不正确。

            伪素数：满足2^(n-1)%n==1的合数 n 。

                          满足a^(n-1)%n==1的合数n叫做以 a 为底的伪素数。

            Carmichael数：对于所有小于 n 的底数 a，都满足a^(n-1)%n==1 的数。（前10亿个数中，仅有600多个这样的数存在）。

            

            在具体实现用fermat素性测试中，我们一般是随机选取有限个数的底数 a ，对 n 进行素性测试。若全部满足，则认为数字 n 是质数，若有一个不满足，则认为数字 n 不是质数。

三、Miller-Rabin素性测试。

原理：如果p是素数，x是小于p的正整数，且x^2 % p==1，那么x==1 或者 x==p-1。

思路：Miller-Rabin素性测试是对fermat素性测试的优化，极大地提高了正确性，虽然仍是一个不确定算法。

          我们以 a==2，n==341为例，演示一下该测试是如何进行的。（2^340%341==1，但是341并不是一个质数）

          2^340%341==1  ==> (2^170%341)^2 %341==1

==>    2^170%341==1    或者    2^170%341==340，而2^170%341==1，定理继续适用

         2^170%341==1  ==>  (2^85%341)^2  %341 ==1

==>   2^85%341==1   或者  2^85%341==340 ，很遗憾的是，两个都不成立，与上述所提到的原理相悖，所以341不是素数。

         

         Millar-Rabin素性测试： n-1=d*2^r，如果n是一个素数，那么a^d%n==1，或者是存在某个 i 使得 a^(d*2^i)%n==n-1（0<=i<r）。

        强伪素数：通过以 a 为底的Millar-Rabin测试的合数称为以 a 为底的强伪素数。

代码： 

对底数 a 的选择，我是直接求10^5以内的素数，并选取 1 ——> maxn3 的素数为底，当然，你也可以随机选取。

在使用的时候，要注意 n*n 不能超过 n 的数据类型范围，比如int n，那么 n*n就不能超过 int 范围。

#include<cstdio>
#define maxn1 50000
#define maxn2 10000
#define maxn3 5000
using namespace std;

bool f[maxn1+100];
int p[maxn2];

void pre_a()
{
  int i,j,k,x,y,z;
  p[++p[0]]=2;
  for(i=1;i<=maxn1;i++)
    {
      k=i*2+1;
      if(!f[i])p[++p[0]]=k;
      for(j=2;j<=p[0] && k*p[j]/2<=maxn1;j++)
        {
          f[k*p[j]/2]=1;
          if(k%p[j]==0)break;
        }
    }
}

int pp(int a,long long d,long long n)
{
  long long ans=1,x=a;
  while(d>0)
    {
      if(d&1)ans=(ans*x)%n;
      x=(x*x)%n,d>>=1;
    }
  return ans;  
}

bool ok(int a,long long n)
{
  long long d,t;
  d=n-1;
  while((d&1)==0)d>>=1;
  t=pp(a,d,n);
  while(d!=n-1 && t!=1 && t!=n-1)
    t=(t*t)%n,d<<=1;
  return t==n-1 || d&1;  
}

bool isprime(long long n)
{
  if(n<2)return 0;
  if(n==2)return 1;
  if((n&1)==0)return 0;
  if(n/2<=maxn1)return !f[n/2];
  
  for(int i=1;i<=maxn3;i++)
    if(!ok(p[i],n))return 0;
  return 1;
} 

int main()
{
  pre_a();
  long long n;
  scanf("%I64d",&n);
  if(isprime(n))printf("n是质数\n");
  else printf("n是合数\n");
  return 0;
}