﻿ /*******************************
```/*****************************************************************************/
/*									     */
/*  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[1]=1;
for (i=3; i<6000003; i+=6) {   // odd n values divisible by 3
iters=iters+1;	       // increment number of starting values
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];
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
}
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");
flag=0;		       // clear k iteration count
count=1;		       // initialize k value
for (k=1; k<m-1; k++) {
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[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);
}
```