/******************************************************************************/
/*									      */
/*  FIND PARITY VECTOR							      */
/*  11/26/10 (dkc)							      */
/*									      */
/*  This C program finds the parity vector required when twice the minimum    */
/*  element in a 3n+c cycle is larger than the maximum odd element (and the   */
/*  cycle has large elements).						      */
/*									      */
/******************************************************************************/
#include <math.h>
#include <stdio.h>
void shiftn(unsigned int *a, unsigned int *b, unsigned int shift, unsigned int n);
void addn(unsigned int *a, unsigned int *b, unsigned int n);
void copyn(unsigned int *a, unsigned int *b, unsigned int n);
void setn(unsigned int *a, unsigned int b, unsigned int n);
void subn(unsigned int *a, unsigned int *b, unsigned int n);
//
unsigned int error[6];	      // error array
unsigned int P[8192],Q[8192],O[8192],T[8192];
unsigned int sv[65536];      // parity vector
unsigned int ln[13000];       // l and n values

int main () {
unsigned int i,count,p,n,index;
FILE *Outfp;
Outfp = fopen("out14uf.dat","w");
p=2048;
for (i=0; i<65536; i++)
   sv[i]=0;
setn(P, 0, 8192);		// clear P array
error[0]=0;			// clear error array
error[1]=0;
error[2]=0;
error[3]=0;
error[4]=0;
error[5]=0;
for (i=0; i<13000; i++) 	 // clear output array
   ln[i]=0;
index=0;
n=1;
sv[0]=1;
count=1;		   // set index
setn(O, 0, p);		   // load 1
O[0]=0x20000000;
setn(P, 0, p);
P[0]=0x20000000;	   // set P
for (i=0; i<30000; i++) {
   sv[count]=1;
   count=count+1;
   if (count>65535) {
      error[2]=2;
      goto zskip;
      }
   n=n+1;
   copyn(P, Q, p);
   addn(Q, Q, p);
   addn(P, Q, p);
   shiftn(Q, P, 1, p);
   if (P[p-1]!=0) {	   // check for overflow of P array
      error[2]=3;
      goto zskip;
      }
   copyn(P, T, p);
   subn(O, T, p);
   if ((T[0]&0x80000000)!=0) {
      sv[count]=0;
      count=count+1;
      if (count>65535) {
	 error[2]=2;
	 goto zskip;
	 }
      copyn(P, Q, p);
      shiftn(Q, P, 1, p); // P=0.5*P
      ln[index]=(count<<16)|n;
      index=index+1;
      if (index>12999) {
	 error[2]=4;
	 goto zskip;
	 }
      }
   }
zskip:
error[5]=count;
for (i=0; i<6; i++)
   printf("%d \n",error[i]);
if ((error[2]==0)||(error[2]==4)) {
   for (i=0; i<index; i++)
      fprintf(Outfp,"%#010x, \n",ln[i]);
   }
fclose(Outfp);
return(0);
}