poj 2191Mersenne Composite Numbers 大整数因式分解+素数判断

Mersenne Composite Numbers

Description

One of the world-wide cooperative computing tasks is the "Grand Internet Mersenne Prime Search" -- GIMPS -- striving to find ever-larger prime numbers by examining a particular category of such numbers. 
A Mersenne number is defined as a number of the form (2p–1), where p is a prime number -- a number divisible only by one and itself. (A number that can be divided by numbers other than itself and one are called "composite" numbers, and each of these can be uniquely represented by the prime numbers that can be multiplied together to generate the composite number — referred to as its prime factors.) 
Initially it looks as though the Mersenne numbers are all primes. 

If, however, we are having a "Grand Internet" search, that must not be the case. 
Where k is an input parameter, compute all the Mersenne composite numbers less than 2k -- where k <= 63 (that is, it will fit in a 64-bit signed integer on the computer). In Java, the "long" data type is a signed 64 bit integer. Under gcc and g++ (C and C++ in the programming contest environment), the "long long" data type is a signed 64 bit integer. 

Input

Input contains a single number, without leading or trailing blanks, giving the value of k. As promised, k <= 63.

Output

One line per Mersenne composite number giving first the prime factors (in increasing order) separate by asterisks, an equal sign, the Mersenne number itself, an equal sign, and then the explicit statement of the Mersenne number, as shown in the sample output. Use exactly this format. Note that all separating white space fields consist of one blank.

Sample Input

31

Sample Output

23 * 89 = 2047 = ( 2 ^ 11 ) - 1

47 * 178481 = 8388607 = ( 2 ^ 23 ) - 1

233 * 1103 * 2089 = 536870911 = ( 2 ^ 29 ) - 1

Source

​​Pacific Northwest 2004​​

算法分析：

题意：

给你一个n，从1~n，满足2^i-1为合数且i为素数，则满足

实现：

大整数因式分解Pollard_rho算法和大整数素数判断miller_rabin。但pow函数不知道为什么会编译错误，手写一个快速幂

代码实现：

#include<cstdio>  
#include<cstring>  
#include<cstdlib>  
#include<cctype>  
#include<cmath>  
#include<iostream>  
#include<sstream>  
#include<iterator>  
#include<algorithm>  
#include<string>  
#include<vector>  
#include<set>  
#include<map>  
#include<stack>  
#include<deque>  
#include<queue>  
#include<list> 
using namespace std;
 
 
//****************************************************************
// Miller_Rabin 算法进行素数测试
//速度快，而且可以判断 <2^63的数
//****************************************************************
const int S=20;//随机算法判定次数，S越大，判错概率越小
 
 
//计算 (a*b)%c.   a,b都是long long的数，直接相乘可能溢出的
//  a,b,c <2^63
long long mult_mod(long long a,long long b,long long c)       //快速乘计算   (a*b)%c
{
    a%=c;
    b%=c;
    long long ret=0;
    while(b)
    {
        if(b&1){ret+=a;ret%=c;}
        a<<=1;
        if(a>=c)a%=c;
        b>>=1;
    }
    return ret;
}
 
 
 
//计算  x^n %c  把快速幂分成这两部分写的好处，可以防止数越界
long long pow_mod(long long x,long long n,long long mod)//快速幂 计算  x^n%c
{
    if(n==1)return x%mod;
    x%=mod;
    long long tmp=x;
    long long ret=1;
    while(n)
    {
        if(n&1) ret=mult_mod(ret,tmp,mod);
        tmp=mult_mod(tmp,tmp,mod);
        n>>=1;
    }
    return ret;
}
 
 
 
 
 
//以a为基,n-1=x*2^t      a^(n-1)=1(mod n)  验证n是不是合数
//一定是合数返回true,不一定返回false
bool check(long long a,long long n,long long x,long long t)
{
    long long ret=pow_mod(a,x,n);
    long long last=ret;
    for(int i=1;i<=t;i++)
    {
        ret=mult_mod(ret,ret,n);
        if(ret==1&&last!=1&&last!=n-1) return true;//合数
        last=ret;
    }
    if(ret!=1) return true;
    return false;
}
 
// Miller_Rabin()算法素数判定
//是素数返回true.(可能是伪素数，但概率极小)
//合数返回false;
 
bool Miller_Rabin(long long n)
{
    if(n<2)return false;
    if(n==2)return true;
    if((n&1)==0) return false;//偶数
    long long x=n-1;
    long long t=0;
    while((x&1)==0){x>>=1;t++;}
    for(int i=0;i<S;i++)
    {
        long long a=rand()%(n-1)+1;//rand()需要stdlib.h头文件
        if(check(a,n,x,t))
            return false;//合数
    }
    return true;
}
 
 
//************************************************
//pollard_rho 算法进行质因数分解
//************************************************
long long factor[100];//质因数分解结果（刚返回时是无序的）
int tol;//质因数的个数。数组小标从0开始
 
long long gcd(long long a,long long b)
{
    if(a==0)return 1;//???????
    if(a<0) return gcd(-a,b);
    while(b)
    {
        long long t=a%b;
        a=b;
        b=t;
    }
    return a;
}
 
long long Pollard_rho(long long x,long long c)
{
    long long i=1,k=2;
    long long x0=rand()%x;
    long long y=x0;
    while(1)
    {
        i++;
        x0=(mult_mod(x0,x0,x)+c)%x;
        long long d=gcd(y-x0,x);
        if(d!=1&&d!=x) return d;
        if(y==x0) return x;
        if(i==k){y=x0;k+=k;}
    }
    
}
//对n进行素因子分解
void findfac(long long n)
{
    if(Miller_Rabin(n))//素数
    {
        factor[tol++]=n;
        return;
    }
    long long p=n;
    while(p>=n)p=Pollard_rho(p,rand()%(n-1)+1);
    findfac(p);
    findfac(n/p);
}
long long qpow(long long a,long long b)
{
    long long  res=1;
    while(b>0)
    {
        if(b%2) res=(res*a);
        a=(a*a);
        b/=2;
    }
    return res;
}

int main()
{
    //srand(time(NULL));//需要time.h头文件//POJ上G++不能加这句话
    long long n;
    while(scanf("%I64d",&n)!=EOF)
    {
        for(int i=2;i<=n;i++)
        {   
        long long m=qpow(2,i)-1;
    
        if(Miller_Rabin(m)||!Miller_Rabin(i))
            continue;
        
        tol=0;
        findfac(m);
        sort(factor,factor+tol);
        for(int i=0;i<tol-1;i++)
        printf("%I64d * ",factor[i]);
         printf("%I64d = %lld = ( 2 ^ %d ) - 1",factor[tol-1],m,i);
        printf("\n");
        }
    }
    return 0;
}