﻿ /*******************************
```/******************************************************************************/
/*									      */
/*  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 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
}
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
}
}
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
if ((s[t]&mask)==0)		 // check if bit not set
s[t]=s[t]|mask;		 // set bit in array
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);
}
```