POJ 1811 Prime Test

链接： http://poj.org/problem?id=1811

判断一个数是否为素数，如果不是的话，输出最小的素因子  范围 n <=2^54

从大小上来看，普通方法一定搞不定!故只能参考大牛blog 及代码 搞的

采用Miller-Rabin随机素数测试法 + Pollard_rho启发式分解因子法（算法导论有介绍，具体东东木看很懂也）

其中随机素数测试法基于两个定理和一个推论：

1.如果n是素数那么对于a<n，有a^(n - 1) mod n = 1. 该定理被称为费马小定理。

2.如果n是奇素数，那么x^2 mod n =1.仅有两个根，x = 1,和x = -1.

3.如果模n存在1的非平凡平方根，那么n是合数。

根据上述条件就可以判断一个数是否为素数了！

 Pollard_rho 启发式分解因子法表示木看懂，看导论去吧自己！为什么每次 c--  呢，求解释！

#include<stdio.h>
#include<stdlib.h>
#include<time.h>
typedef __int64 int64 ;
#define INF 99999999999999
#define Rand_time 5					//随机次数与正确率有关
#define C 16381
int64 MIN;
int64 Gcd(int64 a,int64 b)			//最大公约数，不解释
{
	return a%b?Gcd(b,a%b):b;
}
int64 Multi_mod(int64 a,int64 b,int64 n)	//(a*b)%n 直接计算有超出int64表示范围，故采取加法变形来计算值
{											
	int64 s = 0;
	while(b){
		if(b&1)s = (s+a)%n;
		a = (a+a)%n;
		b >>= 1;
	}
	return s;
}
int64 Pow_mod(int64 a,int64 b,int64 n)		//标准的快速求幂，不解释
{
	int64 s=1;
	while(b){
		if(b&1) s = Multi_mod(s,a,n);
		a = Multi_mod(a,a,n);
		b >>= 1;
	}
	return s;
}

int Prime_test(int64 n)
{
	int64 u = n-1,pre,x;
	int i,j,k=0;
	if( n == 2 || n == 3 || n==5 || n==7 || n==11)return 1;
	if( n == 1 || (!(n%2)) ||(!(n%3)) || (!(n%5)) || (!(n%7)) || (!(n%11)))return 0;
	for(;!(u&1);k++,u>>=1);
	srand((int64)time(0));		//随机参数
	for(i=0;i<Rand_time;i++){
		x = rand()%(n-2) + 2;	//保证数据大于等于2
		if(x%n==0)continue;
		x = Pow_mod(x,u,n);		//计算x^u
		pre = x;				//保留x,既x^2 mod n =1 的根
		for(j=0;j<k;j++){		//如果n是素数那么x^2mod n = 1,仅有两个根（不同余），+1和-1（n-1）
			x = Multi_mod(x,x,n);
			if(x==1 && pre != 1 &&pre != (n-1))
				return 0;
			pre = x;
		}
		if( x != 1)return 0;	//费马小定理 a^(n-1) mod n 恒等于 1 则可以认为 n 是素数
	}
	return 1;
}
int64 rho(int64 n ,int c)
{
	int64 x,y,d;
	int i = 1,j =2;
	srand((int64)time(0));
	x = rand()%(n-1)+1;
	y = x;
	while(1){
		i ++;
		x = (Multi_mod(x,x,n)+c)%n;
		d = Gcd(y-x+n,n);
		if(d!=1 && d!=n)return d;
		if(x==y)return n;
		if(i==j){
			y = x;
			j<<=1;
		}
	}
}
void find_factor(int64 n ,int c)
{
	int64 r,d;
	if( n <=1)return ;
	if(Prime_test(n)){
		if(MIN > n) MIN = n;
		return ;
	}
	r = rho(n,c--);			//不明白为什么 c-- ，求包养
	d = n/r;
	find_factor(d,c);
	find_factor(r,c);
}
int main()
{
	int T;
	int64 n;
	scanf("%d",&T);
	while(T--){
		scanf("%I64d",&n);
		if(Prime_test(n))puts("Prime");
		else {
			MIN  = INF;
			find_factor(n,C);
			printf("%I64d\n",MIN);
		}
	}
	return 0;
}