﻿ /*******************************
/******************************************************************************/
/*									      */
/*  COMPUTE 3N+1 OR 3N-1 SEQUENCE					      */
/*  02/20/10 (dkc)							      */
/*									      */
/*  This C program generates least-residue trees and verifies that all least- */
/*  residues are present.  The trees for orders from k=2 to 16 are computed.  */
/*  That x<0 for limbs with negative limits and x>0 for limbs with positive   */
/*  limits is verified. 						      */
/*									      */
/******************************************************************************/
#include <stdio.h>
#include <math.h>
int main () {
//
unsigned int c=-1;			    // c value
unsigned int l=17;			    // limb length
//
// negative limits for limb lengths of 4,9,14,17,22,27,30
// positive limits for limb lengths of 2,5,7,10,12,15,18,20,23,25,28
//
unsigned int h,i,j,k,n,offset,order,flag;   // indices
unsigned int first,oldn,total,count,shift;  // flags, counters
unsigned int v[512],w[512];		    // sub-sequence arrays
unsigned int g[196608]; 		    // least-residue array
unsigned short ho[512],he[512]; 	    // histograms
unsigned int y[40]={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
FILE *Outfp;
Outfp = fopen("out8.dat","w");
if (l>30) {
printf("limb length too large \n");
goto zskip;
}
if (l==1) {
printf("one-element limbs are not processed \n");
goto zskip;
}
if ((l==3)||(l==6)||(l==8)||(l==11)||(l==13)||(l==16)||(l==19)||(l==21)||(l==24)||(l==26)||(l==29)) {
printf("invalid limb length \n");
goto zskip;
}
if ((c!=1)&&(c!=-1)) {
printf("invalid c value \n");
goto zskip;
}
//
// begin order loop
//
for (shift=2; shift<17; shift++) {
order=(1<<shift)*3;
if (order>196608)			       // check if order in range
goto zskip;
printf("order=%d \n",order);
//
// set flag
//
if ((l==4)||(l==9)||(l==14)||(l==17)||(l==22)||(l==27)||(l==30))
flag=0;
else
flag=1;
for (i=0; i<512; i++) {		       // clear histograms of lengths
ho[i]=0;				       // for limbs in S(L)
he[i]=0;				       // for limbs in E
}
for (i=0; i<order; i++)		       // clear least-residue array
g[i]=0;
g[1]=1;				       // unreachable sub-sequence
g[2]=2;				       //   "
if (c==1)
g[4]=4;				       //   "
else {
g[10]=10; 			       // on the end of the dead limb
g[5]=5;				       // {54,27,80,40} when order>48
g[14]=14;
g[7]=7;
g[20]=20;
}
if (c==1)
offset=8;
else
offset=4;
for (i=0; i<(order/24); i++) 	       // set U values
g[12*i+(order/2)+offset]=12*i+(order/2)+offset;
for (i=0; i<(order/4); i++)		       // set T values
g[2*i+(order/2)+1]=2*i+(order/2)+1;
count=0;
//
// left (limbs in A)
//
if (c==1)
offset=1;
else
offset=7;
for (h=(order/3)+offset; h<order/2; h+=8) {
first=1;				       // set first sub-sequence flag
oldn=0;				       // clear length
for (i=order/2; i<order; i+=2) {	       // possible starting elements
k=i;				       // set starting element of limb
v[0]=k;			       // save starting element
n=1;				       // set index
while (k==(k/2)*2) {		       // check for even element
k=k/2;			       // next element of limb
v[n]=k;			       // save element
n=n+1;			       // increment index
}
if (c==-1) {
if ((k==5)||(k==7))
goto bskip;
}
for (j=1; j<100000; j++) {	       // iterate until end of limb
if (k==h)			       // end of limb
if (k==1)
goto bskip;
k=3*k+c;			       // next element
if (k>order)		       // not in tree, start over
goto bskip;
if (k==(k/8)*8)		       // not valid limb, start over
goto bskip;
v[n]=k;			       // save element
n=n+1;			       // increment index
while (k==(k/2)*2) {	       // check for even element
k=k/2;			       // next element
v[n]=k;			       // save element
n=n+1;			       // increment index
}
}
goto zskip;
if (first==1) {		       // check first-time flag
for (j=0; j<n; j++) 	       // save sub-sequence
w[j]=v[j];
oldn=n;			       // save length
first=0;			       // clear first-time flag
}
else {
if (n>oldn) {		       // check for longer limb
for (j=0; j<n; j++)	       // replace sub-sequence
w[j]=v[j];
oldn=n;			       // save length
}
}
bskip:
n=0;
}
if (oldn<512)			       // update length histogram
ho[oldn]=ho[oldn]+1;
else
printf("error A: ho array not big enough \n");
if (oldn>512) {
printf("error: w array not big enough \n");
goto zskip;
}
for (j=0; j<oldn; j++) {		       // update least-residues array
if (g[w[j]]==0)
g[w[j]]=w[j];
else {
printf("error A: redundant value in sequence=%d \n",w[j]);
goto zskip;
}
}
if (oldn==l) {
count=count+1;
j=(3*w[oldn-1]+c)/2;
if (flag==1) {
if (w[0]<=j) {
printf("error: less than or equal \n");
printf("limb=(%d, %d)=>%d \n",w[0],w[oldn-1],j);
goto zskip;
}
}
else {
if (w[0]>=j) {
printf("error: greater than or equal \n");
printf("limb=(%d, %d)=>%d \n",w[0],w[oldn-1],j);
goto zskip;
}
}
}
}
printf("A count=%d \n",count);
//
// center (limbs in C)
//
if (c==1)
offset=5;
else
offset=3;
for (h=(order/3)+offset; h<order/2; h+=8) {
first=1;
oldn=0;
for (i=order/2; i<order; i+=2) {
k=i;
v[0]=k;
n=1;
while (k==(k/2)*2) {
k=k/2;
v[n]=k;
n=n+1;
}
if (c==-1) {
if ((k==5)||(k==7))
goto dskip;
}
for (j=1; j<100000; j++) {
if (k==h)
goto cskip;
if (k==1)
goto dskip;
k=3*k+c;
if (k>order)
goto dskip;
if (k==(k/8)*8)
goto dskip;
v[n]=k;
n=n+1;
while (k==(k/2)*2) {
k=k/2;
v[n]=k;
n=n+1;
}
}
printf("error C \n");
goto zskip;
cskip:
if (first==1) {
for (j=0; j<n; j++)
w[j]=v[j];
oldn=n;
first=0;
}
else {
if (n>oldn) {
for (j=0; j<n; j++)
w[j]=v[j];
oldn=n;
}
}
dskip:
n=0;
}
if (oldn<512)
ho[oldn]=ho[oldn]+1;
else
printf("error C: ho array not big enough \n");
if (oldn>512) {
printf("error: w array not big enough \n");
goto zskip;
}
for (j=0; j<oldn; j++) {
if (g[w[j]]==0)
g[w[j]]=w[j];
else {
printf("error B: redundant value in sequence=%d \n",w[j]);
goto zskip;
}
}
if (oldn==l) {
count=count+1;
j=(3*w[oldn-1]+c)/2;
if (flag==1) {
if (w[0]<=j) {
printf("error: less than or equal \n");
printf("limb=(%d, %d)=>%d \n",w[0],w[oldn-1],j);
goto zskip;
}
}
else {
if (w[0]>=j) {
printf("error: greater than or equal \n");
printf("limb=(%d, %d)=>%d \n",w[0],w[oldn-1],j);
goto zskip;
}
}
}
}
printf("C count=%d \n",count);
//
// left of center (limbs in B)
//
if (c==1)
offset=3;
else
offset=5;
for (h=(order/3)+offset; h<order/2; h+=8) {
first=1;
oldn=0;
for (i=order/2; i<order; i+=2) {
k=i;
v[0]=k;
n=1;
while (k==(k/2)*2) {
k=k/2;
v[n]=k;
n=n+1;
}
if (c==-1) {
if ((k==5)||(k==7))
goto fskip;
}
for (j=1; j<100000; j++) {
if (k==h)
goto eskip;
if (k==1)
goto fskip;
k=3*k+c;
if (k>order)
goto fskip;
if (k==(k/8)*8)
goto fskip;
v[n]=k;
n=n+1;
while (k==(k/2)*2) {
k=k/2;
v[n]=k;
n=n+1;
}
}
printf("error B \n");
goto zskip;
eskip:
if (first==1) {
for (j=0; j<n; j++)
w[j]=v[j];
oldn=n;
first=0;
}
else {
if (n>oldn) {
for (j=0; j<n; j++)
w[j]=v[j];
oldn=n;
}
}
fskip:
n=0;
}
if (oldn<512)
ho[oldn]=ho[oldn]+1;
else
printf("error B: ho array not big enough \n");
if (oldn>512) {
printf("error: w array not big enough \n");
goto zskip;
}
for (j=0; j<oldn; j++) {
if (g[w[j]]==0)
g[w[j]]=w[j];
else {
printf("error C: redundant value in sequence=%d \n",w[j]);
goto zskip;
}
}
if (oldn==l) {
count=count+1;
j=(3*w[oldn-1]+c)/2;
if (flag==1) {
if (w[0]<=j) {
printf("error: less than or equal \n");
printf("limb=(%d, %d)=>%d \n",w[0],w[oldn-1],j);
goto zskip;
}
}
else {
if (w[0]>=j) {
printf("error: greater than or equal \n");
printf("limb=(%d, %d)=>%d \n",w[0],w[oldn-1],j);
goto zskip;
}
}
}
}
printf("B count=%d \n",count);
//
//  right of center (limbs in D)
//
if (c==1)
offset=7;
else
offset=1;
for (h=(order/3)+offset; h<order/2; h+=8) {
first=1;
oldn=0;
for (i=order/2; i<order; i+=2) {
k=i;
v[0]=k;
n=1;
while (k==(k/2)*2) {
k=k/2;
v[n]=k;
n=n+1;
}
if (c==-1) {
if ((k==5)||(k==7))
goto hskip;
}
for (j=1; j<100000; j++) {
if (k==h)
goto gskip;
if (k==1)
goto hskip;
k=3*k+c;
if (k>order)
goto hskip;
if (k==(k/8)*8)
goto hskip;
v[n]=k;
n=n+1;
while (k==(k/2)*2) {
k=k/2;
v[n]=k;
n=n+1;
}
}
printf("error D \n");
goto zskip;
gskip:
if (first==1) {
for (j=0; j<n; j++)
w[j]=v[j];
oldn=n;
first=0;
}
else {
if (n>oldn) {
for (j=0; j<n; j++)
w[j]=v[j];
oldn=n;
}
}
hskip:
n=0;
}
if (oldn<512)
ho[oldn]=ho[oldn]+1;
else
printf("error D: ho array not big enough \n");
if (oldn>512) {
printf("error: w array not big enough \n");
goto zskip;
}
for (j=0; j<oldn; j++) {
if (g[w[j]]==0)
g[w[j]]=w[j];
else {
printf("error D: redundant value in sequence=%d \n",w[j]);
goto zskip;
}
}
if (oldn==l) {
count=count+1;
j=(3*w[oldn-1]+c)/2;
if (flag==1) {
if (w[0]<=j) {
printf("error: less than or equal \n");
printf("limb=(%d, %d)=>%d \n",w[0],w[oldn-1],j);
goto zskip;
}
}
else {
if (w[0]>=j) {
printf("error: greater than or equal \n");
printf("limb=(%d, %d)=>%d \n",w[0],w[oldn-1],j);
goto zskip;
}
}
}
}
printf("D count=%d \n",count);
//
// down (elements of E)
//
total=0;					// clear counter
if (c==1)
offset=2;
else
offset=22;
for (h=(order/2)+offset; h<3*(order/4); h+=24) {  // possible starting elements(not divisible by 3)
k=h;					// set starting element
v[0]=k;					// save starting element
n=1;					// set index
while (k==(k/2)*2) {			// check for even element
k=k/2; 				// next element
v[n]=k;				// save element
n=n+1; 				// increment index
}
if (c==-1) {
if ((k==5)||(k==7))
goto jskip;
}
for (j=1; j<100000; j++) {		// iterate until end of limb
k=3*k+c;				// next element
if (k>order)				// check if in tree
goto jskip; 			// invalid limb, start over
v[n]=k;				// save element
n=n+1; 				// increment index
if (k==(k/8)*8) {			// check if 8 divides element
v[n]=k/2;				// save element
n=n+1;				// increment index
goto iskip; 			// valid limb
}
while (k==(k/2)*2) {			// check for even element
k=k/2;				// next element
v[n]=k;				// save element
n=n+1;				// increment index
}
}
goto zskip;
iskip:
if (n>512) {
printf("error: v array not big enough \n");
goto zskip;
}
if (n<512)				// update histogram of lengths
he[n]=he[n]+1;
else
printf("error: he array not big enough \n");
total=total+1;				// increment number of limbs in E
for (j=0; j<n; j++) {			// update least-residues array
if (g[v[j]]==0)
g[v[j]]=v[j];
else {
printf("error E: redundant value in sequence=%d \n",v[j]);
goto zskip;
}
}
jskip:
n=0;
}
if ((he[2]!=0)||(he[3]!=0)||(he[5]!=0)||(he[6]!=0)||(he[8]!=0)||(he[11]!=0)) {
printf("E error: invalid limb length \n");
goto zskip;
}
//
// DEAD LIMBS (ending in an even natural number)
//
for (i=3; i<(order/2); i+=6) {
k=i;
while (k<(order/2))
k=k*2;
v[0]=k;
n=1;
while (k==(k/2)*2) {
k=k/2;
v[n]=k;
n=n+1;
}
if (c==-1) {
if ((k==5)||(k==7))
goto mskip;
}
for (j=1; j<100000; j++) {
k=3*k+c;
if (k>order)
goto mskip;
v[n]=k;
n=n+1;
if (k==(k/8)*8) {
v[n]=k/2;
n=n+1;
goto lskip;
}
while (k==(k/2)*2) {
k=k/2;
v[n]=k;
n=n+1;
}
}
printf("error \n");
goto zskip;
lskip:
if (n>512) {
printf("error: v array not big enough \n");
goto zskip;
}
total=total+1;
for (j=0; j<n; j++) {
if (g[v[j]]==0)
g[v[j]]=v[j];
else {
printf("error F: redundant value in sequence=%d \n",v[j]);
goto zskip;
}
}
mskip:
n=0;
}
if (total!=(order/24))
printf("error: E and F count=%d \n",total);
//
// CHECK IF ALL ELEMENTS PRESENT
//
for (h=1; h<order; h++) {
if (g[h]!=h) {
fprintf(Outfp,"error: not covered=%d \n",h);
printf("error: not covered=%d \n",h);
}
}
//
// CHECK LIMB LENGTHS
//
for (h=0; h<512; h++) {
if (ho[h]==0)
continue;
for (i=1; i<40; i++) {
if (h==y[i]) {
printf("error: invalid limb length=%d \n",y[i]);
goto zskip;
}
}
}
printf("\n");
}
zskip:
fclose(Outfp);
return(0);
}
