/*****************************************************************************/
/*									     */
/*  COMPUTE 3N+C SEQUENCE						     */
/*  02/22/13 (dkc)							     */
/*									     */
/*  This C program computes 3n+c sequences starting with odd integers	     */
/*  divisible by 3.  When an attachment point is encountered, the sequence   */
/*  is extended backwards ("odd" paths are taken) to another odd integer     */
/*  divisible by 3.  The differences in the absolute values of the odd	     */
/*  integers divisible by 3 are histogrammed.				     */
/*									     */
/*****************************************************************************/
#include <stdio.h>
#include <math.h>
void add64(unsigned int *a, unsigned int *b);
void sub64(unsigned int *a, unsigned int *b);
void div64_32(unsigned int *A, unsigned int *B, unsigned int C);
int main () {
unsigned int c=1;	    // c value
unsigned int limit=10;   // 10, 100, 1000, 10000, 100000, 1000000, 10000000
			    // When "limit" is less than or equal to 100, the
			    // individual t1 and t2 values are output.
unsigned int relprime=1;    // If set, exclude (t2,c)!=1.
int scale=1536; 	    // Divide |t2|-|t1| by this amount.
			    // This scaling factor is for "limit" equal to
			    // 10000.
int offset=550; 	    // This offset (to make the histogram bins
			    // non-negative) is for "limit" equal to 10000.
unsigned int j,k,m,n,out[2000],flag,count,total,iters,kount,toobig,even1,outrng;
unsigned int L[2],S[2],U[2],V[2],W[2],Z[2],g[1000],h[1000],T[2],K[2],Q[2];
unsigned int histo[10000];
unsigned int THI[2],M[2];
int i,f,temp,thi,diff;
double r;
FILE *Outfp;
Outfp = fopen("out1c.dat","w");
if (c>11) {
   printf("c too large \n");
   goto zskip;
   }
for (i=0; i<10000; i++)       // clear histogram of |t2|-|t1| values
   histo[i]=0;
Z[0]=0;
Z[1]=0;
for (i=0; i<1000; i++) {
   h[i]=0;		       // clear histogram of k values
   g[i]=0;		       // clear histogram of iteration values
   }
outrng=0;		       // big sequence value counter
even1=0;		       // number of failures when t2 attaches to the
			       // first even element after t1
THI[0]=0;		       // set maximum t2 value
THI[1]=0;
toobig=0;		       // number of values not histogrammed
kount=0;		       // number of failures (|t2|<|t1|)
iters=0;		       // number of starting values
total=0;		       // number of k values
W[0]=0; 		       // load c
W[1]=c;
for (f=-(int)limit; f<=(int)limit; f++) {
//for (f=0; f<=(int)limit; f++) {  // natural numbers
   i=3+f*6;		       // t2, odd n values divisible by 3
   if (c!=1) {
      if (relprime!=0) {
	 if (i==(i/(int)c)*(int)c)
	    continue;
	 }
      }
   iters=iters+1;	       // increment number of starting values
   if (i>0) {
      V[0]=0;		       // load n
      V[1]=i;
      }
   else {
      V[0]=0xffffffff;
      V[1]=(unsigned int)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++) {
      T[0]=V[0];
      T[1]=V[1];
      add64(V,V);
      add64(T,V);
      add64(W, V);	       // 3n+c
      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
      if ((V[1]&7)==0) {
	 V[1]=V[1]>>1;
	 V[1]=V[1]|(V[0]<<31);
	 V[0]=(int)V[0]>>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
	 V[1]=V[1]>>1;
	 V[1]=V[1]|(V[0]<<31);
	 V[0]=(int)V[0]>>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
	 goto askip;
	 }
      while ((V[1]&1)==0) {
	 V[1]=V[1]>>1;
	 V[1]=V[1]|(V[0]<<31);
	 V[0]=(int)V[0]>>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
	 }
      }
   fprintf(Outfp,"error \n");
   goto zskip;
//
// attachment point found, extend sequence backwards using 0-1 sequence vector,
//
//    k=s[i];
//    while (k!=(k/3)*3) {
//	 if (k==(k/2)*2) {
//	    if ((k-c)==((k-c)/3)*3)
//	       k=(k-c)/3;
//	    else
//	       k=k*2;
//	    }
//	 else
//	    k=k*2;
//	 }
askip:
      K[0]=out[2*m-2];
      K[1]=out[2*m-3];
      S[0]=K[0];
      S[1]=K[1];
//
// do adjustment so that the |t2|>|t1| test will pass for cycles
//
      if (c==1) {
	 if ((i!=3)&&(i!=-27)&&(i!=-15)&&(i!=-9)&&(i!=-3)) {
	    K[1]=K[1]>>1;
	    K[1]=K[1]|(K[0]<<31);
	    K[0]=(int)K[0]>>1;
	    }
	 }
      if (c==5) {
	 if ((i!=51)&&(i!=63)&&(i!=153)&&(i!=-45)&&(i!=-15)&&(i!=15)&&(i!=-135)) {
	    K[1]=K[1]>>1;
	    K[1]=K[1]|(K[0]<<31);
	    K[0]=(int)K[0]>>1;
	    }
	 }
      if (c==7) {
	 if ((i!=-21)&&(i!=-189)&&(i!=-63)&&(i!=21)) {
	    K[1]=K[1]>>1;
	    K[1]=K[1]|(K[0]<<31);
	    K[0]=(int)K[0]>>1;
	    }
	 }
      if (c==11) {
	 if ((i!=-33)&&(i!=-27)&&(i!=33)&&(i!=-297)&&(i!=-105)&&(i!=-99)&&(i!=-81)) {
	    K[1]=K[1]>>1;
	    K[1]=K[1]|(K[0]<<31);
	    K[0]=(int)K[0]>>1;
	    }
	 }
//
// find odd integer divisible by 3
//
      if ((K[0]&0x80000000)!=0) {
	 T[0]=K[0];
	 T[1]=K[1];
	 sub64(Z,T);
	 div64_32(T,Q,3);
	 sub64(Z,Q);
	 }
      else
	 div64_32(K,Q,3);
      T[0]=Q[0];
      T[1]=Q[1];
      add64(Q,Q);
      add64(T,Q);
      sub64(K,Q);
      while ((Q[0]!=0)||(Q[1]!=0)) {
	 if ((K[1]&1)==0) {
	    U[0]=W[0];
	    U[1]=W[1];
	    sub64(K,U);
	    if ((U[0]&0x80000000)!=0) {
	       T[0]=U[0];
	       T[1]=U[1];
	       sub64(Z,T);
	       div64_32(T,Q,3);
	       sub64(Z,Q);
	       }
	    else
	       div64_32(U,Q,3);
	    T[0]=Q[0];
	    T[1]=Q[1];
	    add64(T,T);
	    add64(Q,T);
	    sub64(U,T);
	    if ((T[0]==0)&&(T[1]==0)) {
	       K[0]=Q[0];
	       K[1]=Q[1];
	       }
	    else
	       add64(K,K);
	    }
	 else
	    add64(K,K);
	 if ((K[0]&0x80000000)!=0) {
	    T[0]=K[0];
	    T[1]=K[1];
	    sub64(Z,T);
	    div64_32(T,Q,3);
	    sub64(Z,Q);
	    }
	 else
	    div64_32(K,Q,3);
	 T[0]=Q[0];
	 T[1]=Q[1];
	 add64(Q,Q);
	 add64(T,Q);
	 sub64(K,Q);
	 }
//
// histogram scaled |t2|-|t1| values
//
   T[0]=K[0];
   T[1]=K[1];
   if ((T[0]&0x80000000)!=0)
      sub64(Z,T);
   if ((T[0]!=0)||((T[1]&0x80000000)!=0))
      outrng=outrng+1;
   L[0]=T[0];
   L[1]=T[1];
   sub64(THI,L);
   if ((L[0]&0x80000000)!=0) {
      THI[0]=T[0];
      THI[1]=T[1];
      }
   temp=i;
   if (temp<0)
      temp=-temp;
   if (limit==10000) {
      diff=(temp-(int)T[1])/scale;
      diff=diff+offset;
      if (diff<0) {
//	 printf("offset too small, delta=%d \n",diff);
	 toobig=toobig+1;
	 goto ajump;
	 }
      if (diff>9999) {
	 printf("array not big enough, temp=%d \n",diff);
	 goto zskip;
	 }
      histo[diff]=histo[diff]+1;
      }
//
// compare t values
//
ajump:
   if ((T[1]>(unsigned int)temp)||(T[0]!=0)) {
//
      L[0]=K[0];
      L[1]=K[1];
      add64(W, L);
      n=0;			 // clear count
      while ((L[1]&1)==0) {
	 L[1]=(L[1]>>1)|(L[0]<<31);
	 L[0]=(int)L[0]>>1;
	 n=n+1;
	 }
      L[0]=K[0];
      L[1]=K[1];
      for (j=0; j<n; j++) {
	 M[0]=L[0];
	 M[1]=L[1];
	 add64(L, L);
	 add64(M, L);
	 add64(W, L);
	 if (limit<=100)
	    printf("even element in first jump from t1=%d \n",L[1]);
	 M[0]=L[0];
	 M[1]=L[1];
	 add64(M, M);
	 if ((M[0]==S[0])&&(M[1]==S[1])) {
	    even1=even1+1;
	    if (limit<=100) {
	       printf("attaches to jth even: j=%d, S=%d %d, L=%d %d \n",j,S[0],S[1],K[0],K[1]);
	       }
	    }
	 L[1]=(L[1]>>1)|(L[0]<<31);
	 L[0]=(int)L[0]>>1;
	 }
      kount=kount+1;
      if (limit<=100) {
	 printf("failure: t2=%d, t1=%d \n",i,K[1]);
	 fprintf(Outfp,"failure: t2=%d, t1=%d \n",i,K[1]);
	 }
      }
   if (limit<=100) {
      printf("f=%d, attacher=%d, k=%d, %d, attached to=%d, %d \n",f,i,out[2*m-2],out[2*m-3],K[0],K[1]);
      fprintf(Outfp,"f=%d, attacher=%d, k=%d, %d, attached to=%d, %d \n",f,i,out[2*m-2],out[2*m-3],K[0],K[1]);
      }
//
// find number of odd elements in each jump to the first attachment point
//
   flag=0;		       // clear k iteration count
   count=1;		       // initialize k value
   for (k=0; k<(m-1); k++) {
      V[0]=out[2*k];	       // load n'
      V[1]=out[2*k+1];
      if ((V[1]&1)==0) {
	 V[0]=out[2*k+2];      // load n
	 V[1]=out[2*k+3];
	 if ((V[1]&1)==0) {
	    count=count/2;     // k=k/2
	    if (count>999) {
	       printf("array not big enough \n");
	       goto zskip;
	       }
	    h[count]=h[count]+1;  // histogram k value
	    if (limit<100)
	       printf("count=%d \n",count);
	    count=0;	       // clear k value
	    flag=flag+1;       // increment k iteration count
	    }
	 }
      count=count+1;	       // k=k+1
      }
   if (flag>999) {
      printf("array not big enough \n");
      goto zskip;
      }
   g[flag]=g[flag]+1;	       // histogram k iteration count
   total=total+flag;	       // increment total number of k values
   }
printf("\n");
printf("histogram of number of odd elements in jumps to attachment point, total=%d \n",total);
fprintf(Outfp,"histogram of number of odd elements in jumps to attachment point, total=%d \n",total);
for (i=1; i<20; i++) {
   r=(double)h[i];
   r=r/(double)total;
   printf("i=%d, h[i]=%d, %f \n",i,h[i],r);
   fprintf(Outfp," %d     %f \n",h[i],r);
   }
printf("\n");
fprintf(Outfp,"\n");
printf("histogram of number of jumps to attachment point, total=%d \n",iters);
fprintf(Outfp,"histogram of jumps to attachment point, total=%d \n",iters);
for (i=2; i<20; i++) {
   r=(double)g[i];
   r=r/(double)iters;
   printf("i=%d, g[i]=%d, %f \n",i-1,g[i],r);
   fprintf(Outfp," %d     %f \n",g[i],r);
   }
printf("\n");
fprintf(Outfp,"\n");
if (((THI[0]==0)&&(THI[1]&0x80000000)==0)) {
   thi=(int)THI[1];
   printf("largest t2 value=%d, portion of t1 range=%e \n",thi,(double)thi/(double)(6*limit+3));
   fprintf(Outfp,"largest t2 value=%d, portion of t1 range=%e \n",thi,(double)thi/(double)(6*limit+3));
   }
else {
   printf("largest t2 value=%#010x %#010x \n",THI[0],THI[1]);
   fprintf(Outfp,"largest t2 value=%#010x %#010x \n",THI[0],THI[1]);
   }
printf("failure count=%d, ratio=%f \n",kount,(double)(iters-kount)/(double)iters);
fprintf(Outfp,"failure count=%d, ratio=%f \n",kount,(double)(iters-kount)/(double)iters);
printf("failures due to attachments to first even=%d \n",even1);
fprintf(Outfp,"failures due to attachments to first even=%d \n",even1);
printf("big sequence values=%d \n",outrng);
fprintf(Outfp,"big sequence values=%d \n",outrng);
//
if (limit==10000) {
   printf("\n");
   fprintf(Outfp,"\n");
   printf("histogram of scaled differences \n");
   fprintf(Outfp,"histogram of scaled differences \n");
   for (i=0; i<600; i++) {
      printf("i=%d, n=%d \n",i-offset,histo[i]);
      fprintf(Outfp,"i=%d, n=%d \n",i-offset,histo[i]);
      }
   printf("number of values not histogrammed=%d \n",toobig);
   fprintf(Outfp,"number of values not histogrammed=%d \n",toobig);
   }
zskip:
fclose(Outfp);
return(0);
}