/*****************************************************************************/
/* */
/* FIND CYCLES IN THE 3N+C SEQUENCE (cycles with attachment points only) */
/* 11/17/13 (dkc) */
/* */
/* This C program finds cycles in the 3n+c sequence. Attachment points 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 attachment points. The values in */
/* the output array "cflag" indicate which attachment points are associated */
/* with a cycle. The output array "size" gives the total number of */
/* attachment points 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, "K" is the number of odd */
/* elements, and "a" is the number of primary attachment points. "count" */
/* gives the number of cycles found for a given c value. */
/* */
/*****************************************************************************/
#include <stdio.h>
#include <math.h>
#include "table4.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=1; // beginning c value
unsigned int cend=499; // ending c value
unsigned int mode=0; // double-word flag. When the mode is set to 0,
// only the second word of attachment points is
// output (to save memory).
unsigned int iters=1; // 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). A good c value
// to tune the parameters with is 17021 (where
// there should be 258 cycles). For some larger
// c values, "iters" must be set to 2, 1, or 0 to
// find all cycles. When "iters" is set to 0, every
// "twig" of the tree in the range determined by
// "limit" is checked. Apparently, every cycle
// with attachment points has a no-jump or one-
// jump attachment point. It's a good indication
// that all the cycles have been found if the
// cycles given for "iters" equal to 1 (and say
// "limit" equal to 18000000) match those given
// for "iters" equal to 0 (and say "limit" equal
// to 24000000). Usually, a jump results in gain
// of the initial value; cycles with large elements
// can then be found starting with a relatively
// small initial value.
//
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. Setting
// "iters" to 1 and "limit" to 18000000 is not adequate
// for finding all the cycles for c=31639, 39409, 71515,
// 85105, 158195, 169495, 178825, and 185357 (c
// values where (K+L,K) values of some cycles equal
// generalized continued-fraction convergents of
// log(3)/log(2)). Apparently, setting "iters" to 0
// and "limit" to 60000000 is required to find all
// the cycles for c=185357 (setting "limit" to
// 66000000 doesn't produce any more cycles).
//
unsigned int max=0x4000000; // maximum allowable value of first word of double-
// word
int i;
unsigned int first,acount,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[10000*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>5*199999) {
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;
attind=0;
for (c=cstart; c<=cend; c+=2) {
if (c==(c/3)*3)
continue;
printf("c=%d \n",c);
kount=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<17983; 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);
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;
}
//
// 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==0)
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>9999) {
printf("array not big enough (gom) \n");
goto zskip;
}
//
// new cycle
//
kount=kount+1;
oldlop=lop;
lodds=0;
acount=0;
flag=0;
T[0]=Y[0];
T[1]=Y[1];
V[0]=T[0];
V[1]=T[1];
//
// begin loop
//
aloop:first=0;
while ((T[1]&3)==0) {
T[1]=(T[1]>>2)|(T[0]<<30);
T[0]=(int)T[0]>>2;
if ((T[1]&1)==1)
goto qskip;
if (first==0) {
acount=acount+1;
first=1;
}
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]=T[0];
bigatt[3*attind+2]=T[1];
attind=attind+1;
if (attind>399) {
printf("array too small (bigatt) \n");
goto zskip;
}
}
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;
}
lop=lop+1;
if ((T[0]==V[0])&&(T[1]==V[1])) {
acount=acount-1;
goto rskip;
}
}
if ((T[1]&1)==0) {
T[1]=(T[1]>>1)|(T[0]<<31);
T[0]=(int)T[0]>>1;
}
qskip:lodds=lodds+1;
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 attachment points
//
rskip:if (flag!=0) {
printf("big attachment point: 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]=acount;
bigind=bigind+1;
if (bigind>199) {
printf("array not big enough \n");
goto zskip;
}
}
printf("L=%d, K=%d, a=%d, count=%d \n",g-2*lodds,lodds,acount,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 attachment points \n");
fprintf(Outfp,"// big attachment points \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);
}