miller_rabin素数判断和pollard_rho的素数因子分解算法

这类属于模板题吧。  感觉严格的数学证明 很是 NB，这两个都是 概率算法，不过准确率还是很高的。

下面是这两个模板。。留着以后用。。

Miller_rabin素数判断算法，一般的算法的复杂度是O(sqrt(n)) ，对于大数无能为力。

应用了费马大定理 和 二次探测。  

LL Mulmod(LL a,LL b,LL mod)//(a*b)%mod,防止溢出，把b看成二进制{LL ret=0;while(b){if(b&1){ret+=a; ret%=mod;}b>>=1;a+=a; a%=mod;}return ret;}LL Powmod(LL a,LL b,LL mod)//(a^b)%mod, 快速幂，同样把b看成二进制{LL ret=1;while(b){if(b&1) ret=Mulmod(ret,a,mod);b>>=1;a=Mulmod(a,a,mod);}return ret;}bool Witness(LL a,LL n)//a is a rand  number , n is test number;{LL t=0,u=n-1;while(!(u&1)){t++; u>>=1;}LL x0=Powmod(a,u,n);for(LL i=1;i<=t;i++){LL x1=Mulmod(x0,x0,n);if(x1==1 && x0!=1 && x0!=(n-1) ) return true;//二次探测的检测x0=x1;}if(x0!=1) return true;//费马大定理的检测return false;}bool Miller_rabin(LL n,int t) // n is the test number , t is the test times;{if(n<2) return false;if(n==2) return true;if( !(n&1) ) return false;srand(time(NULL));for(int i=0;i<t;i++)//检测t次{LL a=rand()%(n-1)+1;if(Witness(a,n)) return false ;}return true;}

Pollard_rho整数分解为素数因子 算法，一般的算法是 试除法，效率低。

//要用到Miller_rabin();LL gcd(LL a,LL b)// 求gcd(a,b);{if(!b) return a;return gcd(b,a%b);}LL pollard(LL n, int c) // n is the test number , c is a constant{LL x,y,d,i=1,k=2;x=rand()%(n-1)+1;y=x;while(true){i++;x=(Mulmod(x,x,n)+c)%n;d=gcd(abs(y-x),n);//随机出y x ,使得gcd(y-x,n) 在区间(1,n)间为 因数,y-x不一定大于0if(d>1 && d<n) return d;if(x==y) return n;if(i==k) y=x,k<<=1;}}void findfac(LL n) {if(n==1) return ;if(Miller_rabin(n,4)){save.push_back(n); return;//存入vector<LL > save中去得到的是无序的因子}LL p=n; while(p>=n) p=pollard(p,rand()%(n-1)+1);findfac(p);//递归找出所有因数findfac(n/p);}

POJ 2191：  http://poj.org/problem?id=2191

就是这两个算法的应用，题意是 梅森素数是 2^p-1的形式且这个数是素数，如果满足这个形式但是不是素数，称为梅森合数。给你k(k<=63)， 让你输出p<=k的所有梅森合数，并且 把这些梅森合数分解成质数因数的形式。

当然可以全部算出之后打表。。。。 还是应用一下这两个算法吧。。

CODE：

#include<cstdio>#include<cstring>#include<algorithm>#include<iostream>#include<string>#include<vector>#include<cmath>#include<ctime>using namespace std;#define INF 1000000000typedef unsigned long long LL;#define N 1010vector<LL > save;LL Mulmod(LL a,LL b,LL mod){LL ret=0;while(b){if(b&1){ret+=a; ret%=mod;}b>>=1;a+=a; a%=mod;}return ret;}LL Powmod(LL a,LL b,LL mod){LL ret=1;while(b){if(b&1) ret=Mulmod(ret,a,mod);b>>=1;a=Mulmod(a,a,mod);}return ret;}bool Witness(LL a,LL n)//a is a rand test number , n is test number;{LL t=0,u=n-1;while(!(u&1)){t++; u>>=1;}LL x0=Powmod(a,u,n);for(LL i=1;i<=t;i++){LL x1=Mulmod(x0,x0,n);if(x1==1 && x0!=1 && x0!=(n-1) ) return true;x0=x1;}if(x0!=1) return true;return false;}bool Miller_rabin(LL n,int t) // n is the test number , t is the test times 3?;{if(n<2) return false;if(n==2) return true;if( !(n&1) ) return false;for(int i=0;i<t;i++){LL a=rand()%(n-1)+1;if(Witness(a,n)) return false ;}return true;}LL gcd(LL a,LL b){if(!b) return a;return gcd(b,a%b);}LL pollard(LL n, int c) // n is the test number , c is a constant{LL x,y,d,i=1,k=2;x=rand()%(n-1)+1;y=x;while(true){i++;x=(Mulmod(x,x,n)+c)%n;d=gcd(y-x,n);if(d>1 && d<n) return d;if(x==y) return n;if(i==k) y=x,k<<=1;}}void findfac(LL n) // k is a constant;{if(n==1) return ;if(Miller_rabin(n,4)){save.push_back(n); return;}LL p=n;while(p>=n) p=pollard(p,rand()%(n-1)+1);findfac(p);findfac(n/p);}vector<LL > rr,num[100];int main(){#ifndef ONLINE_JUDGE    freopen ("in.txt" , "r" , stdin);#endifsrand(time(NULL));//POJ c++ 编译LL n,b;rr.clear();for(int i=0;i<100;i++) num[i].clear();for(int i=2;i<=63;i++)    {        if(!Miller_rabin((LL)i,4)) continue;        LL tmp=((LL)1<<i)-1;        if(!Miller_rabin(tmp,4)) rr.push_back((LL)i);    }    for(int i=0;i<rr.size();i++ )    {        save.clear();        LL now=rr[i];        LL tmp=((LL)1<<now)-1;        findfac(tmp);        num[now]=save;        sort(num[now].begin(),num[now].end());    }while(scanf("%I64d",&n)!=EOF){    for(int i=0;i<rr.size();i++)        {            LL now=rr[i];            if(now>n) break;            printf("%I64d",num[now][0]);            for(int j=1;j<num[now].size();j++)                printf(" * %I64d",num[now][j]);            printf(" = %I64d = ( 2 ^ %I64d ) - 1\n",((LL)1<<now)-1,now);        }}    return 0;}