/*****************************************************************************/
/*									     */
/*  FIND LOOPS IN 3N+C SEQUENCE 					     */
/*  05/14/09 (dkc)							     */
/*									     */
/*  This C program finds loops in the 3n+c sequence.  Connection points are  */
/*  output.								     */
/*									     */
/*****************************************************************************/
#include <math.h>
#include <stdio.h>
int main () {
unsigned int e,f,g,h,j,m,n,bsave;
int limit,temp,max,i,k,save[32],v[65536];
int out[2000*2];
int c=-163;			  // c value
unsigned int iters=2;		 // number of "jumps" from the initial value
				 // (an odd natural number divisible by 3)
FILE *Outfp;
Outfp = fopen("out0a.dat","w");
limit=1000000;
max=1<<29;
m=0;
for (i=3; i<limit; i+=6) {
   k=(int)i;			 // initial sequence value
   for (e=0; e<iters; e++) {
      k=k+c;
      n=0;			 // clear count
      while (k==(k/2)*2) {
	 k=k/2;
	 n=n+1;
	 }
      temp=k;
      if (temp<0)
	 temp=-temp;
      for (j=0; j<n; j++) {
	 if (temp>max)		 // check for overflow
	    goto askip;
	 temp=temp*3;
	 k=k*3;
	 }
      k=(k-c)/2;		 // sequence value after a "jump"
      if (k==(k/2)*2)
	 goto bskip;
      }
   temp=k;
   if (temp<0)
      temp=-temp;
   if (temp>max)		 // check for overflow
      goto askip;
   k=3*k+c;
bskip:
   h=0; 			 // clear index
   n=0; 			 // clear index
   while (k==(k/2)*2) { 	 // find next odd sequence value
      v[n]=k;			 // save sequence value
      n=n+1;			 // increment index
      if (h<32) {		 // check index
	 save[h]=k;		 // save even sequence value
	 h=h+1; 		 // increment index
	 }
      else {			 // array too small
	 printf("error: array too small \n");
	 goto zskip;
	 }
      k=k/2;
      }
   bsave=k;			 // save odd sequence value
   for (j=1; j<10000000; j++) {
      v[n]=k;			 // save sequence value
      n=n+1;			 // increment index
      if (n>65535)		 // array full
	 goto askip;
      temp=k;
      if (temp<0)
	 temp=-temp;
      if (temp>max)
	 goto askip;
      k=3*k+c;
      while (k==(k/2)*2) {
	 for (g=0; g<h; g++) {
	    if (k==save[g]) {
	       for (f=g; f<n; f++) {
		  if (v[f]==(v[f]/8)*8) {
		     for (e=0; e<m; e++) {
			if (((int)(n-g)==out[2*e])&&((v[f]/4)==out[2*e+1])) {
			   goto askip;
			   }
			}
		     printf("n=%d  k=%d  i=%d \n",n-g,v[f]/4,i);
		     fprintf(Outfp,"n=%d  k=%d \n",n-g,v[f]/4);
		     out[2*m]=n-g;
		     out[2*m+1]=v[f]/4;
		     m=m+1;
		     if (m>=2000) {
			printf("output array full \n");
			goto zskip;
			}
		     goto askip;
		     }
		  }
		  for (e=0; e<m; e++) {
		     if (((int)(n-g)==out[2*e])&&((int)bsave==out[2*e+1])) {
			goto askip;
			}
		     }
		  printf("n=%d  k=%d \n",n-g,bsave);
		  fprintf(Outfp,"n=%d  k=%d \n",n-g,bsave);
		  out[2*m]=n-g;
		  out[2*m+1]=(int)bsave;
		  m=m+1;
		  if (m>=2000) {
		     printf("output array full \n");
		     goto zskip;
		     }
		  goto askip;
	       }
	    }
	 v[n]=k;
	 n=n+1;
	 if (n>65535)
	    goto askip;
	 k=k/2;
	 }
      }
askip:
   n=0;
   }
zskip:
fclose(Outfp);
return(0);
}