/* 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); }