/*****************************************************************************/
/*									     */
/*  COMPUTE 3N+1 SEQUENCE (probabilities)				     */
/*  12/26/08 (dkc)							     */
/*									     */
/*  This C program determines the number of jumps required to reach an even  */
/*  natural number.							     */
/*									     */
/*****************************************************************************/
#include <stdio.h>
#include <math.h>
unsigned int lmbd(unsigned int mode, unsigned int a);
void mul64_32(unsigned int a, unsigned int b, unsigned int c, unsigned int *p);
void div64_32(unsigned int *a, unsigned int *b, unsigned int);
void add64(unsigned int *a, unsigned int *b);
void sub64(unsigned int *c, unsigned int *d);
int main () {
unsigned int i,j,k,m,n,out[2000],flag,count,total,iters;
unsigned int V[2],W[2],X[3],Y[2],g[1000],h[1000];
double r;
FILE *Outfp;
Outfp = fopen("out1.dat","w");
for (i=0; i<1000; i++) {
   h[i]=0;		       // clear histogram of k values
   g[i]=0;		       // clear histogram of iteration values
   }
iters=0;		       // number of starting values
total=0;		       // number of k values
W[0]=0; 		       // load 1
W[1]=1;
for (i=3; i<6000003; i+=6) {   // odd n values divisible by 3
   iters=iters+1;	       // increment number of starting values
   V[0]=0;		       // load n
   V[1]=i;
   m=0; 		       // clear output array pointer
   out[0]=V[0]; 	       // save n value
   out[1]=V[1];
   m=1; 		       // increment output array pointer
   n=1; 		       // number of terms in sequence (not used)
   for (j=1; j<100000000; j++) {
      mul64_32(V[0], V[1], 3, X);
      V[0]=X[1];
      V[1]=X[2];
      add64(W, V);	       // 3n+1
      if (m<1000) {
	 out[2*m]=V[0];        // save n value
	 out[2*m+1]=V[1];
	 m=m+1; 	       // increment output array pointer
	 }
      else {
	 printf("error: output array not big enough \n");
	 goto zskip;
	 }
      n=n+1;		       // increment number of terms in sequence
      div64_32(V, Y, 8);       // (3n+1)/8
      mul64_32(Y[0], Y[1], 8, X);
      Y[0]=X[1];	       // [(3n+1)/8]*8
      Y[1]=X[2];
      sub64(V, Y);
      if ((Y[0]==0)&&(Y[1]==0)) {  // 8 divides 3n+1
	 div64_32(V, V, 2);    // n=n/2
	 if (m<1000) {
	    out[2*m]=V[0];     // save n value
	    out[2*m+1]=V[1];
	    m=m+1;	       // increment output array pointer
	    }
	 n=n+1; 	       // increment number of terms in sequence
	 div64_32(V, V, 2);    // n=n/2
	 if (m<1000) {
	    out[2*m]=V[0];     // save n value
	    out[2*m+1]=V[1];
	    m=m+1;	       // increment output array pointer
	    }
	 else {
	    printf("error: output array not big enough \n");
	    goto zskip;
	    }
	 n=n+1; 	       // increment number of terms in sequence
	 goto askip;
	 }
      div64_32(V, Y, 2);       // n/2
      mul64_32(Y[0], Y[1], 2, X);  // (n/2)*2
      Y[0]=X[1];
      Y[1]=X[2];
      sub64(V, Y);	       // n-(n/2)*2
      while ((Y[0]==0)&&(Y[1]==0)) {  // 2 divides n
	 div64_32(V, V, 2);    // n=n/2
	 if (m<1000) {
	    out[2*m]=V[0];     // save n value
	    out[2*m+1]=V[1];
	    m=m+1;	       // increment output array pointer
	    }
	 else {
	    printf("error: output array not big enough \n");
	    goto zskip;
	    }
	 n=n+1; 	       // increment number of terms in sequence
	 div64_32(V, Y, 2);    // n/2
	 mul64_32(Y[0], Y[1], 2, X);  // (n/2)*2
	 Y[0]=X[1];
	 Y[1]=X[2];
	 sub64(V, Y);	       // n-(n/2)*2
	 }
      }
   fprintf(Outfp,"error \n");
askip:
   flag=0;		       // clear k iteration count
   count=1;		       // initialize k value
   for (k=1; k<m-1; k++) {
      V[0]=out[2*k-2];	       // load n'
      V[1]=out[2*k-1];
      div64_32(V, Y, 2);       // n'/2
      mul64_32(Y[0], Y[1], 2, X);  // (n'/2)*2
      Y[0]=X[1];
      Y[1]=X[2];
      sub64(V, Y);	       //  n'-(n'/2)*2
      if ((Y[0]==0)&&(Y[1]==0)) {   // 2 divides n'
	 V[0]=out[2*k];        // load n
	 V[1]=out[2*k+1];
	 div64_32(V, Y, 2);    // n/2
	 mul64_32(Y[0], Y[1], 2, X);  // (n/2)*2
	 Y[0]=X[1];
	 Y[1]=X[2];
	 sub64(V, Y);	       // n-(n/2)*2
	 if ((Y[0]==0)&&(Y[1]==0)) {  // 2 divides n
	    count=count/2;     // k=k/2
//	    fprintf(Outfp," %d \n",count);
	    h[count]=h[count]+1;  // histogram k value
	    count=0;	       // clear k value
	    flag=flag+1;       // increment k iteration count
	    }
	 }
      count=count+1;	       // k=k+1
      }
// fprintf(Outfp,"\n");
   g[flag]=g[flag]+1;	       // histogram k iteration count
   total=total+flag;	       // increment total number of k values
   }
fprintf(Outfp,"histogram of k values, total=%d \n",total);
for (i=0; i<100; i++) {
   r=(double)h[i];
   r=r/(double)total;
   fprintf(Outfp," %d     %f \n",h[i],r);
   }
fprintf(Outfp,"\n");
fprintf(Outfp,"histogram of iteration values, total=%d \n",iters);
for (i=0; i<100; i++) {
   r=(double)g[i];
   r=r/(double)iters;
   fprintf(Outfp," %d     %f \n",g[i],r);
   }
zskip:
fclose(Outfp);
return(0);
}