POJ2191

POJ2191

#include <iostream>#include <cstdlib>#include <cstring>#include <algorithm>#include <cstdio>#include<cmath>const int Times = 10;const int N = 5500;using namespace std;typedef long long LL;LL ct, cnt;LL fac[N], num[N];LL gcd(LL a, LL b){    return b? gcd(b, a % b) : a;}LL multi(LL a, LL b, LL m){    LL ans = 0;    a %= m;    while(b)    {        if(b & 1)        {            ans = (ans + a) % m;            b--;        }        b >>= 1;        a = (a + a) % m;    }    return ans;}LL quick_mod(LL a, LL b, LL m){    LL ans = 1;    a %= m;    while(b)    {        if(b & 1)        {            ans = multi(ans, a, m);            b--;        }        b >>= 1;        a = multi(a, a, m);    }    return ans;}bool Miller_Rabin(LL n){    if(n == 2) return true;    if(n < 2 || !(n & 1)) return false;    LL m = n - 1;    int k = 0;    while((m & 1) == 0)    {        k++;        m >>= 1;    }    for(int i=0; i<Times; i++)    {        LL a = rand() % (n - 1) + 1;        LL x = quick_mod(a, m, n);        LL y = 0;        for(int j=0; j<k; j++)        {            y = multi(x, x, n);            if(y == 1 && x != 1 && x != n - 1) return false;            x = y;        }        if(y != 1) return false;    }    return true;}LL pollard_rho(LL n, LL c){    LL i = 1, k = 2;    LL x = rand() % (n - 1) + 1;    LL y = x;    while(true)    {        i++;        x = (multi(x, x, n) + c) % n;        LL d = gcd((y - x + n) % n, n);        if(1 < d && d < n) return d;        if(y == x) return n;        if(i == k)        {            y = x;            k <<= 1;        }    }}void find(LL n, int c){    if(n == 1) return;    if(Miller_Rabin(n))    {        fac[ct++] = n;        return ;    }    LL p = n;    LL k = c;    while(p >= n) p = pollard_rho(p, c--);    find(p, k);    find(n / p, k);}bool witness(LL a,LL n){    LL t,d,x;    d=1;    int i=ceil(log(n-1.0)/log(2.0)) - 1;    for(; i>=0; i--)    {        x=d;        d=(d*d)%n;        if(d==1 && x!=1 && x!=n-1) return true;        if( ((n-1) & (1<<i)) > 0)            d=(d*a)%n;    }    return d==1? false : true;}bool miller_rabin(LL n){    if(n==2)    return true;    if(n==1 || ((n&1)==0))    return false;    for(int i=0; i<50; i++)    {        LL a=rand()*(n-2)/RAND_MAX +1;        if(witness(a, n))    return false;    }    return true;}bool isprime(int x){    for(int i=2; i<=sqrt(x)+1; i++)    {        if(x%i==0)        {            return false;        }    }    return true;}int main(){    int k;    while(scanf("%d",&k)!=EOF)    {        for(int i=2; i<=k; i++)        {            if(i>=59&&i<=63)            {                printf("179951 * 3203431780337 = 576460752303423487 = ( 2 ^ 59 ) - 1\n");                break;            }            else            {                if(isprime(i))                {                    LL t=pow(2,i)-1;                    //printf("%lld \n",t);                    if(!miller_rabin(t))                    {                        ct = 0;                        find(t, 120);                        sort(fac, fac + ct);                        num[0] = 1;                        int k = 1;                        for(int i=1; i<ct; i++)                        {                            if(fac[i] == fac[i-1])                                ++num[k-1];                            else                            {                                num[k] = 1;                                fac[k++] = fac[i];                            }                        }                        cnt = k;                        for(int i=0; i<cnt-1; i++)                        {                            printf("%lld * ",fac[i]);                        }                        printf("%lld = %lld = ( 2 ^ %d ) - 1\n",fac[cnt-1],t,i);                    }                }            }        }    }    return 0;}