POJ 1811 *** Prime Test（详解Miiler_Rabin算法与Pollard_Rho算法）

题意:对于一个给定的数n，判断n是否为质数，n如果为质数，输出“Prime”，如果为合数，则输出其最小的素因子。

分析：其实这道题就是对于Miller_Rabin算法和Pollard_Rho算法的应用，具体的详解如下（刚买的椰子味身体乳香香的，感觉自己变成了一颗大椰子哈哈）

Miller_Rabin随机性素数测试方法：

先给出费马定理的定义：

如果p是素数，则a^(p-1)==1 mod p 对所有a∈Zp*都成立；Zp*为1到p-1中与p互质的数的集合。

算 法利用费马定理，试验了多个随机选取的基值a，来判断a^(p-1) ?= 1 mod p。

算法在实现的过程中会用到下面一些辅助过程，贴一下伪代码吧。//伪代码基本来自于算导

辅助过程1（求a^b mod n）

pow_mod(a,b,n)
1   d=1
2   while(b)
3      if b&1 
4         d=d*a mod n
5      a<<=1
5      b>>=1
6   return d

辅助 过程2（求a^(n-1) mod n）

辅助过程2与辅助过程1不同在于，过程2中间在寻找一个模n为1的非平凡平方根。贴了伪代码之后细说。

Witness(a,n)
1   let t and u be such that t>=1,u is odd,and n-1=u*2^t
2   Xo=pow_mod(a,u,n)
3   for i=1 to t
4      Xi=Xi-1^2 mod n
5      if Xi==1 and Xi-1 != 1 and Xi-1 != n-1
6         return true //一定是合数
7   if Xi != 1
8      return true //一定是合数
9   return false  //不一定是合数</span>

上面代码是将a^(n-1) mod n 转换为 (a^u)^(2^t) mod n（第一行），然后第二行求出 Xo=a^u mod n的值之后继续进行计算。但是在计算过程中，5、6行代码利用了下面的定理：

如果p是一个奇素数且e>=1，则方程
      x^2==1 mod p^e
仅有两个解，即 x==1和x==-1

上述定理的证明如下：

x^2==1 mod p^e 等价于 p^e|(x-1)(x+1)，但是他们不同时成立，否则p也能整除他们的差(x+1)-(x-1)=2。
如果gcd(p,x-1)==1,则 p^e|(x+1)，则x==-1 mod p^e。
如果gcd(p,x+1)==1,则p^e|(x-1)，则x=1 mod p^e。

但是如果一个数x，满足方程x^2==1 mod n，然而x却不等于1 或者 -1 ，则称x是一个以模n为1的非平凡平方根。

上述定理的逆否命题为：

如果存在模n为1的非平凡平方根，那么n不可能是奇素数或者奇素数的幂。即n为合数。

过程2中的5、6行代码，就是在判断X i 是否是模n为1的非平凡平方根，如果是，那么n为合数。

有了上面的基础之后再看Miller_Rabin的算法就比较轻松了，伪代码如下:

Miller_Rabin(n,s)
1   for i=1 to s
2      a=random(1,n-1)
3      if Witness(a,n)
4         return COMPOSITE //   definitely
5   return PRIME //   almost surely

Pollard_Rho算法：

这个算法主要是基于Birthday Trick来提高概率的。从头讲起。

对于一个很大的奇数n，如果想要知道他的因子，那么我们可以进行试除，从3一直试除到n-1，传统的试除法很暴力。

那么我们可以变得不一样，那就是从我们random(3,n-1)，这样的试除法只做一次的话，成功找到n的因子的概率为1/(n-3)，概率仍然很小，但是如果通过Birthday Trick，那么概率会大大提高。

Birthday Trick（生日悖论）

如果一个房间里有23个或23以上的人，那么至少有两个人的生日相同的概率大于50%。

在[1,10000]中选一个数，选中50的概率为 1/10000
在[1,10000]中选两个数i、j，i-j==50的概率约为1/5000
那么在[1,10000]中选k个数，这k个数两两相减得到一个数为50的概率是多少呢？

嗯。。。事实上这个概率是随着k的增加而增加的，总之最后会近乎为1。

于是我们可以这样子想，我们随机的选取k个数，然后判断Xi-Xj是否为n的因子，但是其实这样做并没有什么用。。

但是如果我们把要求降低，比如，判断gcd(Xi-Xj,n)是否等于1，这样的话如果n=p*q(p、q都为素数)那么如果Xi-Xj=2p(3p,4p...)都是可以判断的了。

所以我们可以在[2,n-1]范围内选取k个数，来进行判断，但是这样做需要的数据量可能很大，不方便进行存储。于是这个时候Circle Detection就上场了。

我们并不需要一次性生成k个数，我们可以利用一个函数来依次生成一个一个数，并且进行相减的检查。这个函数生成的数看上去就像随机的一样。。

