www.pudn.com > calc.rar > collatz.c

/* collatz.c */ 
#ifdef _WIN32 
#include "unistd_DOS.h" 
#else 
#include  
#endif 
#include  
#include  
#include  
#include "integer.h" 
#include "fun.h" 
#define UPPER_BOUND 4294967294UL /* 2^32 - 2 */ 
#define CYCLE_MAX 150 
#define BRANCH_MAX 20 
 
MPI *GEN_COLLATZ(MPI *Aptr, USL d, MPI *m[], MPI *X[]) 
/* 
 * The generalized Collatz function. 
 * Here d is the number of branches. 
 * m[0],...,m[d-1] are non-zero integers (the multipliers). 
 * X[0],...,X[d-1] are the shifts. 
 * T(Aptr)=INT(m[r]*Aptr/d)+X[r], if Aptr(MOD(d))=r. 
 */ 
 
{ 
	USL r; 
	MPI *T1, *T2, *T3; 
 
	r = MOD_(Aptr, d); 
	T1 = MULTI(Aptr, m[r]); 
	T2 = INT_(T1, d); 
	FREEMPI(T1); 
	T3 = ADDI(T2, X[r]); 
	FREEMPI(T2); 
	return (T3); 
} 
 
void P_CYCLE(MPI *Aptr, USL d, MPI *m[], MPI *X[], USI c) 
 
/* the cycle starting with *Aptr is printed */ 
{ 
	MPI *B, *Temp, *MIN_ELT, *TEMP1; 
	unsigned int i = 0; 
	FILE *outfile; 
	char buff[20]; 
 
	strcpy(buff, "cycle.out"); 
	outfile = fopen(buff, "a"); 
	B = COPYI(Aptr); 
	MIN_ELT = COPYI(Aptr); 
	do{ 
		Temp = B; 
		B = GEN_COLLATZ(B, d, m, X); 
		FREEMPI(Temp); 
		if(RSV(B, MIN_ELT) < 0){ 
			TEMP1 = MIN_ELT; 
			MIN_ELT = COPYI(B); 
			FREEMPI(TEMP1); 
		} 
        }while (!EQUALI(Aptr, B)); 
 
	Temp = B; 
	B = COPYI(MIN_ELT); 
        FREEMPI(Temp); 
	do { 
		PRINTI(B); 
		FPRINTI(outfile, B); 
		printf("\n"); 
		fprintf(outfile, "\n"); 
		Temp = B; 
		B = GEN_COLLATZ(B, d, m, X); 
		FREEMPI(Temp); 
		i++; 
	} 
	while (!EQUALI(MIN_ELT, B)); 
	printf("length of cycle %u is %u\n", c, i); 
	fprintf(outfile, "length of cycle %u is %u\n", c, i); 
	fclose(outfile); 
	FREEMPI(B); 
	FREEMPI(MIN_ELT); 
	return; 
} 
 
unsigned int IS_CYCLE(MPI *Aptr, USL d, MPI *m[], MPI *X[], MPI *Z[], USI c) 
/* the cycle starting from Aptr is distinct from those starting from  
 * Z[0],..., Z[c]  if and only if IS_CYCLE() = 1. 
 */ 
{ 
	USL k; 
	MPI *B, *Temp; 
	unsigned int i, t; 
 
	B = COPYI(Aptr); 
	for (k = 0; 1; k++) 
        { 
		for (i = 0; i < c; i++) 
	        { 
			if (EQUALI(B, Z[i]))  
			{ 
				FREEMPI(B); 
				return (0); 
			} 
		} 
		Temp = B; 
		B = GEN_COLLATZ(B, d, m, X); 
		FREEMPI(Temp); 
		t = EQUALI(Aptr, B); 
		if (t) 			/* found a new cycle */ 
		{ 
			FREEMPI(B); 
			break; 
		} 
	} 
	return (1); 
} 
 
