﻿ /*******************************
```/******************************************************************************/
/*									      */
/*  COMPUTE 3N+1 OR 3N-1 SEQUENCE					      */
/*  01/03/09 (dkc)							      */
/*									      */
/*  This C program generates a least-residue tree for k<=15.  Numerous	      */
/*  properties are verified.						      */
/*									      */
/******************************************************************************/
#include <stdio.h>
#include <math.h>
int main () {
//
unsigned int order=98304;  // 3*2**k
unsigned int l=51;			    // length
int e=-1;				    // c
//
unsigned int h,i,j,k,m,n;		    // indices
unsigned int first,oldn,count,total;	    // flags, counters
unsigned int v[8192],w[8192];		    // sub-sequence arrays
unsigned int a[6144],b[6144],c[6144],d[6144];  // S(L) arrays
unsigned int g[98304];			    // least-residue array
unsigned int offset,delta,min;
unsigned short ho[512],he[512]; 	    // histograms
double r;
FILE *Outfp;
Outfp = fopen("out10.dat","w");
if (order>98304)			    // check if order in range
goto zskip;
min=0x7fffffff;
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<98304; i++) 		    // clear least-residue array
g[i]=0;
g[1]=1; 				    // unreachable sub-sequence
g[2]=2; 				    //	 "
if (e==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 (e==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;
fprintf(Outfp," ORDER=%d \n",order);
//
// left (limbs in A)
//
if (order<=24576)
fprintf(Outfp," LIMBS IN A \n");
m=0;					    // clear index for A array
if (e==1)
delta=1;
else
delta=7;
for (h=(order/3)+delta; 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 (e==-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+e;			    // next element
if (k>order)			    // not in tree, start over
goto bskip;
if (k==(k/8)*8)		    // not valid limb, start over
goto bskip; 		    // (prevents taking even path at nodes)
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 (order<=24576) {
if (w[0]==(w[0]/3)*3) {		    // check for dead limb
fprintf(Outfp,"D");                // output limb
fprintf(Outfp," %d (%d,%d)=>%d \n",oldn,w[0],w[oldn-1],(3*w[oldn-1]+e)/2);
}
else{
fprintf(Outfp,"  %d (%d,%d)=>%d \n",oldn,w[0],w[oldn-1],(3*w[oldn-1]+e)/2 );
if (w[0]==(w[0]/4)*4)		    // check if second element is odd
printf("error: second element not odd \n");
}
}
r=(float)(w[0])-(float)((3*w[oldn-1]+e)/2);
r=r/(float)(w[0]);
if (oldn==l)
printf("A=%d %d %d %f \n",w[0],w[oldn-1],oldn,r);
a[3*m]=w[0]; 			    // save limb
a[3*m+1]=w[oldn-1];
a[3*m+2]=oldn;
m=m+1;
if (oldn<512)			    // update length histogram
ho[oldn]=ho[oldn]+1;
else
printf("error A: ho array not big enough \n");
if (oldn>8192) {
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 (w[0]!=(w[0]/3)*3) {
for (j=0; j<oldn; j++) {
if (w[j]<min)
min=w[j];
}
}
if (oldn==5) {			    // check for various lengths
if (w[0]!=(w[0]/3)*3)
printf("error: length=5, but not dead \n");
}
if (oldn==10) {			     // check for various lengths
if (w[0]!=(w[0]/3)*3)
printf("error: length=10, but not dead \n");
}
if (oldn==15) {			     // check for various lengths
if (w[0]!=(w[0]/3)*3)
printf("error: length=15, but not dead \n");
}
if (oldn==18) {			     // check for various lengths
if (w[0]!=(w[0]/3)*3)
printf("error: length=18, but not dead \n");
}
if (oldn==23) {			     // check for various lengths
if (w[0]!=(w[0]/3)*3)
printf("error: length=23, but not dead \n");
}
if (oldn==28) {			     // check for various lengths
if (w[0]!=(w[0]/3)*3)
printf("error: length=28, but not dead \n");
}
if (oldn==6)
printf("error: length=6 \n");
if (oldn==8)
printf("error: length=8 \n");
if (oldn==11)
printf("error: length=11 \n");
if (oldn==13)
printf("error: length=13 \n");
if (oldn==16)
printf("error: length=16 \n");
if (oldn==19)
printf("error: length=19 \n");
if (oldn==21)
printf("error: length=21 \n");
}
printf("minimum=%d \n",min);
//
// center (limbs in C)
//
min=0x7fffffff;
if (order<=24576) {
fprintf(Outfp,"\n");
fprintf(Outfp," LIMBS IN C \n");
}
m=0;
if (e==1)
delta=5;
else
delta=3;
for (h=(order/3)+delta; 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 (e==-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+e;
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: h=%d \n",h);
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 (order<=24576) {
if (w[0]==(w[0]/3)*3) {
fprintf(Outfp,"D");
fprintf(Outfp," %d (%d,%d)=>%d \n",oldn,w[0],w[oldn-1],(3*w[oldn-1]+e)/2);
}
else {
fprintf(Outfp,"  %d (%d,%d)=>%d \n",oldn,w[0],w[oldn-1],(3*w[oldn-1]+e)/2 );
if (w[0]==(w[0]/4)*4)
printf("error: second element not odd \n");
}
}
r=(float)(w[0])-(float)((3*w[oldn-1]+e)/2);
r=r/(float)(w[0]);
if (oldn==l)
printf("C=%d %d %d %f \n",w[0],w[oldn-1],oldn,r);
c[3*m]=w[0];
c[3*m+1]=w[oldn-1];
c[3*m+2]=oldn;
m=m+1;
if (oldn<512)
ho[oldn]=ho[oldn]+1;
else
printf("error C: ho array not big enough \n");
if (oldn>8192) {
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 (w[0]!=(w[0]/3)*3) {
for (j=0; j<oldn; j++) {
if (w[j]<min)
min=w[j];
}
}
if (oldn==5) {
if (w[0]!=(w[0]/3)*3)
printf("error: length=5, but not dead \n");
}
if (oldn==10) {			     // check for various lengths
if (w[0]!=(w[0]/3)*3)
printf("error: length=10, but not dead \n");
}
if (oldn==15) {			     // check for various lengths
if (w[0]!=(w[0]/3)*3)
printf("error: length=15, but not dead \n");
}
if (oldn==18) {			     // check for various lengths
if (w[0]!=(w[0]/3)*3)
printf("error: length=18, but not dead \n");
}
if (oldn==23) {			     // check for various lengths
if (w[0]!=(w[0]/3)*3)
printf("error: length=23, but not dead \n");
}
if (oldn==28) {			     // check for various lengths
if (w[0]!=(w[0]/3)*3)
printf("error: length=28, but not dead \n");
}
if (oldn==6)
printf("error: length=6 \n");
if (oldn==8)
printf("error: length=8 \n");
if (oldn==11)
printf("error: length=11 \n");
if (oldn==13)
printf("error: length=13 \n");
if (oldn==16)
printf("error: length=16 \n");
if (oldn==19)
printf("error: length=19 \n");
if (oldn==21)
printf("error: length=21 \n");
}
printf("minimum=%d \n",min);
//
// left of center (limbs in B)
//
min=0x7fffffff;
if (order<=24576) {
fprintf(Outfp,"\n");
fprintf(Outfp," LIMBS IN B \n");
}
m=0;
if (e==1)
delta=3;
else
delta=1;
for (h=(order/3)+delta; 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 (e==-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+e;
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 (order<=24576) {
if (w[0]==(w[0]/3)*3) {
fprintf(Outfp,"D");
fprintf(Outfp," %d (%d,%d)=>%d \n",oldn,w[0],w[oldn-1],(3*w[oldn-1]+e)/2);
}
else {
fprintf(Outfp,"  %d (%d,%d)=>%d \n",oldn,w[0],w[oldn-1],(3*w[oldn-1]+e)/2 );
if (w[0]==(w[0]/4)*4)
printf("error: second element not odd \n");
}
}
r=(float)(w[0])-(float)((3*w[oldn-1]+e)/2);
r=r/(float)(w[0]);
if (oldn==l)
printf("B=%d %d %d %f \n",w[0],w[oldn-1],oldn,r);
b[3*m]=w[0];
b[3*m+1]=w[oldn-1];
b[3*m+2]=oldn;
m=m+1;
if (oldn<512)
ho[oldn]=ho[oldn]+1;
else
printf("error B: ho array not big enough \n");
if (oldn>8192) {
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 (w[0]!=(w[0]/3)*3) {
for (j=0; j<oldn; j++) {
if (w[j]<min)
min=w[j];
}
}
if (oldn==5) {
if (w[0]!=(w[0]/3)*3)
printf("error: length=5, but not dead \n");
}
if (oldn==10) {			     // check for various lengths
if (w[0]!=(w[0]/3)*3)
printf("error: length=10, but not dead \n");
}
if (oldn==15) {			     // check for various lengths
if (w[0]!=(w[0]/3)*3)
printf("error: length=15, but not dead \n");
}
if (oldn==18) {			     // check for various lengths
if (w[0]!=(w[0]/3)*3)
printf("error: length=18, but not dead \n");
}
if (oldn==23) {			     // check for various lengths
if (w[0]!=(w[0]/3)*3)
printf("error: length=23, but not dead \n");
}
if (oldn==28) {			     // check for various lengths
if (w[0]!=(w[0]/3)*3)
printf("error: length=28, but not dead \n");
}
if (oldn==6)
printf("error: length=6 \n");
if (oldn==8)
printf("error: length=8 \n");
if (oldn==11)
printf("error: length=11 \n");
if (oldn==13)
printf("error: length=13 \n");
if (oldn==16)
printf("error: length=16 \n");
if (oldn==19)
printf("error: length=19 \n");
if (oldn==21)
printf("error: length=21 \n");
}
printf("minimum=%d \n",min);
//
//  right of center (limbs in D)
//
min=0x7fffffff;
if (order<=24576) {
fprintf(Outfp,"\n");
fprintf(Outfp," LIMBS IN D \n");
}
m=0;
if (e==1)
delta=7;
else
delta=5;
for (h=(order/3)+delta; 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 (e==-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+e;
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 (order<=24576) {
if (w[0]==(w[0]/3)*3) {
fprintf(Outfp,"D");
fprintf(Outfp," %d (%d,%d)=>%d \n",oldn,w[0],w[oldn-1],(3*w[oldn-1]+e)/2);
}
else {
fprintf(Outfp,"  %d (%d,%d)=>%d \n",oldn,w[0],w[oldn-1],(3*w[oldn-1]+e)/2 );
if (w[0]==(w[0]/4)*4)
printf("error: second element not odd \n");
}
}
r=(float)(w[0])-(float)((3*w[oldn-1]+e)/2);
r=r/(float)(w[0]);
if (oldn==l)
printf("D=%d %d %d %f \n",w[0],w[oldn-1],oldn,r);
d[3*m]=w[0];
d[3*m+1]=w[oldn-1];
d[3*m+2]=oldn;
m=m+1;
if (oldn<512)
ho[oldn]=ho[oldn]+1;
else
printf("error D: ho array not big enough \n");
if (oldn>8192) {
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 (w[0]!=(w[0]/3)*3) {
for (j=0; j<oldn; j++) {
if (w[j]<min)
min=w[j];
}
}
if (oldn==5) {
if (w[0]!=(w[0]/3)*3)
printf("error: length=5, but not dead \n");
}
if (oldn==10) {			     // check for various lengths
if (w[0]!=(w[0]/3)*3)
printf("error: length=10, but not dead \n");
}
if (oldn==15) {			     // check for various lengths
if (w[0]!=(w[0]/3)*3)
printf("error: length=15, but not dead \n");
}
if (oldn==18) {			     // check for various lengths
if (w[0]!=(w[0]/3)*3)
printf("error: length=18, but not dead \n");
}
if (oldn==23) {			     // check for various lengths
if (w[0]!=(w[0]/3)*3)
printf("error: length=23, but not dead \n");
}
if (oldn==28) {			     // check for various lengths
if (w[0]!=(w[0]/3)*3)
printf("error: length=28, but not dead \n");
}
if (oldn==6)
printf("error: length=6 \n");
if (oldn==8)
printf("error: length=8 \n");
if (oldn==11)
printf("error: length=11 \n");
if (oldn==13)
printf("error: length=13 \n");
if (oldn==16)
printf("error: length=16 \n");
if (oldn==19)
printf("error: length=19 \n");
if (oldn==21)
printf("error: length=21 \n");
}
printf("minimum=%d \n",min);
//
// down (elements of E)
//
if (order<=24576) {
fprintf(Outfp,"\n");
fprintf(Outfp," LIMBS IN E \n");
}
else
fprintf(Outfp,"LIMBS OF EQUAL LENGTH (limbs in A attaching to limbs in E) \n");
total=0;				     // clear counter
oldn=0; 				     // clear counter
if (e==1)
delta=2;
else
delta=22;
for (h=(order/2)+delta; 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 (e==-1) {
if ((k==5)||(k==7))
goto jskip;
}
for (j=1; j<100000; j++) {		     // iterate until end of limb
k=3*k+e;				     // next element
if (k>order) {			     // check if in tree
w[oldn]=h;			     // save starting element
oldn=oldn+1;			     // increment index
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>8192) {
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
if (n==5)
printf("error: length=5 \n");
if (n==6)
printf("error: length=6 \n");
if (n==8)
printf("error: length=8 \n");
if (n==11)
printf("error: length=11 \n");
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;
}
}
if (order<=24576) {
if (v[0]==(v[0]/3)*3) {		     // check for dead limb
fprintf(Outfp,"D");                 // output limb
fprintf(Outfp," %d (%d,%d) \n",n,v[0],v[n-1]);
}
else {
fprintf(Outfp,"  %d (%d,%d) \n",n,v[0],v[n-1]);
if (v[0]==(v[0]/4)*4)		     // check if second element is odd
printf("error: second element not odd \n");
}
}
else {
for (i=0; i<(order/48); i+=2) {	     // check if limb in A attaches
if ((3*a[3*i+1]+e)/2==v[0]) {	     // compare last to first element
if (n==a[3*i+2]) {		     // check if equal lengths
fprintf(Outfp," E=%d %d %d   A=%d %d %d \n",v[0],v[n-1],n,a[3*i],a[3*i+1],a[3*i+2]);
if ((n==4)||(n==17)) {
if (v[0]<a[3*i])
printf("error: length=4 or 17 \n");
}
else {
if (v[0]>a[3*i])
printf("error: length!=4 or 17 \n");
}
r=(float)(a[3*i])-(float)(v[0]);
r=r/(float)(a[3*i]);
goto jskip;
}
}
}
}
jskip:
n=0;
}
//
// CHECK FOR LIMBS OF EQUAL LENGTH
// (limbs in A(even) attaching to limbs in A, B, C, or D)
//
if (oldn>=8192) {
printf("error: w array not big enough \n");
goto zskip;
}
printf("\n");
printf("LIMBS OF EQUAL LENGTH \n");
if (order<=24576) {
fprintf(Outfp,"\n");
fprintf(Outfp,"ATTACHMENTS (limbs in A(even) to limbs in A, B, C, D) \n");
}
else {
fprintf(Outfp,"\n");
fprintf(Outfp,"LIMBS OF EQUAL LENGTH (limbs in A(even) attaching to limbs in A, B, C, D) \n");
}
n=0;
for (h=0; h<oldn; h++) {		     // limbs not dead, not in E
k=w[h];				     // starting element
if (k<(order/3)*2)			     // limit for attaching to 2-element limbs
n=n+1;				     // count
for (i=0; i<(order/48); i++) {	     // limbs ending in A
if ((3*a[3*i+1]+e)/2==k) {	     // check if attaches to starting element
count=a[3*i+2];		     // length of limb
for (j=0; j<(order/48); j++) {      // limbs ending in A
if (a[3*j]==k) {		     // check if equal starting elements
if (k<(order/3)*2) {	     // check limit
if (a[3*j+2]==count) {     // check for equal lengths
printf(" A=%d %d %d    A=%d %d %d \n",a[3*j],a[3*j+1],a[3*j+2],a[3*i],a[3*i+1],count);
if (order>24576) {
if (a[3*i]==(a[3*i]/3)*3)
fprintf(Outfp," A=%d %d %d    A=%d %d %d \n",a[3*j],a[3*j+1],a[3*j+2],a[3*i],a[3*i+1],count);
else
fprintf(Outfp," A=%d %d %d    A=%d %d %d   not dead \n",a[3*j],a[3*j+1],a[3*j+2],a[3*i],a[3*i+1],count);
}
if (count==17) {
if (a[3*j]<=a[3*i])
printf("error: length=17 \n");
}
else {
if (a[3*j]>=a[3*i])
printf("error: length!=17 \n");
}
}
}
else {
if (a[3*j+2]!=2)
printf("error: length not equal to 2 \n");
if (count==2) {
if (a[3*j]>=a[3*i])
printf("error: length==2 \n");
}
}
if (order<=24576)
fprintf(Outfp," A=%d %d %d   A=%d %d %d \n",a[3*j],a[3*j+1],a[3*j+2],a[3*i],a[3*i+1],count);
goto kskip;
}
}
for (j=0; j<(order/48); j++) {
if (b[3*j]==k) {
if (k<(order/3)*2) {
if (b[3*j+2]==count) {
printf(" B=%d %d %d    A=%d %d %d \n",b[3*j],b[3*j+1],b[3*j+2],a[3*i],a[3*i+1],count);
if (order>24576) {
if (a[3*i]==(a[3*i]/3)*3)
fprintf(Outfp," B=%d %d %d    A=%d %d %d \n",b[3*j],b[3*j+1],b[3*j+2],a[3*i],a[3*i+1],count);
else
fprintf(Outfp," B=%d %d %d    A=%d %d %d   not dead \n",b[3*j],b[3*j+1],b[3*j+2],a[3*i],a[3*i+1],count);
}
if (count==17) {
if (b[3*j]<=a[3*i])
printf("error: length=17 \n");
}
else {
if (b[3*j]>=a[3*i])
printf("error: length!=17 \n");
}
}
}
else {
if (b[3*j+2]!=2)
printf("error: length not equal to 2 \n");
if (count==2) {
if (b[3*j]>=a[3*i])
printf("error: length==2 \n");
}
}
if (order<=24576)
fprintf(Outfp," B=%d %d %d   A=%d %d %d \n",b[3*j],b[3*j+1],b[3*j+2],a[3*i],a[3*i+1],count);
goto kskip;
}
}
for (j=0; j<(order/48); j++) {
if (c[3*j]==k) {
if (k<(order/3)*2) {
if (c[3*j+2]==count) {
printf(" C=%d %d %d    A=%d %d %d \n",c[3*j],c[3*j+1],c[3*j+2],a[3*i],a[3*i+1],count);
if (order>24576) {
if (a[3*i]==(a[3*i]/3)*3)
fprintf(Outfp," C=%d %d %d    A=%d %d %d \n",c[3*j],c[3*j+1],c[3*j+2],a[3*i],a[3*i+1],count);
else
fprintf(Outfp," C=%d %d %d    A=%d %d %d   not dead \n",c[3*j],c[3*j+1],c[3*j+2],a[3*i],a[3*i+1],count);
}
if (count==17) {
if (c[3*j]<=a[3*i])
printf("error: length=17 \n");
}
else {
if (c[3*j]>=a[3*i])
printf("error: length!=17 \n");
}
}
}
else {
if (c[3*j+2]!=2)
printf("error: length not equal to 2 \n");
if (count==2) {
if (c[3*j]>=a[3*i])
printf("error: length==2 \n");
}
}
if (order<=24576)
fprintf(Outfp," C=%d %d %d   A=%d %d %d \n",c[3*j],c[3*j+1],c[3*j+2],a[3*i],a[3*i+1],count);
goto kskip;
}
}
for (j=0; j<(order/48); j++) {
if (d[3*j]==k) {
if (k<(order/3)*2) {
if (d[3*j+2]==count) {
printf(" D=%d %d %d    A=%d %d %d \n",d[3*j],d[3*j+1],d[3*j+2],a[3*i],a[3*i+1],count);
if (order>24576) {
if (a[3*i]==(a[3*i]/3)*3)
fprintf(Outfp," D=%d %d %d    A=%d %d %d \n",d[3*j],d[3*j+1],d[3*j+2],a[3*i],a[3*i+1],count);
else
fprintf(Outfp," D=%d %d %d    A=%d %d %d  not dead \n",d[3*j],d[3*j+1],d[3*j+2],a[3*i],a[3*i+1],count);
}
if (count==17) {
if (d[3*j]<=a[3*i])
printf("error: length=17 \n");
}
else {
if (d[3*j]>=a[3*i])
printf("error: length!=17 \n");
}
}
}
else {
if (d[3*j+2]!=2)
printf("error: length not equal to 2 \n");
if (count==2) {
if (d[3*j]>=a[3*i])
printf("error: length==2 \n");
}
}
if (order<=24576)
fprintf(Outfp," D=%d %d %d   A=%d %d %d \n",d[3*j],d[3*j+1],d[3*j+2],a[3*i],a[3*i+1],count);
goto kskip;
}
}
printf("error: k=%d \n",k);
goto zskip;
}
}
printf("error: k=%d \n",k);
goto zskip;
kskip:
k=0;
}
printf("count=%d \n",n);
printf("\n");
fprintf(Outfp,"count=%d \n",n);
printf("number of limbs in E=%d \n",total);
//
// 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 (e==-1) {
if ((k==5)||(k==7))
goto mskip;
}
for (j=1; j<100000; j++) {
k=3*k+e;
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>8192) {
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;
}
printf("number of limbs in E and F=%d \n",total);
if (total!=(order/24))
printf("error: correct count=%d \n",order/24);
//
// CHECK IF ALL ELEMENTS PRESENT
//
fprintf(Outfp,"\n");
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);
}
}
//
// HISTOGRAMS OF LENGTHS
//
fprintf(Outfp,"\n");
fprintf(Outfp,"HISTOGRAM (odds) \n");
for (h=0; h<64; h++)
fprintf(Outfp," %d %d %d %d %d %d %d %d \n",ho[8*h],ho[8*h+1],ho[8*h+2],
ho[8*h+3],ho[8*h+4],ho[8*h+5],ho[8*h+6],ho[8*h+7]);
fprintf(Outfp,"\n");
fprintf(Outfp,"HISTOGRAM (evens) \n");
for (h=0; h<64; h++)
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]);
//
// LIMBS IN A(odd)
//
if (e==1)
delta=14;
else
delta=10;
for (h=(order/2)+delta; h<3*(order/4); h+=24) { // possible starting elements (not divisible by 3)
k=h; 				     // set starting element
for (i=0; i<(order/48); i++) {	     // limbs ending in A
if ((3*a[3*i+1]+e)/2==k) {	     // check if attaches to starting element
count=a[3*i+2];		     // length of limb
for (j=0; j<(order/48); j++) {      // limbs ending in A
if (a[3*j]==k) {		     // check if equal starting elements
if ((a[3*j+2]!=2)&&(a[3*j+2]!=4))
printf("error A: length=%d \n",a[3*j+2]);
if ((count==2)&&(a[3*j+2]==2)) {
if (a[3*j]>=a[3*i])
printf("error: length==2 \n");
r=(float)(a[3*i])-(float)(a[3*j]);
r=r/(float)(a[3*i]);
printf("A: length=2, ratio=%f %d %d %d %d \n",r,a[3*i],a[3*i+1],a[3*j],a[3*j+1]);
}
if ((count==4)&&(a[3*j+2]==4)) {
if (a[3*j]<=a[3*i])
printf("error: length==4 \n");
r=(float)(a[3*j])-(float)(a[3*i]);
r=r/(float)(a[3*i]);
}
goto nskip;
}
}
for (j=0; j<(order/48); j++) {
if (b[3*j]==k) {
if ((b[3*j+2]!=2)&&(b[3*j+2]!=4))
printf("error B: length=%d \n",b[3*j+2]);
if ((count==2)&&(b[3*j+2]==2)) {
if (b[3*j]>=a[3*i])
printf("error: length==2 \n");
r=(float)(a[3*i])-(float)(b[3*j]);
r=r/(float)(a[3*i]);
}
if ((count==4)&&(b[3*j+2]==4)) {
if (b[3*j]<=a[3*i])
printf("error: length==4 \n");
r=(float)(b[3*j])-(float)(a[3*i]);
r=r/(float)(a[3*i]);
}
goto nskip;
}
}
for (j=0; j<(order/48); j++) {
if (c[3*j]==k) {
if ((c[3*j+2]!=2)&&(c[3*j+2]!=4))
printf("error C: length=%d \n",c[3*j+2]);
if ((count==2)&&(c[3*j+2]==2)) {
if (c[3*j]>=a[3*i])
printf("error: length==2 \n");
r=(float)(a[3*i])-(float)(c[3*j]);
r=r/(float)(a[3*i]);
printf("C: length=2, ratio=%f %d %d %d %d \n",r,a[3*i],a[3*i+1],c[3*j],c[3*j+1]);
}
if ((count==4)&&(c[3*j+2]==4)) {
if (c[3*j]<=a[3*i])
printf("error: length==4 \n");
r=(float)(c[3*j])-(float)(a[3*i]);
r=r/(float)(a[3*i]);
}
goto nskip;
}
}
for (j=0; j<(order/48); j++) {
if (d[3*j]==k) {
if ((d[3*j+2]!=2)&&(d[3*j+2]!=4))
printf("error D: length=%d \n",d[3*j+2]);
if ((count==2)&&(d[3*j+2]==2)) {
if (d[3*j]>=a[3*i])
printf("error: length==2 \n");
r=(float)(a[3*i])-(float)(d[3*j]);
r=r/(float)(a[3*i]);
}
if ((count==4)&&(d[3*j+2]==4)) {
if (d[3*j]<=a[3*i])
printf("error: length==4 \n");
r=(float)(d[3*j])-(float)(a[3*i]);
r=r/(float)(a[3*i]);
}
goto nskip;
}
}
printf("error: A[odd}, k=%d \n",k);
goto zskip;
}
}
printf("error: A[odd], k=%d \n",k);
goto zskip;
nskip:
k=0;
}
zskip:
fclose(Outfp);
return(0);
}
```