/******************************************************************************/
/* */
/* COMPUTE 3N+1 SEQUENCE */
/* 11/23/09 (dkc) */
/* */
/* This C program computes limbs in S from a sequence vector consisting */
/* of three peaks and two troughs. */
/* */
/* Limit1 and limit2 should be set to 0 to allow "tumbles". I should be */
/* set to order/2+2 and J should be set to order. */
/* */
/******************************************************************************/
#include <stdio.h>
#include <math.h>
unsigned int o1c=5; // usually 5 // count to first trough
unsigned int p1c=4; // usually 4 // count to first peak
unsigned int o2c=5; // usually 5 // count to second trough
unsigned int p2c=3; // usually 3 // count to second peak
unsigned int limit1=0; // usually 5 (oc1)
unsigned int limit2=0; // usually 5 (oc2)
unsigned int test=1; // test flag (for small solutions)
unsigned int I[8]={0x00000000, 0x00000000, 0x00000000, 0x00000000,
0x00000000, 0x00000000, 0x00000000, 0x0c000002}; // min i
unsigned int J[8]={0x00000000, 0x00000000, 0x00000000, 0x00000000,
0x00000000, 0x00000000, 0x00000000, 0x18000000}; // max i
//
void add256(unsigned int *a, unsigned int *b);
void sub256(unsigned int *c, unsigned int *d);
void shift256(unsigned int *a, unsigned int *b, unsigned int shift);
void lshift256(unsigned int *a, unsigned int *b, unsigned int shift);
void copy256(unsigned int *a, unsigned int *b);
void set256(unsigned int *a, unsigned int b);
int main () {
//
unsigned int N[8],O[8],P[8],Q[8],R[8],S[8],T[8],U[8],V[8],X[8],W[8],Y[8],Z[8];
unsigned int count,iters,roll,mask,H[8],M[8];
unsigned int z,o1,p1,o2,p2,saveo1,saveo2;
FILE *Outfp;
Outfp = fopen("out5a.dat","w");
iters=0;
if ((o1c>6)||(o1c<1)) {
printf("error: o1c=%d \n",o1c);
goto zskip;
}
if ((p1c>5)||(p1c<1)) {
printf("error: p1c=%d \n",p1c);
goto zskip;
}
if ((o2c>7)||(o2c<1)) {
printf("error: o2c=%d \n",o2c);
goto zskip;
}
if ((p2c>5)||(p2c<1)) {
printf("error: p2c=%d \n",p2c);
goto zskip;
}
set256(T, 1); // t=1
set256(H, 4); // h=4
//
// begin loop
//
aloop:
if ((I[7]&0x7ffffff)==0)
printf("i=%#10x \n",I[7]);
copy256(I, O); // y=i
shift256(O, P, 1); // y/2
//
// compute o1 from y
//
if (o1c==6) {
lshift256(P, Q, 9);
add256(P, Q);
lshift256(P, R, 7);
add256(R, Q);
lshift256(P, R, 6);
add256(R, Q);
lshift256(P, R, 4);
add256(R, Q);
lshift256(P, R, 3);
add256(R, Q); // (y/2)*729
shift256(Q, Y, 12); // (((y/2)*729)/4096
}
if (o1c==5) {
lshift256(P, Q, 7);
add256(P, Q);
lshift256(P, R, 6);
add256(R, Q);
lshift256(P, R, 5);
add256(R, Q);
lshift256(P, R, 4);
add256(R, Q);
lshift256(P, R, 1);
add256(R, Q); // (y/2)*243
shift256(Q, Y, 10); // ((y/2)*243)/1024
}
if (o1c==4) {
lshift256(P, Q, 6);
add256(P, Q);
lshift256(P, R, 4);
add256(R, Q); // (y/2)*81
shift256(Q, Y, 8); // ((y/2)*81)/256
}
if (o1c==3) {
lshift256(P, Q, 4);
add256(P, Q);
lshift256(P, R, 3);
add256(R, Q);
lshift256(P, R, 1);
add256(R, Q); // (y/2)*27
shift256(Q, Y, 6); // ((y/2)*27)/64
}
if (o1c==2) {
lshift256(P, Q, 3);
add256(P, Q); // (y/2)*9
shift256(Q, Y, 4); // ((y/2)*9)/16
}
if (o1c==1) {
lshift256(P, Q, 1);
add256(P, Q); // (y/2)*3
shift256(Q, Y, 2); // ((y/2)*3)/4
}
add256(T, Y); // o1
if ((Y[7]&1)==0) // if o1 is even, continue
goto bskip;
//
// check if 2**x divides o1+1 (for p1)
//
//
set256(P, 1);
add256(Y, P); // o1+1
if (p1c==5)
mask=0x1f;
if (p1c==4)
mask=0xf;
if (p1c==3)
mask=0x7;
if (p1c==2)
mask=0x3;
if (p1c==1)
mask=0x1;
if ((P[7]&mask)!=0) // if ((o1+1)!=((o1+1)/*2**x)*2**x) continue
goto bskip;
//
// compute p1 from o1
//
set256(P, 1);
add256(Y, P); // o1+1
if (p1c==5) {
shift256(P, P, 5); // (o1+1)/32
lshift256(P, Q, 7);
add256(P, Q);
lshift256(P, R, 6);
add256(R, Q);
lshift256(P, R, 5);
add256(R, Q);
lshift256(P, R, 4);
add256(R, Q);
lshift256(P, R, 1);
add256(R, Q); // ((o1+1)/32)*243
}
if (p1c==4) {
shift256(P, P, 4); // (o1+1)/16
lshift256(P, Q, 6);
add256(P, Q);
lshift256(P, R, 4);
add256(R, Q); // ((o1+1)/16)*81
}
if (p1c==3) {
shift256(P, P, 3); // (o1+1)/8
lshift256(P, Q, 4);
add256(P, Q);
lshift256(P, R, 3);
add256(R, Q);
lshift256(P, R, 1);
add256(R, Q); // ((o1+1)/8)*27
}
if (p1c==2) {
shift256(P, P, 2); // (o1+1)/4
lshift256(P, Q, 3);
add256(P, Q); // ((o1+1)/4)*9
}
if (p1c==1) {
shift256(P, P, 1); // (o1+1)/2
lshift256(P, Q, 1);
add256(P, Q); // ((o1+1)/2)*3
}
set256(X, 1);
sub256(Q, X);
if ((X[7]&1)==0) // if p1 is even, continue
goto bskip;
//
// compute o2 from p1
//
copy256(X, N);
if (o2c==7) {
lshift256(N, P, 11);
add256(N, P);
lshift256(N, Q, 7);
add256(Q, P);
lshift256(N, Q, 3);
add256(Q, P);
lshift256(N, Q, 1);
add256(Q, P); // j=j*2187
shift256(P, Q, 14); // j=j/16384
}
if (o2c==6) {
lshift256(N, P, 9);
add256(N, P);
lshift256(N, Q, 7);
add256(Q, P);
lshift256(N, Q, 6);
add256(Q, P);
lshift256(N, Q, 4);
add256(Q, P);
lshift256(N, Q, 3);
add256(Q, P); // j=j*729
shift256(P, Q, 12); // j=j/4096
}
if (o2c==5) {
lshift256(N, P, 7);
add256(N, P);
lshift256(N, Q, 6);
add256(Q, P);
lshift256(N, Q, 5);
add256(Q, P);
lshift256(N, Q, 4);
add256(Q, P);
lshift256(N, Q, 1);
add256(Q, P); // j=j*243
shift256(P, Q, 10); // j=j/1024
}
if (o2c==4) {
lshift256(N, P, 6);
add256(N, P);
lshift256(N, Q, 4);
add256(Q, P); // j=j*81
shift256(P, Q, 8); // j=j/256
}
if (o2c==3) {
lshift256(N, P, 4);
add256(N, P);
lshift256(N, Q, 3);
add256(Q, P);
lshift256(N, Q, 1);
add256(Q, P); // j=j*27
shift256(P, Q, 6); // j=j/64
}
if (o2c==2) {
lshift256(N, P, 3);
add256(N, P); // j=j*9
shift256(P, Q, 4); // j=j/16
}
if (o2c==1) {
lshift256(N, P, 1);
add256(N, P); // j=j*3
shift256(P, Q, 2); // j=j/4
}
set256(Z, 1);
add256(Q, Z); // o2=j+1
if ((Z[7]&1)==0) // if o2 is even, continue
goto bskip;
//
// check if 2**y divides o2+1 (for p2)
//
set256(N, 1);
add256(Z, N);
if (p2c==5)
mask=0x1f;
if (p2c==4)
mask=0xf;
if (p2c==3)
mask=0x7;
if (p2c==2)
mask=0x3;
if (p2c==1)
mask=0x1;
if ((N[7]&mask)!=0) // if (j!=(j/2**y)*2**y) continue
goto bskip;
//
// compute p2 from o2
//
set256(N, 1);
add256(Z, N); // o2+1
if (p2c==5) {
shift256(N, N, 5); // (o2+1)/32
lshift256(N, P, 7);
add256(N, P);
lshift256(N, Q, 6);
add256(Q, P);
lshift256(N, Q, 5);
add256(Q, P);
lshift256(N, Q, 4);
add256(Q, P);
lshift256(N, Q, 1);
add256(Q, P); // p2=((o2+1)/32)*243-1
}
if (p2c==4) {
shift256(N, N, 4); // (o2+1)/16
lshift256(N, P, 6);
add256(N, P);
lshift256(N, Q, 4);
add256(Q, P); // p2=((o2+1)/16)*81-1
}
if (p2c==3) {
shift256(N, N, 3); // (o2+1)/8
lshift256(N, P, 4);
add256(N, P);
lshift256(N, Q, 3);
add256(Q, P);
lshift256(N, Q, 1);
add256(Q, P); // p2=((o2+1)/8)*27-1
}
if (p2c==2) {
shift256(N, N, 2); // (o2+1)/4
lshift256(N, P, 3);
add256(N, P); // p2=((o2+1)/4)*9-1
}
if (p2c==1) {
shift256(N, N, 1); // (o2+1)/2
lshift256(N, P, 1);
add256(N, P); // p2=((o2+1)/2)*3-1
}
set256(W, 1);
sub256(P, W);
if ((W[7]&1)==0) // if p2 is even, continue
goto bskip;
//
// compute order
//
set256(P, 6); // order=6
set256(M, 1); // order'=1
copy256(P, Q);
sub256(O, Q); // y-order
while ((Q[0]&0x80000000)==0) { // while (order<y) order=order*2
add256(P, P); // order=order*2
add256(M, M); // order'=order'*2
copy256(P, Q);
sub256(O, Q); // y-order
}
//
shift256(M, Q, 1); // order/12
copy256(Y, R);
sub256(Q, R); // (order/12)-o1
// if ((R[0]&0x80000000)==0) // if (o1<(order/12)) continue
// goto bskip;
//
copy256(Z, R);
sub256(Q, R); // (order/12)-o2
// if ((R[0]&0x80000000)==0) // if (o2<(order/12)) continue
// goto bskip;
//
lshift256(Q, Q, 2); // order/3
copy256(W, R);
sub256(Q, R); // (order/3)-p2
if ((R[0]&0x80000000)==0) // if (p2<(order/3)) continue
goto bskip;
//
shift256(O, U, 1); // y=y/2
copy256(U, V);
sub256(Y, V); // o1-y
count=0;
while ((V[0]|V[1]|V[2]|V[3]|V[4]|V[5]|V[6]|V[7])!=0) {
copy256(U, S);
add256(U, U);
add256(S, U);
add256(T, U); // y=3*y+1
if ((U[7]&7)==0) { // if 8 divides y, continue
goto bskip;
}
if ((U[7]&3)!=0) { // if 4 doesn't divide y, continue
goto bskip;
}
count=count+1;
if (count>10) {
printf("error: count too big \n");
goto zskip;
}
if (count>limit1)
goto askip;
shift256(U, U, 2); // y=y/4
copy256(U, V);
sub256(Y, V); // o1-y
}
askip:
copy256(X, U); // p1
copy256(U, V);
sub256(Z, V); // o2-p1
count=0;
while ((V[0]|V[1]|V[2]|V[3]|V[4]|V[5]|V[6]|V[7])!=0) {
copy256(U, S);
add256(U, U);
add256(S, U);
add256(T, U); // p1=3*p1+1
if ((U[7]&7)==0) { // if 8 divides p1, continue
goto bskip;
}
if ((U[7]&3)!=0) { // if 4 doesn't divide p1, continue
goto bskip;
}
count=count+1;
if (count>10) {
printf("error: count too big \n");
goto zskip;
}
if (count>limit2)
goto cskip;
shift256(U, U, 2); // p1=p1/4
copy256(U, V);
sub256(Z, V); // o2-p1
}
cskip:
/* printf("\n");
fprintf(Outfp,"\n");
printf("order=%#10x %#10x %#10x %#10x \n",P[0],P[1],P[2],P[3]);
printf(" %#10x %#10x %#10x %#10x \n",P[4],P[5],P[6],P[7]);
printf(" y=%#10x %#10x %#10x %#10x \n",O[0],O[1],O[2],O[3]);
printf(" %#10x %#10x %#10x %#10x \n",O[4],O[5],O[6],O[7]);
printf(" o1=%#10x %#10x %#10x %#10x \n",Y[0],Y[1],Y[2],Y[3]);
printf(" %#10x %#10x %#10x %#10x \n",Y[4],Y[5],Y[6],Y[7]);
printf(" p1=%#10x %#10x %#10x %#10x \n",X[0],X[1],X[2],X[3]);
printf(" %#10x %#10x %#10x %#10x \n",X[4],X[5],X[6],X[7]);
printf(" o2=%#10x %#10x %#10x %#10x \n",Z[0],Z[1],Z[2],Z[3]);
printf(" %#10x %#10x %#10x %#10x \n",Z[4],Z[5],Z[6],Z[7]);
printf(" p2=%#10x %#10x %#10x %#10x \n",W[0],W[1],W[2],W[3]);
printf(" %#10x %#10x %#10x %#10x \n",W[4],W[5],W[6],W[7]);
printf("\n");
fprintf(Outfp,"order=%#10x %#10x %#10x %#10x \n",P[0],P[1],P[2],P[3]);
fprintf(Outfp," %#10x %#10x %#10x %#10x \n",P[4],P[5],P[6],P[7]);
fprintf(Outfp," y=%#10x %#10x %#10x %#10x \n",O[0],O[1],O[2],O[3]);
fprintf(Outfp," %#10x %#10x %#10x %#10x \n",O[4],O[5],O[6],O[7]);
fprintf(Outfp," o1=%#10x %#10x %#10x %#10x \n",Y[0],Y[1],Y[2],Y[3]);
fprintf(Outfp," %#10x %#10x %#10x %#10x \n",Y[4],Y[5],Y[6],Y[7]);
fprintf(Outfp," p1=%#10x %#10x %#10x %#10x \n",X[0],X[1],X[2],X[3]);
fprintf(Outfp," %#10x %#10x %#10x %#10x \n",X[4],X[5],X[6],X[7]);
fprintf(Outfp," o2=%#10x %#10x %#10x %#10x \n",Z[0],Z[1],Z[2],Z[3]);
fprintf(Outfp," %#10x %#10x %#10x %#10x \n",Z[4],Z[5],Z[6],Z[7]);
fprintf(Outfp," p2=%#10x %#10x %#10x %#10x \n",W[0],W[1],W[2],W[3]);
fprintf(Outfp," %#10x %#10x %#10x %#10x \n",W[4],W[5],W[6],W[7]);
fprintf(Outfp,"\n"); */
if (test==0)
goto dskip;
printf("\n");
fprintf(Outfp,"\n");
printf("\n");
fprintf(Outfp,"\n");
o1=Y[7];
printf("o1=%d \n",o1);
fprintf(Outfp,"o1=%d \n",o1);
z=O[7]/2;
printf(" %d",z);
fprintf(Outfp," %d",z);
count=0;
roll=1;
while (z!=o1) { // find o1
z=3*z+1;
count=count+1;
if (count<=3*o1c) {
printf(" %d",z);
fprintf(Outfp," %d",z);
roll=roll+1;
if ((roll&7)==0)
fprintf(Outfp,"\n");
saveo1=z;
}
else
break;
while (z==(z/2)*2) {
z=z/2;
count=count+1;
if (count<=3*o1c) {
printf(" %d",z);
fprintf(Outfp," %d",z);
roll=roll+1;
if ((roll&7)==0)
fprintf(Outfp,"\n");
saveo1=z;
}
else
break;
}
}
printf("\n");
fprintf(Outfp,"\n");
//
printf("\n");
fprintf(Outfp,"\n");
count=0;
p1=X[7];
printf("p1=%d \n",p1);
fprintf(Outfp,"p1=%d \n",p1);
z=Y[7];
printf(" %d",z);
fprintf(Outfp," %d",z);
roll=1;
while (z!=p1) { // find p1
z=3*z+1;
count=count+1;
if (count<=2*p1c) {
printf(" %d",z);
fprintf(Outfp," %d",z);
roll=roll+1;
if ((roll&7)==0)
fprintf(Outfp,"\n");
}
else
break;
while (z==(z/2)*2) {
z=z/2;
count=count+1;
if (count<=2*p1c) {
printf(" %d",z);
fprintf(Outfp," %d",z);
roll=roll+1;
if ((roll&7)==0)
fprintf(Outfp,"\n");
}
else
break;
}
}
printf("\n");
fprintf(Outfp,"\n");
//
printf("\n");
fprintf(Outfp,"\n");
o2=Z[7];
printf("o2=%d \n",o2);
fprintf(Outfp,"o2=%d \n",o2);
z=p1;
printf(" %d",z);
fprintf(Outfp," %d",z);
roll=1;
count=0;
while (z!=o2) { // find o2
z=3*z+1;
count=count+1;
if (count<=3*o2c) {
printf(" %d",z);
fprintf(Outfp," %d",z);
roll=roll+1;
if ((roll&7)==0)
fprintf(Outfp,"\n");
saveo2=z;
}
else
break;
while (z==(z/2)*2) {
z=z/2;
count=count+1;
if (count<=3*o2c) {
printf(" %d",z);
fprintf(Outfp," %d",z);
roll=roll+1;
if ((roll&7)==0)
fprintf(Outfp,"\n");
saveo2=z;
}
else
break;
}
}
printf("\n");
fprintf(Outfp,"\n");
//
printf("\n");
fprintf(Outfp,"\n");
count=0;
p2=W[7];
printf("p2=%d \n",p2);
fprintf(Outfp,"p2=%d \n",p2);
z=Z[7];
printf(" %d",z);
fprintf(Outfp," %d",z);
roll=1;
while (z!=p2) { // find p2
z=3*z+1;
count=count+1;
if (count<=2*p2c) {
printf(" %d",z);
fprintf(Outfp," %d",z);
roll=roll+1;
if ((roll&7)==0)
fprintf(Outfp,"\n");
}
else
break;
while (z==(z/2)*2) {
z=z/2;
count=count+1;
if (count<=2*p2c) {
printf(" %d",z);
fprintf(Outfp," %d",z);
roll=roll+1;
if ((roll&7)==0)
fprintf(Outfp,"\n");
}
else
break;
}
}
printf("\n");
fprintf(Outfp,"\n");
if ((saveo1-o1)<10) {
printf("case1: o1=%d, delta=%d \n",o1,saveo1-o1);
fprintf(Outfp,"case1: o1=%d, delta=%d \n",o1,saveo1-o1);
}
else {
if ((o1-(saveo1/6))<10) {
printf("case2: o1=%d, delta=%d \n",o1,o1-(saveo1/6));
fprintf(Outfp,"case2: o1=%d, delta=%d \n",o1,o1-(saveo1/6));
}
else {
if ((saveo1*6)-o1<10) {
printf("case3: o1=%d, delta=%d \n",o1,saveo1*6-o1);
fprintf(Outfp,"case3: o1=%d, delta=%d \n",o1,saveo1*6-o1);
}
else {
if (((3*o1+1)-(saveo1/12))<10) {
printf("case4: o1=%d, delta=%d \n",o1,(3*o1+1-saveo1/12));
fprintf(Outfp,"case4: o1=%d, delta=%d \n",o1,(3*o1+1-saveo1/12));
}
else {
printf("case not covered: o1=%d, saveo1=%d \n",o1,saveo1);
fprintf(Outfp,"case not covered: o1=%d, saveo1=%d \n",o1,saveo1);
}
}
}
}
if ((saveo2-o2)<10) {
printf("case1: o2=%d, delta=%d \n",o2,saveo2-o2);
fprintf(Outfp,"case1: o2=%d, delta=%d \n",o2,saveo2-o2);
}
else {
if ((o2-(saveo2/6))<10) {
printf("case2: o2=%d, delta=%d \n",o2,o2-(saveo2/6));
fprintf(Outfp,"case2: o2=%d, delta=%d \n",o2,o2-(saveo2/6));
}
else {
if ((saveo2*6)-o2<10) {
printf("case3: o2=%d, delta=%d \n",o2,saveo2*6-o2);
fprintf(Outfp,"case3: o2=%d, delta=%d \n",o2,saveo2*6-o2);
}
else {
if (((3*o2+1)-(saveo2/12))<10) {
printf("case4: o2=%d, delta=%d \n",o2,3*o2+1-(saveo2/12));
fprintf(Outfp,"case4: o2=%d, delta=%d \n",o2,3*o2+1-(saveo2/12));
}
else {
printf("case not covered: o2=%d, saveo2=%d \n",o2,saveo2);
fprintf(Outfp,"case not covered: o2=%d, saveo2=%d \n",o2,saveo2);
}
}
}
}
dskip:
iters=iters+1; // check for enough solutions
if (iters>100)
goto zskip;
bskip:
copy256(I, P);
add256(H, I);
sub256(J, P);
if ((P[0]&0x80000000)==0)
goto aloop;
zskip:
fclose(Outfp);
return(0);
}