/*  FIND CYCLES IN THE 3N+C SEQUENCE (cycles without attachment points)      */
/*  03/20/13 (dkc)							     */
/*									     */
/*  This C program finds cycles in the 3n+c sequence.  Odd elements are      */
/*  output.  Double-precision (64-bit) words are used.	Robert Floyd's       */
/*  cycle finding algorithm (Knuth 1969, "The Art of Computer Programming,   */
/*  Vol. 1: Fundamental Algorithms", pp. 4-7, exercise 7) is used.           */
/*									     */
/*  The output array "sinp" contains the odd elements.  The values in        */
/*  the output array "cflag" indicate which odd elements are associated      */
/*  with a cycle.  The output array "size" gives the total number of         */
/*  odd elements for each c value.  The output array "cval" gives the        */
/*  corresponding c values.  The output variable "numbc" gives the number    */
/*  of c values.  After some minor editing to fill in array sizes, the	     */
/*  output can be made into an "include" file for other programs.            */
/*									     */
/*  Output "big" cycles indicate whether parameters (in particular "limit")  */
/*  need to be modified to find cycles that may have been missed.  "L" is    */
/*  is the number of even elements in the cycles and "K" is the number of    */
/*  odd elements.  "count" gives the number of cycles found for a given c    */
/*  value.								     */
/*									     */
/*****************************************************************************/
#include <stdio.h>
#include <math.h>
#include "table.h"
unsigned int carry(unsigned int a, unsigned int b, unsigned int sum);
unsigned int lmbd(unsigned int mode, unsigned int a);
void add64(unsigned int *a, unsigned int *b);
void sub64(unsigned int *a, unsigned int *b);
void div64_32(unsigned int *dividend, unsigned int *quotient,
	      unsigned int divisor);
void mul64_32(unsigned int a0, unsigned int a2, unsigned int m0,
	      unsigned int *product);
