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