﻿ /*******************************
```/******************************************************************************/
/*									      */
/*  COMPUTE 3N+1 OR 3N-1 SEQUENCE					      */
/*  02/13/10 (dkc)							      */
/*									      */
/*  This C program finds long limbs in S.  Sequence vectors are generated.    */
/*									      */
/******************************************************************************/
#include <stdio.h>
#include <math.h>
void div256_32(unsigned int a0, unsigned int a1, unsigned int a2,
unsigned int a3, unsigned int a4, unsigned int a5,
unsigned int a6, unsigned int a7, unsigned int *quotient,
unsigned int d7);
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 () {
//
int c=-1;				  // c
//
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};	  // I
unsigned int g,h,i,count,scount,temp;  // indices, etc.
unsigned int O[8],P[8],Q[8],R[8],U[8],V[8],X[8],W[8],Y[8],Z[8];
unsigned int first,sum,save,len,sv[2048];
FILE *Outfp;
Outfp = fopen("out11.dat","w");
if (c==1)
set256(Z, 1);
else {
Z[0]=0xffffffff;
Z[1]=0xffffffff;
Z[2]=0xffffffff;
Z[3]=0xffffffff;
Z[4]=0xffffffff;
Z[5]=0xffffffff;
Z[6]=0xffffffff;
Z[7]=0xffffffff;
}
//
// begin loop
//
for (h=2; h<256; h++) {
set256(Y, 1);
for (i=0; i<h; i++) {	    // for (i=0; i<h; i++)
copy256(Y, X);
if (Y[0]>0x40000000) {
printf("h=%d, order too big \n",h);
goto zskip;
}
}
//
copy256(Z, V);
sub256(Y, V);
shift256(V, Y, 1);		    // n=(n-c)/2;
if ((Y[7]&1)==0)
continue;
//
count=2*(h+1);
g=h; 			    // power of 3 count
lshift256(Y, U, 2);
set256(O, 3);
temp=0;
copy256(U, W);
sub256(O, W);
while ((W[0]&0x80000000)!=0) {   // while (order<(4*n))
temp=temp+1;
copy256(U, W);
sub256(O, W);
}
//
div256_32(O[0],O[1],O[2],O[3],O[4],O[5],O[6],O[7],W, 3);   // order/3
copy256(W, U);
sub256(Y, U);
while ((U[0]&0x80000000)!=0) {   // while (n<(order/3)) {
g=g+1;
set256(U, 1);
sub256(Y, U);
if ((U[0]|U[1]|U[2]|U[3]|U[4]|U[5]|U[6]|U[7])==0) // if (n==1)
if (c==-1) {
set256(U, 5);
sub256(Y, U);
if ((U[0]|U[1]|U[2]|U[3]|U[4]|U[5]|U[6]|U[7])==0)
set256(U, 7);
sub256(Y, U);
if ((U[0]|U[1]|U[2]|U[3]|U[4]|U[5]|U[6]|U[7])==0)
set256(U, 17);
sub256(Y, U);
if ((U[0]|U[1]|U[2]|U[3]|U[4]|U[5]|U[6]|U[7])==0)
set256(U, 25);
sub256(Y, U);
if ((U[0]|U[1]|U[2]|U[3]|U[4]|U[5]|U[6]|U[7])==0)
set256(U, 37);
sub256(Y, U);
if ((U[0]|U[1]|U[2]|U[3]|U[4]|U[5]|U[6]|U[7])==0)
set256(U, 55);
sub256(Y, U);
if ((U[0]|U[1]|U[2]|U[3]|U[4]|U[5]|U[6]|U[7])==0)
set256(U, 41);
sub256(Y, U);
if ((U[0]|U[1]|U[2]|U[3]|U[4]|U[5]|U[6]|U[7])==0)
set256(U, 61);
sub256(Y, U);
if ((U[0]|U[1]|U[2]|U[3]|U[4]|U[5]|U[6]|U[7])==0)
set256(U, 91);
sub256(Y, U);
if ((U[0]|U[1]|U[2]|U[3]|U[4]|U[5]|U[6]|U[7])==0)
}
copy256(Y, U);
count=count+1;
if ((Y[7]&0x7)==0)
while ((Y[7]&1)==0) {
shift256(Y, Y, 1);	    //	     n=n/2;
count=count+1;
}			     //       }
copy256(W, U);
sub256(Y, U);
} 			     //    }
//
g=g+1;
copy256(Y, P);
set256(Y, 1);
lshift256(Y, Y, h);
copy256(Z, W);
sub256(Y, W);
copy256(W, Y);
count=count+1;
shift256(O, W, 1);		     //  order/2
copy256(Y, V);
sub256(W, V);
while ((V[0]&0x80000000)==0) {    //  while (n<(order/2)) {
copy256(Z, U);
sub256(Y, U);		      //     n-c
div256_32(U[0],U[1],U[2],U[3],U[4],U[5],U[6],U[7],Q,3);	//     (n-c)/3
copy256(Q, R);
sub256(U, R);
if ((R[0]|R[1]|R[2]|R[3]|R[4]|R[5]|R[6]|R[7])==0) {
g=g+1;
copy256(Q, Y);
count=count+2;
}			     //        }
else {
count=count+1;
}
copy256(Y, V);
sub256(W, V);
}
printf("order=(%#10x %#10x %#10x %#10x \n",O[0],O[1],O[2],O[3]);
printf("       %#10x %#10x %#10x %#10x), length=%d, h=%d, g=%d \n",O[4],O[5],O[6],O[7],count,h,g);
printf("limb= (%#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("       %#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]);
fprintf(Outfp,"order=(%#10x %#10x %#10x %#10x \n",O[0],O[1],O[2],O[3]);
fprintf(Outfp,"       %#10x %#10x %#10x %#10x), length=%d, h=%d, g=%d \n",O[4],O[5],O[6],O[7],count,h,g);
fprintf(Outfp,"limb= (%#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,"       %#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]);
for (i=0; i<178; i++) {
if (count==y[i]) {
printf("invalid length=%d \n",count);
fprintf(Outfp,"invalid length=%d \n",count);
goto zskip;
}
}
scount=count+1;
//
// COMPUTE SEQUENCE VECTOR
//
i=0;
if ((Y[7]&1)==0) {
count=1;
while ((Y[7]&1)==0) {
count=count+1;
shift256(Y, Y, 1);
}
sv[0]=count;
i=1;
}
first=1;
copy256(P, W);
sub256(Y, W);
while ((W[0]|W[1]|W[2]|W[3]|W[4]|W[5]|W[6]|W[7])!=0) {
if ((Y[7]&1)==0) {
shift256(Y, Y, 1);
count=count+1;
}
else {
copy256(Y, W);
if (first==0) {
sv[i]=count;
i=i+1;
if (i>2047) {
printf("error: sequence vector array not big enough \n");
goto zskip;
}
}
else
first=0;
count=0;
}
copy256(P, W);
sub256(Y, W);
}
sv[i]=count;
i=i+1;
//
// COMPUTE X, Y, AND Z
//
len=i;
if (len!=g) {
printf("error: len=%d, g=%d \n",len,g);
goto zskip;
}
sum=0;
for (i=0; i<len; i++)
sum=sum+sv[i];
save=sum;
if ((sum+len)!=scount) {
printf("error: length=%d, sum=%d \n",scount-1,sum+len);
goto zskip;
}
set256(X, 1);
for (i=0; i<sum; i++)
set256(Y, 1);
for (i=0; i<len; i++) {
copy256(Y, W);
}
set256(U, 0);
for (h=0; h<len; h++) {
set256(W, 1);
for (i=0; i<h; i++) {
copy256(W, V);
}
sum=0;
for (i=0; i<(len-h); i++)
sum=sum+sv[i];
lshift256(W, W, sum);
}
shift256(U, U, 1);
printf("     X=%#10x %#10x %#10x %#10x \n",X[0],X[1],X[2],X[3]);
printf("       %#10x %#10x %#10x %#10x, exp=%d \n",X[4],X[5],X[6],X[7],save);
printf("     Y=%#10x %#10x %#10x %#10x \n",Y[0],Y[1],Y[2],Y[3]);
printf("       %#10x %#10x %#10x %#10x, exp=%d \n",Y[4],Y[5],Y[6],Y[7],len);
printf("     Z=%#10x %#10x %#10x %#10x \n",U[0],U[1],U[2],U[3]);
printf("       %#10x %#10x %#10x %#10x \n",U[4],U[5],U[6],U[7]);
fprintf(Outfp,"     X=%#10x %#10x %#10x %#10x \n",X[0],X[1],X[2],X[3]);
fprintf(Outfp,"       %#10x %#10x %#10x %#10x,exp=%d \n",X[4],X[5],X[6],X[7],save);
fprintf(Outfp,"     Y=%#10x %#10x %#10x %#10x \n",Y[0],Y[1],Y[2],Y[3]);
fprintf(Outfp,"       %#10x %#10x %#10x %#10x, exp=%d \n",Y[4],Y[5],Y[6],Y[7],len);
fprintf(Outfp,"     Z=%#10x %#10x %#10x %#10x \n",U[0],U[1],U[2],U[3]);
fprintf(Outfp,"       %#10x %#10x %#10x %#10x \n",U[4],U[5],U[6],U[7]);
sub256(X, Y);
printf(" numer=%#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]);
fprintf(Outfp," numer=%#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]);
temp=temp-1;
shift256(U, U, temp);    // Z/(order/6)
div256_32(U[0],U[1],U[2],U[3],U[4],U[5],U[6],U[7],W,3);  // Z/(order/2)
if ((Y[0]&0x80000000)!=0) {
set256(X, 0);
sub256(X, Y);
}
sub256(Y, W);
printf(" delta=%#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]);
fprintf(Outfp," delta=%#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]);
printf("\n");
fprintf(Outfp,"\n");