/*****************************************************************************/
/* */
/* COMPUTE 3N+C SEQUENCE */
/* 02/22/13 (dkc) */
/* */
/* This C program computes 3n+c sequences starting with odd integers */
/* divisible by 3. When an attachment point is encountered, the sequence */
/* is extended backwards ("odd" paths are taken) to another odd integer */
/* divisible by 3. The differences in the absolute values of the odd */
/* integers divisible by 3 are histogrammed. */
/* */
/*****************************************************************************/
#include <stdio.h>
#include <math.h>
void add64(unsigned int *a, unsigned int *b);
void sub64(unsigned int *a, unsigned int *b);
void div64_32(unsigned int *A, unsigned int *B, unsigned int C);
int main () {
unsigned int c=1; // c value
unsigned int limit=10; // 10, 100, 1000, 10000, 100000, 1000000, 10000000
// When "limit" is less than or equal to 100, the
// individual t1 and t2 values are output.
unsigned int relprime=1; // If set, exclude (t2,c)!=1.
int scale=1536; // Divide |t2|-|t1| by this amount.
// This scaling factor is for "limit" equal to
// 10000.
int offset=550; // This offset (to make the histogram bins
// non-negative) is for "limit" equal to 10000.
unsigned int j,k,m,n,out[2000],flag,count,total,iters,kount,toobig,even1,outrng;
unsigned int L[2],S[2],U[2],V[2],W[2],Z[2],g[1000],h[1000],T[2],K[2],Q[2];
unsigned int histo[10000];
unsigned int THI[2],M[2];
int i,f,temp,thi,diff;
double r;
FILE *Outfp;
Outfp = fopen("out1c.dat","w");
if (c>11) {
printf("c too large \n");
goto zskip;
}
for (i=0; i<10000; i++) // clear histogram of |t2|-|t1| values
histo[i]=0;
Z[0]=0;
Z[1]=0;
for (i=0; i<1000; i++) {
h[i]=0; // clear histogram of k values
g[i]=0; // clear histogram of iteration values
}
outrng=0; // big sequence value counter
even1=0; // number of failures when t2 attaches to the
// first even element after t1
THI[0]=0; // set maximum t2 value
THI[1]=0;
toobig=0; // number of values not histogrammed
kount=0; // number of failures (|t2|<|t1|)
iters=0; // number of starting values
total=0; // number of k values
W[0]=0; // load c
W[1]=c;
for (f=-(int)limit; f<=(int)limit; f++) {
//for (f=0; f<=(int)limit; f++) { // natural numbers
i=3+f*6; // t2, odd n values divisible by 3
if (c!=1) {
if (relprime!=0) {
if (i==(i/(int)c)*(int)c)
continue;
}
}
iters=iters+1; // increment number of starting values
if (i>0) {
V[0]=0; // load n
V[1]=i;
}
else {
V[0]=0xffffffff;
V[1]=(unsigned int)i;
}
m=0; // clear output array pointer
out[0]=V[0]; // save n value
out[1]=V[1];
m=1; // increment output array pointer
n=1; // number of terms in sequence (not used)
for (j=1; j<100000000; j++) {
T[0]=V[0];
T[1]=V[1];
add64(V,V);
add64(T,V);
add64(W, V); // 3n+c
if (m<1000) {
out[2*m]=V[0]; // save n value
out[2*m+1]=V[1];
m=m+1; // increment output array pointer
}
else {
printf("error: output array not big enough \n");
goto zskip;
}
n=n+1; // increment number of terms in sequence
if ((V[1]&7)==0) {
V[1]=V[1]>>1;
V[1]=V[1]|(V[0]<<31);
V[0]=(int)V[0]>>1;
if (m<1000) {
out[2*m]=V[0]; // save n value
out[2*m+1]=V[1];
m=m+1; // increment output array pointer
}
else {
printf("error: output array not big enough \n");
goto zskip;
}
n=n+1; // increment number of terms in sequence
V[1]=V[1]>>1;
V[1]=V[1]|(V[0]<<31);
V[0]=(int)V[0]>>1;
if (m<1000) {
out[2*m]=V[0]; // save n value
out[2*m+1]=V[1];
m=m+1; // increment output array pointer
}
else {
printf("error: output array not big enough \n");
goto zskip;
}
n=n+1; // increment number of terms in sequence
goto askip;
}
while ((V[1]&1)==0) {
V[1]=V[1]>>1;
V[1]=V[1]|(V[0]<<31);
V[0]=(int)V[0]>>1;
if (m<1000) {
out[2*m]=V[0]; // save n value
out[2*m+1]=V[1];
m=m+1; // increment output array pointer
}
else {
printf("error: output array not big enough \n");
goto zskip;
}
n=n+1; // increment number of terms in sequence
}
}
fprintf(Outfp,"error \n");
goto zskip;
//
// attachment point found, extend sequence backwards using 0-1 sequence vector,
//
// k=s[i];
// while (k!=(k/3)*3) {
// if (k==(k/2)*2) {
// if ((k-c)==((k-c)/3)*3)
// k=(k-c)/3;
// else
// k=k*2;
// }
// else
// k=k*2;
// }
askip:
K[0]=out[2*m-2];
K[1]=out[2*m-3];
S[0]=K[0];
S[1]=K[1];
//
// do adjustment so that the |t2|>|t1| test will pass for cycles
//
if (c==1) {
if ((i!=3)&&(i!=-27)&&(i!=-15)&&(i!=-9)&&(i!=-3)) {
K[1]=K[1]>>1;
K[1]=K[1]|(K[0]<<31);
K[0]=(int)K[0]>>1;
}
}
if (c==5) {
if ((i!=51)&&(i!=63)&&(i!=153)&&(i!=-45)&&(i!=-15)&&(i!=15)&&(i!=-135)) {
K[1]=K[1]>>1;
K[1]=K[1]|(K[0]<<31);
K[0]=(int)K[0]>>1;
}
}
if (c==7) {
if ((i!=-21)&&(i!=-189)&&(i!=-63)&&(i!=21)) {
K[1]=K[1]>>1;
K[1]=K[1]|(K[0]<<31);
K[0]=(int)K[0]>>1;
}
}
if (c==11) {
if ((i!=-33)&&(i!=-27)&&(i!=33)&&(i!=-297)&&(i!=-105)&&(i!=-99)&&(i!=-81)) {
K[1]=K[1]>>1;
K[1]=K[1]|(K[0]<<31);
K[0]=(int)K[0]>>1;
}
}
//
// find odd integer divisible by 3
//
if ((K[0]&0x80000000)!=0) {
T[0]=K[0];
T[1]=K[1];
sub64(Z,T);
div64_32(T,Q,3);
sub64(Z,Q);
}
else
div64_32(K,Q,3);
T[0]=Q[0];
T[1]=Q[1];
add64(Q,Q);
add64(T,Q);
sub64(K,Q);
while ((Q[0]!=0)||(Q[1]!=0)) {
if ((K[1]&1)==0) {
U[0]=W[0];
U[1]=W[1];
sub64(K,U);
if ((U[0]&0x80000000)!=0) {
T[0]=U[0];
T[1]=U[1];
sub64(Z,T);
div64_32(T,Q,3);
sub64(Z,Q);
}
else
div64_32(U,Q,3);
T[0]=Q[0];
T[1]=Q[1];
add64(T,T);
add64(Q,T);
sub64(U,T);
if ((T[0]==0)&&(T[1]==0)) {
K[0]=Q[0];
K[1]=Q[1];
}
else
add64(K,K);
}
else
add64(K,K);
if ((K[0]&0x80000000)!=0) {
T[0]=K[0];
T[1]=K[1];
sub64(Z,T);
div64_32(T,Q,3);
sub64(Z,Q);
}
else
div64_32(K,Q,3);
T[0]=Q[0];
T[1]=Q[1];
add64(Q,Q);
add64(T,Q);
sub64(K,Q);
}
//
// histogram scaled |t2|-|t1| values
//
T[0]=K[0];
T[1]=K[1];
if ((T[0]&0x80000000)!=0)
sub64(Z,T);
if ((T[0]!=0)||((T[1]&0x80000000)!=0))
outrng=outrng+1;
L[0]=T[0];
L[1]=T[1];
sub64(THI,L);
if ((L[0]&0x80000000)!=0) {
THI[0]=T[0];
THI[1]=T[1];
}
temp=i;
if (temp<0)
temp=-temp;
if (limit==10000) {
diff=(temp-(int)T[1])/scale;
diff=diff+offset;
if (diff<0) {
// printf("offset too small, delta=%d \n",diff);
toobig=toobig+1;
goto ajump;
}
if (diff>9999) {
printf("array not big enough, temp=%d \n",diff);
goto zskip;
}
histo[diff]=histo[diff]+1;
}
//
// compare t values
//
ajump:
if ((T[1]>(unsigned int)temp)||(T[0]!=0)) {
//
L[0]=K[0];
L[1]=K[1];
add64(W, L);
n=0; // clear count
while ((L[1]&1)==0) {
L[1]=(L[1]>>1)|(L[0]<<31);
L[0]=(int)L[0]>>1;
n=n+1;
}
L[0]=K[0];
L[1]=K[1];
for (j=0; j<n; j++) {
M[0]=L[0];
M[1]=L[1];
add64(L, L);
add64(M, L);
add64(W, L);
if (limit<=100)
printf("even element in first jump from t1=%d \n",L[1]);
M[0]=L[0];
M[1]=L[1];
add64(M, M);
if ((M[0]==S[0])&&(M[1]==S[1])) {
even1=even1+1;
if (limit<=100) {
printf("attaches to jth even: j=%d, S=%d %d, L=%d %d \n",j,S[0],S[1],K[0],K[1]);
}
}
L[1]=(L[1]>>1)|(L[0]<<31);
L[0]=(int)L[0]>>1;
}
kount=kount+1;
if (limit<=100) {
printf("failure: t2=%d, t1=%d \n",i,K[1]);
fprintf(Outfp,"failure: t2=%d, t1=%d \n",i,K[1]);
}
}
if (limit<=100) {
printf("f=%d, attacher=%d, k=%d, %d, attached to=%d, %d \n",f,i,out[2*m-2],out[2*m-3],K[0],K[1]);
fprintf(Outfp,"f=%d, attacher=%d, k=%d, %d, attached to=%d, %d \n",f,i,out[2*m-2],out[2*m-3],K[0],K[1]);
}
//
// find number of odd elements in each jump to the first attachment point
//
flag=0; // clear k iteration count
count=1; // initialize k value
for (k=0; k<(m-1); k++) {
V[0]=out[2*k]; // load n'
V[1]=out[2*k+1];
if ((V[1]&1)==0) {
V[0]=out[2*k+2]; // load n
V[1]=out[2*k+3];
if ((V[1]&1)==0) {
count=count/2; // k=k/2
if (count>999) {
printf("array not big enough \n");
goto zskip;
}
h[count]=h[count]+1; // histogram k value
if (limit<100)
printf("count=%d \n",count);
count=0; // clear k value
flag=flag+1; // increment k iteration count
}
}
count=count+1; // k=k+1
}
if (flag>999) {
printf("array not big enough \n");
goto zskip;
}
g[flag]=g[flag]+1; // histogram k iteration count
total=total+flag; // increment total number of k values
}
printf("\n");
printf("histogram of number of odd elements in jumps to attachment point, total=%d \n",total);
fprintf(Outfp,"histogram of number of odd elements in jumps to attachment point, total=%d \n",total);
for (i=1; i<20; i++) {
r=(double)h[i];
r=r/(double)total;
printf("i=%d, h[i]=%d, %f \n",i,h[i],r);
fprintf(Outfp," %d %f \n",h[i],r);
}
printf("\n");
fprintf(Outfp,"\n");
printf("histogram of number of jumps to attachment point, total=%d \n",iters);
fprintf(Outfp,"histogram of jumps to attachment point, total=%d \n",iters);
for (i=2; i<20; i++) {
r=(double)g[i];
r=r/(double)iters;
printf("i=%d, g[i]=%d, %f \n",i-1,g[i],r);
fprintf(Outfp," %d %f \n",g[i],r);
}
printf("\n");
fprintf(Outfp,"\n");
if (((THI[0]==0)&&(THI[1]&0x80000000)==0)) {
thi=(int)THI[1];
printf("largest t2 value=%d, portion of t1 range=%e \n",thi,(double)thi/(double)(6*limit+3));
fprintf(Outfp,"largest t2 value=%d, portion of t1 range=%e \n",thi,(double)thi/(double)(6*limit+3));
}
else {
printf("largest t2 value=%#010x %#010x \n",THI[0],THI[1]);
fprintf(Outfp,"largest t2 value=%#010x %#010x \n",THI[0],THI[1]);
}
printf("failure count=%d, ratio=%f \n",kount,(double)(iters-kount)/(double)iters);
fprintf(Outfp,"failure count=%d, ratio=%f \n",kount,(double)(iters-kount)/(double)iters);
printf("failures due to attachments to first even=%d \n",even1);
fprintf(Outfp,"failures due to attachments to first even=%d \n",even1);
printf("big sequence values=%d \n",outrng);
fprintf(Outfp,"big sequence values=%d \n",outrng);
//
if (limit==10000) {
printf("\n");
fprintf(Outfp,"\n");
printf("histogram of scaled differences \n");
fprintf(Outfp,"histogram of scaled differences \n");
for (i=0; i<600; i++) {
printf("i=%d, n=%d \n",i-offset,histo[i]);
fprintf(Outfp,"i=%d, n=%d \n",i-offset,histo[i]);
}
printf("number of values not histogrammed=%d \n",toobig);
fprintf(Outfp,"number of values not histogrammed=%d \n",toobig);
}
zskip:
fclose(Outfp);
return(0);
}