/******************************************************************************/
/*									      */
/*  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);
      add256(Y, X);
      add256(X, Y);		    //	  n=n*3;
      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))
      add256(O, O);		    //	  order=order*2;
      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)
	 goto askip;		    //	     goto askip;
      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)
	    goto askip;
	 set256(U, 7);
	 sub256(Y, U);
	 if ((U[0]|U[1]|U[2]|U[3]|U[4]|U[5]|U[6]|U[7])==0)
	    goto askip;
	 set256(U, 17);
	 sub256(Y, U);
	 if ((U[0]|U[1]|U[2]|U[3]|U[4]|U[5]|U[6]|U[7])==0)
	    goto askip;
	 set256(U, 25);
	 sub256(Y, U);
	 if ((U[0]|U[1]|U[2]|U[3]|U[4]|U[5]|U[6]|U[7])==0)
	    goto askip;
	 set256(U, 37);
	 sub256(Y, U);
	 if ((U[0]|U[1]|U[2]|U[3]|U[4]|U[5]|U[6]|U[7])==0)
	    goto askip;
	 set256(U, 55);
	 sub256(Y, U);
	 if ((U[0]|U[1]|U[2]|U[3]|U[4]|U[5]|U[6]|U[7])==0)
	    goto askip;
	 set256(U, 41);
	 sub256(Y, U);
	 if ((U[0]|U[1]|U[2]|U[3]|U[4]|U[5]|U[6]|U[7])==0)
	    goto askip;
	 set256(U, 61);
	 sub256(Y, U);
	 if ((U[0]|U[1]|U[2]|U[3]|U[4]|U[5]|U[6]|U[7])==0)
	    goto askip;
	 set256(U, 91);
	 sub256(Y, U);
	 if ((U[0]|U[1]|U[2]|U[3]|U[4]|U[5]|U[6]|U[7])==0)
	    goto askip;
	 }
      copy256(Y, U);
      add256(Y, Y);		    //	   2*n
      add256(U, Y);		    //	   3*n
      add256(Z, Y);		    //	   n=3*n+c
      count=count+1;
      if ((Y[7]&0x7)==0)
	 goto askip;		    //	      goto askip;
      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);
   add256(Y, Y);		     //  n=n*2;
   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);
      add256(Q, R);		     //     ((n-c)/3)*2
      add256(Q, R);		     //     ((n-c)/3)*3
      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);
	 add256(Y, Y);		     //        n=n*2;
	 count=count+2;
	 }			     //        }
      else {
	 add256(Y, Y);		     //        n=n*2;
	 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);
	 add256(Y, W);
	 add256(W, Y);
	 add256(Z, Y);
	 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++)
      add256(X, X);
   set256(Y, 1);
   for (i=0; i<len; i++) {
      copy256(Y, W);
      add256(Y, W);
      add256(W, Y);
      }
   set256(U, 0);
   for (h=0; h<len; h++) {
      set256(W, 1);
      for (i=0; i<h; i++) {
	 copy256(W, V);
	 add256(W,V);
	 add256(V,W);
	 }
      sum=0;
      for (i=0; i<(len-h); i++)
	 sum=sum+sv[i];
      lshift256(W, W, sum);
      add256(W, U);
      }
   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");
askip:
   count=0;
   }
zskip:
fclose(Outfp);
return(0);
}