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 200000using 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;}