/******************************************************************************/
/* */
/* 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
goto askip;
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
}
}
printf("error A \n"); // limb not found
goto zskip;
askip:
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
}
}
printf("error E \n"); // limb not found
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);
}