/*****************************************************************************/
/*									     */
/*  EUCLIDEAN G.C.D.							     */
/*  01/12/07 (dkc)							     */
/*									     */
/*****************************************************************************/
unsigned int euclid(unsigned int d, unsigned int e) {
unsigned int a,b,temp;
a=d;
b=e;
if (b>a) {
   temp=a;
   a=b;
   b=temp;
   }
loop:
   temp=a-(a/b)*b;
   a=b;
   b=temp;
   if (b!=0) goto loop;
return(a);
}
/******************************************************************************/
/*									      */
/*  FIND APPROXIMATIONS TO THE PARITY VECTOR REQUIRED FOR A CYCLE	      */
/*  09/23/10 (dkc)							      */
/*									      */
/*  This C program finds parity vectors for live limbs in S of a least-       */
/*  residue tree.							      */
/*									      */
/******************************************************************************/
#include <math.h>
unsigned int euclid(unsigned int a, unsigned int b);
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 maxn=65536;	  // set maximum n value
//
far unsigned int error[4];	  // error array
far unsigned int P[4096],MAX[4096],MIN[4096],T[4096];  // P, max, min arrays
far unsigned int sv[131072];	  // parity vector
far unsigned int output[200];	  // limb lengths (with exceptions)
far unsigned int output1[200];	  // limb lengths
far unsigned int ln[200];	  // l and n values
far unsigned int X[4096];	  // X array
int main () {
unsigned int c,d,length,offset,temp,flag1,flag2,flag3,flag4,count;
unsigned int f,g,h,i,j,k,l,n,p,a,b,oldg,oldh,iters,factor,oldoff,deloff;
int del,delo,delg;
if (maxn>65536) {		  // set size of X array
   error[0]=0xffffffff;
   goto zskip;
   }
if (maxn<=1024)
   p=64;
else {
   if (maxn<=2048)
      p=128;
   else {
      if (maxn<=4096)
	 p=256;
      else {
	 if (maxn<=8192)
	    p=512;
	 else {
	    if (maxn<=16384)
	       p=1024;
	    else {
	       if (maxn<=32768)
		  p=2048;
	       else
		  p=4096;
	       }
	    }
	 }
      }
   }
setn(sv, 0, 131072);		// clear s.v. array
setn(P, 0, 4096);		// clear P array
setn(X, 0, 4096);		// clear X array
error[0]=0;			// clear error array
error[1]=0;
error[2]=0;
error[3]=0;
iters=0;			// clear exception count
for (i=0; i<200; i++) { 	// clear output arrays
   output[i]=0;
   output1[i]=0;
   ln[i]=0;
   }
count=0;			// clear limb count
c=0;				// clear c
d=0;				// clear d
setn(X, 0, p);
X[0]=0x10000000;		// x=1/2
for (i=0; i<maxn; i++) {	// maximum n value
   if (((X[0]+0xe0000000)&0x80000000)==0) {  // x>1
      mul3shft(X, 2, p);	// x=(3/4)x
      c=c+1;			// increment c
      }
   else {			// x<1
      mul3shft(X, 1, p);	// x=(3/2)x
      d=d+1;			// increment d
      }
   if (X[p-1]!=0) {		// check for overflow of X array
      error[2]=1;
      goto zskip;
      }
   length=3*c+2*d;		// length of limb
   l=2*c+d+1;			// length of cycle
   if (l>=131072) {		// check for overflow of s.v. array
      error[2]=2;
      goto zskip;
      }
   n=c+d;			// number of odd elements in cycle
   error[3]=n;
   factor=euclid(l, n); 	// find G.C.D. of l and n
   f=l/factor;			// length of non-repeating portion of s.v.
//
//  compute parity vector
//
   flag1=0;			// clear offset count
   flag2=0;			// clear exception flag
   flag3=0;			// clear first-time flag
   flag4=0;			// clear exception flag
   oldg=0;			// clear old index
   oldh=0;			// clear old index
   oldoff=0;			// clear old offset
   for (offset=0; offset<l; offset++) {
      for (j=1; j<=l; j++) {
	 a=j*n;
	 if (a==((a/l)*l))	// a=(j*n)/l
	    a=a/l;
	 else
	    a=(a/l)+1;
	 b=(j-1)*n;
	 if (b==((b/l)*l))	// b=((j-1)*n)/l
	    b=b/l;
	 else
	    b=(b/l)+1;
	 k=a-b; 		// parity value
	 temp=j+offset; 	// store value using circular address		ing
	 if (temp>l)
	    temp=temp-l;
	 sv[temp]=k;
	 }
//
// check parity vector
//
      if (sv[1]!=0)
	 continue;
      if (sv[2]!=1)
	 continue;
      if (sv[3]!=0)
	 continue;
      flag1=flag1+1;		// increment offset count
      setn(MAX, 0, p);		// set maximum
      setn(MIN, 0, p);		// set minimum
      MIN[0]=0x7fffffff;
      g=0;
      h=0;
      setn(P, 0, p);
      P[0]=0x20000000;		// set P
      for (j=1; j<=l; j++) {
	 if (sv[j]==1) {
	    copyn(MAX, T, p);	// find maximum
	    negn(T, p);
	    addn(P, T, p);
	    if ((T[0]&0x80000000)==0) {
	       copyn(P, MAX, p);
	       g=j;
	       }
	    copyn(P, T, p);	// find minimum
	    negn(T, p);
	    addn(MIN, T, p);
	    if ((T[0]&0x80000000)==0) {
	       copyn(P, MIN, p);
	       h=j;
	       }
	    mul3shft(P, 1, p);	// P=1.5*P
	    }
	 else
	    rshiftn(P, P, 1, p); // P=0.5*P
	 }
      if (P[p-1]!=0) {		// check for overflow of P array
	 error[2]=3;
	 goto zskip;
	 }
      copyn(MIN, T, p); 	// compare twice the minimum to maximum
      addn(T, T, p);
      negn(T, p);
      addn(MAX, T, p);
      if ((T[0]&0x80000000)==0)
	 goto askip;
      del=(int)(g-h);		// difference in indices
      delo=(int)(oldg-oldh);	// difference in old indices
      if (del<0)
	 del=del+l;
      del=del-(del/f)*f;
      if (delo<0)
	 delo=delo+l;
      delo=delo-(delo/f)*f;
      if ((del!=delo)&&(flag3!=0))
	 flag2=length;		// set exception flag
      if ((flag3!=0)&&(flag1>factor)) {
	 deloff=offset-oldoff;	// difference in offsets
	 delg=g-oldg;		// difference in indices
	 if (delg<0)
	    delg=delg+l;
	 delg=delg-(delg/f)*f;
	 if ((int)deloff!=delg)
	    flag4=length;	// set exception flag
	 }
      oldg=g;			// set old index
      oldh=h;			// set old index
      oldoff=offset;		// set old offset
      flag3=1;			// set first-time flag
      }
   if (flag1!=0) {
      error[0]=error[0]+1;	// increment limb count
      output1[count]=(length<<16)|factor;    // store length and factor
      ln[count]=(l<<16)|n;	// store l and n values
      count=count+1;
      if (flag2==length) {	// save exception
	 output[iters]=length;
	 iters=iters+1;
	 error[1]=iters;
	 if (iters>200)
	    goto zskip;
	 }
      if (flag4==length) {	// save exception
	 output[iters]=length|0x80000000;
	 iters=iters+1;
	 error[1]=iters;
	 if (iters>200)
	    goto zskip;
	 }
      }
askip:
   flag1=0;
   }
zskip:
return(0);
}