/*****************************************************************************/
/* */
/* FIND CYCLES IN THE 3N+C SEQUENCE (cycles with attachment points only) */
/* 02/19/13 (dkc) */
/* */
/* This C program finds cycles in the 3n+c sequence. Attachment points are */
/* output. Double-precision (64-bit) words are used. */
/* */
/*****************************************************************************/
#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=48001; // beginning c value
unsigned int cend=48097; // ending c value
unsigned int mode=0; // double-word flag
unsigned int iters=1; // number of "jumps" from the initial odd integer
// divisible by 3, usually set to 1 (every cycle having
// an attachment point has at least one one-jump or
// no-jump attachment point). A jump usually results
// in gain of the initial value. Cycles with large
// elements can then be found starting with relatively
// small initial values.
unsigned int limit=1500000; // must be a multiple of 6, should be about half
// of "max" to allow for cycles having a constrained
// dynamic range (a preferable value is 3000000, but
// execution time is too long [cycles with (K+L,K)
// values equal to generalized continued-fraction
// convergents of log(3)/log(2) are most likely to
// be missed])
unsigned int max=0x80000000; // maximum allowable value of second word of double-
// word (when "bigmax" is set to 0), apparently adequate
// for c values of up to about 20000. For larger c
// values, both words are required for some cycles.
// When the mode is set to 0, only the second word
// of attachment points is output. (When "iters" is
// greater than 1, "max" usually must be increased to
// accommodate the gain.)
unsigned int bigmax=0; // if set, "max" is the upper bound of the first word of
// the double word (a typical value of "max" would be 1)
int i;
unsigned int first,acount;
unsigned int c,t,e,f,g,h,j,m,n,count,total,flag,wrap,oldlop,lodds,bigind;
unsigned int save[32*2],v[32768*2];
unsigned char cyc[60000];
unsigned int lopcnt[500];
unsigned int out[15000*2];
unsigned int outn[15000];
unsigned int factor[400];
unsigned int bigcyc[100*5];
unsigned int index,z,lopind,lop,K[2],T[2],U[2],L[2],C[2],V[2],Y[2],Z[2],P[3];
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>49865)
printf("warning: cycle may not be primitive \n");
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);
C[0]=0;
C[1]=c;
count=0;
z=0;
if (c!=1) {
factor[0]=c;
count=1;
}
for (i=0; i<1228; i++) {
t=(unsigned int)table[i];
if (c==(c/t)*t) {
factor[count]=t;
count=count+1;
}
}
m=0;
lop=0;
for (i=3-(int)limit; i<(int)limit; i+=6) {
K[0]=0;
if (i<0) { // initial sequence value
K[1]=-i;
sub64(Z, K);
}
else
K[1]=i;
for (e=0; e<iters; e++) {
add64(C, K);
n=0; // clear count
while ((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];
if ((T[0]&0x80000000)!=0)
sub64(Z, T);
for (j=0; j<n; j++) {
if (bigmax==0) {
if ((T[0]!=0)||(T[1]>max))
goto askip;
}
else {
if (T[0]>max)
goto askip;
}
L[0]=T[0];
L[1]=T[1];
add64(T, T);
add64(L, T);
L[0]=K[0];
L[1]=K[1];
add64(K, K);
add64(L, 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 bskip;
}
T[0]=K[0];
T[1]=K[1];
if ((T[0]&0x80000000)!=0)
sub64(Z, T);
if (bigmax==0) {
if ((T[0]!=0)||(T[1]>max))
goto askip;
}
else {
if (T[0]>max)
goto askip;
}
L[0]=K[0];
L[1]=K[1];
add64(K, K);
add64(L, K);
add64(C, K);
bskip:
h=0; // clear index
n=0; // clear index
while ((K[1]&1)==0) {
v[2*n]=K[0];
v[2*n+1]=K[1];
n=n+1; // increment index
if (h<32) { // check index
save[2*h]=K[0];
save[2*h+1]=K[1];
h=h+1; // increment index
}
else { // array too small
printf("error: array too small \n");
goto zskip;
}
K[1]=(K[1]>>1)|(K[0]<<31);
K[0]=(int)K[0]>>1;
}
for (j=1; j<1000000; j++) {
v[n*2]=K[0];
v[n*2+1]=K[1];
n=n+1; // increment index
if (n>32767) // array full
goto askip;
T[0]=K[0];
T[1]=K[1];
if ((T[0]&0x80000000)!=0)
sub64(Z, T);
if (bigmax==0) {
if ((T[0]!=0)||(T[1]>max))
goto askip;
}
else {
if (T[0]>max)
goto askip;
}
L[0]=K[0];
L[1]=K[1];
add64(K, K);
add64(L, K);
add64(C, K);
while ((K[1]&1)==0) {
for (g=0; g<h; g++) {
if ((K[0]==save[g*2])&&(K[1]==save[g*2+1])) {
for (f=g; f<n; f++) {
if ((v[f*2+1]&7)==0) {
Y[1]=(v[f*2+1]>>2)|(v[f*2]<<30);
Y[0]=(int)v[f*2]>>2;
for (e=0; e<m; e++) {
if (((n-g)==outn[e])&&(Y[0]==out[2*e])&&(Y[1]==out[2*e+1])) {
goto askip;
}
}
// check if primitive
T[0]=Y[0];
T[1]=Y[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;
}
// attachment point
oldlop=lop;
lodds=0;
acount=0;
flag=0;
T[0]=v[f*2];
T[1]=v[f*2+1];
V[0]=T[0];
V[1]=T[1];
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]&0xe0000000)!=0))
flag=1;
printf("n=%d, k=%d \n",n-g,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;
outn[m]=(n-g);
out[2*m]=T[0];
out[2*m+1]=T[1];
m=m+1;
if (m>=15000) {
printf("output array full \n");
goto zskip;
}
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];
L[1]=T[1];
add64(T, T);
add64(L, T);
add64(C, T);
if ((T[0]!=V[0])||(T[1]!=V[1]))
goto aloop;
// new cycle
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]=n-g-2*lodds;
bigcyc[bigind*5+3]=lodds;
bigcyc[bigind*5+4]=acount;
bigind=bigind+1;
if (bigind>99) {
printf("array not big enough \n");
goto zskip;
}
}
printf("L=%d, K=%d, a=%d, big=%d \n",n-g-2*lodds,lodds,acount,bigind);
z=z+1;
goto askip;
}
}
}
}
v[n*2]=K[0];
v[n*2+1]=K[1];
n=n+1;
if (n>32767)
goto askip;
K[1]=(K[1]>>1)|(K[0]<<31);
K[0]=(int)K[0]>>1;
}
}
askip:n=0;
}
lopcnt[lopind]=lop;
lopind=lopind+1;
printf("\n");
}
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);
printf("\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]);
}
zskip:
fclose(Outfp);
return(0);
}