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