当前位置: 代码迷 >> 综合 >> poj 2191 Mersenne Composite Numbers 大数分解
  详细解决方案

poj 2191 Mersenne Composite Numbers 大数分解

热度:56   发布时间:2024-01-19 05:45:51.0

题意:

给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;}elseprintf("* %I64d ",it->first);--it->second; }printf("= %I64d = ( 2 ^ %d ) - 1\n",p-1,i);	}		}	return 0;	
}