/******************************************************************************/ /* */ /* 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>
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); void mul32_32(unsigned int a, unsigned int m, unsigned int *p); void div64_32(unsigned int *a, unsigned int *q, unsigned int d); void sub64(unsigned int *a, unsigned int *b); far unsigned int error[6]; // error array far unsigned int P[12288],MAX[12288],MIN[12288],T[12288]; // P, max, min arrays far unsigned int X[12288]; // X array far unsigned char sv[262144]; // parity vector far unsigned int output[200]; // limb lengths far unsigned int ln[200]; // l and n values far unsigned int M[2],Q[2],U[2]; int main () { // unsigned int maxn=65536*2; // set maximum n value // unsigned int c,d,length,offset,flag1,count,index; unsigned int i,j,k,l,n,p,a,b; if (maxn>131072) { // 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 { if (maxn<=65536) p=4096; else p=8192; } } } } } } for (i=0; i<262144; i++) sv[i]=0; setn(P, 0, 8192); // clear P array setn(X, 0, 8192); // clear X array error[0]=0; // clear error array error[1]=0; error[2]=0; error[3]=0; error[4]=0; error[5]=0; for (i=0; i<200; i++) { // clear output arrays output[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>=262144) { // check for overflow of s.v. array error[2]=2; goto zskip; } n=c+d; // number of odd elements in cycle error[3]=n; // // compute parity vector // flag1=0; // clear offset count for (j=1; j<=l; j++) { mul32_32(j, n, M); div64_32(M, Q, l); mul32_32(Q[1], l, U); sub64(M, U); if ((U[0]==0)&&(U[1]==0)) a=Q[1]; else a=Q[1]+1; mul32_32(j-1, n, M); div64_32(M, Q, l); mul32_32(Q[1], l, U); sub64(M, U); if ((U[0]==0)&&(U[1]==0)) b=Q[1]; else b=Q[1]+1; k=a-b; // parity value if (k>0x20000000) { error[0]=10; goto zskip; } sv[j]=(unsigned char)k; } // // check parity vector // for (offset=1; offset<l; offset++) { if (sv[offset]!=0) continue; index=offset+1; if (index>l) index=index-l; if (sv[index]!=1) continue; index=offset+2; if (index>l) index=index-l; if (sv[index]!=0) continue; flag1=flag1+1; // increment offset count setn(MAX, 0, p); // set maximum setn(MIN, 0, p); // set minimum MIN[0]=0x7fffffff; setn(P, 0, p); P[0]=0x20000000; // set P index=offset; for (j=1; j<=l; j++) { if (sv[index]==1) { copyn(MAX, T, p); // find maximum negn(T, p); addn(P, T, p); if ((T[0]&0x80000000)==0) { copyn(P, MAX, p); } copyn(P, T, p); // find minimum negn(T, p); addn(MIN, T, p); if ((T[0]&0x80000000)==0) { copyn(P, MIN, p); } mul3shft(P, 1, p); // P=1.5*P } else rshiftn(P, P, 1, p); // P=0.5*P index=index+1; if (index>l) index=index-l; } 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; error[4]=error[4]+1; } if (flag1!=0) { error[0]=error[0]+1; // increment limb count output[count]=length; // store length and factor ln[count]=(l<<16)|n; // store l and n values count=count+1; } askip: flag1=0; } zskip: return(0); }