这个函数是：

f(x)=(x^2+a) mod n

这个函数最后生成的数据形状大概类似罗马字ρ(rho)，对于这个函数，重要的一点探测环的出现。

探测环出现的算法大概是这个意思：如果A在一个圆圈上面行走，那么A怎么知道自己已经走了一圈了呢？那么这个时候需要B了， 如果B的速度是A的两倍，那么当B第一次追赶上A的时候，A一定已经走完一圈了。

这大概就是Pollar_Rho我能理解到的程度了，以后如果有更深入的理解那么继续补充。

Pollar_Rho的伪代码如下:

Pollar_Rho(n)
1   i=1
2   X1=random(0,n-1)
3   Y=X1
4   k=2
5   while(true)
7      i=i+1
8      Xi=(Xi-1^2-1) mod n
6      d=gcd(Y-X,n)
7      if d != 1 and d != n
8         return d
9      if i==k
10       Y=Xi
11       k+=k

poj1811

题目的代码下也有一些解释，模板基本来自于kuangbin。

代码如下：

#pragma warning(disable:4996)
#include<iostream>
#include<cstdio>
#include<cmath>
#include<stack>
#include<queue>
#include<cstring>
#include<sstream>
#include<set>
#include<string>
#include<iterator>
#include<vector>
#include<map>
#include<algorithm>
using namespace std;
typedef long long ll;



//计算 (a*b)%c.   a,b都是long long的数，直接相乘可能溢出的
//这里相乘是利用b=2^n1+2^n2+...+2^n,然后代入(a*b)%c中展开相乘
ll multi_mod(ll a, ll b, ll n) {
	a %= n;
	b %= n;
	ll res = 0;
	while (b) {
		if (b & 1) {
			res += a;
			res %= n;
		}
		a <<= 1;
		if (a >= n)a %= n;
		b >>= 1;
	}
	return res;
}


//计算  n^a % mod
ll pow_mod(ll n, ll a, ll mod)
{
	if (a == 1)return n%mod;
	n %= mod;
	ll res = 1;
	while (a) {
		if (a & 1)res=multi_mod(res, n, mod);
		n = multi_mod(n, n, mod);
		a >>= 1;
	}
	return res;
}

//以a为基,n-1=x*2^t      a^(n-1)==1(mod n)  验证n是不是合数
//一定是合数返回true,不一定返回false
bool check(ll a, ll n,ll x, ll t)
{
	ll res = pow_mod(a, x, n);
	ll last = res;
	for (int i = 1; i <= t; i++)
	{
		res = multi_mod(res, res, n);
		if (res == 1 && last != 1 && last != n - 1)
			return 1;
		last = res;
	}
	if (res != 1)return 1;
	return false;
}

// Miller_Rabin()算法素数判定
//是素数返回true.(可能是伪素数，但概率极小)
//合数返回false;

bool Miller_Rabin(long long n)
{
	if (n<2)return false;
	if (n == 2)return true;
	if ((n & 1) == 0) return false;//偶数
	ll x = n - 1, t=0;
	while ((x&1)==0) { t++, x /= 2; }
	for (int i = 0; i < 20; ++i) {
		ll a = rand() % (n - 1) + 1;
		if (check(a, n, x, t))
			return 0;//合数
	}
	return 1;
	return true;
}



//最大公约数的判断
ll gcd(ll a,ll b){
	if (a < 0)return gcd(-a, b);
	if (a == 0)return 1;
	while (b) {
		ll t = a%b;
		a = b;
		b = t;
	}
	return a;
}

//寻找其中一个因子
ll pollard_rho(ll n,ll c) {
	ll i = 1, k = 2, x0, y;
	x0 = rand() % (n - 1) + 1;
	y = x0;
	while (1) {
		i++;
		x0 = (multi_mod(x0, x0, n) + c) % n;
		ll d = gcd(y - x0, n);
		if (d != 1 && d != n)return d;
		if (x0 == y)return n;
		if (i == k) {
			y = x0;
			k += k;
		}
	}
}

//寻找素因子
ll factor[100];
int tot;
void findp(ll n) {
	if (Miller_Rabin(n)) {
		factor[tot++] = n;
		return;
	}
	ll p=n;
	while (p >= n)p = pollard_rho(n, rand() % (n - 1) + 1);
	findp(p);
	findp(n / p);
}


int main(void) {
	int t;
	ll n;
	cin >> t;
	while (t--) {
		cin >> n;
		if (Miller_Rabin(n)) {
			cout << "Prime" << endl;
			continue;
		}
		tot = 0;
		findp(n);
		ll ans = factor[0];
		for (int i = 1; i < tot; ++i)
			if (factor[i]<ans)
				ans = factor[i];
		cout << ans << endl;
	}
	return 0;
}