/*****************************************************************************/
/* */
/* FIND CYCLES IN THE 3N+C SEQUENCE */
/* 03/22/11 (dkc) */
/* */
/* This C program finds cycles in the 3n+c sequence. Attachment points in */
/* cycles and t values (odd integers divisible by 3) of the extended */
/* sequences are output. Cycles not having attachment points are also */
/* found (in this case, odd elements of the cycles are output). Lengths of */
/* the cycles (where the element after an odd element i in the 3n+c */
/* sequence is defined to be 3i+c) are also output. Domains of the */
/* exponential functions for the primary attachment points and t values are */
/* also output. The number of even elements in the cycle (denoted by L) */
/* and the number of odd elements in the cycle (denoted by K) are also */
/* output. (In this case, the element after an odd element i is defined */
/* to be (3i+c)/2.) */
/* */
/*****************************************************************************/
#include <stdio.h>
#include <math.h>
/*****************************************************************************/
/* */
/* regenerate cycles */
/* */
/*****************************************************************************/
int regen(int c, int s, int flag, int *w) {
int k,tres;
unsigned int j,count;
//
// find odd natural number divisible by 3
//
k=s;
while (k!=(k/3)*3) {
if ((k&1)==0) {
if ((k-c)==((k-c)/3)*3)
k=(k-c)/3;
else
k=k*2;
}
else
k=k*2;
}
tres=k;
//
k=s;
count=2;
if (flag==0)
goto zskip;
w[1]=k;
while ((k&1)==0) {
k=k/2;
w[count]=k;
count=count+1;
}
for (j=1; j<1000000; j++) {
k=3*k+c;
w[count]=k;
count=count+1;
if (count>32767) {
tres=0;
goto zskip;
}
while ((k&1)==0) {
if (k==s)
goto zskip;
k=k/2;
w[count]=k;
count=count+1;
if (count>32767) {
tres=0;
goto zskip;
}
}
}
tres=0;
zskip:
w[0]=count-1;
return(tres);
}
int main () {
//
int cstart=1103; // beginning c value
int cend=1199; // ending c value
unsigned int iters=2; // number of "jumps" from the initial value
// // usually set to 1 or 2
//
unsigned int d,e,f,g,h,j,m,n,o,count,sume,sumo,ne,no,ma,mi,delta,numb;
int limit,temp,max,i,k,save[32],v[32768],w[32768],t,u,tres,saves,oldk,bsave;
int out[4000*2];
int factor[200],t1,t2,c,maxodd,minodd,tmpodd;
unsigned int x[2000],kx[2000],levens,lodds,first,cmin,cmax;
unsigned int table[167]={5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,
73,79,83,89,97,101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,
179,181,191,193,197,199,211,223,227,229,233,239,241,251,257,263,269,271,277,
281,283,293,307,311,313,317,331,337,347,349,353,359,367,373,379,383,389,397,
401,409,419,421,431,433,439,443,449,457,461,463,467,479,487,491,499,503,509,
521,523,541,547,557,563,569,571,577,587,593,599,601,607,613,617,619,631,641,
643,647,653,659,661,673,677,683,691,701,709,719,727,733,739,743,751,757,761,
769,773,787,797,809,811,821,823,827,829,839,853,857,859,863,877,881,883,887,
907,911,919,929,937,941,947,953,967,971,977,983,991,997,1009};
double lambdae,lambdao,y,tmax,tmin,kmax,kmin,maxdel,mindel;
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 (cstart<0) {
printf("negative c value \n");
goto zskip;
}
if (cstart>5045) {
printf("warning: Cycle may not be primitive. \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<0) {
printf("negative c value \n");
goto zskip;
}
if (cend>5045) {
printf("warning: Cycle may not be primitive. \n");
goto zskip;
}
maxdel=-1000000.0;
mindel=1000000.0;
for (c=cstart; c<=cend; c+=2) {
if (c==(c/3)*3)
continue;
printf("c=%d \n",c);
printf("\n");
fprintf(Outfp,"c=%d \n",c);
fprintf(Outfp,"\n");
count=0;
if (c!=1) {
factor[0]=c;
count=1;
}
for (i=0; i<167; i++) {
t=(int)table[i];
if (c==(c/t)*t) {
factor[count]=t;
count=count+1;
}
}
limit=1500000; // must be a multiple of 6
max=1<<29;
m=0;
for (i=3-limit; i<limit; i+=6) {
k=(int)i; // initial sequence value
for (e=0; e<iters; e++) {
k=k+c;
n=0; // clear count
while ((k&1)==0) {
k=k/2;
n=n+1;
}
temp=k;
if (temp<0)
temp=-temp;
for (j=0; j<n; j++) {
if (temp>max) // check for overflow
goto askip;
temp=temp*3;
k=k*3;
}
k=(k-c)/2; // sequence value after a "jump"
if ((k&1)==0)
goto bskip;
}
temp=k;
if (temp<0)
temp=-temp;
if (temp>max) // check for overflow
goto askip;
k=3*k+c;
bskip:
h=0; // clear index
n=0; // clear index
while ((k&1)==0) {
v[n]=k; // save sequence value
n=n+1; // increment index
if (h<32) { // check index
save[h]=k; // save even sequence value
h=h+1; // increment index
}
else { // array too small
printf("error: array too small \n");
goto zskip;
}
k=k/2;
}
bsave=k; // save odd sequence value
for (j=1; j<1000000; j++) {
v[n]=k; // save sequence value
n=n+1; // increment index
if (n>32767) // array full
goto askip;
temp=k;
if (temp<0)
temp=-temp;
if (temp>max)
goto askip;
k=3*k+c;
while ((k&1)==0) {
for (g=0; g<h; g++) {
if (k==save[g]) {
for (f=g; f<n; f++) {
if ((v[f]&7)==0) {
for (e=0; e<m; e++) {
if (((int)(n-g)==out[2*e])&&((v[f]/4)==out[2*e+1])) {
goto askip;
}
}
t=v[f]/4;
if (t<0)
t=-t;
for (e=0; e<count; e++) {
u=factor[e];
if (t==(t/u)*u)
goto askip;
}
tres=regen(c, v[f]/4, 1, w);
if (tres==0) {
printf("return value equals 0 \n");
goto zskip;
}
t1=tres;
first=1;
temp=v[f]/4;
if (temp<0)
temp=-temp;
kx[0]=temp;
sume=temp;
temp=tres;
if (temp<0)
temp=-temp;
x[0]=temp;
sumo=temp;
ne=1;
no=1;
printf("n=%d k=%d t=%d \n",n-g,v[f]/4,tres);
fprintf(Outfp,"n=%d k=%d t=%d \n",n-g,v[f]/4,tres);
out[2*m]=n-g;
out[2*m+1]=v[f]/4;
m=m+1;
if (m>=4000) {
printf("output array full \n");
goto zskip;
}
o=w[0];
saves=v[f]/4;
delta=0;
while ((saves&3)==0) {
saves=saves/4;
delta=delta+2;
if ((saves&1)==1)
goto qskip;
tres=regen(c, saves, 0, w);
if (tres==0) {
printf("return value equals 0 \n");
goto zskip;
}
if (first==1) {
t2=tres;
first=0;
}
temp=tres;
if (temp<0)
temp=-temp;
x[no]=temp;
sumo=sumo+temp;
no=no+1;
if (no>999) {
printf("array not big enough \n");
goto zskip;
}
printf(" k=%d t=%d \n",saves,tres);
fprintf(Outfp," k=%d t=%d \n",saves,tres);
out[2*m]=n-g;
out[2*m+1]=saves;
m=m+1;
if (m>=4000) {
printf("output array full \n");
goto zskip;
}
e=0;
}
if ((saves&1)==0) {
saves=saves/2;
delta=delta+1;
}
qskip: saves=3*saves+c;
delta=delta+1;
oldk=0;
for (d=2+delta; d<o; d++) {
if (((w[d]&1)==0)&&(w[d]!=saves)&&((w[d]-c)==((w[d]-c)/3)*3)) {
tres=regen(c, w[d], 0, w);
if (tres==0) {
printf("return value equals 0 \n");
goto zskip;
}
if (first==1) {
t2=tres;
first=0;
}
temp=tres;
if (temp<0)
temp=-temp;
x[no]=temp;
sumo=sumo+temp;
no=no+1;
if (no>999) {
printf("array not big enough \n");
goto zskip;
}
if ((w[d]*4)!=oldk) {
temp=w[d];
if (temp<0)
temp=-temp;
kx[ne]=temp;
sume=sume+temp;
ne=ne+1;
}
oldk=w[d];
printf(" k=%d t=%d \n",w[d],tres);
fprintf(Outfp," k=%d t=%d \n",w[d],tres);
out[2*m]=n-g;
out[2*m+1]=w[d];
m=m+1;
if (m>=4000) {
printf("output array full \n");
goto zskip;
}
}
if ((w[d]&1)==1)
saves=3*w[d]+c;
e=0;
}
if (o>3) {
if (w[1]==(w[o-2]*4)) {
printf("unexpected duplicate: w=%d %d %d %d \n",w[1],w[2],w[o-2],w[o-1]);
temp=w[o-2];
if (temp<0)
temp=-temp;
sume=sume-temp;
ne=ne-1;
}
if (((w[1]*4)==w[o-2])&&((w[o-3]&1)!=1)) {
temp=w[1];
if (temp<0)
temp=-temp;
sume=sume-temp;
ne=ne-1;
}
}
lambdae=(double)sume/ne;
lambdao=(double)sumo/no;
if (ne==1)
goto mskip;
printf("lambda (k)=%e n=%d, lambda (t)=%e n=%d \n",lambdae,ne,lambdao,no);
fprintf(Outfp,"lambda (k)=%e n=%d, lambda (t)=%e n=%d\n",lambdae,ne,lambdao,no);
ma=0;
mi=0x7fffffff;
for (e=0; e<(int)ne; e++) {
y=-log((double)kx[e]/lambdae)/lambdae;
if (kx[e]>ma) {
ma=kx[e];
kmax=y;
}
if (kx[e]<mi) {
mi=kx[e];
kmin=y;
}
}
printf("max/min (k)=(%d, %d), domain=(%e, %e) \n",ma,mi,kmax,kmin);
fprintf(Outfp,"max/min (k)=(%d, %d), domain=(%e, %e) \n",ma,mi,kmax,kmin);
if (kmin<-kmax) {
printf("warning: Maximum x not greater than absolute value of minimum. \n");
}
mskip: if (no==1)
goto nskip;
if ((no==2)&&(t1==-t2)) {
printf("Count equals 2 and absolute values of t values same.\n");
goto nskip;
}
ma=0;
mi=0x7fffffff;
for (e=0; e<(int)no; e++) {
y=-log((double)x[e]/lambdao)/lambdao;
if (x[e]>ma) {
ma=x[e];
tmax=y;
}
if (x[e]<mi) {
mi=x[e];
tmin=y;
}
}
y=tmin-tmax;
if (y>maxdel) {
maxdel=y;
cmax=c;
}
if (y<mindel) {
mindel=y;
cmin=c;
}
printf("max/min (t)=(%d, %d), domain=(%e, %e), d=%e \n",ma,mi,tmax,tmin,y);
fprintf(Outfp,"max/min (t)=(%d %d), domain=(%e, %e), d=%e \n",ma,mi,tmax,tmin,y);
if (tmin<-tmax) {
printf("warning: Maximum x not greater than absolute value of minimum. \n");
}
if (ne!=1) {
if (tmax>kmax) {
printf("warning: Not within t domain. \n");
}
if (tmin<kmin) {
printf("warning: Not within t domain. \n");
}
}
nskip: k=w[1];
levens=0;
while ((k&1)==0) {
k=k/2;
levens=levens+1;
}
lodds=1;
maxodd=0;
minodd=0x7fffffff;
for (j=1; j<1000000; j++) {
tmpodd=k;
if (tmpodd<0)
tmpodd=-tmpodd;
if (tmpodd>maxodd)
maxodd=tmpodd;
if (tmpodd<minodd)
minodd=tmpodd;
k=3*k+c;
while ((k&1)==0) {
if (k==w[1]) {
levens=levens-1;
goto oskip;
}
k=k/2;
levens=levens+1;
}
levens=levens-1;
lodds=lodds+1;
}
printf("error \n");
goto zskip;
oskip: printf("L=%d, K=%d, min=%d, max=%d \n",levens,lodds,minodd,maxodd);
fprintf(Outfp,"L=%d, K=%d, min=%d, max=%d \n",levens,lodds,minodd,maxodd);
if ((2*minodd>maxodd)&&(lodds!=1)) {
printf("warning: Twice mininum odd element is greater than maximum odd element. \n");
fprintf(Outfp,"warning: Twice mininum odd element is greater than maximum odd element. \n");
}
printf("\n");
fprintf(Outfp,"\n");
goto askip;
}
}
for (e=0; e<m; e++) {
if (((int)(n-g)==out[2*e])&&(bsave==out[2*e+1])) {
goto askip;
}
}
t=bsave;
if (t<0)
t=-t;
for (e=0; e<count; e++) {
u=factor[e];
if (t==(t/u)*u)
goto askip;
}
printf("n=%d o=%d",n-g,bsave);
fprintf(Outfp,"n=%d o=%d",n-g,bsave);
out[2*m]=n-g;
out[2*m+1]=(int)bsave;
m=m+1;
if (m>=4000) {
printf("output array full \n");
goto zskip;
}
numb=2;
levens=0;
lodds=1;
maxodd=bsave;
if (maxodd<0)
maxodd=-maxodd;
minodd=bsave;
if (minodd<0)
minodd=-minodd;
k=bsave;
for (j=1; j<1000000; j++) {
k=3*k+c;
while ((k&1)==0) {
k=k/2;
levens=levens+1;
}
levens=levens-1;
lodds=lodds+1;
if (k==bsave) {
lodds=lodds-1;
goto pskip;
}
out[2*m]=n-g;
out[2*m+1]=k;
m=m+1;
if (m>=4000) {
printf("output array full \n");
goto zskip;
}
tmpodd=k;
if (tmpodd<0)
tmpodd=-tmpodd;
if (maxodd<tmpodd)
maxodd=tmpodd;
if (minodd>tmpodd)
minodd=tmpodd;
printf(" %d",k);
fprintf(Outfp," %d",k);
numb=numb+1;
if (numb>8) {
printf("\n");
fprintf(Outfp,"\n");
numb=0;
}
}
printf("error \n");
goto zskip;
pskip: if (numb!=0) {
printf("\n");
fprintf(Outfp,"\n");
}
printf("L=%d, K=%d \n",levens,lodds);
fprintf(Outfp,"L=%d, K=%d \n",levens,lodds);
if (minodd!=maxodd) {
if ((2*minodd)>maxodd) {
printf("minodd=%d, maxodd=%d \n",minodd,maxodd);
fprintf(Outfp,"minodd=%d, maxodd=%d \n",minodd,maxodd);
}
}
printf("\n");
fprintf(Outfp,"\n");
goto askip;
}
}
v[n]=k;
n=n+1;
if (n>32767)
goto askip;
k=k/2;
}
}
askip:
n=0;
}
printf("\n");
fprintf(Outfp,"\n");
}
printf("\n");
fprintf(Outfp,"\n");
printf("maximum interval=%e, minimum interval=%e \n",maxdel,mindel);
fprintf(Outfp,"maximum interval=%e, minimum interval=%e \n",maxdel,mindel);
printf("c (for maximum)=%d, c (for minimum)=%d \n",cmax,cmin);
fprintf(Outfp,"c (for maximum)=%d, c (for minimum)=%d \n",cmax,cmin);
zskip:
fclose(Outfp);
return(0);
}