﻿ /*******************************
/******************************************************************************/
/*									      */
/*  COMPUTE 3N+1 OR 3N-1 SEQUENCE					      */
/*  02/13/10 (dkc)							      */
/*									      */
/*  This C program computes general 3n+1 or 3n-1 sequences.  The third-to-last*/
/*  element of a sub-sequence must be odd and the second-to-last element of   */
/*  a sub-sequence must be divisible by 8.  The third-to-last element must be */
/*  larger than the order divided by 6.  Certain lengths of such sequences    */
/*  are not permissible.						      */
/*									      */
/******************************************************************************/
#include <stdio.h>
#include <math.h>
void div64_32(unsigned int *a, unsigned int *b, unsigned int);
void add64(unsigned int *a, unsigned int *b);
void sub64(unsigned int *c, unsigned int *d);
void shift(unsigned int *a, unsigned int *b, unsigned int shift);
int main () {
//
int c=1;			      // c
//
unsigned int h,i,j,n,np,order,count;   // indices
unsigned int he[2048];			    // histogram
unsigned int y[178]={0,3,6,8,11,	  // I	(invalid lengths)
13,16,19,21,24,	  // I
26,29,32,34,37,	  // I
39,42,44,47,50,	  // J
52,55,57,60,63,	  // J
65,68,70,73,75,	  // K
78,81,83,86,88,	  // K
91,94,96,99,101,	  // K
104,106,109,112,114,  // L
117,119,122,125,127,  // L
130,132,135,138,140,  // L
143,145,148,150,153,  // M
156,158,161,163,166,  // M
171,174,176,179,	  // M
184,187,189,192,194,  // I'
197,200,202,205,207,  // I'
210,212,215,218,220,  // J'
223,225,228,231,233,  // J'
236,238,241,243,246,  // K'
249,251,254,256,259,  // K'
262,264,267,269,272,  // K'
275,277,280,282,285,  // K'
287,290,293,295,298,  // L'
300,303,306,308,311,  // L'
313,316,318,321,324,  // M'
326,329,331,334,337,  // M'
339,342,344,347,349,  // N
352,355,357,360,362,  // N
365,368,370,373,375,  // N
378,380,383,386,388,  // O
391,393,396,399,401,  // O
404,406,409,412,414,  // O
417,419,422,424,427,  // P
430,432,435,437,440,  // P
443,445,448,450,453,  // P
458,461,463,466}; //
unsigned int z[75]={1,9,	 // I  (first elements of length pairs which
14,22,	 // I	when decremented by 3 give invalid
27,35,	 // I	lengths)
40,45,	 // J
53,58,	 // J
66,71,76,	 // K
84,89,	 // K
97,102,	 // K
107,115,	 // L
120,128,	 // L
133,141,	 // L
146,151,	 // M
159,164,	 // M
172,177,	 // M
182,190,	 // I'
195,203,	 // I'
208,213,	 // J'
221,226,	 // J'
234,239,244, // K'
252,257,	 // K'
265,270,	 // K'
278,283,	 // K'
288,296,	 // L'
301,309,	 // L'
314,319,	 // M'
327,332,	 // M'
340,345,350, // N
358,363,	 // N
371,376,	 // N
384,389,	 // O
394,402,	 // O
407,415,	 // O
420,428,	 // P
433,438,	 // P
446,451,	 // P
459,464};	 //
unsigned int V[2],W[2],Y[2],M[2],N[2],O[2],P[2];
FILE *Outfp;
Outfp = fopen("out2.dat","w");
if (c==1) {
W[1]=1;
}
else {
W[1]=0xffffffff;
}
for (i=0; i<2048; i++)			    // clear histogram of lengths
he[i]=0;				    // for limbs ending in E
//
// DEAD LIMBS (ending in an even natural number)
//
count=0;
for (i=9999999; i>3; i-=6) {
if (count<200) {
printf("\n");
fprintf(Outfp,"\n");
}
V[1]=i;
n=1;			 // element count
//    if (count<200)
//	 printf("k=%d %d,  n=%d \n",V[0],V[1],n);
order=3;			 // set order
while (order<i)
order=order*2;
O[0]=0;
O[1]=order;
//
// begin loop
//
aloop:
if (V[0]>0x10000000)
goto cskip;
P[0]=V[0];
P[1]=V[1];
Y[0]=V[0];
Y[1]=V[1];
n=n+1;			 // increment element count
//    if (count<200)
//	 printf("k=%d %d, n=%d \n",V[0],V[1],n);
//
Y[0]=O[0];
Y[1]=O[1];
sub64(V, Y);		 // n-order
while ((Y[0]&0x80000000)==0) {
Y[0]=O[0];
Y[1]=O[1];
sub64(V, Y);		 // n-order
}
//
shift(V, Y, 3);
sub64(V, Y);
if ((Y[0]==0)&&(Y[1]==0)) {  // 8 divides 3n+c
np=n+1;		  // increment element count
M[0]=0;
M[1]=i;
shift(O, Y, 1);
N[0]=Y[0];
N[1]=Y[1];
sub64(M, N);
while ((N[0]&0x80000000)!=0) {
np=np+1;
N[0]=Y[0];
N[1]=Y[1];
sub64(M, N);
}
if (count<200)  {
shift(V, Y, 1);
printf("limb=(%#4x %#10x, %#4x %#10x), order=%#4x %#10x, n=%d",M[0],M[1],Y[0],Y[1],O[0],O[1],np);
fprintf(Outfp,"limb=(%#4x %#10x, %#4x %#10x), order=%#4x %#10x, n=%d",M[0],M[1],Y[0],Y[1],O[0],O[1],np);
}
div64_32(O, Y, 6);	  // order/6
sub64(P, Y);
if ((Y[0]&0x80000000)==0) {
Y[0]=V[0];
Y[1]=V[1];
sub64(O, Y);
if ((Y[0]&0x80000000)==0) {
printf("error: Next-to-last element does not determine order.\n");
goto zskip;
}
if (count<200) {
printf(", gt \n");
fprintf(Outfp,", gt \n");
}
for (j=0; j<75; j++) {
if (z[j]==np)
goto bskip;
}
if (np<471)
np=np-3;
bskip:	    if (np<2048)
he[np]=he[np]+1;
else {
printf("length too large to histogram \n");
goto zskip;
}
}
else {
if (count<200) {
printf("\n");
fprintf(Outfp,"\n");
}
}
count=count+1;
}
shift(V, Y, 1);
sub64(V, Y);	       // n-(n/2)*2
while ((Y[0]==0)&&(Y[1]==0)) {  // 2 divides n
shift(V, V, 1);
n=n+1; 	       // increment number of terms in sequence
//	 if (count<200)
//	    printf("k=%d %d, n=%d \n",V[0],V[1],n);
shift(V, Y, 1);
sub64(V, Y);	       // n-(n/2)*2
}
if (c==1) {
Y[0]=0;
Y[1]=1;
sub64(V, Y);
if ((Y[0]!=0)||(Y[1]!=0))
goto aloop;
}
if (c==-1) {
Y[0]=0;
Y[1]=1;
sub64(V,Y);
if ((Y[0]==0)&&(Y[1]==0))
goto cskip;
Y[0]=0;
Y[1]=5;
sub64(V,Y);
if ((Y[0]==0)&&(Y[1]==0))
goto cskip;
Y[0]=0;
Y[1]=7;
sub64(V,Y);
if ((Y[0]==0)&&(Y[1]==0))
goto cskip;
Y[0]=0;
Y[1]=17;
sub64(V,Y);
if ((Y[0]==0)&&(Y[1]==0))
goto cskip;
Y[0]=0;
Y[1]=25;
sub64(V,Y);
if ((Y[0]==0)&&(Y[1]==0))
goto cskip;
Y[0]=0;
Y[1]=37;
sub64(V,Y);
if ((Y[0]==0)&&(Y[1]==0))
goto cskip;
Y[0]=0;
Y[1]=55;
sub64(V,Y);
if ((Y[0]==0)&&(Y[1]==0))
goto cskip;
Y[0]=0;
Y[1]=41;
sub64(V,Y);
if ((Y[0]==0)&&(Y[1]==0))
goto cskip;
Y[0]=0;
Y[1]=61;
sub64(V,Y);
if ((Y[0]==0)&&(Y[1]==0))
goto cskip;
Y[0]=0;
Y[1]=91;
sub64(V,Y);
if ((Y[0]==0)&&(Y[1]==0))
goto cskip;
goto aloop;
}
cskip:n=0;
}
//
// check for valid lengths
//
for (i=0; i<178; i++) {
if (he[y[i]]!=0) {
printf("error: n=%d \n",y[i]);
goto zskip;
}
}
//
// HISTOGRAM OF LENGTHS
//
fprintf(Outfp,"\n");
fprintf(Outfp,"HISTOGRAM \n");
printf("\n");
for (h=0; h<100; h++) {
if (h==58)
printf("\n");
fprintf(Outfp," %d %d %d %d %d %d %d %d \n",he[8*h],he[8*h+1],he[8*h+2],
he[8*h+3],he[8*h+4],he[8*h+5],he[8*h+6],he[8*h+7]);
printf(" %d %d %d %d %d %d %d %d \n",he[8*h],he[8*h+1],he[8*h+2],
he[8*h+3],he[8*h+4],he[8*h+5],he[8*h+6],he[8*h+7]);
}
zskip:
fclose(Outfp);
return(0);
}
