/******************************************************************************/
/* */
/* COMPUTE 3N+1 OR 3N-1 SEQUENCE */
/* 02/13/10 (dkc) */
/* */
/* This C program computes general 3n+1 or 3n-1 sequences. The third-to-last*/
/* element of a sub-sequence must be odd and the second-to-last element of */
/* a sub-sequence must be divisible by 8. The third-to-last element must be */
/* larger than the order divided by 6. Certain lengths of such sequences */
/* are not permissible. */
/* */
/******************************************************************************/
#include <stdio.h>
#include <math.h>
void div64_32(unsigned int *a, unsigned int *b, unsigned int);
void add64(unsigned int *a, unsigned int *b);
void sub64(unsigned int *c, unsigned int *d);
void shift(unsigned int *a, unsigned int *b, unsigned int shift);
int main () {
//
int c=1; // c
//
unsigned int h,i,j,n,np,order,count; // indices
unsigned int he[2048]; // histogram
unsigned int y[178]={0,3,6,8,11, // I (invalid lengths)
13,16,19,21,24, // I
26,29,32,34,37, // I
39,42,44,47,50, // J
52,55,57,60,63, // J
65,68,70,73,75, // K
78,81,83,86,88, // K
91,94,96,99,101, // K
104,106,109,112,114, // L
117,119,122,125,127, // L
130,132,135,138,140, // L
143,145,148,150,153, // M
156,158,161,163,166, // M
171,174,176,179, // M
184,187,189,192,194, // I'
197,200,202,205,207, // I'
210,212,215,218,220, // J'
223,225,228,231,233, // J'
236,238,241,243,246, // K'
249,251,254,256,259, // K'
262,264,267,269,272, // K'
275,277,280,282,285, // K'
287,290,293,295,298, // L'
300,303,306,308,311, // L'
313,316,318,321,324, // M'
326,329,331,334,337, // M'
339,342,344,347,349, // N
352,355,357,360,362, // N
365,368,370,373,375, // N
378,380,383,386,388, // O
391,393,396,399,401, // O
404,406,409,412,414, // O
417,419,422,424,427, // P
430,432,435,437,440, // P
443,445,448,450,453, // P
458,461,463,466}; //
unsigned int z[75]={1,9, // I (first elements of length pairs which
14,22, // I when decremented by 3 give invalid
27,35, // I lengths)
40,45, // J
53,58, // J
66,71,76, // K
84,89, // K
97,102, // K
107,115, // L
120,128, // L
133,141, // L
146,151, // M
159,164, // M
172,177, // M
182,190, // I'
195,203, // I'
208,213, // J'
221,226, // J'
234,239,244, // K'
252,257, // K'
265,270, // K'
278,283, // K'
288,296, // L'
301,309, // L'
314,319, // M'
327,332, // M'
340,345,350, // N
358,363, // N
371,376, // N
384,389, // O
394,402, // O
407,415, // O
420,428, // P
433,438, // P
446,451, // P
459,464}; //
unsigned int V[2],W[2],Y[2],M[2],N[2],O[2],P[2];
FILE *Outfp;
Outfp = fopen("out2.dat","w");
if (c==1) {
W[0]=0; // load 1
W[1]=1;
}
else {
W[0]=0xffffffff; // load -1
W[1]=0xffffffff;
}
for (i=0; i<2048; i++) // clear histogram of lengths
he[i]=0; // for limbs ending in E
//
// DEAD LIMBS (ending in an even natural number)
//
count=0;
for (i=9999999; i>3; i-=6) {
if (count<200) {
printf("\n");
fprintf(Outfp,"\n");
}
V[0]=0; // load k
V[1]=i;
n=1; // element count
// if (count<200)
// printf("k=%d %d, n=%d \n",V[0],V[1],n);
order=3; // set order
while (order<i)
order=order*2;
O[0]=0;
O[1]=order;
//
// begin loop
//
aloop:
if (V[0]>0x10000000)
goto cskip;
P[0]=V[0];
P[1]=V[1];
Y[0]=V[0];
Y[1]=V[1];
add64(V, Y);
add64(Y, V);
add64(W, V); // 3n+c
n=n+1; // increment element count
// if (count<200)
// printf("k=%d %d, n=%d \n",V[0],V[1],n);
//
Y[0]=O[0];
Y[1]=O[1];
sub64(V, Y); // n-order
while ((Y[0]&0x80000000)==0) {
add64(O, O);
Y[0]=O[0];
Y[1]=O[1];
sub64(V, Y); // n-order
}
//
shift(V, Y, 3);
add64(Y, Y);
add64(Y, Y);
add64(Y, Y);
sub64(V, Y);
if ((Y[0]==0)&&(Y[1]==0)) { // 8 divides 3n+c
np=n+1; // increment element count
M[0]=0;
M[1]=i;
shift(O, Y, 1);
N[0]=Y[0];
N[1]=Y[1];
sub64(M, N);
while ((N[0]&0x80000000)!=0) {
add64(M, M);
np=np+1;
N[0]=Y[0];
N[1]=Y[1];
sub64(M, N);
}
if (count<200) {
shift(V, Y, 1);
printf("limb=(%#4x %#10x, %#4x %#10x), order=%#4x %#10x, n=%d",M[0],M[1],Y[0],Y[1],O[0],O[1],np);
fprintf(Outfp,"limb=(%#4x %#10x, %#4x %#10x), order=%#4x %#10x, n=%d",M[0],M[1],Y[0],Y[1],O[0],O[1],np);
}
div64_32(O, Y, 6); // order/6
sub64(P, Y);
if ((Y[0]&0x80000000)==0) {
Y[0]=V[0];
Y[1]=V[1];
add64(V, Y);
sub64(O, Y);
if ((Y[0]&0x80000000)==0) {
printf("error: Next-to-last element does not determine order.\n");
goto zskip;
}
if (count<200) {
printf(", gt \n");
fprintf(Outfp,", gt \n");
}
for (j=0; j<75; j++) {
if (z[j]==np)
goto bskip;
}
if (np<471)
np=np-3;
bskip: if (np<2048)
he[np]=he[np]+1;
else {
printf("length too large to histogram \n");
goto zskip;
}
}
else {
if (count<200) {
printf("\n");
fprintf(Outfp,"\n");
}
}
count=count+1;
}
shift(V, Y, 1);
add64(Y, Y);
sub64(V, Y); // n-(n/2)*2
while ((Y[0]==0)&&(Y[1]==0)) { // 2 divides n
shift(V, V, 1);
n=n+1; // increment number of terms in sequence
// if (count<200)
// printf("k=%d %d, n=%d \n",V[0],V[1],n);
shift(V, Y, 1);
add64(Y, Y);
sub64(V, Y); // n-(n/2)*2
}
if (c==1) {
Y[0]=0;
Y[1]=1;
sub64(V, Y);
if ((Y[0]!=0)||(Y[1]!=0))
goto aloop;
}
if (c==-1) {
Y[0]=0;
Y[1]=1;
sub64(V,Y);
if ((Y[0]==0)&&(Y[1]==0))
goto cskip;
Y[0]=0;
Y[1]=5;
sub64(V,Y);
if ((Y[0]==0)&&(Y[1]==0))
goto cskip;
Y[0]=0;
Y[1]=7;
sub64(V,Y);
if ((Y[0]==0)&&(Y[1]==0))
goto cskip;
Y[0]=0;
Y[1]=17;
sub64(V,Y);
if ((Y[0]==0)&&(Y[1]==0))
goto cskip;
Y[0]=0;
Y[1]=25;
sub64(V,Y);
if ((Y[0]==0)&&(Y[1]==0))
goto cskip;
Y[0]=0;
Y[1]=37;
sub64(V,Y);
if ((Y[0]==0)&&(Y[1]==0))
goto cskip;
Y[0]=0;
Y[1]=55;
sub64(V,Y);
if ((Y[0]==0)&&(Y[1]==0))
goto cskip;
Y[0]=0;
Y[1]=41;
sub64(V,Y);
if ((Y[0]==0)&&(Y[1]==0))
goto cskip;
Y[0]=0;
Y[1]=61;
sub64(V,Y);
if ((Y[0]==0)&&(Y[1]==0))
goto cskip;
Y[0]=0;
Y[1]=91;
sub64(V,Y);
if ((Y[0]==0)&&(Y[1]==0))
goto cskip;
goto aloop;
}
cskip:n=0;
}
//
// check for valid lengths
//
for (i=0; i<178; i++) {
if (he[y[i]]!=0) {
printf("error: n=%d \n",y[i]);
goto zskip;
}
}
//
// HISTOGRAM OF LENGTHS
//
fprintf(Outfp,"\n");
fprintf(Outfp,"HISTOGRAM \n");
printf("\n");
for (h=0; h<100; h++) {
if (h==58)
printf("\n");
fprintf(Outfp," %d %d %d %d %d %d %d %d \n",he[8*h],he[8*h+1],he[8*h+2],
he[8*h+3],he[8*h+4],he[8*h+5],he[8*h+6],he[8*h+7]);
printf(" %d %d %d %d %d %d %d %d \n",he[8*h],he[8*h+1],he[8*h+2],
he[8*h+3],he[8*h+4],he[8*h+5],he[8*h+6],he[8*h+7]);
}
zskip:
fclose(Outfp);
return(0);
}