﻿ /*******************************
```/******************************************************************************/
/*									      */
/*  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);
if ((T[0]&0x80000000)==0) {
copyn(P, MAX, p);
}
copyn(P, T, p);	// find minimum
negn(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
negn(T, p);
if ((T[0]&0x80000000)==0)
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;
}