[POJ 1811]Prime Test---Miller-Rabin算法&Pollard-rho算法

Prime Test

Time Limit: 6000MS		Memory Limit: 65536K
Total Submissions: 33257		Accepted: 8787
Case Time Limit: 4000MS

Description
Given a big integer number, you are required to find out whether it's a prime number.
Input
The first line contains the number of test cases T (1 <= T <= 20 ), then the following T lines each contains an integer number N (2 <= N < 254).
Output
For each test case, if N is a prime number, output a line containing the word "Prime", otherwise, output a line containing the smallest prime factor of N.
Sample Input

2
5
10

Sample Output

Prime
2

Source
POJ Monthly

题意：给你T和T个数据，问数据是否为素数，若是，输出“Prime”，否则输出其最小质因数。。。

显然N迷之大，朴素的只能卡一点点。。。

这就要用莫名其妙异常难懂的Miller-Rabin算法和Pollard-rho算法了。。。

那么。。。先是Miller-Rabin。。。

Miller-Rabin算法:

这是一种神奇的随机的素数判定算法，也就是说，答案可能出错。。。但是可能性很小。。。

先是两个神奇的定理：

费马小定理：百度百科

对于一个质数p，若取任意整数a，满足gcd(p,a)=1，则有

ap−1≡1(modp)
二次探测定理： （这个没百科）

对于0<x<p，若p是素数，则方程：

x2≡1(modp)

的解为：

x1=1,x2=p−1
费马小定理的逆命题不成立，但否逆命题成立，所以可以利用以下一点：

对于任意整数a<p，不满足ap−1≡1(modp)，则p为合数。

所以我们可以不断在区间[2,p−1] 
范围内随机取a，并进行判定。在s次（好像20就差不多了。。。）判定不为合数之后，我们就可以说这个数是质数。

这显然还不够精确，so我们可以先把p−1分解成2t×u(p-1不停除2直到u为奇数为止)     
的形式，然后令x[0]=aumodp,，那么将x[0]平方t
次就是(au)2tmodp    
的值，我们设x[i]为x[0]平方i次的值，根据二次探测定理，若x[i]等于1，则x[i−1]等于1或p-1，不满足则p为合数。

显然类似abmodp的式子都要用快速幂运算，当n太大时，形如a×bmodp   
的式子都要用分治的思想来计算，不然还不如暴枚。。。

时间复杂度：带常数为                          ，误差率为1/2^s 
O(slog3n)

代码如下：

const int MAXN = 65;
long long n, x[MAXN];

<pre name="code" class="cpp">long long multi(long long a,long long b,long long p){//快速乘。。。
long long ans=0;
while(b){
if(b&1LL)
ans=(ans+a)%p;
a=(a+a)%p;
b>>=1;
}
return ans;//快速幂乘法版,顺手模上p 。。。
}

long long qpow(long long a,long long b,long long p){//啊哈和上面迷之相似 。。。
long long ans=1;
while(b){
if(b&1LL)
ans=multi(ans,a,p);
a=multi(a,a,p);
b>>=1;
}
return ans;
}

bool miller_rabin(long long n){
if(n==2)
return true;//显然是特判。。。
int s=20,t=0;//s为判断次数（精度），t为2的次方数
long long u=n-1;//u初始化为n/2^1
while(!(u&1)){//当u还为偶数
t++;
u>>=1;//+1次方并让u除以2
}
while(s--){//判断s次
long long a=rand()%(n-2)+2;//随机0+2~n-3+2区间取a
x[0]=qpow(a,u,n);//初始化为a的u次方模n
for(int i=1;i<=t;i++){//一次次乘上2
x[i]=multi(x[i-1],x[i-1],n);
if(x[i]==1 && x[i-1]!=1 && x[i-1]!=n-1)
return false;//二次判断区
}
if(x[t]!=1)//同上
return false;
}
return true;
}

Pollard-Rho算法： 

