/*****************************************************************************/
/* */
/* EUCLIDEAN G.C.D. */
/* 01/12/07 (dkc) */
/* */
/*****************************************************************************/
unsigned int euclid(unsigned int d, unsigned int e) {
unsigned int a,b,temp;
a=d;
b=e;
if (b>a) {
temp=a;
a=b;
b=temp;
}
loop:
temp=a-(a/b)*b;
a=b;
b=temp;
if (b!=0) goto loop;
return(a);
}
/******************************************************************************/
/* */
/* FIND APPROXIMATIONS TO THE PARITY VECTOR REQUIRED FOR A CYCLE */
/* 09/23/10 (dkc) */
/* */
/* This C program finds parity vectors for live limbs in S of a least- */
/* residue tree. */
/* */
/******************************************************************************/
#include <math.h>
unsigned int euclid(unsigned int a, unsigned int b);
void mul3shft(unsigned int *c, unsigned int b, unsigned int n);
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);
//
far unsigned int maxn=65536; // set maximum n value
//
far unsigned int error[4]; // error array
far unsigned int P[4096],MAX[4096],MIN[4096],T[4096]; // P, max, min arrays
far unsigned int sv[131072]; // parity vector
far unsigned int output[200]; // limb lengths (with exceptions)
far unsigned int output1[200]; // limb lengths
far unsigned int ln[200]; // l and n values
far unsigned int X[4096]; // X array
int main () {
unsigned int c,d,length,offset,temp,flag1,flag2,flag3,flag4,count;
unsigned int f,g,h,i,j,k,l,n,p,a,b,oldg,oldh,iters,factor,oldoff,deloff;
int del,delo,delg;
if (maxn>65536) { // set size of X array
error[0]=0xffffffff;
goto zskip;
}
if (maxn<=1024)
p=64;
else {
if (maxn<=2048)
p=128;
else {
if (maxn<=4096)
p=256;
else {
if (maxn<=8192)
p=512;
else {
if (maxn<=16384)
p=1024;
else {
if (maxn<=32768)
p=2048;
else
p=4096;
}
}
}
}
}
setn(sv, 0, 131072); // clear s.v. array
setn(P, 0, 4096); // clear P array
setn(X, 0, 4096); // clear X array
error[0]=0; // clear error array
error[1]=0;
error[2]=0;
error[3]=0;
iters=0; // clear exception count
for (i=0; i<200; i++) { // clear output arrays
output[i]=0;
output1[i]=0;
ln[i]=0;
}
count=0; // clear limb count
c=0; // clear c
d=0; // clear d
setn(X, 0, p);
X[0]=0x10000000; // x=1/2
for (i=0; i<maxn; i++) { // maximum n value
if (((X[0]+0xe0000000)&0x80000000)==0) { // x>1
mul3shft(X, 2, p); // x=(3/4)x
c=c+1; // increment c
}
else { // x<1
mul3shft(X, 1, p); // x=(3/2)x
d=d+1; // increment d
}
if (X[p-1]!=0) { // check for overflow of X array
error[2]=1;
goto zskip;
}
length=3*c+2*d; // length of limb
l=2*c+d+1; // length of cycle
if (l>=131072) { // check for overflow of s.v. array
error[2]=2;
goto zskip;
}
n=c+d; // number of odd elements in cycle
error[3]=n;
factor=euclid(l, n); // find G.C.D. of l and n
f=l/factor; // length of non-repeating portion of s.v.
//
// compute parity vector
//
flag1=0; // clear offset count
flag2=0; // clear exception flag
flag3=0; // clear first-time flag
flag4=0; // clear exception flag
oldg=0; // clear old index
oldh=0; // clear old index
oldoff=0; // clear old offset
for (offset=0; offset<l; offset++) {
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 address ing
if (temp>l)
temp=temp-l;
sv[temp]=k;
}
//
// check parity vector
//
if (sv[1]!=0)
continue;
if (sv[2]!=1)
continue;
if (sv[3]!=0)
continue;
flag1=flag1+1; // increment offset count
setn(MAX, 0, p); // set maximum
setn(MIN, 0, p); // set minimum
MIN[0]=0x7fffffff;
g=0;
h=0;
setn(P, 0, p);
P[0]=0x20000000; // set P
for (j=1; j<=l; j++) {
if (sv[j]==1) {
copyn(MAX, T, p); // find maximum
negn(T, p);
addn(P, T, p);
if ((T[0]&0x80000000)==0) {
copyn(P, MAX, p);
g=j;
}
copyn(P, T, p); // find minimum
negn(T, p);
addn(MIN, T, p);
if ((T[0]&0x80000000)==0) {
copyn(P, MIN, p);
h=j;
}
mul3shft(P, 1, p); // P=1.5*P
}
else
rshiftn(P, P, 1, p); // P=0.5*P
}
if (P[p-1]!=0) { // check for overflow of P array
error[2]=3;
goto zskip;
}
copyn(MIN, T, p); // compare twice the minimum to maximum
addn(T, T, p);
negn(T, p);
addn(MAX, T, p);
if ((T[0]&0x80000000)==0)
goto askip;
del=(int)(g-h); // difference in indices
delo=(int)(oldg-oldh); // difference in old indices
if (del<0)
del=del+l;
del=del-(del/f)*f;
if (delo<0)
delo=delo+l;
delo=delo-(delo/f)*f;
if ((del!=delo)&&(flag3!=0))
flag2=length; // set exception flag
if ((flag3!=0)&&(flag1>factor)) {
deloff=offset-oldoff; // difference in offsets
delg=g-oldg; // difference in indices
if (delg<0)
delg=delg+l;
delg=delg-(delg/f)*f;
if ((int)deloff!=delg)
flag4=length; // set exception flag
}
oldg=g; // set old index
oldh=h; // set old index
oldoff=offset; // set old offset
flag3=1; // set first-time flag
}
if (flag1!=0) {
error[0]=error[0]+1; // increment limb count
output1[count]=(length<<16)|factor; // store length and factor
ln[count]=(l<<16)|n; // store l and n values
count=count+1;
if (flag2==length) { // save exception
output[iters]=length;
iters=iters+1;
error[1]=iters;
if (iters>200)
goto zskip;
}
if (flag4==length) { // save exception
output[iters]=length|0x80000000;
iters=iters+1;
error[1]=iters;
if (iters>200)
goto zskip;
}
}
askip:
flag1=0;
}
zskip:
return(0);
}