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