对于一个大整数n，若我们取任意一个数x使得x是n
的质因数的几率很小，但如果取两个数x1以及x2使得它们的差是n的因数的几率就提高了（原因未知。。。），如果取x1以及x2使得gcd(abs(x1−x2),n)>1 
的概率就更高了（依旧未知。。。）。这就是Pollard-Rho算法的主要思想。。。

对于满足gcd(abs(x1−x2),n)>1的x1和x2，gcd(abs(x1−x2),n)就是n的一个因数，只需要判断它是否为素数，若为素数，则是n的质因数，否则递归继续。。。

其中判断素数就要用上面的Miller-Rabin
算法。

那么我们怎样不断取得x1和x2呢？x1在区间[1,n]中随机出来，而x2则由

x[i]=(x[i-1]*x[i-1]%n+c)%n

推算出来，其中c为任意给定值，事实（才怪）证明，这样就是比较(kan)优(lian)的。

这个。。。就不写代码了。。。反正和正文毫无差别。。。（偷懒）

正文代码：

#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<iostream>
using namespace std;

long long x[65],ans;

long long multi(long long a,long long b,long long p){//又是快速乘。。。
long long ans=0;
while(b){
if(b&1LL)
ans=(ans+a)%p;
a=(a+a)%p;
b>>=1;
}
return ans;
}

long long qpow(long long a,long long b,long long p){//再次迷之相似 。。。
long long ans=1;
while(b){
if(b&1LL)
ans=multi(ans,a,p);
a=multi(a,a,p);
b>>=1;
}
return ans;
}

bool miller_rabin(long long n){
if(n==2)
return true;//。。。
int s=20,t=0;//s为判断次数（精度），t为2的次方数
long long u=n-1;//u初始化为n/2^1
while(!(u&1)){//当u还为偶数
t++;
u>>=1;//+1次方并让u除以2
}
while(s--){//判断s次
long long a=rand()%(n-2)+2;//随机0+2~n-3+2区间取a
x[0]=qpow(a,u,n);//初始化为a的u次方模n
for(int i=1;i<=t;i++){//一次次乘上2
x[i]=multi(x[i-1],x[i-1],n);
if(x[i]==1 && x[i-1]!=1 && x[i-1]!=n-1)
return false;//二次判断区
}
if(x[t]!=1)//同上
return false;
}
return true;
}

long long gcd(long long a,long long b){
return b?gcd(b,a%b):a;//辗转相除法求最大公约数
}

long long pollard_rho(long long n,int c){
long long i=1,k=2,x=rand()%(n-1)+1,y=x;
//i:走了几步，x:当前位置，y:起点
while(true){
i++;
x=(multi(x,x,n)+c)%n;//迈一步,x^2+c%n;
long long p=gcd((y-x+n)%n,n);//取公约数
if(p!=1 && p!=n)
return p;
if(y==x)
return n;//迈不动了就跳出
if(i==k){
y=x;
k<<=1;
}//更新起点，让踩到的x值更随机。。。
}
}

void find(long long n,int c){
if(n==1)
return;
if(miller_rabin(n)){//若为素数
ans=min(ans,n);
return;
}
long long p=n,k=c;
while(p>=n)//查找因子
p=pollard_rho(p,c--);
find(p,k);//递归查找最小质因子
find(n/p,k);
}

int main(){
long long n,i;
int t;
scanf("%d",&t);
for(int v=1;v<=t;v++){
ans=1LL<<62;
scanf("%lld",&n);
find(n,120);
if(ans==n)
printf("Prime\n");
else printf("%lld\n",ans);
}
return 0;
}

以下为浪费时间：

由于该算法在推算x[i]
时，最后必定会出现x[i+j]=x[i]，然后把x[i]按照一种神奇的方式写下来，就会长得像希腊字母ρ
。。。

如下图：



然后发明者叫做Pollard，所以就叫做了这样一个神奇的名字。。。