/******************************************************************************/
/*									      */
/*  FIND LIMBS IN S							      */
/*  09/02/10 (dkc)							      */
/*									      */
/*  This C program finds a limb in S having a given parity vector.  n is the .*/
/*  number of odd elements in the limb and l is the length of the limb.       */
/*  "Offset" is the rotation amount of Halbeisen and Hungerbuhler's parity    */
/*  vector.								      */
/*									      */
/*  This program is for use on the C64. 				      */
/*									      */
/******************************************************************************/
#include <math.h>
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);
void mult3(unsigned int *a, unsigned int n);
unsigned int limb(unsigned int A, unsigned int B, unsigned int shift,
		  unsigned int j, unsigned int delta, unsigned int *O);
far unsigned int error[2];	// flags
far unsigned int G[4];
far unsigned int sv[158];	// parity vector
far unsigned int ab[100*2]={	// a and b values
 0, 1, 0, 2, 1, 2, 1, 3, 2, 3,
 2, 4, 3, 4, 4, 4, 4, 5, 5, 5,
 5, 6, 6, 6, 7, 6, 7, 7, 8, 7,
 8, 8, 9, 8, 9, 9, 10, 9, 11, 9,
 11, 10, 12, 10, 12, 11, 13, 11, 14, 11,
 14, 12, 15, 12, 15, 13, 16, 13, 16, 14,
 17, 14, 18, 14, 18, 15, 19, 15, 19, 16,
 20, 16, 21, 16, 21, 17, 22, 17, 22, 18,
 23, 18, 23, 19, 24, 19, 25, 19, 25, 20,
 26, 20, 26, 21, 27, 21, 28, 21, 28, 22,
 29, 22, 29, 23, 30, 23, 31, 23, 31, 24,
 32, 24, 32, 25, 33, 25, 33, 26, 34, 26,
 35, 26, 35, 27, 36, 27, 36, 28, 37, 28,
 38, 28, 38, 29, 39, 29, 39, 30, 40, 30,
 40, 31, 41, 31, 42, 31, 42, 32, 43, 32,
 43, 33, 44, 33, 45, 33, 45, 34, 46, 34,
 46, 35, 47, 35, 47, 36, 48, 36, 49, 36,
 49, 37, 50, 37, 50, 38, 51, 38, 52, 38,
 52, 39, 53, 39, 53, 40, 54, 40, 54, 41,
 55, 41, 56, 41, 56, 42, 57, 42, 57, 43};
int main () {
//
unsigned int len=58;			// specified limb length
unsigned int shift=5;			// right-shift amount of order
unsigned int delta=4;			// S increment (usually 2)
unsigned int J[2]={0x00006000, 0x00000000}; // order
unsigned int I[2]={0x00003000, 0x00000000}; // order/2
unsigned int K[2]={0x00002000, 0x00000000}; // order/3
unsigned int m=2;			// number of words
//
unsigned int T[2],U[2],length,offset,temp;
unsigned int i,j,k,l,n,a,b;
error[0]=0;				// clear error and flag array
error[1]=0;
rshiftn(J, J, shift, m);		// right-shift order
if ((J[0]==0)||(J[0]==1)) {
   error[1]=0xffffffff;
   goto zskip;
   }
rshiftn(I, I, shift, m);		// right-shift order/2
rshiftn(K, K, shift, m);		// right-shift order/3
setn(T, 2, m);				// load 2
addn(T, I, m);				// load order/2+2
setn(G, 0, m*2);
for (i=0; i<158; i++)			// clear parity vector
   sv[i]=0xffffffff;
for (i=0; i<100; i++) {
   length=3*ab[2*i]+2*ab[2*i+1];	// length of limb
   if (length!=len)			// check for specified length
      continue;
   l=2*ab[2*i]+ab[2*i+1]+1;		// length of cycle
   n=ab[2*i]+ab[2*i+1]; 		// number of odd elements in cycle
   offset=0;				// set index offsets
   if (length==4)			// l=3, n=2, shift=18, delta=4
      offset=0;
   if (length==7)			// l=5, n=3, shift=18, delta=4
      offset=0;
   if (length==9)			// l=6, n=4, shift=18, delta=4
      offset=0;
   if (length==12)			// l=8, n=5, shift=18, delta=4
      offset=0; 			// or 5
   if (length==14)			// l=9, n=6, shift=18, delta=4
      offset=0; 			// or 3,6
   if (length==17)			// l=11, n=7, shift=18, delta=4
      offset=0; 			// or 5,8
   if (length==20)			// l=13, n=8, shift=18, delta=4
      offset=5; 			// or 10, not 0
   if (length==22)			// l=14, n=9, shift=18, delta=4
      offset=11;			// or 5,8, not 0
   if (length==25)			// l=16, n=10, shift=18, delta=4
      offset=0; 			// or 5,8,13
   if (length==27)			// l=17, n=11, shift=18, delta=4
      offset=11;			// or 5,8, not 0,14
   if (length==30)			// l=19, n=12, shift=16, delta=4
      offset=0; 			// or 5,8,13,16
   if (length==33)			// l=21, n=13, shift=14, delta=4
      offset=5; 			// or 10,18 not 0,13
//
// Note: Sequence vector giving smallest e value unknown beyond this point.
//	 Sequence vector giving smallest e value assumed for 38 and 43 since
//	 all possible rotates give solutions.
//
   if (length==35)			// l=22, n=14, shift=13, delta=4
      offset=5; 			// or 8,11,16,19, not 0
   if (length==38)			// l=24, n=15, shift=13, delta=4
      offset=0; 			// or 5,8,13,16,21
   if (length==40)			// l=25, n=16, shift=12, delta=4
      offset=5; 			// or 16,19 not 0,8,11,22
   if (length==43)			// l=27, n=17, shift=12, delta=4
      offset=0; 			// or 5,8,13,16,21,24
   if (length==45)			// l=28, n=18, shift=12, delta=4
      offset=5; 			// 5 or 19, not 0,8,11,14,22,25
   if (length==48)			// l=30, n=19, shift=11, delta=4
      offset=5; 			// or 8,13,16,27 not 0,19
   if (length==51)			// l=32, n=20, shift=8, delta=4
      offset=5; 			// or 13,21,29, not 0,8,16,24
   if (length==53)			// l=33, n=21, shift=8, delta=4
      offset=5; 			// or 5,8,11,16,19,22,27,30, not 0
   if (length==56)			// l=35, n=22, shift=8, delta=4
      offset=5; 			// or 8,13,16,21,24,29,32, not 0
   if (length==58)			// l=36, n=23, shift=5, delta=4
      offset=11;			// or 5, ??16,19,22,27,30,33, not 0,8,11
//
// Note: Offsets haven't been determined for lengths greater than 58.
//
//
//  compute parity vector
//
   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 addressing
      if (temp>l)
	 temp=temp-l;
      sv[temp-1]=k;
      }
//
// find limb in S having given sequence vector
//
   U[0]=I[0];				// starting value
   U[1]=I[1];
aloop:
   j=limb(U[0], U[1], 13-shift, length-1, delta, G);  // find limb
   if (j==0)				// check if finished
      goto zskip;
   U[0]=G[2];				// load ending value
   U[1]=G[3];
   copyn(K, T, m);			// load order/3
   negn(T, m);				// load -order/3
   addn(G, T, m);			// odd element minus order/3
   if ((T[0]&0x80000000)!=0) {		// if negative, continue
      error[1]=error[1]+1;
      goto aloop;
      }
   error[0]=1;				// save limb count
   }
zskip:
return(0);
}