/******************************************************************************/
/* */
/* COMPUTE 3N+C SEQUENCE (tree for k<=24, more efficient algorithm) */
/* 02/14/09 (dkc) */
/* */
/* This C program finds the minimum order for which a given limb length */
/* occurs. The sequence vector of the limb is computed. The maximum */
/* and minimum x/y ratios are computed for every order. */
/* */
/******************************************************************************/
#include <stdio.h>
#include <math.h>
int main () {
//
int c=-1; // c
unsigned int initshift=1; // initial shift
//
int m,k,del,r,zsave,z[162];
unsigned int a,b,d,e,f,g,h,i,j,l,n,t,u,mask,flag,offset,shift,order,count,first;
unsigned int asave,bsave,p,maxset;
unsigned int s[131072]; // duplicate array, minimum size=(order/12)/32
unsigned int sv[2048]; // sequence vector
unsigned int dist[200]; // histogram of z values
double delta,x,maxx,minx,oldmax,oldmin,rat;
int num[11]={1,-1,5,-17,13,-217,-139,1631,-3299,6487,-46075};
int den[11]={4,8,32,64,256,512,2048,8192,16384,65536,131072};
unsigned int len[11]={2,4,7,9,12,14,17,20,22,25,27};
unsigned int pi[11]={1,1,2,2,6,7,20,40,58,151,162};
//
int z2[1]={2}; // length=2
int z4[1]={10}; // length=4
int z7[2]={76,58}; // length=7
int z9[2]={260,206}; // length=9
int z12[6]={1688,986,1148,1004,842,1364}; // length=12
int z14[7]={5320,3700,2782,3214,3268,4348}; // length=14
int z17[20]={ // length=17
32944,23224,11290,20308,12586,17752,12748,11434,14836,14314,10138,
12892,18904,15988,11596,20632,14044,27112,17716,15772};
int z20[40]={ // length=20
61630,54718,50110,49534,44926,69064,59740,93112,68092,59092,
81448,55132,117520,88504,68488,54484,76840,71836,127888,63880,
125944,110608,80584,53908,64924,79612,89980,106000,201760,49300,
60316,100024,73672,66004,88360,166768,110392,72700,98728,143440};
int z22[58]={ // length=22
120938,131306,132602,134060,142970,144428,145130,156092,156794,
158252,158522,169916,172346,193082,199832,202964,206204,212468,
226292,238712,247028,249944,252536,278132,304376,308264,326192,
340016,360752,386024,438512,508496,
145724,163220,171644,173588,185468,187412,189140,189464,213656,
215384,223700,229208,273272,273704,287528,339368,391856,613472};
int z25[151]={ // length=25
527726,673454,574382,728750,523118,569774,636590,564590,619886,811694,
611246,666542,486254,532910,492734,488126,2383904,1374968,1062632,693704,
532100,675956,1269776,609140,1757936,891740,1157240,956792,1325072,656444,
807176,740360,1215416,2068976,1012088,711740,878456,1619696,773480,1019000,
1314704,1674992,1074296,1114256,1372880,736220,765308,1145576,1169552,1035920,
567092,820604,1112312,529598,698312,1672400,1408016,1167608,835292,613748,
3713600,584894,983900,1582832,1547984,982136,899336,915320,1307576,614972,
1077392,1859024,914024,3083744,728444,661628,1882352,969320,1075448,1176464,
2072864,798428,731612,1231760,451262,578108,903548,703100,1191260,758396,
877160,1514936,810344,624764,1250552,973532,1409744,773084,1139600,1072784,
1465040,828380,694748,603956,970472,659252,735176,525620,851060,650612,
1066844,790472,656840,2663840,982280,1390520,705908,1232912,781832,572276,
1934624,1269992,919928,837128,619580,1989920,703496,1701560,1052264,1532432,
851816,527492,666236,2348912,907112,890588,744968,1897760,568964,712820,
928604,624260,814952,490628,768116,2197280,844040,1252280,1007336,562484};
int z27[162]={ // length=27
11206336,9316768,8057056,7217248,7112272,6657376,6284128,6272464,6035296,
5869408,5758816,5712592,5642608,5339344,5170216,5090512,5082736,4924624,
4814032,4709488,4662832,4610344,4460656,4294768,4289584,4237096,4190440,
4184176,4040752,4009648,3988264,3875512,3874864,3822376,3817192,3764272,
3760816,3711784,3639316,3594928,3574192,3568360,3537256,3502264,3484336,
3408304,3402472,3297712,3291880,3288424,3283888,3266068,3253432,3222328,
3173296,3122536,3101800,3087544,3017236,3012376,3011944,2986132,2976952,
2973496,2935912,2851348,2825320,2811496,2807608,2786872,2776180,2763544,
2740756,2737300,2700904,2697016,2620984,2618716,2597656,2576920,2571412,
2550676,2527348,2510392,2500618,2496568,2487064,2460820,2436952,2411032,
2385976,2384788,2369884,2361460,2340724,2300440,2286616,2274196,2271064,
2260372,2251786,2250868,2203996,2200756,2183260,2176024,2174836,2160472,
2149780,2146648,2093404,2085898,2065162,2064244,2053336,2050420,2043292,
2036056,2034868,2017372,1975306,1942744,1939828,1938316,1925194,1924276,
1910452,1906780,1899274,1892956,1877404,1820218,1817140,1799860,1788682,
1782364,1774858,1772428,1766812,1759306,1752988,1706548,1664266,1661836,
1659676,1654330,1648714,1648012,1642396,1634890,1554700,1549084,1543738,
1541578,1537420,1529914,1524298,1444108,1436602,1430986,1419322,1326010};
FILE *Outfp;
Outfp = fopen("out12.dat","w");
if (c<-1) {
printf("c too small \n");
goto zskip;
}
if (c>101) {
printf("c too big \n");
goto zskip;
}
//
// begin limb length loop
//
for (g=0; g<11; g++) {
rat=(double)num[g];
rat=rat/(double)den[g];
for (i=0; i<200; i++)
dist[i]=0;
l=len[g];
p=pi[g];
if (l==2) {
for (i=0; i<p; i++)
z[i]=z2[i];
}
if (l==4) {
for (i=0; i<p; i++)
z[i]=z4[i];
}
if (l==7) {
for (i=0; i<p; i++)
z[i]=z7[i];
}
if (l==9) {
for (i=0; i<p; i++)
z[i]=z9[i];
}
if (l==12) {
for (i=0; i<p; i++)
z[i]=z12[i];
}
if (l==14) {
for (i=0; i<p; i++)
z[i]=z14[i];
}
if (l==17) {
for (i=0; i<p; i++)
z[i]=z17[i];
}
if (l==20) {
for (i=0; i<p; i++)
z[i]=z20[i];
}
if (l==22) {
for (i=0; i<p; i++)
z[i]=z22[i];
}
if (l==25) {
for (i=0; i<p; i++)
z[i]=z25[i];
}
if (l==27) {
for (i=0; i<p; i++)
z[i]=z27[i];
}
flag=0;
oldmax=-1.0;
oldmin=-1.0;
//
// begin order loop
//
for (shift=initshift; shift<25; shift++) {
maxset=0;
maxx=-1.0;
minx=1000000;
order=(1<<shift)*3;
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;
if (order>50331648) {
printf("order too big \n");
goto zskip;
}
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<131072; i++) // clear duplicate array
s[i]=0;
//
// limbs in A, B, C, and D
//
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;
if (n!=l) // continue if not specified length
goto bskip;
if (m!=(m/3)*3) {
r=(int)(m)-(int)(h/2); // difference of first elements
x=(double)r;
if (x<0.0)
x=-x;
x=x/(double)(m);
if (x>maxx)
maxx=x;
if (x<minx)
minx=x;
maxset=1;
for (f=0; f<16; f++) {
if (n==len[f]) {
del=(int)(m)*num[f]-r*den[f];
zsave=0;
for (e=0; e<p; e++) {
if (del==c*z[e]) {
zsave=z[e];
if (dist[e]==0) {
if ((num[f]<0)&&(c>0))
goto eskip;
if ((num[f]>0)&&(c<0))
goto eskip;
delta=(double)num[f];
delta=delta*(double)(order/2);
delta=delta-(double)(c*z[e]);
fprintf(Outfp,"order=%d, length=%d, delta=%e, e=%d \n",order,l,delta,z[e]);
}
eskip: dist[e]=dist[e]+1;
goto fskip;
}
}
printf("error: No solution of Diophantine equation. \n");
goto zskip;
}
}
//
// COMPUTE SEQUENCE VECTOR
//
fskip: if (l==2) {
sv[0]=1;
f=1;
goto dskip;
}
f=0;
if ((m&1)==0) {
count=0;
while ((m&1)==0) {
count=count+1;
m=m/2;
}
sv[0]=count;
f=1;
}
first=1;
while (m!=k) {
if ((m&1)==0) {
m=m/2;
count=count+1;
}
else {
m=3*m+c;
if (first==0) {
sv[f]=count;
f=f+1;
if (f>2047) {
printf("error: sequence vector array not big enough \n");
goto zskip;
}
}
else
first=0;
count=0;
}
}
sv[f]=count;
f=f+1;
dskip: a=0;
b=0;
for (d=0; d<f; d++) {
if (sv[d]==1)
a=a+1;
else
b=b+1;
}
if (flag==0) {
asave=a;
bsave=b;
}
if ((a!=asave)||(b!=bsave)) {
printf("error: a=%d, b=%d \n",a,b);
goto zskip;
}
if (l==4) {
if (f!=2) {
printf("error: l=4, f=%d \n",f);
goto zskip;
}
if ((sv[0]!=1)||(sv[1]!=1)) {
printf("error: l=4, sv[1]=%d \n",sv[1]);
goto zskip;
}
}
if (l==7) {
if (f!=3) {
printf("error: l=7, f=%d \n",f);
goto zskip;
}
if ((sv[0]!=1)||(sv[1]!=2)||(sv[2]!=1)) {
printf("error: l=7, sv[1]=%d, sv[2]=%d \n",sv[1],sv[2]);
goto zskip;
}
}
if (l==9) {
if (f!=4) {
printf("error: l=9, f=%d \n",f);
goto zskip;
}
if ((sv[0]!=1)||(sv[1]!=2)||(sv[2]!=1)||(sv[3]!=1)) {
printf("error: l=9, sv[1]=%d, sv[2]=%d, sv[3]=%d \n",sv[1],sv[2],sv[3]);
goto zskip;
}
}
if (l==12) {
if (f!=5) {
printf("error: l=12, f=%d \n",f);
goto zskip;
}
if ((sv[0]!=1)||(sv[1]!=2)||(sv[4]!=1)) {
printf("error: l=12, sv[1]=%d, sv[2]=%d, sv[3]=%d, sv[4]=%d \n",sv[1],sv[2],sv[3],sv[4]);
goto zskip;
}
if ((sv[2]==1)&&(sv[3]==1)) {
printf("error: l=12, sv[1]=%d, sv[2]=%d, sv[3]=%d, sv[4]=%d \n",sv[1],sv[2],sv[3],sv[4]);
goto zskip;
}
if ((sv[2]==2)&&(sv[3]==2)) {
printf("error: l=12, sv[1]=%d, sv[2]=%d, sv[3]=%d, sv[4]=%d \n",sv[1],sv[2],sv[3],sv[4]);
goto zskip;
}
}
if (l==14) {
if (f!=6) {
printf("error: l=14, f=%d \n",f);
goto zskip;
}
if ((sv[0]!=1)||(sv[1]!=2)||(sv[4]!=1)||(sv[5]!=1)) {
printf("error: l=14, sv[1]=%d, sv[2]=%d, sv[3]=%d, sv[4]=%d, sv[5]=%d \n",sv[1],sv[2],sv[3],sv[4],sv[5]);
goto zskip;
}
if ((sv[2]==1)&&(sv[3]==1)) {
printf("error: l=14, sv[1]=%d, sv[2]=%d, sv[3]=%d, sv[4]=%d, sv[5]=%d \n",sv[1],sv[2],sv[3],sv[4],sv[5]);
goto zskip;
}
if ((sv[2]==2)&&(sv[3]==2)) {
printf("error: l=14, sv[1]=%d, sv[2]=%d, sv[3]=%d, sv[4]=%d, sv[5]=%d \n",sv[1],sv[2],sv[3],sv[4],sv[5]);
goto zskip;
}
}
if (l==17) {
if (f!=7) {
printf("error: l=17, f=%d \n",f);
goto zskip;
}
if ((sv[0]!=1)||(sv[1]!=2)||(sv[6]!=1)) {
printf("error: l=17, sv[1]=%d, sv[2]=%d, sv[3]=%d, sv[4]=%d, sv[5]=%d, sv[6]=%d \n",sv[1],sv[2],sv[3],sv[4],sv[5],sv[6]);
goto zskip;
}
if ((sv[4]==2)&&(sv[5]==2)) {
printf("error: l=17, sv[1]=%d, sv[2]=%d, sv[3]=%d, sv[4]=%d, sv[5]=%d, sv[6]=%d \n",sv[1],sv[2],sv[3],sv[4],sv[5],sv[6]);
goto zskip;
}
}
if (l==20) {
if (f!=8) {
printf("error: l=20, f=%d \n",f);
goto zskip;
}
if ((sv[0]!=1)||(sv[1]!=2)||(sv[2]!=2)||(sv[7]!=1)) {
printf("error: l=20, sv[1]=%d, sv[2]=%d, sv[3]=%d, sv[4]=%d, sv[5]=%d, sv[6]=%d, sv[7]=%d \n",sv[1],sv[2],sv[3],sv[4],sv[5],sv[6],sv[7]);
goto zskip;
}
if ((sv[5]==2)&&(sv[6]==2)) {
printf("error: l=20, sv[1]=%d, sv[2]=%d, sv[3]=%d, sv[4]=%d, sv[5]=%d, sv[6]=%d, sv[7]=%d \n",sv[1],sv[2],sv[3],sv[4],sv[5],sv[6],sv[7]);
goto zskip;
}
if ((sv[3]==2)&&(sv[4]==2)&&(sv[5]==2)) {
printf("error: l=20, sv[1]=%d, sv[2]=%d, sv[3]=%d, sv[4]=%d, sv[5]=%d, sv[6]=%d, sv[7]=%d \n",sv[1],sv[2],sv[3],sv[4],sv[5],sv[6],sv[7]);
goto zskip;
}
if (sv[2]==1) {
printf("error: l=20, sv[1]=%d, sv[2]=%d, sv[3]=%d, sv[4]=%d, sv[5]=%d, sv[6]=%d, sv[7]=%d \n",sv[1],sv[2],sv[3],sv[4],sv[5],sv[6],sv[7]);
goto zskip;
}
}
if (l==22) {
if (f!=9) {
printf("error: l=22, f=%d \n",f);
goto zskip;
}
if ((sv[0]!=1)||(sv[1]!=2)||(sv[7]!=1)||(sv[8]!=1)) {
printf("error: l=22, sv[1]=%d, sv[2]=%d, sv[3]=%d, sv[4]=%d, sv[5]=%d, sv[6]=%d, sv[7]=%d, sv[8]=%d \n",sv[1],sv[2],sv[3],sv[4],sv[5],sv[6],sv[7],sv[8]);
goto zskip;
}
if ((sv[4]==2)&&(sv[5]==2)&&(sv[6]==2)) {
printf("error: l=22, sv[1]=%d, sv[2]=%d, sv[3]=%d, sv[4]=%d, sv[5]=%d, sv[6]=%d, sv[7]=%d, sv[8]=%d \n",sv[1],sv[2],sv[3],sv[4],sv[5],sv[6],sv[7],sv[8]);
goto zskip;
}
if ((sv[2]==1)&&(sv[3]==1)) {
printf("error: l=22, sv[1]=%d, sv[2]=%d, sv[3]=%d, sv[4]=%d, sv[5]=%d, sv[6]=%d, sv[7]=%d, sv[8]=%d \n",sv[1],sv[2],sv[3],sv[4],sv[5],sv[6],sv[7],sv[8]);
goto zskip;
}
}
if (l==25) {
if (f!=10) {
printf("error: l=25, f=%d \n",f);
goto zskip;
}
if ((sv[0]!=1)||(sv[1]!=2)||(sv[9]!=1)) {
printf("error: l=25, sv[1]=%d, sv[2]=%d, sv[3]=%d, sv[4]=%d, sv[5]=%d, sv[6]=%d, sv[7]=%d, sv[8]=%d, sv[9]=%d \n",sv[1],sv[2],sv[3],sv[4],sv[5],sv[6],sv[7],sv[8],sv[9]);
goto zskip;
}
if ((sv[2]==1)&&(sv[3]==1)) {
printf("error: l=25, sv[1]=%d, sv[2]=%d, sv[3]=%d, sv[4]=%d, sv[5]=%d, sv[6]=%d, sv[7]=%d, sv[8]=%d, sv[9]=%d \n",sv[1],sv[2],sv[3],sv[4],sv[5],sv[6],sv[7],sv[8],sv[9]);
goto zskip;
}
if ((sv[7]==2)&&(sv[8]==2)) {
printf("error: l=25, sv[1]=%d, sv[2]=%d, sv[3]=%d, sv[4]=%d, sv[5]=%d, sv[6]=%d, sv[7]=%d, sv[8]=%d, sv[9]=%d \n",sv[1],sv[2],sv[3],sv[4],sv[5],sv[6],sv[7],sv[8],sv[9]);
goto zskip;
}
if ((sv[6]==2)&&(sv[7]==2)&&(sv[8]==2)) {
printf("error: l=25, sv[1]=%d, sv[2]=%d, sv[3]=%d, sv[4]=%d, sv[5]=%d, sv[6]=%d, sv[7]=%d, sv[8]=%d, sv[9]=%d \n",sv[1],sv[2],sv[3],sv[4],sv[5],sv[6],sv[7],sv[8],sv[9]);
goto zskip;
}
if ((sv[2]==1)&&(sv[3]==2)&&(sv[4]==1)&&(sv[5]==2)&&(sv[6]==2)&&(sv[7]==1)&&(sv[8]==2))
printf("error: l=25, sv[1]=%d, sv[8]=%d \n",sv[1],sv[8]);
if ((sv[2]==1)&&(sv[3]==2)&&(sv[4]==1)&&(sv[5]==2)&&(sv[6]==2)&&(sv[7]==2)&&(sv[8]==1))
printf("error: l=25, sv[1]=%d, sv[7]=%d \n",sv[1],sv[7]);
if ((sv[2]==2)&&(sv[3]==1)&&(sv[4]==1)&&(sv[5]==2)&&(sv[6]==2)&&(sv[7]==1)&&(sv[8]==2))
printf("error: l=25, sv[2]=%d, sv[8]=%d \n",sv[2],sv[8]);
if ((sv[2]==2)&&(sv[3]==1)&&(sv[4]==1)&&(sv[5]==2)&&(sv[6]==2)&&(sv[7]==2)&&(sv[8]==1))
printf("error: l=25, sv[2]=%d, sv[7]=%d \n",sv[2],sv[7]);
}
flag=flag+1;
}
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 (maxset!=0) {
printf(" l=%d, o=%d, max=%e, min=%e, r=%e \n",l,order,maxx,minx,rat);
if (((num[g]>0)&&(c>0))||((num[g]<0)&&(c<0))) {
if (oldmax>maxx) {
printf("error: not monotonic, old maximum=%e \n",oldmax);
goto zskip;
}
if (oldmin>minx) {
printf("error: not monotonic, old minimum=%e \n",oldmin);
goto zskip;
}
}
oldmax=maxx;
oldmin=minx;
}
}
fprintf(Outfp,"\n");
}
zskip:
fclose(Outfp);
return(0);
}