/******************************************************************************/
/* */
/* COMPUTE 3N+1 OR 3N-1 SEQUENCE (tree for k<=24, more efficient algorithm) */
/* 01/29/10 (dkc) */
/* */
/* This C program shows that |n|(order/2) is much greater than |c|e for */
/* orders from 3*2 to 3*2**24. */
/* */
/******************************************************************************/
#include <stdio.h>
#include <math.h>
int main () {
//
int c=1; // c
unsigned int l=23; // specified limb length
unsigned int p=75; // count
unsigned int initshift=1; // initial shift
//
unsigned int e,f,h,i,j,m,n,t,u,mask,offset,shift,order; // indices
int k,temp1;
unsigned int s[131072]; // duplicate array, minimum size=(order/12)/32
unsigned int diste[200];
int num[16]={1,-1,7,5,-17,47,13,-217,295,-139,1909,1631,-3299,
13085,6487,-46075};
int den[16]={4,8,16,32,64,128,256,512,1024,2048,4096,8192,16384,
32768,65536,131072};
unsigned int len[16]={2,4,5,7,9,10,12,14,15,17,18,20,22,23,25,27};
//
//int z[1]={2}; // length=2
//int z[1]={10}; // length=4
//int z[1]={20}; // length=5
//int z[2]={76,58}; // length=7
//int z[2]={260,206}; // length=9
//int z[3]={520,412,340}; // length=10
//int z[6]={1688,986,1148,1004,842,1364}; // length=12
//int z[7]={5320,3700,2782,3214,3268,4348}; // length=14
//int z[10]={10640,7400,5960,6428,4988,6536,4340,8696,5564,4916}; // length=15
//int z[20]={ // length=17
// 32944,23224,11290,20308,12586,17752,12748,11434,14836,14314,10138,
// 12892,18904,15988,11596,20632,14044,27112,17716,15772};
//int z[20]={ // length=18
// 65888,46448,22580,40616,25172,35504,25496,22868,29672,28628,20276,
// 25784,37808,31976,23192,41264,28088,54224,35432,31544};
//int z[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 z[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 z[75]={ // length=23
205948,223444,224380,226684,241876,244180,245116,246772,249688,262612,265204,
267508,268120,270424,272764,273016,285940,288856,290260,291448,293752,298612,
308008,312184,313588,316504,317044,324856,326440,328744,339832,343288,344692,
347176,359848,360496,370936,374824,378280,378928,381232,386164,399664,405928,
406504,412336,412408,424936,427312,430768,452584,458416,458992,447400,477424,
494056,499888,505072,528976,546544,547408,556264,575056,608752,616528,633952,
652384,678736,680032,721504,772048,783712,877024,1016992,1226944};
/* int z[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 z[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}; */
//
int r,del;
double temp;
FILE *Outfp;
Outfp = fopen("out7.dat","w");
if ((c!=1)&&(c!=-1)) {
printf("invalid c value \n");
goto zskip;
}
if (l>27) {
printf("limb length too large \n");
goto zskip;
}
if (c==1) {
if ((l==4)||(l==9)||(l==14)||(l==17)||(l==22)||(l==27)) {
printf("input error: c=1 but limit is negative \n");
goto zskip;
}
}
else {
if ((l!=4)&&(l!=9)&&(l!=14)&&(l!=17)&&(l!=22)&&(l!=27)) {
printf("input error: c=-1 but limit is positive \n");
goto zskip;
}
}
if (c==1)
offset=8;
else
offset=4;
for (i=0; i<200; i++) // clear distribution of e values
diste[i]=0;
//
// begin order loop
//
for (shift=initshift; shift<=24; shift++) {
if ((c==-1)&&(shift==3)) // hung up in a loop
continue;
order=(1<<shift)*3;
if (order>50331648)
goto zskip;
fprintf(Outfp,"ORDER=%d c=%d \n",order,c);
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;
r=(int)(m)-(int)(h/2); // difference of first elements
for (f=0; f<16; f++) {
if (n==len[f]) {
temp=(double)(m)*num[f]-(double)(r)*den[f];
del=(int)temp;
for (e=0; e<p; e++) {
if (del==c*z[e]) {
temp=(double)num[f];
if (temp<0.0)
temp=-temp;
temp1=c;
if (temp1<0)
temp1=-temp1;
temp=temp*(order/2)-temp1*z[e];
if (temp<0.0) {
printf("error: negative value \n");
goto zskip;
}
if (diste[e]==0)
printf("order=%d, e=%d, delta=%f \n",order,z[e],temp);
diste[e]=diste[e]+1;
goto bskip;
}
}
printf("error: No solution of Diophantine equation found.\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;
}
}
zskip:
fclose(Outfp);
return(0);
}