#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>

bool isPrime(__int64 n);
bool miller_rabin(__int64 n,int s);
bool witness(__int64 a,__int64 n);
__int64 modular_exponentiation(__int64 a,__int64 u,__int64 n);
void pollard_rho(__int64 n);
void factorizer(__int64 n); // for n <= 1024*1024
int binary_length(__int64 x);
__int64 euclid(__int64 a,__int64 b);
__int64 euclid_stein(__int64 a,__int64 b);
__int64 rand64();


__int64 min;

int primes[173] = {2, 3, 4, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797, 809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887, 907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013, 1019, 1021};

int main()
{
	srand((unsigned)time(NULL));		
    int t,i;	// the number of test cases, 1 <= t <= 20	
	__int64 n;
    scanf("%d",&t);	
	for (i=0;i<t;i++) {		
		scanf("%I64d",&n);		
		if (isPrime(n)) 
			printf("%s\n","Prime");		
		else
		{			
			min = n;
			pollard_rho(n);			
			printf("%I64d\n",min);
		}			
	}  	
    return 0;
}
//////////////////////////////////////////////////////////
/*
int main()
{	
	__int64 i=1050028,j,p;
	while(1)
	{
		next:
		i++;
		printf("%I64d\n",i);
		min = i;
		if (isPrime(i)) 
			{printf("%s\n","Prime");continue;	}
		else 
			printf("%s\n","Composite");	

		pollard_rho(i);
		p = min;
		//is prime?
		
		for(j=2;j<=sqrt(i);j++)
		{
			if (i%j==0){

				if (j!=p) {printf("pollard WRONG!!!! %I64d j:%I64d p:%I64d\n ",i,j,p);return 1;}
				goto next;
			}
			
		}
		if (!isPrime(i)) {printf("isPrime WRONG!!!! %I64d\n",i);return 1;}		
		
		
	}
	return 0;
	
	
	//pollard_rho(4);
	//printf("%I64d\n",min);
}
*/
////////////////////////////////////////////////
// 2 < N < 2^54
bool 
isPrime(__int64 n) 
{
	int j;
	if (n==1) 
		return false;
	for (j=0;j<173;j++)
	{
		if (n==primes[j]) return true;
		if (n%primes[j]==0) return false; // n is composite.
	}	
	
	return miller_rabin(n,50);	// s取50应该差不多了
}

bool 
miller_rabin(__int64 n,int s)
{
	int j;
	__int64 a;
	for (j=1;j<=s;j++) {
		// Make sure that a is in [ 1, n - 1 ].
		a = rand64()%(n-1)+1;
		if (witness(a,n))
			return false;
	}
	return true;
}


bool 
witness(__int64 a,__int64 n)
{
	__int64 u=n-1,z=0,j;	
	while (u%2==0) 
	{
		u=u/2;
		z++;
	}
	__int64 x0 = modular_exponentiation(a,u,n),x1;
	for (j=1;j<=z;j++) 
	{
		x1 = x0*x0%n;
		if (x1==1&&x0!=1&&x0!=n-1)
			return true;
		x0=x1;
	}
	if (x1!=1) return true;
	return false;
}

// a^b mod n, 算法到论 Corollary 31.35
__int64 
modular_exponentiation(__int64 a,__int64 b,__int64 n) 
{	
	__int64 c=0, d=1;
	int i;	
	for (i=binary_length(b)-1;i>=0;i--)
	{
		c=c*2;
		d=(d*d)%n;
		if (b&(1<<i))
		{
			c = c + 1;
			d = ( d * a ) % n;
		}
	}
	return d;
}

// 返回x用二进制表示的长度
int 
binary_length(__int64 x)
{
	int result = 0;
	while (x>0)
	{
		x=x>>1;
	}
	return result;
}

__int64
rand64()
{
	__int64 result = 0;	
	int i;
	
	for (i=0;i<=62;i++)
	{
		if (rand()%2==0) result = result|1LL<<i;		
	}
	result = result|1;		// 保证为奇数
	result = result|1LL<<62;	// 保证足够的长	
	return result;
}


void 
pollard_rho(__int64 n)
{	
	if (n<=1048576/* 1024*1024 */) {factorizer(n);return;}
		
	__int64 i = 1;
	__int64 x1 = rand64() % n;	
	__int64 y = x1;
	__int64 d;
	__int64 k = 2;	

	while(1)
	{
		i++;
///		printf("n:%I64d x1:%I64d\n",n,x1);
		x1 = (x1 * x1 -1) % n;
		//d = euclid_stein(abs(y-x1), n);
		d = euclid(abs(y-x1), n);
		if (d!=1 && d!=n)
		{
			if (d < min) min = d;
			if (!isPrime(d)) 
			{
				pollard_rho(d);
			}
			if (n/d < min) min = n/d;			
			if (!isPrime(n/d))
			{
				pollard_rho(n/d);
			}		 
///			printf("d:%I64d\n",d);
			return;
		}
        if (i==k) {
         	y=x1;
            k*=2;
        }		
	}	
}




void 
factorizer(__int64 n)
{
	int j;
	for (j=0;j<173;j++)
	{
		if (n%primes[j]==0) {if (min>primes[j]) min=primes[j];return;}
	}
}


// 欧几里德，求最大公约数
__int64 
euclid(__int64 a,__int64 b)
{
	if (b==0) return a;
	else return euclid(b,a%b);
}


/*
Euclid & stein 算法
与传统的Euclid算法不同的是它只有移位和加减法.
*/
/*
__int64
euclid_stein(__int64 a,__int64 b)
{
	__int64 temp;
    if(a < b)   //arrange so that a>b
	{
		temp = a; 
        a = b; 
        b = temp; 
    } 
    if ( !b )    // the base case, b == 0 
        return a; 
    if (((a&1) == 0 ) && ( (b&1) == 0 ))	//a and b are even 
        return euclid_stein( a>>1, b>>1 )<<1; 
    else if ( (b&1) )	// only a is even 
        return euclid_stein( a>>1, b ); 
    else if ( (a&1) )	// only b is even 
        return euclid_stein( a, b>>1 ); 
	else
		return euclid_stein( b, (a-b)>>1 );	// a and b are odd 
}
*/

• 乘法溢出 (5939) frkstyc 2005-12-24 16:28:39
  • Re:乘法溢出 (66) liuyuyangfoam 2005-12-25 12:37:36
    • Re:乘法溢出 (147) liuyuyangfoam 2005-12-25 13:32:17
      • Re:乘法溢出 (36) liuyuyangfoam 2005-12-25 13:37:04