collatz.c

calc大数库

/* collatz.c */
#ifdef _WIN32
#include "unistd_DOS.h"
#else
#include <unistd.h>
#endif
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#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;
}