poj 2191 Mersenne Composite Numbers 大数分解

题意：

给k，求i是素数且在1~k内并且2^i-1是合数的情况，并将它分解。

分析：

枚举1至i然后用miller_rabin素数判定和pollard_rho因数分解即可。

代码：

//poj 2191
//sep9
#include <iostream>
#include <map>
#include <vector>
#define gcc 10007
#define max_prime 200000
using namespace std;
typedef long long ll;
vector<int> primes;
vector<bool> is_prime;
ll gcd(ll a,ll b)
{
	ll m=1;
	if(!b)	return a;
	while(m){
		m=a%b;
		a=b;
		b=m;
	}
	return a;
}

ll multi_mod(ll a,ll b,ll mod)
{
	ll sum=0;
	while(b){
		if(b&1)	sum=(sum+a)%mod;
		a=(a+a)%mod;
		b>>=1;
	}	
	return sum;
}

ll pow_mod(ll a,ll b,ll mod)
{
	ll sum=1;
	while(b)
	{
		if(b&1) sum=multi_mod(sum,a,mod);
		a=multi_mod(a,a,mod);
		b>>=1;
	}	
	return sum;
} 

bool miller_rabin(ll n,int times)
{
	if(n<2) return false;
	if(n==2) return true;
	if(!(n&1)) return false;
	ll q=n-1;
	int k=0;
	while(!(q&1)){
		++k;
		q>>=1;
	}
	for(int i=0;i<times;++i){
		ll a=rand()%(n-1)+1;
		ll x=pow_mod(a,q,n);
		if(x==1) continue;
		for(int j=0;j<k;++j){
			ll buf=multi_mod(x,x,n);
			if(buf==1&&x!=1&&x!=n-1)
				return false;
			x=buf;
		}
		if(x!=1)
			return false;
	}		
	return true;
}

ll pollard_rho(ll n)
{
	ll d,i,k,x,y;
	x=rand()%(n-1)+1;
	y=x;
	i=1;
	k=2;
	do
	{
		++i;
		d=gcd(n+y-x,n);
		if(d>1&&d<n)
			return d;
		if(i==k)
			y=x,k*=2;
		x=((multi_mod(x,x,n)-gcc)%n+n)%n;
	}while(y!=x);
	return n;
}


bool isPrime(ll n)
{
	if(n<=max_prime) return is_prime[n];
	return miller_rabin(n,20);
}

void factorize(ll n,map<ll,int>& factors)
{
	if(isPrime(n))
		++factors[n];
	else{
		for(int i=0;i<primes.size();++i){
			int p=primes[i];
			while(n%p==0){
				++factors[p];
				n/=p;	
			}	
		}
		if(n!=1){
			if(isPrime(n))
				++factors[n];
			else{
				ll d=pollard_rho(n);
				factorize(d,factors);
				factorize(n/d,factors);
			}
		}
	}
}

void ini_primes()
{
	is_prime=vector<bool>(max_prime+1,true);
	is_prime[0]=is_prime[1]=false;
	for(int i=2;i<=max_prime;++i){
		if(is_prime[i]){
			primes.push_back(i);
			for(int j=i*2;j<=max_prime;j+=i)
				is_prime[j]=false;
		}
	}			
}

int main()
{
	ini_primes();
	int i,k;
	scanf("%d",&k);
	ll p=1;
	for(i=1;i<=k;++i){
		p=2*p;	
		if(is_prime[i])
			if(isPrime(p-1)==0){
				map<ll,int> factor;
				factorize(p-1,factor);
				int first=0;
				for(map<ll,int>::iterator it=factor.begin();it!=factor.end();++it)	
					while(it->second){
						if(first==0){
							printf("%I64d ",it->first);
							first=1;
						}
						else
							printf("* %I64d ",it->first);
						--it->second; 
					}
				printf("= %I64d = ( 2 ^ %d ) - 1\n",p-1,i);	
			}		
	}	
	return 0;	
}