void CYCLE(USL d, MPI *m[], MPI *X[], USL INFINITY, USL RANGE)  
/* 
 * This function searches all trajectories of the generalized Collatz 
 * function which start from p, |p| <= RANGE/2 (RANGE an even integer). 
 */ 
{	 
	MPI *A, *B, **Z, *Temp, *B1; 
	unsigned int c = 0, n, r, end, *FLAG, t = 0, i, s = 0; 
	USI diverge_flag=0; 
	long int p, q, LIMIT, k; 
	FILE *outfile; 
	char buff[20]; 
 
	strcpy(buff, "cycle.out"); 
	LIMIT = CYCLE_MAX; 
	FLAG = (unsigned int *)ccalloc(RANGE + 1, sizeof(unsigned int)); 
	Z = (MPI **)mmalloc(LIMIT * sizeof(MPI *)); 
	INTSETI((int *)FLAG, RANGE + 1, RANGE + 1); 
	for (n = 0; n <= RANGE; n++) { 
		if (FLAG[n] > n) { 
			if ((n & 1) == 0) 
				p = (n >> 1); 
			else 
				p = -((n + 1) >> 1); 
			q = p; 
			printf("p = %ld\n", p); 
			A = CHANGEI(p); 
			B =  COPYI(A); 
			end = 0; 
			for (k = 0; k <= UPPER_BOUND; k++)  
			{				 
				Temp = A; 
				A = GEN_COLLATZ(A, d, m, X); 
				FREEMPI(Temp); 
				if (A->D >= INFINITY) { 
				if(diverge_flag == 0){  
					outfile = fopen(buff, "a"); 
		 
				        fprintf(outfile, "divergent trajectory starting from q = %ld\n", q); 
				        printf("divergent trajectory starting from q = %ld\n", q); 
					fclose(outfile); 
				  diverge_flag = 1; 
				} 
				  FREEMPI(A); 
				  FREEMPI(B); 
				  break; 
				} 
				Temp = B; 
				B1 = GEN_COLLATZ(B, d, m, X); 
				B = GEN_COLLATZ(B1, d, m, X); 
				FREEMPI(B1); 
				FREEMPI(Temp); 
				t = EQUALI(A, B); 
				if (t)  
				{ 
					 /* cycle detected starting at p */ 
					FREEMPI(B); 
					break; 
				} 
				if (A->D == 0 && A->V[0] <= RANGE/2)  
				{ 
					p = A->V[0]; 
					r = (A->S >= 0 ? 2 * p : 2 * p - 1); 
					if (FLAG[r] < n) { 
						end = 1; 
						FREEMPI(A); 
						FREEMPI(B); 
						break; 
					} 
					else 
						FLAG[r] = n; 
				} 
			} 
			if (end == 1) 
				continue; 
			if (t) { 
				if (c == 0) { 
					/* store A, the start of first cycle */ 
					Z[0] =  A; 
					printf("\ncycle %u:\n", c + 1); 
					outfile = fopen(buff, "a"); 
					fprintf(outfile, "\ncycle %u:\n", c + 1); 
					fclose(outfile); 
					P_CYCLE(A, d, m, X, c + 1); 
					c++; 
				} 
				else { 
					/* store start of subsequent cycles */  
					s = IS_CYCLE(A, d, m, X, Z, c); 
					if (s == 0) 
						FREEMPI(A); 
					else 
					{ 
						if (c >= LIMIT) 
						{ 
							LIMIT++; 
					/*		ptr = Z;*/ 
							Z = (MPI **)rrealloc((char *)Z, (c + 1) * sizeof(MPI *), sizeof(MPI *)); 
							/*if (ptr != Z) 
								free((char *) ptr);*/ 
						} 
						Z[c] = A; 
						printf("\ncycle %u:\n", c + 1); 
						outfile = fopen(buff, "a"); 
						fprintf(outfile, "\ncycle %u:\n", c + 1); 
						fclose(outfile); 
						P_CYCLE(A, d, m, X, c + 1); 
						c++; 
					}	 
				} 
			} 
		} 
	}								 
	for (i = 0; i < c; i++) 
		FREEMPI(Z[i]); 
	ffree((char *)Z, LIMIT * sizeof(MPI *)); 
	ffree((char *)FLAG, (RANGE + 1) * sizeof(USI)); 
	return; 
}						 
 
void CYCLEX() 
{ 
	unsigned int i, RANGE, u; 
	USL d, INFINITY; 
	MPI **m, **X; 
	FILE *outfile; 
	char buff[20]; 
 
	strcpy(buff, "cycle.out"); 
	if(access(buff, R_OK) == 0) 
	  unlink(buff); 
	outfile = fopen(buff, "w"); 
	printf("enter 'INFINITY', (the cut-off number of allowable digits (base %u) for an iterate - say <20):", R0); 
	scanf("%lu", &INFINITY); 
	printf("INFINITY = %lu\n", INFINITY); 
	fprintf(outfile, "INFINITY = %lu\n", INFINITY); 
	printf("enter an even integer 'RANGE' (<= (say)500,000): (searches all trajectories of the generalized Collatz function which start from p, |p| <= RANGE)/2\n"); 
	scanf("%u", &RANGE); 
	printf("RANGE = %u\n", RANGE); 
	fprintf(outfile, "RANGE = %u\n", RANGE); 
	printf("enter the number d of branches (<= %u)\n", BRANCH_MAX); 
	scanf("%lu", &d); 
        printf("the number of branches is d = %lu\n", d); 
        fprintf(outfile, "the number of branches is d = %lu\n", d); 
 
	m = (MPI **)mmalloc(d * sizeof(MPI *)); 
	X = (MPI **)mmalloc(d * sizeof(MPI *)); 
	for (i = 0; i < d; i++) { 
		printf("enter the nonzero multiplier m[%u]\n", i); 
		m[i] = INPUTI(&u);; 
		printf("m[%u] = ", i);  
		fprintf(outfile, "m[%u] = ", i);  
		PRINTI(m[i]);  
		FPRINTI(outfile, m[i]);  
		printf("\n"); 
		fprintf(outfile, "\n"); 
	} 
 
	for (i = 0; i < d; i++) { 
		printf("enter the shift X[%u]\n", i); 
		X[i] = INPUTI(&u);; 
		printf("X[%u] = ", i);  
		fprintf(outfile, "X[%u] = ", i);  
		PRINTI(X[i]);  
		FPRINTI(outfile, X[i]);  
		printf("\n"); 
		fprintf(outfile, "\n"); 
	} 
	fclose(outfile); 
	Flush(); 
printf("\n\n"); 
	CYCLE(d, m, X, INFINITY, RANGE); 
	for (i = 0; i < d; i++) 
		FREEMPI(m[i]); 
	ffree((char *)m, d * sizeof(MPI *)); 
	for (i = 0; i < d; i++) 
		FREEMPI(X[i]); 
	ffree((char *)X, d * sizeof(MPI *)); 
	return; 
}