/******************************************************************************/
/* */
/* 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);
}