/******************************************************************************/
/*									      */
/*  FIND THE MINIMUM ELEMENT IN A CYCLE 				      */
/*  09/02/10 (dkc)							      */
/*									      */
/*  This C program finds the minimum element in a cycle.  c is set to 1 or -1.*/
/*  Potential cycles are in limbs in S of a least-residue tree.  n is the     */
/*  number of odd elements in the cycle and l is the length of the cycle.     */
/*									      */
/*  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);
far unsigned int error[2];	// flags
far unsigned int G[2];
far unsigned int H[2];		// greater than order/2, less than order
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=48;			// specified limb length
unsigned int shift=11;			// right-shift amount of order
unsigned int delta=2;			// 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 C[2]={0x00000000, 0x00000001}; // c
//unsigned int C[2]={0xffffffff, 0xffffffff}; // c
unsigned int maxiters=1;		// maximum number of limbs in S
unsigned int m=2;			// number of words
//
unsigned int F[2],T[2],U[2],length,offset,temp;
unsigned int i,j,k,l,n,a,b,oldg,index,iters;
rshiftn(J, J, shift, m);		// right-shift order
rshiftn(I, I, shift, m);		// right-shift order/2
rshiftn(K, K, shift, m);		// right-shift order/3
copyn(J, F, m); 			// load order
negn(F, m);				// load -order
setn(T, 2, m);				// load 2
addn(T, I, m);				// load order/2+2
setn(U, delta, m);			// load delta
error[0]=0;				// clear error and flag array
error[1]=0;
for (i=0; i<158; i++)			// clear parity vector
   sv[i]=0xffffffff;
iters=0;				// clear limb count
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=2
      offset=0;
   if (length==7)			// l=5, n=3, shift=18, delta=2
      offset=0;
   if (length==9)			// l=6, n=4, shift=18, delta=2
      offset=0;
   if (length==12)			// l=8, n=5, shift=18, delta=2
      offset=0; 			// or 5
   if (length==14)			// l=9, n=6, shift=18, delta=2
      offset=0; 			// or 3,6
   if (length==17)			// l=11, n=7, shift=18, delta=2
      offset=0; 			// or 5,8
   if (length==20)			// l=13, n=8, shift=18, delta=2
      offset=5; 			// or 10, not 0
   if (length==22)			// l=14, n=9, shift=18, delta=2
      offset=11;			// or 5,8, not 0
   if (length==25)			// l=16, n=10, shift=18, delta=2
      offset=0; 			// or 5,8,13
   if (length==27)			// l=17, n=11, shift=18, delta=2
      offset=11;			// or 5,8, not 0,14
   if (length==30)			// l=19, n=12, shift=16, delta=2
      offset=0; 			// or 5,8,13,16
   if (length==33)			// l=21, n=13, shift=14, delta=2
      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=2
      offset=5; 			// or 8,11,16,19, not 0
   if (length==38)			// l=24, n=15, shift=13, delta=2
      offset=0; 			// or 5,8,13,16,21
   if (length==40)			// l=25, n=16, shift=12, delta=2
      offset=5; 			// or 16,19 not 0,8,11,22
   if (length==43)			// l=27, n=17, shift=12, delta=2
      offset=0; 			// or 5,8,13,16,21,24
   if (length==45)			// l=28, n=18, shift=12, delta=2
      offset=5; 			// 5 or 19, not 0,8,11,14,22,25
   if (length==48)			// l=30, n=19, shift=11, delta=2
      offset=5; 			// or 8,13,16,27 not 0,19
   if (length==51)			// l=32, n=20, shift=8, delta=2
      offset=5; 			// or 13,21,29, not 0,8,16,24
   if (length==53)			// l=33, n=21, shift=8, delta=2
      offset=5; 			// or 5,8,11,16,19,22,27,30, not 0
   if (length==56)			// l=35, n=22, shift=8, delta=2
      offset=5; 			// or 8,13,16,21,24,29,32, not 0
   if (length==58)			// l=36, n=23, shift=7, delta=2
      offset=5; 			// or ??? 8,11,16,19,22,27,30,33, not ?
//
// 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
//
   copyn(I, H, m);			// save starting element of limb
aloop:
      copyn(H, G, m);			// load starting element of limb
      index=0;				// clear index
      oldg=0;
      for (j=0; j<(length-1); j++) {
	 if ((G[m-1]&1)==0) {		// check for even element
	    if (oldg==0) {		// check sequence
	       if (sv[index]!=1)
		  goto askip;
	       }
	    else {
	       if (sv[index]!=0)
		  goto askip;
	       }
	    oldg=1;
	    rshiftn(G, G, 1, m);	// halve element
	    index=index+1;		// increment index
	    }
	 else { 			// odd element
	    mult3(G, m);		// compute next even element
	    addn(C, G, m);		//    "
	    copyn(F, T, m);		// load -order
	    addn(G, T, m);		// even element minus order
	    if ((T[0]&0x80000000)==0)	// if non-negative, continue
	       goto askip;
	    oldg=0;
	    }
	 }
      if ((G[m-1]&1)==1) {		// check for odd element
	 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 askip;
	    }
	 iters=iters+1; 		// increment limb count
	 error[0]=iters;		// save limb count
	 if (iters==maxiters)
	    goto zskip;
	 }
askip:
   addn(U, H, m);			// increment starting element
   copyn(H, T, m);			// load starting element
   negn(T, m);				// -starting element
   addn(J, T, m);			// order minus starting element
   if ((T[0]&0x80000000)==0)		// if negative, exit loop
      goto aloop;
   }
zskip:
return(0);
}