/*****************************************************************************/
/*									     */
/*  REGENERATE LOOPS							     */
/*  12/28/08 (dkc)							     */
/*									     */
/*  This C program regenerates loops for the 3n+c sequence given entry	     */
/*  points ("list.dat" is verified).  A check for 1-2 sequence vectors is    */
/*  made.  "jumping over" is verified to be a 1-2 sequence vector.           */
/*									     */
/*****************************************************************************/
#include <stdio.h>
#include <math.h>
int main () {
int c=125;	// c value
int iters=22;	 // number of entry points
//int s[1]={-68}; // c=1
//int s[7]={2,38,1562,1496,374,1436,1334}; // c=5
//int s[1]={10}; // c=7
//int s[4]={8,2,62,26}; // c=11
//int s[9]={4,454,502,574,682,844,1048,262,358}; // c=13
//int s[5]={140,62,50,8,2}; // c=17
//int s[2]={28,10}; // c=19
//int s[18]={14,20,1244,392,98,110,
//	     -23704,-5926,-12526,-19330,-8134,-14686,-12364,-11878,-7486,-9436,
//	     -9052,-7330}; // c=23
//int s[3]={34,88,22}; // c=25
//int s[19]={8,2,-226,98,44,788552,197138,194366,79082,73934,41606,18794,1707830,
//	     1441016,360254,577232,144308,36902,7622}; // c=29
//int s[5]={-380,106,58,52,34}; // c=31
//int s[3]={68,104,26}; // c=35
//int s[3]={46,58,76}; // c=37
//int s[5]={128,32,8,2,38}; // c=41
//int s[3]={28,16,4}; // c=43
//int s[12]={146,170,206,260,44,38,26,20,494,200,50,98}; // c=47
//int s[4]={814,400,100,178}; // c=49
//int s[3]={350,230,206}; // c=53
//int s[7]={76,634,244,82,64,16,4}; // c=55
//int s[14]={104,26,92,38,32,8,2,1064,266,740,434,596,578,362}; // c=59
//int s[8]={16,4,24958,5974,3850,2878,1258,586}; // c=61
//int s[5]={332,134,116,62,38}; // c=65
//int s[5]={1060,340,184,46,34}; // c=67
//int s[14]={116,248,62,41360,10340,17042,48488,12122,18692,7970,22472,5618,
//	     13874,22724}; // c=71
//int s[15]={58,40,10,22,202,166,94,1216,304,76,958,328,82,100,46}; // c=73
//int s[13]={-18448,-4612,-1978,-14254,-5326,-3874,284,128,32,8,2,74,20}; // c=77
//int s[20]={1864,466,1060,658,268,136,34,88,22,70,46,28,
//	     460,172,142,106,64,16,4,52}; // c=79
//int s[7]={1052,218,566,314,1646,638,260}; // c=83
//int s[13]={2194,2146,1174,844,826,736,184,46,364,214,142,112,28}; // c=85
//int s[3]={266,122,68}; // c=89
//int s[10]={400,100,118,1018,748,220,64,16,4,58}; // c=91
//int s[21]={236,68,182,92,140,50,-1870,-1150,-1426,2288,572,416,104,26,
//	     254,128,32,8,2,122,74}; // c=95
//int s[10]={52,34,106,64,16,4,-5516,-1010,-2222,-1250}; // c=97
//int s[13]={152,38,92,116,74,98,62,200,50,44,224,56,14}; // c=101
//int s[12]={160,40,10,70,52,286,214,184,46,124,106,58}; // c=103
//int s[17]={9104,2276,1862,554,512,128,32,8,2,230,194,176,44,86,80,20,
//	     68}; // c=107
//int s[5]={2170,658,274,130,76}; // c=109
//int s[12]={1064,266,128,32,8,2,50,-13678,-13138,-5668,-5098,-3076}; // c=113
//int s[12]={1036,196,124,52,964,442,286,136,34,46,100,-422}; // c=115
//int s[27]={80,20,56,14,44,38,1196,254,296,74,218,116,110,92,
//	     -622,4922,3698,3350,2978,2654,2162,1286,512,128,32,8,2}; // c=119
//int s[29]={628,196,148,124,58,52,40,10,34,790,670,562,520,130,304,76,238,
//	     202,184,46,106,70,-10244,-4178,-5042,-55928,-13982,-7166,
//	     -4340}; // c=121
  int s[22]={32,8,3218,1238,926,764,752,188,734,572,136328,34082,61538,42764,
	     42584,10646,29462,23108,12812,8852,4364,3002}; // c=125
//int s[15]={382,274,220,166,154,94,82,1030,418,622,598,256,64,16,4}; // c=127
//int s[37]={380,266,104,26,338,326,272,68,1694,668,542,236,218,158,110,92,74,
//	     50,-756598,-638224,-159556,-29884,-170044,-31750,-134566,
//	     -155020,-41182,-21766,-146716,-41182,-21766,-155020,-133108,
//	     -70342,-108628,-43018,-60694}; // c=131
//int s[13]={460,358,352,88,22,298,142,2728,682,1942,250,226,118}; // c=133
//int s[20]={80000,20000,5000,1250,13418,5066,2726,1934,3806,3374,2972,
//	     3212,1916,1178,776,194,692,476,206,164}; // c=137
//int s[58]={1510,916,868,862,706,538,490,484,358,220,202,142,88,22,76,
//	     -16472,-4118,-13556,-12098,-4502,-11612,-10316,-10154,-9182,
//	     -5078,-8876,-8858,-8492,-7994,-7886,-7418,-7238,-5942,-7034,
//	     -7022,-6446,-6374,-6062,-5798,-5414,-5366,-4982,-4694,-17624,
//	     -4406,-5294,-17240,-4310,-7724,-14546,-5420,-10328,-2582,
//	     -35102,-19658,-21602,-8066,-2990}; // c=139
//int s[26]={158,122,116,33224,8306,24476,18032,4508,7178,2384,596,1148,1022,
//	     698,482,392,98,350,224,56,14,536,134,86,68,44}; // c=143
//int s[20]={328,82,184,46,166,802,616,154,226,94,64,16,4,
//	     -1130,490,220,-5486,-3758,-1454,-482}; // c=145
//int s[12]={5798,4826,4412,1334,29450,8348,1868,608,152,38,194,110}; // c=149
//int s[14]={502,370,226,2854,1108,1018,856,214,730,538,406,190,118,82}; // c=151
//
int jump[48*2]={
   13,262, 17,2, 23,98, 25,22, 41,8, 43,4, 47,50, 49,100, 55,16,
   59,26, 59,8, 61,4, 67,46, 71,5618, 73,304, 77,8, 79,4,
   79,22, 85,46, 85,28, 91,16, 95,8, 95,26, 97,4, 101,50,
   101,14, 103,10, 103,46, 107,128, 107,2, 107,20, 113,8,
   119,32, 119,14, 121,10, 121,46, 125,2, 131,26, 133,22, 137,32,
   137,5000, 137,194, 139,-4310, 139,-2582, 143,596, 143,56, 145,16,
   149,152};
int dest[48]={358,5,110,28,2,1,98,25,4,92,2,1,34,32734,76,2,1,28,79,7,4,2,
   74,1,44,71,70,43,32,68,5,2,8,35,34,304,32,299,298,8,1250,107,-7724,-9911,
   149,14,4,38};
int i,k,n,max,order,temp,flag,flag1,flag2,flag3,endjmp;
unsigned int j,count,first;
FILE *Outfp;
Outfp = fopen("out0e.dat","w");
//
// compute order (of loop)
//
for (i=0; i<iters; i++) {
   k=s[i];
   flag1=0;
   for (j=0; j<48; j++) {
      if ((c==jump[2*j])&&(k==jump[2*j+1])) {
	 flag1=1;
	 endjmp=dest[j];
	 }
      }
   max=k;
   if (max<0)
      max=-max;
   while (k==(k/2)*2)
      k=k/2;
   for (j=1; j<10000; j++) {
      k=3*k+c;
      temp=k;
      if (temp<0)
	 temp=-temp;
      if (temp>max)
	 max=temp;
      while (k==(k/2)*2) {
	 if (k==s[i])
	    goto bskip;
	 k=k/2;
	 }
      }
   printf("error: s[i]=%d \n",s[i]);
   goto zskip;
bskip:
   order=3;
   while (order<max)
      order=order*2;
//
// find odd natural number divisible by 3
//
   k=s[i];
   max=k;
   if (max<0)
      max=-max;
   while (k!=(k/3)*3) {
      if (k==(k/2)*2) {
	 if ((k-c)==((k-c)/3)*3) {
	    k=(k-c)/3;
	    temp=k;
	    if (temp<0)
	       temp=-temp;
	    if (temp>max)
	       max=temp;
	    }
	 else {
	    k=k*2;
	    temp=k;
	    if (temp<0)
	       temp=-temp;
	    if (temp>max)
	       max=temp;
	    }
	 }
      else {
	 k=k*2;
	 temp=k;
	 if (temp<0)
	    temp=-temp;
	 if (temp>max)
	    max=temp;
	 }
      }
//
// include even natural numbers to the left of the odd natural number divisible
// by 3
//
   temp=k;
   if (temp<0)
      temp=-temp;
   while (temp<(order/2)) {
      temp=temp*2;
      k=k*2;
      }
//
// compute sequence
//
   count=1;
   first=1;
   while (k==(k/2)*2) {
      k=k/2;
      count=count+1;
      }
   flag2=0;
   if (flag1==1) {
      printf(" %d",k);
      fprintf(Outfp," %d",k);
      flag3=1;
      }
   flag=1;
   for (j=1; j<10000; j++) {
      k=3*k+c;
      count=count+1;
      n=1;
      if (flag1==1) {
	 if (flag3==1) {
	    printf(" %d",k);
	    fprintf(Outfp," %d",k);
	    }
	 if (k==endjmp) {
	    flag2=1;
	    flag3=0;
	    }
	 }
      while (k==(k/2)*2) {
	 if (k==s[i]) {
	    flag=0;
	    if (first==1)
	       first=0;
	    else
	       goto askip;
	    }
	 k=k/2;
	 count=count+1;
	 n=n+1;
	 if (flag1==1) {
	    if (flag3==1) {
	       printf(" %d",k);
	       fprintf(Outfp," %d",k);
	       }
	    if (k==endjmp) {
	       flag2=1;
	       flag3=0;
	       }
	    }
	 if ((n>3)&&(flag==1)) {
	    printf("error: not a 1-2 sequence vector, c=%d, s=%d \n",c,s[i]);
	    fprintf(Outfp,"error: not a 1-2 sequence vector, c=%d, s=%d \n",c,s[i]);
	    goto zskip;
	    }
	 }
      }
   printf("error \n");
   fprintf(Outfp,"error \n");
askip:
   if (flag1==1) {
      printf("\n");
      fprintf(Outfp,"\n");
      if (flag2==0)
	 printf("error:  end jump=%d \n",endjmp);
      }
   if (flag1==1)
      printf("order=%d, count=%d, c=%d, d=%d \n",order,count-1,c,s[i]);
   else
      printf("order=%d, count=%d, c=%d \n",order,count-1,c);
   fprintf(Outfp,"order=%d, count=%d \n",order,count-1);
   if (max>order) {
      printf("  error:  order=%d, maximum=%d \n",order,max);
      fprintf(Outfp,"  error:  order=%d, maximum=%d \n",order,max);
      }
   }
zskip:
fclose(Outfp);
return(0);
}