/*****************************************************************************/
/* */
/* FIND CYCLES IN THE 3N+C SEQUENCE */
/* 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. Cycles without */
/* attachment points are also found. */
/* */
/*****************************************************************************/
#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);
unsigned int jumps(unsigned int K0, unsigned int K1, unsigned int iters,
unsigned int c, unsigned int max, unsigned int *output);
unsigned int cycle(unsigned int K0, unsigned int K1, unsigned int max,
unsigned int c, unsigned int *output);
unsigned int regen(unsigned int K0, unsigned int K1, unsigned int *output,
unsigned int c);
unsigned int check(unsigned int gop, unsigned int MIN0, unsigned int MIN1,
unsigned int *gom, unsigned int gomind);
void dummy(unsigned int A, unsigned int B, unsigned int C, unsigned int D,
unsigned int E, unsigned int F, unsigned int G, unsigned int H);
far unsigned int error[10];
far unsigned int clkout[10000*4];
far unsigned int output[5];
far unsigned char cyc[90000];
far unsigned int out[90000];
far unsigned int lopcnt[500];
far unsigned int factor[400];
far unsigned int bigcyc[200*5];
far unsigned int bigatt[400*3];
far unsigned int gom[10000*3];
far unsigned int noatt[1000*3];
far unsigned int oddout[1000*4];
far unsigned int cval[400];
int main () {
//
unsigned int cstart=1; // beginning c value
unsigned int cend=499; // ending c value
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 to find all
// cycles. When "iters" is set to 0, every "twig"
// of the tree in the range determined by "limit"
// is checked.
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
unsigned int outflag=0; // output flag (if set to 0, (L+K, K) values and
// minima are output, otherwise attachment points
// are output)
int i;
unsigned int first,acount,kount,attind,o,gomind,gop,numbc,noattind,oddind;
unsigned int n,c,t,e,g,h,j,count,flag,oldlop,lodds,bigind,clkind;
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],itersp,limitp;
for (j=0; j<10; j++)
error[j]=0;
if (cstart==(cstart/3)*3) {
error[0]=1;
goto zskip;
}
if (cstart==(cstart/2)*2) {
error[0]=2;
goto zskip;
}
if (cend==(cend/3)*3) {
error[0]=3;
goto zskip;
}
if (cend==(cend/2)*2) {
error[0]=4;
goto zskip;
}
if (cend>199999*5) {
error[0]=5;
goto zskip;
}
index=0;
n=0;
for (h=(unsigned int)cstart; h<=(unsigned int)cend; h+=2) {
if (h==(h/3)*3)
continue;
cval[index]=h;
index=index+1;
if (index>=400) {
error[0]=10;
goto zskip;
}
}
numbc=index;
for (h=0; h<100*5; h++)
bigcyc[h]=0;
output[0]=0;
output[1]=0;
bigind=0;
Z[0]=0;
Z[1]=0;
index=0;
oddind=0;
lopind=0;
attind=0;
clkind=0;
for (c=cstart; c<=cend; c+=2) {
if (c==(c/3)*3)
continue;
if ((c==31639)||(c==39409)||(c==71515)||(c==85105)||(c==169495)||(c==178825)) {
itersp=0;
limitp=24000000;
}
else {
if (c==158195) {
itersp=0;
limitp=30000000;
}
else {
if (c==185357) {
itersp=0;
limitp=60000000;
}
else {
itersp=iters;
limitp=limit;
}
}
}
error[1]=c;
kount=0;
gomind=0;
noattind=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)limitp; i<(int)limitp; i+=6) {
K[1]=(unsigned int)i;
error[2]=K[1];
if (i<0) // initial sequence value
K[0]=0xffffffff;
else
K[0]=0;
//
// do jumps
//
/* for (h=0; h<itersp; 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;
} */
//
// equivalent C64 assembly language
//
if (itersp>0) {
t=jumps(K[0],K[1],itersp,c,max,output);
if (t==0)
goto askip;
K[0]=output[0];
K[1]=output[1];
}
//
// 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;
}
//
// find cycle
//
/* 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; */
//
// equivalent C64 assembly language
//
o=cycle(K[0],K[1],max,c,output);
if (o==2) {
error[0]=11;
// goto zskip;
goto askip;
}
if (o==0)
goto askip;
K[0]=output[0];
K[1]=output[1];
KP[0]=K[0];
KP[1]=K[1];
//
// find attachment point
//
/* MIN[0]=K[0];
MIN[1]=K[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; */
//
// equivalent C64 assembly language (with additional check for cycles without
// attachment points)
//
t=regen(KP[0],KP[1],output,c);
MIN[0]=output[0];
MIN[1]=output[1];
gop=output[2];
Y[0]=output[3];
Y[1]=output[4];
g=gop>>16;
o=gop&0xffff;
if (t==0) {
if (noattind==0)
goto cjump;
t=check(gop, MIN[0], MIN[1], noatt, noattind);
if (t==0)
goto askip;
cjump: noatt[3*noattind]=gop;
noatt[3*noattind+1]=MIN[0];
noatt[3*noattind+2]=MIN[1];
noattind=noattind+1;
if (noattind>999) {
error[0]=14;
goto zskip;
}
oddout[4*oddind]=c;
oddout[4*oddind+1]=((g-o)<<16)|o;
oddout[4*oddind+2]=MIN[0];
oddout[4*oddind+3]=MIN[1];
oddind=oddind+1;
if (oddind>999) {
error[0]=13;
goto zskip;
}
goto askip;
}
//
// check for duplicate cycles
//
/* 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;
}
}
} */
//
// equivalent C64 assembly language
//
if (gomind==0)
goto bjump;
t=check(gop, MIN[0], MIN[1], gom, gomind);
if (t==0)
goto askip;
bjump:
gom[3*gomind]=gop;
gom[3*gomind+1]=MIN[0];
gom[3*gomind+2]=MIN[1];
gomind=gomind+1;
if (gomind>9999) {
error[0]=6;
goto zskip;
}
if (outflag==0) {
clkout[4*clkind]=c;
clkout[4*clkind+1]=((g-o)<<16)|o;
clkout[4*clkind+2]=MIN[0];
clkout[4*clkind+3]=MIN[1];
clkind=clkind+1;
if (clkind>9999) {
error[0]=12;
goto zskip;
}
goto askip;
}
//
// 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) {
error[0]=7;
goto zskip;
}
}
out[index]=T[1];
cyc[index]=z;
index=index+1;
if (index>89999) {
error[0]=8;
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) {
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) {
error[0]=9;
goto zskip;
}
}
z=z+1;
askip:n=n+1;
}
if (outflag!=0) {
lopcnt[lopind]=lop;
lopind=lopind+1;
}
}
error[2]=index;
error[3]=lopind;
error[4]=numbc;
error[5]=bigind;
error[6]=attind;
error[7]=clkind;
error[8]=oddind;
zskip:
return(0);
}