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