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

测试地址：Prime Test
题目大意：给出T个整数N（2≤N≤254），判断它们是不是素数，如果不是，输出它的最小质因子。
做法：这一题需要使用Miller-Rabin素数测试和Pollard-rho大数分解算法。
题如其名，思路题目中都告诉你了，先判断是不是素数，如果不是再找它的最小质因子。但是看到这个庞大的数据范围，就连判断素数这一个步骤，平常的筛法复杂度也是O(N−−√)，铁定爆炸，这时候我们就要采用一种虽然不稳定，但大概率能判断一个数是不是素数的方法：Miller-Rabin素数测试。
我们知道费马小定理：如果p是素数，那么对于所有0<a<p，ap−1≡1(modp)。我们可以把它看做判断一个数是不是素数的必要条件。但是如果枚举所有a的话，算法的时间复杂度就很差了，所以我们可以随机取底数a，然后判断以上式子是否成立，如果不成立则表明该数一定是合数，否则说明该数可能是素数。选取的a的个数越多，该数为素数的可能性就越大。事实上，用前10个素数作为a的话，在262范围内判错的概率已经小到可以忽略不计了。
然而存在一种数，即使对于所有的a都满足上述式子，但是它实际上还是一个合数，而且这种数出现的频率还挺高，这种数叫做Camichael数。这时候就要拿出另一个定理——二次探测定理：设有一个素数p，如果一个数x满足x2≡1(modp)，那么要么x=1，要么x=p−1。这个定理拿到这里来有什么用呢？这个实际上就是添加了必要条件。如果我们经过测试得到ap−1≡1(modp)，如果p−1为偶数（因为p是素数所以这是大概率的），那么可以把p−1表示成d×2x这样的形式，那么我们就可以从ad×2x−1=ap−1−−−−√开始，判断这个数是不是等于1或p−1，如果不等于，那么这个数铁定就是个合数了，否则如果这个数还为1，就继续探测ad×2x−2……这样一直下去，直到指数与2互质或者这个数等于p−1为止。这样我们就对Miller-Rabin素数测试进行了强化，称为改进Miller-Rabin素数测试，这种测试对于题目中的这个数据范围已经毫无压力了。
判断素数的问题解决了，接下来要找一个合数的最小质因子。对于这么大的数据范围，我们应该采用一种叫做Pollard-rho的大数分解算法，然后找到一个数的最小质因子（实际上，这个算法理论上可以在O(N−−√4)的时间内找到一个数的所有质因子）。由于篇幅原因，以及这个算法的证明比较复杂，这里就不再写了。而算法中的素数判断就用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\n");
    else
    {
      minfactor=N;
      find_factor(N);
      printf("%lld\n",minfactor);
    }
  }

  return 0;
}