int main () {
//
unsigned int cstart=10601;		 // beginning c value
unsigned int cend=10799;		 // ending c value
unsigned int mode=0;	    // double-word flag.  When the mode is set to 0,
			    // only the second word of odd elements is
			    // output (to save memory).
unsigned int iters=3;	    // number of "jumps" from the initial odd integer
			    // divisible by 3.	The parameters "iters", "limit",
			    // and "max" have to be tuned to find all cycles in
			    // a reasonable amount of time.  Increasing "iters"
			    // usually decreases execution time.  When "iters"
			    // is increased, "max" must usually be increased
			    // (up to at most about 0x4000000).
unsigned int limit=18000000; // must be a multiple of 6, should be as large as
			    // possible relative to "max" to allow for cycles
			    // having a constrained dynamic range.  Increasing
			    // "limit" increases the execution time.  Cycles with
			    // (K+L,K) values equal to continued-fraction convergents
			    // of log(3)/log(2) are most likely to be missed if
			    // the parameters haven't been set properly.
unsigned int max=0x4000000; // maximum allowable value of first word of double-
			    // word
int i;
unsigned int kount,attind,o,gomind,gop;
unsigned int c,t,e,g,h,j,n,count,total,flag,wrap,oldlop,lodds,bigind;
unsigned char cyc[60000];
unsigned int lopcnt[500];
unsigned int factor[400];
unsigned int bigcyc[200*5];
unsigned int bigatt[400*3];
unsigned int gom[5000*3];
unsigned int index,z,lopind,lop,K[2],T[2],U[2],L[2],C[2],V[2],Y[2],Z[2],P[3];
unsigned int KP[2],MIN[2];
FILE *Outfp;
Outfp = fopen("out0a.dat","w");
if (cstart==(cstart/3)*3) {
   printf("c divisible by 3 \n");
   goto zskip;
   }
if (cstart==(cstart/2)*2) {
   printf("c divisible by 2 \n");
   goto zskip;
   }
if (cend==(cend/3)*3) {
   printf("c divisible by 3 \n");
   goto zskip;
   }
if (cend==(cend/2)*2) {
   printf("c divisible by 2 \n");
   goto zskip;
   }
if (cend>50035) {
   printf("Warning: Cycles may not be primitive and may have elements that won't fit in double-words. \n");
   goto zskip;
   }
for (h=0; h<100*5; h++)
   bigcyc[h]=0;
bigind=0;
Z[0]=0;
Z[1]=0;
fprintf(Outfp,"// c=%d \n",cstart);
fprintf(Outfp,"int sinp[     ]={ \n");
if (mode==0)
   wrap=15;
else
   wrap=3;
index=0;
lopind=0;
total=0;
for (c=cstart; c<=cend; c+=2) {
   if (c==(c/3)*3)
      continue;
   printf("c=%d \n",c);
   kount=0;
   attind=0;
   gomind=0;
   C[0]=0;
   C[1]=c;
   count=0;
   z=0;
   if (c!=1) {
      factor[0]=c;
      count=1;
      for (i=1; i<1229; i++) {
	 t=(unsigned int)table[i];
	 if (t>c)
	    break;
	 if (c==(c/t)*t) {
	    factor[count]=t;
	    count=count+1;
	    }
	 }
      }
   lop=0;
   for (i=3-(int)limit; i<(int)limit; i+=6) {
      K[1]=(unsigned int)i;
      if (i<0)			    // initial sequence value
	 K[0]=0xffffffff;
      else
	 K[0]=0;
      for (h=0; h<iters; h++) {
	 add64(C, K);
	 n=0;			    // clear count
	 while ((K[1]&3)==0) {
	    K[1]=(K[1]>>2)|(K[0]<<30);
	    K[0]=(int)K[0]>>2;
	    n=n+2;
	    }
	 if ((K[1]&1)==0) {
	    K[1]=(K[1]>>1)|(K[0]<<31);
	    K[0]=(int)K[0]>>1;
	    n=n+1;
	    }
	 T[0]=K[0];
	 T[1]=K[1];
	 flag=0;
	 if ((K[0]&0x80000000)!=0) {
	    sub64(Z, K);
	    flag=1;
	    }
	 for (j=0; j<n; j++) {
	    if (K[0]>max)
	       goto askip;
	    L[0]=(K[0]<<1)|(K[1]>>31);
	    L[1]=K[1]<<1;
	    add64(L, K);
	    }
	 if (flag!=0)
	    sub64(Z, K);
	 L[0]=C[0];
	 L[1]=C[1];
	 sub64(K, L);
	 K[1]=(L[1]>>1)|(L[0]<<31);
	 K[0]=(int)L[0]>>1;
	 if ((K[1]&1)==0)
	    goto askip;
	 }
//
// check if primitive
//
      T[0]=K[0];
      T[1]=K[1];
      if ((T[0]&0x80000000)!=0)
	 sub64(Z, T);
      for (e=0; e<count; e++) {
	 t=factor[e];
	 div64_32(T, U, t);
	 mul64_32(U[0], U[1], t, P);
	 if ((T[0]==P[1])&&(T[1]==P[2]))
	    goto askip;
	 }
//
      L[0]=(K[0]<<1)|(K[1]>>31);
      L[1]=K[1]<<1;
      add64(L, K);
      add64(C, K);
      T[0]=K[0];
      T[1]=K[1];
      if ((T[0]&0x80000000)!=0)
	 sub64(Z, T);
      if (T[0]>max)
	 goto askip;
      while ((K[1]&3)==0) {
	 K[1]=(K[1]>>2)|(K[0]<<30);
	 K[0]=(int)K[0]>>2;
	 }
      if ((K[1]&1)==0) {
	 K[1]=(K[1]>>1)|(K[0]<<31);
	 K[0]=(int)K[0]>>1;
	 }
//
      KP[0]=K[0];
      KP[1]=K[1];
      L[0]=(KP[0]<<1)|(KP[1]>>31);
      L[1]=KP[1]<<1;
      add64(L, KP);
      add64(C, KP);
      while ((KP[1]&3)==0) {
	 KP[1]=(KP[1]>>2)|(KP[0]<<30);
	 KP[0]=(int)KP[0]>>2;
	 }
      if ((KP[1]&1)==0) {
	 KP[1]=(KP[1]>>1)|(KP[0]<<31);
	 KP[0]=(int)KP[0]>>1;
	 }
//
// begin loop
//
bloop:L[0]=(K[0]<<1)|(K[1]>>31);
      L[1]=K[1]<<1;
      add64(L, K);
      add64(C, K);
      T[0]=K[0];
      T[1]=K[1];
      if ((T[0]&0x80000000)!=0)
	 sub64(Z, T);
      if (T[0]>max)
	 goto askip;
      while ((K[1]&3)==0) {
	 K[1]=(K[1]>>2)|(K[0]<<30);
	 K[0]=(int)K[0]>>2;
	 }
      if ((K[1]&1)==0) {
	 K[1]=(K[1]>>1)|(K[0]<<31);
	 K[0]=(int)K[0]>>1;
	 }
//
      L[0]=(KP[0]<<1)|(KP[1]>>31);
      L[1]=KP[1]<<1;
      add64(L, KP);
      add64(C, KP);
      while ((KP[1]&3)==0) {
	 KP[1]=(KP[1]>>2)|(KP[0]<<30);
	 KP[0]=(int)KP[0]>>2;
	 }
      if ((KP[1]&1)==0) {
	 KP[1]=(KP[1]>>1)|(KP[0]<<31);
	 KP[0]=(int)KP[0]>>1;
	 }
      L[0]=(KP[0]<<1)|(KP[1]>>31);
      L[1]=KP[1]<<1;
      add64(L, KP);
      add64(C, KP);
      T[0]=KP[0];
      T[1]=KP[1];
      if ((T[0]&0x80000000)!=0)
	 sub64(Z, T);
      if (T[0]>max)
	 goto askip;
      while ((KP[1]&3)==0) {
	 KP[1]=(KP[1]>>2)|(KP[0]<<30);
	 KP[0]=(int)KP[0]>>2;
	 }
      if ((KP[1]&1)==0) {
	 KP[1]=(KP[1]>>1)|(KP[0]<<31);
	 KP[0]=(int)KP[0]>>1;
	 }
      if ((K[0]!=KP[0])||(K[1]!=KP[1]))
	 goto bloop;
//
// find attachment point
//
      MIN[0]=KP[0];
      MIN[1]=KP[1];
      flag=0;
      L[0]=(K[0]<<1)|(K[1]>>31);
      L[1]=K[1]<<1;
      add64(L, K);
      add64(C, K);
      g=1;
      if ((K[1]&7)==0) {
	 Y[0]=K[0];
	 Y[1]=K[1];
	 flag=1;
	 }
      while ((K[1]&3)==0) {
	 K[1]=(K[1]>>2)|(K[0]<<30);
	 K[0]=(int)K[0]>>2;
	 g=g+2;
	 }
      if ((K[1]&1)==0) {
	 K[1]=(K[1]>>1)|(K[0]<<31);
	 K[0]=(int)K[0]>>1;
	 g=g+1;
	 }
      T[0]=K[0];
      T[1]=K[1];
      sub64(MIN, T);
      if ((T[0]&0x80000000)==0) {
	 MIN[0]=K[0];
	 MIN[1]=K[1];
	 }
      o=1;
      while ((K[0]!=KP[0])||(K[1]!=KP[1])) {
	 L[0]=(K[0]<<1)|(K[1]>>31);
	 L[1]=K[1]<<1;
	 add64(L, K);
	 add64(C, K);
	 g=g+1;
	 if (((K[1]&7)==0)&&(flag==0)) {
	    Y[0]=K[0];
	    Y[1]=K[1];
	    flag=1;
	    }
	 while ((K[1]&3)==0) {
	    K[1]=(K[1]>>2)|(K[0]<<30);
	    K[0]=(int)K[0]>>2;
	    g=g+2;
	    }
	 if ((K[1]&1)==0) {
	    K[1]=(K[1]>>1)|(K[0]<<31);
	    K[0]=(int)K[0]>>1;
	    g=g+1;
	    }
	 T[0]=K[0];
	 T[1]=K[1];
	 sub64(MIN, T);
	 if ((T[0]&0x80000000)==0) {
	    MIN[0]=K[0];
	    MIN[1]=K[1];
	    }
	 o=o+1;
	 }
      if (flag==1)
	 goto askip;
      gop=(g<<16)|o;
      for (h=0; h<gomind; h++) {
	 if (gop==gom[3*h]) {
	    if ((MIN[0]==gom[3*h+1])&&(MIN[1]==gom[3*h+2]))
	       goto askip;
	    }
	 }
      gom[3*gomind]=gop;
      gom[3*gomind+1]=MIN[0];
      gom[3*gomind+2]=MIN[1];
      gomind=gomind+1;
      if (gomind>4999) {
	 printf("array not big enough (gom) \n");
	 goto zskip;
	 }
//
// new cycle
//
      kount=kount+1;
      oldlop=lop;
      lodds=0;
      flag=0;
      L[0]=(KP[0]<<1)|(KP[1]>>31);
      L[1]=KP[1]<<1;
      add64(L, KP);
      add64(C, KP);
      T[0]=KP[0];
      T[1]=KP[1];
      V[0]=T[0];
      V[1]=T[1];
//
// begin loop
//
aloop:if ((T[1]&3)==0) {
	 T[1]=(T[1]>>2)|(T[0]<<30);
	 T[0]=(int)T[0]>>2;
	 }
      else {
	 T[1]=(T[1]>>1)|(T[0]<<31);
	 T[0]=(int)T[0]>>1;
	 }
      lop=lop+1;
      lodds=lodds+1;
      printf("i=%d, n=%d, k=%#010x %#010x \n",i,g,T[0],T[1]);
      if (mode==0)
	 fprintf(Outfp,"%d,",T[1]);
      else
	 fprintf(Outfp,"%#010x,%#010x,",T[0],T[1]);
      total=total+1;
      if (total>wrap) {
	 fprintf(Outfp,"\n");
	 total=0;
	 }
      cyc[index]=z;
      index=index+1;
      if (index>59999) {
	 printf("error: array not big enough \n");
	 goto zskip;
	 }
      U[0]=T[0];
      U[1]=T[1];
      if ((U[0]&0x80000000)!=0)
	 sub64(Z, U);
      if ((U[0]!=0)||((U[1]&0x80000000)!=0)) {
	 flag=1;
	 bigatt[3*attind]=c;
	 bigatt[3*attind+1]=U[0];
	 bigatt[3*attind+2]=U[1];
	 attind=attind+1;
	 if (attind>399) {
	    printf("array too small (bigatt) \n");
	    goto zskip;
	    }
	 }
      L[0]=(T[0]<<1)|(T[1]>>31);
      L[1]=T[1]<<1;
      add64(L, T);
      add64(C, T);
      if ((T[0]!=V[0])||(T[1]!=V[1]))
	 goto aloop;
//
// save big odd elements
//
      if (flag!=0) {
	 printf("big odd element: c=%d, count=%d \n",c,lop-oldlop);
	 bigcyc[bigind*5]=c;
	 bigcyc[bigind*5+1]=lop-oldlop;
	 bigcyc[bigind*5+2]=g-2*lodds;
	 bigcyc[bigind*5+3]=lodds;
	 bigcyc[bigind*5+4]=0;
	 bigind=bigind+1;
	 if (bigind>199) {
	    printf("array not big enough \n");
	    goto zskip;
	    }
	 }
      printf("L=%d, K=%d, count=%d \n",g-2*lodds,lodds,kount);
      printf("\n");
      z=z+1;
askip:n=0;
      }
   lopcnt[lopind]=lop;
   lopind=lopind+1;
   printf("\n");
   }
fprintf(Outfp,"\n");
fprintf(Outfp,"unsigned char cflag[      +1]={ \n");
total=0;
for (h=0; h<index; h++) {
   fprintf(Outfp,"%d,",cyc[h]);
   total=total+1;
   if (total>20) {
      fprintf(Outfp,"\n");
      total=0;
      }
   }
fprintf(Outfp,"255};");
fprintf(Outfp," //  index=%d \n",index);
fprintf(Outfp,"unsigned int size[   ]={ \n");
for (h=0; h<lopind; h++)
   fprintf(Outfp,"  %d,\n",lopcnt[h]);
fprintf(Outfp,"int cval[   ]={ \n");
index=0;
total=0;
for (h=(unsigned int)cstart; h<=(unsigned int)cend; h+=2) {
   if (h==(h/3)*3)
      continue;
   index=index+1;
   fprintf(Outfp,"%d,",h);
   total=total+1;
   if (total>10) {
      fprintf(Outfp,"\n");
      total=0;
      }
   }
fprintf(Outfp,"\n");
fprintf(Outfp,"unsigned int numbc=%d;   \n",index);
if (bigind==0)
   goto zskip;
printf("big cycles \n");
fprintf(Outfp,"// big cycles \n");
fprintf(Outfp,"\n");
for (h=0; h<bigind; h++) {
   printf("c=%d, count=%d, L=%d, K=%d, a=%d \n",bigcyc[h*5],bigcyc[h*5+1],bigcyc[h*5+2],bigcyc[h*5+3],bigcyc[h*5+4]);
   fprintf(Outfp,"c=%d, count=%d, L=%d, K=%d, a=%d \n",bigcyc[h*5],bigcyc[h*5+1],bigcyc[h*5+2],bigcyc[h*5+3],bigcyc[h*5+4]);
   }
printf("\n");
printf("big odd elements \n");
fprintf(Outfp,"// big odd elements \n");
fprintf(Outfp,"\n");
for (h=0; h<attind; h++) {
   printf("c=%d, k=%#010x %#010x \n",bigatt[3*h],bigatt[3*h+1],bigatt[3*h+2]);
   fprintf(Outfp,"c=%d, k=%#010x %#010x \n",bigatt[3*h],bigatt[3*h+1],bigatt[3*h+2]);
   }
zskip:
fclose(Outfp);
return(0);
}