﻿ /*******************************
/******************************************************************************/
/*									      */
/*  COMPUTE 3N+1 SEQUENCE						      */
/*  11/23/09 (dkc)							      */
/*									      */
/*  This C program computes limbs in S from a sequence vector consisting      */
/*  of four peaks and three troughs.					      */
/*									      */
/*  Limit1, limit2, and limit3 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=4;	// usually 5		// count to second trough
unsigned int p2c=4;	// usually 3		// count to second peak
unsigned int o3c=5;	// usually 5		// count to third trough
unsigned int p3c=4;	// usually 4		// count to third peak
unsigned int limit1=0;	// usually 5 (oc1)
unsigned int limit2=0;	// usually 5 (oc2)
unsigned int limit3=0;	// usually 5 (oc3)
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 z,o1,p1,o2,p2,o3,p3;
FILE *Outfp;
Outfp = fopen("out5b.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;
}
if ((o3c>7)||(o3c<1)) {
printf("error: o3c=%d \n",o3c);
goto zskip;
}
if ((p3c>5)||(p3c<1)) {
printf("error: p3c=%d \n",p3c);
goto zskip;
}
set256(T, 1);		      // t=1
set256(H, 4);		      // h=4
//
// begin loop
//
aloop:
if ((I[7]&0x3fffffff)==2)
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);
lshift256(P, R, 7);
lshift256(P, R, 6);
lshift256(P, R, 4);
lshift256(P, R, 3);
shift256(Q, Y, 12);     // (((y/2)*729)/4096
}
if (o1c==5) {
lshift256(P, Q, 7);
lshift256(P, R, 6);
lshift256(P, R, 5);
lshift256(P, R, 4);
lshift256(P, R, 1);
shift256(Q, Y, 10);     // ((y/2)*243)/1024
}
if (o1c==4) {
lshift256(P, Q, 6);
lshift256(P, R, 4);
shift256(Q, Y, 8);      // ((y/2)*81)/256
}
if (o1c==3) {
lshift256(P, Q, 4);
lshift256(P, R, 3);
lshift256(P, R, 1);
shift256(Q, Y, 6);      // ((y/2)*27)/64
}
if (o1c==2) {
lshift256(P, Q, 3);
shift256(Q, Y, 4);      // ((y/2)*9)/16
}
if (o1c==1) {
lshift256(P, Q, 1);
shift256(Q, Y, 2);      // ((y/2)*3)/4
}
if ((Y[7]&1)==0)	      // if o1 is even, continue
goto bskip;
//
// check if 2**x divides o1+1 (for p1)
//
//
set256(P, 1);
if (p1c==5)
if (p1c==4)
if (p1c==3)
if (p1c==2)
if (p1c==1)
if ((P[7]&mask)!=0)	      // if ((o1+1)!=((o1+1)/*2**x)*2**x) continue
goto bskip;
//
// compute p1 from o1
//
set256(P, 1);
if (p1c==5) {
shift256(P, P, 5);      // (o1+1)/32
lshift256(P, Q, 7);
lshift256(P, R, 6);
lshift256(P, R, 5);
lshift256(P, R, 4);
lshift256(P, R, 1);
}
if (p1c==4) {
shift256(P, P, 4);      // (o1+1)/16
lshift256(P, Q, 6);
lshift256(P, R, 4);
}
if (p1c==3) {
shift256(P, P, 3);      // (o1+1)/8
lshift256(P, Q, 4);
lshift256(P, R, 3);
lshift256(P, R, 1);
}
if (p1c==2) {
shift256(P, P, 2);      // (o1+1)/4
lshift256(P, Q, 3);
}
if (p1c==1) {
shift256(P, P, 1);      // (o1+1)/2
lshift256(P, Q, 1);
}
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);
lshift256(N, Q, 7);
lshift256(N, Q, 3);
lshift256(N, Q, 1);
shift256(P, Q, 14);     // j=j/16384
}
if (o2c==6) {
lshift256(N, P, 9);
lshift256(N, Q, 7);
lshift256(N, Q, 6);
lshift256(N, Q, 4);
lshift256(N, Q, 3);
shift256(P, Q, 12);     // j=j/4096
}
if (o2c==5) {
lshift256(N, P, 7);
lshift256(N, Q, 6);
lshift256(N, Q, 5);
lshift256(N, Q, 4);
lshift256(N, Q, 1);
shift256(P, Q, 10);     // j=j/1024
}
if (o2c==4) {
lshift256(N, P, 6);
lshift256(N, Q, 4);
shift256(P, Q, 8);      // j=j/256
}
if (o2c==3) {
lshift256(N, P, 4);
lshift256(N, Q, 3);
lshift256(N, Q, 1);
shift256(P, Q, 6);      // j=j/64
}
if (o2c==2) {
lshift256(N, P, 3);
shift256(P, Q, 4);      // j=j/16
}
if (o2c==1) {
lshift256(N, P, 1);
shift256(P, Q, 2);      // j=j/4
}
set256(Z, 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);
if (p2c==5)
if (p2c==4)
if (p2c==3)
if (p2c==2)
if (p2c==1)
if ((N[7]&mask)!=0)	   // if (j!=(j/2**y)*2**y) continue
goto bskip;
//
// compute p2 from o2
//
set256(N, 1);
if (p2c==5) {
shift256(N, N, 5);   // (o2+1)/32
lshift256(N, P, 7);
lshift256(N, Q, 6);
lshift256(N, Q, 5);
lshift256(N, Q, 4);
lshift256(N, Q, 1);
}
if (p2c==4) {
shift256(N, N, 4);   // (o2+1)/16
lshift256(N, P, 6);
lshift256(N, Q, 4);
}
if (p2c==3) {
shift256(N, N, 3);   // (o2+1)/8
lshift256(N, P, 4);
lshift256(N, Q, 3);
lshift256(N, Q, 1);
}
if (p2c==2) {
shift256(N, N, 2);   // (o2+1)/4
lshift256(N, P, 3);
}
if (p2c==1) {
shift256(N, N, 1);   // (o2+1)/2
lshift256(N, P, 1);
}
set256(W, 1);
sub256(P, W);
if ((W[7]&1)==0)	   // if p2 is even, continue
goto bskip;
//
// compute o3 from p2
//
copy256(W, N);
if (o3c==7) {
lshift256(N, P, 11);
lshift256(N, Q, 7);
lshift256(N, Q, 3);
lshift256(N, Q, 1);
shift256(P, Q, 14);     // j=j/16384
}
if (o3c==6) {
lshift256(N, P, 9);
lshift256(N, Q, 7);
lshift256(N, Q, 6);
lshift256(N, Q, 4);
lshift256(N, Q, 3);
shift256(P, Q, 12);     // j=j/4096
}
if (o3c==5) {
lshift256(N, P, 7);
lshift256(N, Q, 6);
lshift256(N, Q, 5);
lshift256(N, Q, 4);
lshift256(N, Q, 1);
shift256(P, Q, 10);     // j=j/1024
}
if (o3c==4) {
lshift256(N, P, 6);
lshift256(N, Q, 4);
shift256(P, Q, 8);      // j=j/256
}
if (o3c==3) {
lshift256(N, P, 4);
lshift256(N, Q, 3);
lshift256(N, Q, 1);
shift256(P, Q, 6);      // j=j/64
}
if (o3c==2) {
lshift256(N, P, 3);
shift256(P, Q, 4);      // j=j/16
}
if (o3c==1) {
lshift256(N, P, 1);
shift256(P, Q, 2);      // j=j/4
}
set256(K, 1);
if ((K[7]&1)==0)	      // if o3 is even, continue
goto bskip;
//
// check if 2**z divides o3+1 (for p3)
//
set256(N, 1);
if (p3c==5)
if (p3c==4)
if (p3c==3)
if (p3c==2)
if (p3c==1)
if ((N[7]&mask)!=0)	   // if (j!=(j/2**x)*2**x) continue
goto bskip;
//
// compute p3 from o3
//
set256(N, 1);
if (p3c==5) {
shift256(N, N, 5);   // (o3+1)/32
lshift256(N, P, 7);
lshift256(N, Q, 6);
lshift256(N, Q, 5);
lshift256(N, Q, 4);
lshift256(N, Q, 1);
}
if (p3c==4) {
shift256(N, N, 4);   // (o3+1)/16
lshift256(N, P, 6);
lshift256(N, Q, 4);
}
if (p3c==3) {
shift256(N, N, 3);   // (o3+1)/8
lshift256(N, P, 4);
lshift256(N, Q, 3);
lshift256(N, Q, 1);
}
if (p3c==2) {
shift256(N, N, 2);   // (o3+1)/4
lshift256(N, P, 3);
}
if (p3c==1) {
shift256(N, N, 1);   // (o3+1)/2
lshift256(N, P, 1);
}
set256(L, 1);
sub256(P, L);
if ((L[7]&1)==0)	   // if p3 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
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;
//
copy256(K, R);
sub256(Q, R);	 // (order/12)-o3
// if ((R[0]&0x80000000)==0)  // if (o3<(order/12)) continue
//    goto bskip;
//
lshift256(Q, Q, 2);	 // order/3
copy256(L, R);
sub256(Q, R);	 // (order/3)-p3
if ((R[0]&0x80000000)==0)  // if (p3<(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);
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)
shift256(U, U, 2); // y=y/4
copy256(U, V);
sub256(Y, V);	 // o1-y
}
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);
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 dskip;
shift256(U, U, 2); // p1=p1/4
copy256(U, V);
sub256(Z, V);	 // o2-p1
}
dskip:
copy256(W, U);	 // p2
copy256(U, V);
sub256(K, V);	 // o3-p2
count=0;
while ((V[0]|V[1]|V[2]|V[3]|V[4]|V[5]|V[6]|V[7])!=0) {
copy256(U, S);
if ((U[7]&7)==0) { // if 8 divides p2, continue
goto bskip;
}
if ((U[7]&3)!=0) { // if 4 doesn't divide p2, continue
goto bskip;
}
count=count+1;
if (count>10) {
printf("error: count too big \n");
goto zskip;
}
if (count>limit3)
goto cskip;
shift256(U, U, 2); // p2=p2/4
copy256(U, V);
sub256(K, V);	 // o3-p2
}
cskip:
/* printf("\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("   o3=%#10x %#10x %#10x %#10x \n",K[0],K[1],K[2],K[3]);
printf("      %#10x %#10x %#10x %#10x \n",K[4],K[5],K[6],K[7]);
printf("   p3=%#10x %#10x %#10x %#10x \n",L[0],L[1],L[2],L[3]);
printf("      %#10x %#10x %#10x %#10x \n",L[4],L[5],L[6],L[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,"   o3=%#10x %#10x %#10x %#10x \n",K[0],K[1],K[2],K[3]);
fprintf(Outfp,"      %#10x %#10x %#10x %#10x \n",K[4],K[5],K[6],K[7]);
fprintf(Outfp,"   p3=%#10x %#10x %#10x %#10x \n",L[0],L[1],L[2],L[3]);
fprintf(Outfp,"      %#10x %#10x %#10x %#10x \n",L[4],L[5],L[6],L[7]);
fprintf(Outfp,"\n"); */
if (test==0)
goto eskip;
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");
}
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");
}
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");
}
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");
}
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");
//
printf("\n");
fprintf(Outfp,"\n");
o3=K[7];
printf("o3=%d \n",o3);
fprintf(Outfp,"o3=%d \n",o3);
z=W[7];
printf(" %d",z);
fprintf(Outfp," %d",z);
roll=1;
count=0;
while (z!=o3) {	 // find o3
z=3*z+1;
count=count+1;
if (count<=3*o3c) {
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<=3*o3c) {
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");
count=0;
p3=L[7];
printf("p3=%d \n",p3);
fprintf(Outfp,"p3=%d \n",p3);
z=K[7];
printf(" %d",z);
fprintf(Outfp," %d",z);
roll=1;
while (z!=p3) {	 // find p3
z=3*z+1;
count=count+1;
if (count<=2*p3c) {
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*p3c) {
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");
eskip:
iters=iters+1;	 // check for enough solutions
if (iters>20)
goto zskip;
bskip:
copy256(I, P);