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 16381int64 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;//保证数据大于等于2if(x%n==0)continue;x = Pow_mod(x,u,n);//计算x^upre = 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;}