/******************************************************************************/
/* */
/* 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);
}