/******************************************************************************/
/*									      */
/*  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>

void mul3shft(unsigned int *c, unsigned int b, unsigned int n);
void rshiftn(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 negn(unsigned int *a, unsigned int n);
//
far unsigned int error[6];	  // error array
far unsigned int P[8192],O[8192],T[8192];
far unsigned int sv[65536];	 // parity vector
far unsigned int ln[1000];	 // l and n values

int main () {
unsigned int i,count,p,n,index;
p=256;
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<1000; 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<2000; i++) {
   sv[count]=1;
   count=count+1;
   n=n+1;
   mul3shft(P, 1, p);	   // P=1.5*P
   if (P[p-1]!=0) {	   // check for overflow of P array
      error[2]=3;
      goto zskip;
      }
   copyn(P, T, p);
   negn(T, p);
   addn(O, T, p);
   if ((T[0]&0x80000000)!=0) {
      sv[count]=0;
      count=count+1;
      rshiftn(P, P, 1, p); // P=0.5*P
      ln[index]=(count<<16)|n;
      index=index+1;
      if (index>999) {
	 error[2]=3;
	 goto zskip;
	 }
      }
   }
zskip:
return(0);
}