/******************************************************************************/
/*									      */
/*  FIND THE MINIMUM ELEMENT IN A CYCLE 				      */
/*  07/25/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[6];	// flags
far unsigned int MAX[2];	// maximum odd element
far unsigned int MIN[2];	// minimum odd element
far unsigned int H[2];		// greater than order/2, less than order
far unsigned int sv[158];	// parity vector
far double pv[158];
far unsigned int R[257*2];	// limbs in S
far unsigned int O[100*2];	// odd elements of a limb in S
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=53;			// specified limb length
unsigned int sflag=1;			// limb in S flag
unsigned int shift=1;			// right-shift amount of order
unsigned int delta=20000;		// 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=1000000;		// maximum number of limbs in S
unsigned int m=2;			// number of words
//
unsigned int G[2],T[2],length,offset,temp,saveh;
unsigned int h,i,j,k,l,n,a,b,oldg,index,iters,count;
double max,min;
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
setn(T, 2, m);				// load 2
addn(T, I, m);				// load order/2+2
setn(MAX, 0, m);			// clear maximum odd element
setn(MIN, 0, m);			// clear minimum odd element
error[0]=0;				// clear error and flag array
error[1]=0;
error[2]=0;
error[3]=0;
error[4]=0;
error[5]=0;
setn(O, 0, 100*m);			// clear output array
for (i=0; i<158; i++) { 		// clear parity vector
   sv[i]=0xffffffff;
   pv[i]=0.0;
   }
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;
   if (length==14)			// l=9, n=6, shift=18, delta=2
      offset=0;
   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
   if (length==35)			// l=22, n=14, shift=13, delta=2
      offset=8; 			// or 5,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=19;// for 5,16, 3*MIN!>MAX // or 5,16, 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; // for 5,19 3*MIN!>MAX	// or 19, not 0,8,11,14,22,25
   if (length==48)			// l=30, n=19, shift=11, delta=2
      offset=8; 			// or 5,13,16,27 not 0,19
   if (length==51)			// l=32, n=20, shift=1, delta=20000
      offset=5; 			// or 13,21,29, not 0,8,16,24
   if (length==53)			// l=33, n=21, shift=1, delta=20000
      offset=5; 			// or 8,11,16,19,22,27, not 0,30 ???
   if (length==56)			// l=35, n=22, shift=1, delta=20000
      offset=5; 			// or 8,13,16,21,24,29,32, not 0 ???
//
// Note: Offsets haven't been determined for lengths greater than 56.
//
//
//  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;
      }
//
// check parity vector
//
   pv[1]=1.0;
   for (j=1; j<l; j++) {
      if (sv[j]==1)
	 pv[j+1]=1.5*pv[j];
      else
	 pv[j+1]=0.5*pv[j];
      }
   h=1;
   max=0.0;
   min=1000000;
   for (j=1; j<l; j++) {		// find maximum and minimum
      if (sv[j]==1) {
	 if (pv[j]>max) {
	    max=pv[j];
	    h=j;
	    }
	 if (pv[j]<min)
	    min=pv[j];
	 }
      }
   if (pv[l-1]<max)
      error[0]=h;
   if (pv[l-1]>2.0)
      error[4]=1;
   if (3.0*min<max)
      error[4]=2;
   if (2.0*min<max)
      error[5]=1;			// okay when last two s.v. values are 1, 1
//
// 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;
      copyn(G, R, m);			// store limb element
      count=1;
      error[0]=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
	    R[2*count]=G[0];		// save element
	    R[2*count+1]=G[1];
	    count=count+1;
	    }
	 else { 			// odd element
	    mult3(G, m);		// compute next even element
	    addn(C, G, m);		//    "
	    copyn(J, T, m);		// load order
	    negn(T, m); 		// load -order
	    addn(G, T, m);		// even element minus order
	    if (sflag!=0) {
	       if ((T[0]&0x80000000)==0)   // if non-negative, continue
		  goto askip;
	       }
	    else {
	       if ((T[0]&0x80000000)==0)   // if non-negative, set flag
		  error[0]=1;
	       }
	    oldg=0;
	    R[2*count]=G[0];		// save element of limb
	    R[2*count+1]=G[1];
	    count=count+1;
	    }
	 }
      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
	    goto askip;
	 h=0;
	 setn(MAX, 0, m);
	 setn(MIN, 1, m);
	 negn(MIN, m);
	 MIN[0]=MIN[0]&0x7fffffff;
	 for (j=0; j<count; j++) {	// save odd elements of limb
	    if ((R[m*j+(m-1)]&1)==1) {
	       copyn(R+m*j, O+m*h, m);
	       copyn(R+m*j, T, m);	// find maximum
	       negn(T, m);
	       addn(MAX, T, m);
	       if ((T[0]&0x80000000)!=0) {
		  copyn(R+m*j, MAX, m);
		  saveh=h;
		  }
	       copyn(R+m*j, T, m);	// find minimum
	       negn(T, m);
	       addn(MIN, T, m);
	       if ((T[0]&0x80000000)==0)
		  copyn(R+m*j, MIN, m);
	       h=h+1;
	       }
	    }
	 if ((saveh+1)!=h)		// check if last is largest
	    error[1]=1;
	 copyn(MIN, T, m);		// check if 2*MIN>MAX
	 addn(MIN, T, m);
	 copyn(MAX, G, m);
	 negn(G, m);
	 addn(T, G, m);
	 if ((G[0]&0x80000000)!=0)
	    error[3]=1;
	 copyn(MIN, T, m);		// check if 3*MIN>MAX
	 addn(MIN, T, m);
	 addn(MIN, T, m);
	 copyn(MAX, G, m);
	 negn(G, m);
	 addn(T, G, m);
	 if ((G[0]&0x80000000)!=0)
	    error[3]=0xffffffff;
	 iters=iters+1; 		// increment limb count
	 error[2]=iters;		// save limb count
	 if (iters==maxiters)
	    goto zskip;
	 }
askip:
   setn(T, delta, m);			// load 2
   addn(T, 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);
}