/******************************************************************************/
/*									      */
/*  COMPUTE 3N+C SEQUENCE (tree for k<=24, more efficient algorithm)	      */
/*  02/14/09 (dkc)							      */
/*									      */
/*  This C program histograms the lengths of limbs in A, B, C, and D.  There  */
/*  are outer "c" and "order" loops.  The permutation of A, B, C, and D for   */
/*  different c values is checked; a limb in B is verified to attach to a     */
/*  limb in A or C and a limb in D is verified to attach to a limb in B or D. */
/*  Also, a limb in A[odd] is verified to attach to a 2-element or 4-element  */
/*  limb.  The lengths of the limbs are checked.  The actual number of limbs  */
/*  of a given length is compared to the expected number of limbs for that    */
/*  length.								      */
/*									      */
/******************************************************************************/
#include <stdio.h>
#include <math.h>
int main () {
int c,k;
unsigned int shift,order;
unsigned int g,h,i,j,l,m,n,t,u,mask,count,offset;  // indices
unsigned int s[131072];      // duplicate array, minimum size=(order/12)/32
unsigned int ho[512];			    // histogram of lengths
unsigned int dist[200]; 		    // histogram of z values
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
unsigned int z[10*2]={1,2,	    // 1  (fractions)
		      1,3,	    // 2
		      1,12,	    // 4
		      1,72,	    // 5
		      1,36,	    // 7
		      5,576,	    // 9
		      7,864,	    // 10
		      77,7776,	    // 12
		      13,15552,     // 14
		      815,186624};   // 15
unsigned int v[10]={1,2,4,5,7,9,10,12,14,15};
unsigned int d[4],e,f,b,flag,flag0,flag1,flag2;
double r;
FILE *Outfp;
Outfp = fopen("out9.dat","w");
//
// begin outer loop
//
for (c=-1; c<103; c+=2) {
   if (c==(c/3)*3)
      continue;
//
// begin middle loop
//
   for (shift=2; shift<31; shift++) {
      if ((c==-1)&&(shift==3))	    // hung up in a loop
	 continue;
      if ((c==5)&&(shift==3))
	 continue;
      if ((c==5)&&(shift==6))
	 continue;
      if ((c==7)&&(shift==4))
	 continue;
      if ((c==11)&&(shift==4))
	 continue;
      if ((c==13)&&(shift==5))
	 continue;
      if ((c==13)&&(shift==10))
	 continue;
      if ((c==17)&&(shift==5))
	 continue;
      if ((c==19)&&(shift==5))
	 continue;
      if ((c==23)&&(shift==5))
	 continue;
      if ((c==25)&&(shift==6))
	 continue;
      if ((c==25)&&(shift==8))
	 continue;
      if ((c==29)&&(shift==6))
	 continue;
      if ((c==31)&&(shift==6))
	 continue;
      if ((c==35)&&(shift==6))
	 continue;
      if ((c==35)&&(shift==9))
	 continue;
      if ((c==37)&&(shift==6))
	 continue;
      if ((c==41)&&(shift==6))
	 continue;
      if ((c==43)&&(shift==6))
	 continue;
      if ((c==47)&&(shift==6))
	 continue;
      if ((c==47)&&(shift==8))
	 continue;
      if ((c==49)&&(shift==7))
	 continue;
      if ((c==53)&&(shift==7))
	 continue;
      if ((c==55)&&(shift==7))
	 continue;
      if ((c==55)&&(shift==9))
	 continue;
      if ((c==59)&&(shift==7))
	 continue;
      if ((c==59)&&(shift==10))
	 continue;
      if ((c==61)&&(shift==7))
	 continue;
      if ((c==65)&&(shift==7))
	 continue;
      if ((c==65)&&(shift==9))
	 continue;
      if ((c==65)&&(shift==12))
	 continue;
      if ((c==67)&&(shift==7))
	 continue;
      if ((c==71)&&(shift==7))
	 continue;
      if ((c==73)&&(shift==7))
	 continue;
      if ((c==77)&&(shift==7))
	 continue;
      if ((c==79)&&(shift==7))
	 continue;
      if ((c==83)&&(shift==7))
	 continue;
      if ((c==85)&&(shift==7))
	 continue;
      if ((c==85)&&(shift==10))
	 continue;
      if ((c==89)&&(shift==7))
	 continue;
      if ((c==91)&&(shift==7))
	 continue;
      if ((c==91)&&(shift==12))
	 continue;
      if ((c==91)&&(shift==13))
	 continue;
      if ((c==95)&&(shift==7))
	 continue;
      if ((c==95)&&(shift==10))
	 continue;
      if ((c==97)&&(shift==8))
	 continue;
      if ((c==101)&&(shift==8))
	 continue;
      order=(1<<shift)*3;    // 3*2**k
      if (order>50331648) {
	 break;
	 }
      printf("order=%d c=%d \n",order,c);
      if (c==-1)
	 e=7;
      else
	 e=c-(c/8)*8;
      if (e==1) {
	 d[0]=0;
	 d[1]=1;
	 d[2]=2;
	 d[3]=3;
	 }
      if (e==3) {
	 d[0]=1;
	 d[1]=0;
	 d[2]=3;
	 d[3]=2;
	 }
      if (e==5) {
	 d[0]=2;
	 d[1]=3;
	 d[2]=0;
	 d[3]=1;
	 }
      if (e==7) {
	 d[0]=3;
	 d[1]=2;
	 d[2]=1;
	 d[3]=0;
	 }
      if (c==-1)
	 e=11;
      else
	 e=c-(c/12)*12;
      if ((e==1)||(e==7))
	 offset=8;
      if ((e==5)||(e==11))
	 offset=4;
      for (i=0; i<200; i++)
	 dist[i]=0;
      for (i=0; i<512; i++)		       // clear histogram of lengths
	 ho[i]=0;
      for (i=0; i<131072; i++)		       // clear duplicate array
	 s[i]=0;
//
// limbs in A, B, C, and D
//
g=0;					 // clear counter
count=0;				 // clear counter
for (i=order/2; i<order; i+=2) {	 // possible starting elements
   if ((i-offset)==((i-offset)/12)*12)	 // check for elements of U
      continue;
   k=i; 				 // back-track
   if (k!=(k/3)*3) {			 // check for dead limb
      if ((k-c)!=((k-c)/3)*3) { 	 // check for beginning of limb
	 goto askip;
	 }
      k=(k-c)/3;			 // previous number in sequence
      if (k==(k/3)*3)			 // check for dead limb
	 goto bskip;
      k=k*2;				 // previous number in sequence
      while ((k<(int)(order/2))||((k-c)==((k-c)/3)*3)) { // check for beginning
	 if ((k-c)==((k-c)/3)*3) {
	    k=(k-c)/3;			 // previous number in sequence
	    if (k==(k/3)*3)		 // check for dead limb
	       goto bskip;
	    k=k*2;			 // previous number in sequence
	    }
	 else
	    k=k*2;			 // previous number in sequence
	 }
      }
askip:
   m=k; 				 // save beginning of limb
   n=1; 				 // set length
   while (k==(k/2)*2) { 		 // check for even element
      k=k/2;				 // next element of limb
      n=n+1;				 // increment length
      }
   for (j=1; j<1000000; j++) {		 // iterate until end of limb
      h=3*k+c;				 // next element of limb
      if (h>order) {			 // check for end of limb
	 if (k<(int)(order/3))
	    goto bskip;
	 t=((k-(order/3))-1)/2; 	 // index into array
	 u=t-(t/32)*32; 		 // fractional portion of index
	 t=t/32;			 // integer portion of index
	 mask=1<<u;			 // set mask
	 if ((s[t]&mask)==0)		 // check if bit not set
	    s[t]=s[t]|mask;		 // set bit in array
	 else				 // already set
	    goto bskip;
	 g=g+1; 			 // increment counter
	 e=k-(k/8)*8;
	 e=(e-1)/2;
	 e=d[e];
	 if ((h/2)==(h/4)*2) {
	    f=8;
	    if (e==0) {
	       flag1=0;
	       if ((k/8)==(k/16)*2)
		  flag1=1;
	       flag2=0;
	       if ((c!=-1)&&(c!=11)&&(c!=13)&&(c!=25)&&(c!=29)&&(c!=31)&&(c!=41)&&(c!=43)&&(c!=47)&&(c!=59)&&(c!=61)&&(c!=73)&&(c!=77)&&(c!=79)&&(c!=89)&&(c!=91)&&(c!=95)) {
		  if (flag1==0)
		     flag2=1;
		  }
	       else {
		  if (flag1==1)
		     flag2=1;
		  }
	       if (flag2!=0) {
		  b=h/4;
		  if (b==(b/2)*2) {
		     printf("error: even second element \n");
		     goto zskip;
		     }
		  if (b<(order/3)) {
		     b=3*b+c;
		     b=b/2;
		     if (b==(b/2)*2) {
			printf("error: even fourth element \n");
			goto zskip;
			}
		     if (b<(order/3)) {
			printf("error: not an element of S \n");
			goto zskip;
			}
		     }
		  b=b-(b/8)*8;
		  b=(b-1)/2;
		  f=d[b];
		  }
	       }
	    if (e==2) {
	       if (((h/2)-offset)==(((h/2)-offset)/12)*12)
		  f=6;
	       else {
		  printf("error: C does not attach to U \n");
		  printf("  %d (%d,%d(%d))=>%d(%d) \n",n,m,k,e,h/2,f);
		  goto zskip;
		  }
	       }
	    }
	 else {
	    if ((h/2)>(2*order/3))
	       f=9;
	    else {
	       f=(h/2)-(h/16)*8;
	       f=(f-1)/2;
	       f=d[f];
	       if (e==1) {
		  if ((f!=0)&&(f!=2)) {
		     printf("error: incorrect attachment point \n");
		     goto zskip;
		     }
		  }
	       if (e==3) {
		  if ((f!=1)&&(f!=3)) {
		     printf("error: incorrect attachment point \n");
		     goto zskip;
		     }
		  }
	       }
	    }
//
// check for valid limb lengths
//
	 flag0=0;
	 for (l=0; l<40; l++) {
	    if (n==y[l]) {
	       flag0=1;
	       break;
	       }
	    }
	 if (flag0!=0) {
	    if (m==(m/3)*3) {		 // check for dead limb
	       fprintf(Outfp,"D");
	       fprintf(Outfp," %d (%d,%d(%d))=>%d(%d), c=%d, order=%d \n",n,m,k,e,h/2,f,c,order);
	       printf("unexpected limb length: D");
	       printf(" %d (%d,%d(%d))=>%d(%d), c=%d, order=%d \n",n,m,k,e,h/2,f,c,order);
	       }
	    else{
	       fprintf(Outfp,"  %d (%d,%d(%d))=>%d(%d), c=%d, order=%d \n",n,m,k,e,h/2,f,c,order);
	       printf("unexpected limb length:    %d (%d,%d(%d))=>%d(%d), c=%d, order=%d \n",n,m,k,e,h/2,f,c,order);
	       }
	    }
	 if (n<512)			 // check length
	   ho[n]=ho[n]+1;		 // histogram length
	 else {
	    printf("error: histogram array not big enough \n");
	    goto zskip;
	    }
	 goto bskip;
	 }
      k=h;				 // next element
      if (k==(k/8)*8)			 // not valid limb, start over
	 goto bskip;			 // (prevents taking even path at nodes)
      n=n+1;				 // increment length
      while (k==(k/2)*2) {		 // check for even element
	 k=k/2; 			 // next element
	 n=n+1; 			 // increment length
	 }
      }
bskip:
   n=0;
   }
      if (g!=(order/12))
	 fprintf(Outfp,"error:  incorrect number of entries, delta=%d, c=%d, order=%d \n",g-(order/12),c,order);
//
// check probabilities of occurrence of limb lengths
//
      g=order/6;
      for (i=1; i<10; i++) {
	 r=(double)g;
	 r=r*(double)z[2*i];
	 r=r/(double)z[2*i+1];
	 j=(unsigned int)r;
	 j=j+1;
	 k=(int)j-(int)ho[v[i]];
	 flag=0;
	 if ((i<4)&&((k<0)||(k>1)))
	    flag=1;
	 if ((i==4)&&((k<-1)||(k>2)))
	    flag=1;
	 if ((i==5)&&((k<-1)||(k>2)))
	    flag=1;
	 if ((i==6)&&((k<-1)||(k>2)))
	    flag=1;
	 if ((i==7)&&((k<-2)||(k>3)))
	    flag=1;
	 if ((i==8)&&((k<-2)||(k>3)))
	    flag=1;
	 if ((i==9)&&((k<-2)||(k>3)))
	    flag=1;
	 if (flag!=0) {
	    if ((c==1)||(c==-1)) {
	       printf("error: unexpected number of limbs \n");
	       printf("j=%d, h=%d, v=%d, i=%d, r=%f \n",j,ho[v[i]],v[i],i,r);
	       for (j=0; j<20; j++)
		  printf(" %d %d \n",j,ho[j]);
	       goto zskip;
	       }
	    else
	       printf("warning: unexpected number of limbs, delta=%d \n",k);
	    }
	 }
      }
   }
zskip:
fclose(Outfp);
return(0);
}