【POJ1811】Prime Test-Miller-Rabin素数测试+Pollard-rho大数分解

测试地址:Prime Test
题目大意:给出T个整数N2N254),判断它们是不是素数,如果不是,输出它的最小质因子。
做法:这一题需要使用Miller-Rabin素数测试和Pollard-rho大数分解算法。
题如其名,思路题目中都告诉你了,先判断是不是素数,如果不是再找它的最小质因子。但是看到这个庞大的数据范围,就连判断素数这一个步骤,平常的筛法复杂度也是O(N),铁定爆炸,这时候我们就要采用一种虽然不稳定,但大概率能判断一个数是不是素数的方法:Miller-Rabin素数测试。
我们知道费马小定理:如果p是素数,那么对于所有0<a<pap11(modp)。我们可以把它看做判断一个数是不是素数的必要条件。但是如果枚举所有a的话,算法的时间复杂度就很差了,所以我们可以随机取底数a,然后判断以上式子是否成立,如果不成立则表明该数一定是合数,否则说明该数可能是素数。选取的a的个数越多,该数为素数的可能性就越大。事实上,用前10个素数作为a的话,在262范围内判错的概率已经小到可以忽略不计了。
然而存在一种数,即使对于所有的a都满足上述式子,但是它实际上还是一个合数,而且这种数出现的频率还挺高,这种数叫做Camichael数。这时候就要拿出另一个定理——二次探测定理:设有一个素数p,如果一个数x满足x21(modp),那么要么x=1,要么x=p1。这个定理拿到这里来有什么用呢?这个实际上就是添加了必要条件。如果我们经过测试得到ap11(modp),如果p1为偶数(因为p是素数所以这是大概率的),那么可以把p1表示成d×2x这样的形式,那么我们就可以从ad×2x1=ap1开始,判断这个数是不是等于1p1,如果不等于,那么这个数铁定就是个合数了,否则如果这个数还为1,就继续探测ad×2x2……这样一直下去,直到指数与2互质或者这个数等于p1为止。这样我们就对Miller-Rabin素数测试进行了强化,称为改进Miller-Rabin素数测试,这种测试对于题目中的这个数据范围已经毫无压力了。
判断素数的问题解决了,接下来要找一个合数的最小质因子。对于这么大的数据范围,我们应该采用一种叫做Pollard-rho的大数分解算法,然后找到一个数的最小质因子(实际上,这个算法理论上可以在O(N4)的时间内找到一个数的所有质因子)。由于篇幅原因,以及这个算法的证明比较复杂,这里就不再写了。而算法中的素数判断就用Miller-Rabin即可。综合上述两种方法,在某一些地方适当使用随机化,就可以通过此题了。要注意的是,由于题目中的数字可能很大,两个数相乘可能会爆long long,所以我们就仿照快速幂写一个快速乘来代替原先的乘法即可。
以下是本人代码:

#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <iostream>
#include <algorithm>
#define ll long long
using namespace std;
ll T,N,minfactor,p[11]={2,3,5,7,11,13,17,19,23,29};

ll mult(ll a,ll b,ll mod) //计算a*b%mod
{
  ll s=a,sum=0;
  while(b)
  {
    if (b&1)
    {
      sum+=s;
      if (sum>=mod) sum-=mod; //这里没有必要使用取模,取模运算比减法慢很多
    }
    b>>=1;s<<=1;if (s>=mod) s-=mod;
  }
  return sum;
}

ll power(ll a,ll b,ll mod) //计算a^b%mod
{
  ll s=a,sum=1;
  while(b)
  {
    if (b&1) sum=mult(sum,s,mod);
    b>>=1;s=mult(s,s,mod);
  }
  return sum;
}

bool witness(ll n,ll a) //对整数n,底数a进行测试,返回1表示通过测试
{
  ll p=power(a,n-1,n);
  if (p!=1) return 0;
  else
  {
    ll s=n-1;
    while(!(s%2)&&p==1)
    {
      s>>=1;
      p=power(a,s,n);
    }
    if (p==1||p==n-1) return 1;
    else return 0;
  }
}

bool miller_rabin(ll n) //对整数n进行Miller-Rabin素数测试,返回1表示通过测试
{
  if (n<=29)
  {
    for(int i=0;i<10;i++)
      if (n==p[i]) return 1;
    return 0;
  }
  for(int i=0;i<10;i++)
    if (!witness(n,p[i])) return 0;
  return 1;
}

ll gcd(ll a,ll b) //Euclid算法计算a和b的最大公因数,很经典了
{
  return (b==0)?a:gcd(b,a%b);
}

ll pollard_rho(ll n,ll c) //返回n的一个因子,如果返回n表示失败
{
  ll x=rand()%n,y=x,d,i=1,k=2;
  while(1)
  {
    i++;
    x=(mult(x,x,n)+c)%n;
    d=gcd(y-x+n,n);
    if (d>1&&d<n) return d;
    if (y==x) return n;
    if (i==k) {y=x;k<<=1;}
  }
}

void find_factor(ll n) //对整数n进行分解
{
  if (miller_rabin(n))
  {
    minfactor=min(minfactor,n);
    return;
  }
  ll p=n;
  while(p>=n) p=pollard_rho(n,rand()%(n-1)+1);
  find_factor(p);
  find_factor(n/p);
}

int main()
{
  scanf("%lld",&T);
  while(T--)
  {
    scanf("%lld",&N);
    if (miller_rabin(N)) printf("Prime
");
    else
    {
      minfactor=N;
      find_factor(N);
      printf("%lld
",minfactor);
    }
  }

  return 0;
}
原文地址:https://www.cnblogs.com/Maxwei-wzj/p/9793681.html