/******************************************************************************/
/*									      */
/*  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
	    goto askip;
	 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
	    }
	 }
      printf("error A: h=%d \n",h);        // 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 (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
	 }
      }
   printf("error \n");                       // limb not found
   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);
}