/******************************************************************************/
/*									      */
/*  COMPUTE 3N+1 OR 3N-1 SEQUENCE					      */
/*  02/09/11 (dkc)							      */
/*									      */
/*  This C program computes general 3n+1 or 3n-1 sequences. The fourth-to-last*/
/*  element of a sub-sequence must be odd and the third-to-last element of    */
/*  a sub-sequence must be divisible by 8.  The fourth-to-last element must   */
/*  be larger than the order divided by 6.  Certain lengths of such sequences */
/*  are not permissible.						      */
/*									      */
/*  The number of even and odd elements in the segment are computed (where the*/
/*  element after an odd element i is defined to be (3i+c)/2).		      */
/*									      */
/******************************************************************************/
#include <stdio.h>
#include <math.h>
int main () {
//
int c=-1;			       // c
unsigned int iters=5;	  // segment count
unsigned int iters1=3;	  // segment count (normally set to "iters" value)
//
unsigned int v[8192],itero,itere;
unsigned int a,b,f,g,h,i,j,k,m,n,np,order,oddk,count;	    // indices
unsigned int he[2048],ho[2048*4];		  // histogram
unsigned int y[178]={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
		    104,106,109,112,114,  // L
		    117,119,122,125,127,  // L
		    130,132,135,138,140,  // L
		    143,145,148,150,153,  // M
		    156,158,161,163,166,  // M
		    171,174,176,179,	  // M
		    184,187,189,192,194,  // I'
		    197,200,202,205,207,  // I'
		    210,212,215,218,220,  // J'
		    223,225,228,231,233,  // J'
		    236,238,241,243,246,  // K'
		    249,251,254,256,259,  // K'
		    262,264,267,269,272,  // K'
		    275,277,280,282,285,  // K'
		    287,290,293,295,298,  // L'
		    300,303,306,308,311,  // L'
		    313,316,318,321,324,  // M'
		    326,329,331,334,337,  // M'
		    339,342,344,347,349,  // N
		    352,355,357,360,362,  // N
		    365,368,370,373,375,  // N
		    378,380,383,386,388,  // O
		    391,393,396,399,401,  // O
		    404,406,409,412,414,  // O
		    417,419,422,424,427,  // P
		    430,432,435,437,440,  // P
		    443,445,448,450,453,  // P
		    458,461,463,466}; //
unsigned int z[75]={1,9,	 // I  (first elements of length pairs which
		    14,22,	 // I	when decremented by 3 give invalid
		    27,35,	 // I	lengths)
		    40,45,	 // J
		    53,58,	 // J
		    66,71,76,	 // K
		    84,89,	 // K
		    97,102,	 // K
		    107,115,	 // L
		    120,128,	 // L
		    133,141,	 // L
		    146,151,	 // M
		    159,164,	 // M
		    172,177,	 // M
		    182,190,	 // I'
		    195,203,	 // I'
		    208,213,	 // J'
		    221,226,	 // J'
		    234,239,244, // K'
		    252,257,	 // K'
		    265,270,	 // K'
		    278,283,	 // K'
		    288,296,	 // L'
		    301,309,	 // L'
		    314,319,	 // M'
		    327,332,	 // M'
		    340,345,350, // N
		    358,363,	 // N
		    371,376,	 // N
		    384,389,	 // O
		    394,402,	 // O
		    407,415,	 // O
		    420,428,	 // P
		    433,438,	 // P
		    446,451,	 // P
		    459,464};	 //
unsigned int table[121*2]={
   0, 0,      // 0
   4, 3,      // 1
   5, 4,      // 2   (6, 2) => (3, 2)
   5, 4,      // 3   (7, 3) => (5, 3)
   6, 5,      // 4
   6, 5,      // 5   (10, 5) => (8, 5)
   7, 6,      // 6
   8, 7,      // 7   (14, 7) => (11, 7)
   8, 7,      // 8
   9, 8,      // 9
   9, 8,      // 10
  10, 9,      // 11
  11, 10,     // 12  (22, 12) => (19, 12)
  11, 10,     // 13
  12, 11,     // 14
  12, 11,     // 15
  13, 12,     // 16
  13, 12,     // 17  (29, 17) => (27, 17)
  14, 13,     // 18
  15, 14,     // 19
  15, 14,     // 20
  16, 15,     // 21
  16, 15,     // 22
  17, 16,     // 23
  18, 17,     // 24
  18, 17,     // 25
  19, 18,     // 26
  19, 18,     // 27
  20, 19,     // 28
  20, 19,     // 29  (48, 29) => (46, 29)
  21, 20,     // 30
  22, 21,     // 31
  22, 21,     // 32
  23, 22,     // 33
  23, 22,     // 34
  24, 23,     // 35
  25, 24,     // 36
  25, 24,     // 37
  26, 25,     // 38
  26, 25,     // 39
  27, 26,     // 40
  27, 26,     // 41  (67, 41) => (65, 41)
  28, 27,     // 42
  29, 28,     // 43
  29, 28,     // 44
  30, 29,     // 45
  30, 29,     // 46
  31, 30,     // 47
  32, 31,     // 48
  32, 31,     // 49
  33, 32,     // 50
  33, 32,     // 51
  34, 33,     // 52
  34, 34,     // 53  (87, 53) => (84, 53)     Note:  Even count the same.
  35, 34,     // 54
  36, 35,     // 55
  36, 35,     // 56
  37, 36,     // 57
  37, 36,     // 58
  38, 37,     // 59
  39, 38,     // 60
  39, 38,     // 61
  40, 39,     // 62
  40, 39,     // 63
  41, 40,     // 64
  42, 41,     // 65
  42, 41,     // 66
  43, 42,     // 67
  43, 42,     // 68
  44, 43,     // 69
  44, 43,     // 70
  45, 44,     // 71
  46, 45,     // 72
  46, 45,     // 73
  47, 46,     // 74
  47, 46,     // 75
  48, 47,     // 76
  49, 48,     // 77
  49, 48,     // 78
  50, 49,     // 79
  50, 49,     // 80
  51, 50,     // 81
  51, 50,     // 82
  52, 51,     // 83
  53, 52,     // 84
  53, 52,     // 85
  54, 53,     // 86
  54, 53,     // 87
  55, 54,     // 88
  56, 55,     // 89
  56, 55,     // 90
  57, 56,     // 91
  57, 56,     // 92
  58, 57,     // 93
  58, 57,     // 94
  59, 58,     // 95
  60, 59,     // 96
  60, 59,     // 97
  61, 60,     // 98
  61, 60,     // 99
  62, 61,     // 100
  63, 62,     // 101
  63, 62,     // 102
  64, 63,     // 103
  64, 63,     // 104
  65, 64,     // 105
   0, 65,     // 106
  66, 65,     // 107
  67, 66,     // 108
  67, 66,     // 109
  68, 67,     // 110
  68, 67,     // 111
  69, 68,     // 112
  70, 69,     // 113
  70, 69,     // 114
  71, 70,     // 115
  71, 70,     // 116
  72, 71,     // 117
  72, 0,      // 118
  73, 72,     // 119
  74, 73};    // 120


FILE *Outfp;
Outfp = fopen("out2d.dat","w");
if (iters1>iters) {
   printf("Error: Second number of iterations too large. \n");
   goto zskip;
   }
for (i=0; i<2048; i++) {		    // clear histogram of lengths
   he[i]=0;				    // for limbs ending in E
   }
for (i=0; i<2048*4; i++) {		      // clear histogram of lengths
   ho[i]=0;
   }
//
// DEAD LIMBS (ending in an even natural number)
//
   count=0;
   for (i=99999999; i>3; i-=6) {
      g=0;
      k=i;			 // set k
      v[0]=k;
      n=1;
      order=3;			 // set order
      while (order<k)
	 order=order*2;
aloop:
      if (k>0x30000000) {
//	 printf("error: k=%d \n",k);
	 goto cskip;
	 }
      oddk=k;
      k=3*k+c;			 // next sequence value
      v[n]=k;
      n=n+1;			 // increment element count
      if (n>8191) {
	 printf("error: v array too small \n");
	 goto zskip;
	 }
      while (order<k)
	 order=order*2;
      if (k==(k/8)*8) { 	 // check for connection point
	 g=g+1;
	 if (g<iters)
	    goto dskip;
	 if (g>iters)
	    goto cskip;
	 np=n+1;		  // increment element count
	 m=i;
	 while (m<(order/2)) {
	    m=m*2;
	    np=np+1;
	    }
	 if (oddk>(order/6)) {
	    count=count+1;
	    itero=0;
	    itere=0;
	    f=0;
	    for (h=0; h<n; h++) {
	       if ((v[h]&1)==1)
		  itero=itero+1;
	       else {
		  if ((v[h-1]&1)==0)
		     itere=itere+1;
		  }
	       if (((v[h]&1)==1)&&(v[h+1]==(v[h+1]/8)*8))
		  f=f+1;
	       if (f==iters1) {
//		  printf(" order=%d h=%d, %d %d %d \n",order,h,v[0],v[h],np-n);
//		  printf(" order=%d n=%d, %d %d %d \n",order,n,v[0],v[n-1],np-n);
		  goto wskip;
		  }
	       }
wskip:	    itere=itere+(np-n)+1;
	    if (itero<121) {
	       a=table[2*itero];
	       b=table[2*itero+1];
	       if (iters==iters1) {
		  if ((itere!=a)&&(itere!=b)) {
		     printf("error: odd=%d, even=%d \n",itero,itere);
		     goto zskip;
		     }
		  }
	       if (itere<b) {
		  printf("error: odd=%d, even=%d \n",itero,itere);
		  goto zskip;
		  }
	       if (((itere-b)<16)&&(itero<511))
		  ho[16*itero+itere-b]=ho[16*itero+itere-b]+1;
	       }
	    if (count<400) {
	       printf("l=%d, n=%d \n",itere,itero);
	       fprintf(Outfp,"l=%d, n=%d \n",itere,itero);
	       fprintf(Outfp," order=%d n=%d, %d %d %d \n",order,n,v[0],v[n-1],np-n);
	       }
	    for (j=0; j<75; j++) {
	       if (z[j]==np)
		  goto bskip;
	       }
	    np=np-3;
bskip:	    if (np<2048)
	       he[np]=he[np]+1;
	    }
	 }
dskip:while (k==(k/2)*2) {
	 k=k/2; 		 // next sequence value
	 v[n]=k;
	 n=n+1; 		 // increment element count
	 if (n>8191) {
	    printf("error: v array too small \n");
	    goto zskip;
	    }
	 }
      if (c==1) {
	 if (k!=1)
	    goto aloop;
	 }
      if (c==-1) {
	 if ((k!=1)&&(k!=5)&&(k!=7)&&(k!=17)&&(k!=25)&&(k!=37)&&(k!=55)&&(k!=41)&&(k!=61)&&(k!=91))
	    goto aloop;
	 }
cskip:k=0;
      }
//
// check for valid lengths
//
for (i=0; i<178; i++) {
   if (he[y[i]]!=0) {
      printf("error: n=%d \n",y[i]);
      goto zskip;
      }
   }
//
// HISTOGRAM OF LENGTHS
//
fprintf(Outfp,"\n");
fprintf(Outfp,"HISTOGRAM \n");
printf("\n");
printf("HISTOGRAM OF LENGTHS \n");
for (h=0; h<64; h++) {
   fprintf(Outfp," %d %d %d %d %d %d %d %d \n",he[8*h],he[8*h+1],he[8*h+2],
	   he[8*h+3],he[8*h+4],he[8*h+5],he[8*h+6],he[8*h+7]);
   printf(" %d %d %d %d %d %d %d %d \n",he[8*h],he[8*h+1],he[8*h+2],
	   he[8*h+3],he[8*h+4],he[8*h+5],he[8*h+6],he[8*h+7]);
   }
//
// HISTOGRAM OF LENGTHS
//
fprintf(Outfp,"\n");
fprintf(Outfp,"HISTOGRAM \n");
printf("\n");
printf("HISTOGRAM OF NUMBER OF EVEN ELEMENTS (PER NUMBER OF ODD ELEMENTS)\n");
for (h=0; h<64; h++) {
   fprintf(Outfp," %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d \n",ho[16*h],
	   ho[16*h+1],ho[16*h+2],ho[16*h+3],ho[16*h+4],ho[16*h+5],ho[16*h+6],
	   ho[16*h+7],ho[16*h+8],ho[16*h+9],ho[16*h+10],ho[16*h+11],ho[16*h+12],
	   ho[16*h+13],ho[16*h+14],ho[16*h+15]);
   printf(" %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d \n",ho[16*h],
	   ho[16*h+1],ho[16*h+2],ho[16*h+3],ho[16*h+4],ho[16*h+5],ho[16*h+6],
	   ho[16*h+7],ho[16*h+8],ho[16*h+9],ho[16*h+10],ho[16*h+11],ho[16*h+12],
	   ho[16*h+13],ho[16*h+14],ho[16*h+15]);
   }
zskip:
fclose(Outfp);
return(0);
}