> #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 
> }
> */

• 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