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不一定大于0
		if(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 1000000000
typedef unsigned long long LL;
#define N 1010
vector<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);
#endif
	srand(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;
}