﻿ /*******************************
```/******************************************************************************/
/*									      */
/*  COMPUTE THE MINIMUM ELEMENT IN A LOOP (for M-cycles generated from the    */
/*					   parity vector p)		      */
/*  04/03/10 (dkc)							      */
/*									      */
/******************************************************************************/
#include <stdio.h>
#include <math.h>
#include "ln.h"
#include "sv.h"
void shiftn(unsigned int *a, unsigned int *b, unsigned int shift, unsigned int n);
void lshiftn(unsigned int *a, unsigned int *b, unsigned int shift, unsigned int n);
void div128_64(unsigned int a0, unsigned int a1, unsigned int a2,
unsigned int a3, unsigned int *quotient, unsigned int d2,
unsigned int d3);
void addn(unsigned int *a, unsigned int *b, unsigned int n);
void subn(unsigned int *c, unsigned int *d, 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 shiftn(unsigned int *a, unsigned int *b, unsigned int shift, unsigned int n);
unsigned int normn(unsigned int *a, unsigned int n);
unsigned int orn(unsigned int *a, unsigned int n);

int main () {
unsigned int m=512;			// number of words
unsigned int big=1;			// matrix width flag
int g;
//unsigned int k,savek;
unsigned int pv[2000];
unsigned int h,i,j,l,n,a,b,oldl,oldn,oldsgn,oldm,count3,w;
unsigned int A[1024],B[1024],C[4],S[1024],T[1024],O[1024],flag,flag1,o;
//unsigned int U[1024],V[1024],W[1024];
unsigned char v[79*250];
FILE *Outfp;
Outfp = fopen("out14da.dat","w");
for (i=0; i<79*250; i++)
v[i]=2;
if (big==0)
w=39;
else
w=79;
setn(T, 0, m);
oldl=1;
oldn=1;
oldsgn=0;
oldm=0x7fffffff;
setn(O, 1, m);
for (i=0; i<600; i++) {
l=ln[i]>>16; 	// length of loop
n=ln[i]&0xffff;	    // number of odd elements in loop
setn(S, 0, m);			// clear sum
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;
pv[j-1]=a-b;
if (a!=b) {
setn(A, a-b, m);
for (h=0; h<(n-a); h++) {	   // (a-b)*3**(n-a)
copyn(A, B, m);
}
lshiftn(A, B, j-1, m); 	   // (a-b)*2**(j-1)*3**(n-a)
addn(B, S, m); 		   // sum+=(a-b)*2**(j-1)*3**(n-a)
if ((S[0]&0xc0000000)!=0) {
printf("error A: sum too large \n");
goto zskip;
}
}
}
if (i<w) {
for (j=0; j<l; j++)
v[w*j+i]=(unsigned char)pv[j];
}
//
// check relationships (usually not valid)
//
/* if ((l-oldl)==2) {
copyn(O, U, m);
addn(O, U, m);	 // 3*M
setn(V, 1, m);
lshiftn(V, W, oldl, m);
addn(W, U, m);	 // 3*M+2^l
subn(S, U, m);
o=orn(U, m);
if (o!=0) {
printf("error: l=%d, n=%d \n",l,n);
printf("M=%d %d %d %d \n",O[m-4],O[m-3],O[m-2],O[m-1]);
printf("S=%d %d %d %d \n",S[m-4],S[m-3],S[m-2],S[m-1]);
goto zskip;
}
}
else {
copyn(O, U, m);
addn(O, U, m);	 // 9*M
setn(V, 1, m);
lshiftn(V, W, oldl-1, m);
copyn(W, V, m);
addn(W, V, m);	 // 3*2^(l-1)
addn(V, U, m);	 // 9*M+3*2^(l-1)
setn(V, 1, m);
lshiftn(V, W, oldl+1, m);
addn(W, U, m);	 // 9*M+3*2^(l-1)+2^(l+1)
subn(S, U, m);
o=orn(U, m);
if (o!=0) {
printf("error: l=%d, n=%d \n",l,n);
printf("M=%d %d %d %d \n",O[m-4],O[m-3],O[m-2],O[m-1]);
printf("S=%d %d %d %d \n",S[m-4],S[m-3],S[m-2],S[m-1]);
goto zskip;
}
}
copyn(S, O, m);
*/
//
// check if a rotation of parity vector matches parity vector p (it doesn't)
//
/*
savek=10000;
for (k=0; k<l; k++) {
for (j=0; j<l; j++) {
if (pv[j]!=sv[j])
}
savek=k;
goto bskip;
a=pv[0];
for (j=0; j<(l-1); j++)
pv[j]=pv[j+1];
pv[l-1]=a;
}
printf("error: unexpected parity vector, l=%d, n=%d, k=%d \n",l,n,savek);
for (j=0; j<l; j++)
printf(" %d %d \n",sv[j],pv[j]);
goto zskip;
bskip:
*/
setn(A, 1, m);
for (h=0; h<n; h++) {		// 3**n
copyn(A, B, m);
if ((A[0]&0xc0000000)!=0) {
printf("error: n=%d \n",n);
goto zskip;
}
}
if ((unsigned int)l>=(32*m)) {
printf("error: l=%d \n",l);
goto zskip;
}
setn(B, 1, m);
lshiftn(B, B, l, m); 		// 2**l
subn(B, A, m);			// 2**l-3**n
flag=0;
if ((A[0]&0x80000000)!=0) {		// negate 2**l-3**n
setn(B, 0, m);
subn(B, A, m);
flag=1;
}
copyn(T, B, m);
subn(A, B, m);
if ((B[0]&0x80000000)!=0) {
printf("not monotonic \n");
fprintf(Outfp,"not monotonic \n");
//    goto zskip;
}
copyn(A, T, m);
if ((S[0]&0xff000000)!=0) {
printf("error: sum too big \n");
goto zskip;
}
lshiftn(S, S, 8, m);
flag1=0;
o=normn(A, m);			// count to normalize divisor
g=(int)(32*m)-(int)o-62;		// right shift amount
if (g>0) {
shiftn(S, S, g, m);
shiftn(A, A, g, m);
flag1=(unsigned int)g;
}
o=orn(S, m-4);
if ((o|(S[m-4]&0xc0000000))!=0) {
printf("error B: sum too large \n");
goto zskip;
}
div128_64(S[m-4],S[m-3],S[m-2],S[m-1],C,A[m-2],A[m-1]);    // sum/(2**l-3**n)
shiftn(C, C, 8, 4);
if ((C[0]|C[1]|C[2])!=0) {
printf("warning: x too large \n");
goto zskip;
}
printf("l=%d, n=%d, sign=%d, shift=%d, x=%d %d %d %d \n",l,n,flag,flag1,C[0],C[1],C[2],C[3]);
fprintf(Outfp,"l=%d, n=%d, sign=%d, shift=%d, x=%d %d %d %d \n",l,n,flag,flag1,C[0],C[1],C[2],C[3]);
if ((l-oldl)==3) {
count3=count3+1;
if (count3>3) {
printf("error: l=%d, oldl=%d \n",l,oldl);
goto zskip;
}
if ((count3<2)&&((oldsgn!=0)||(flag!=0)||(C[3]<oldm))) {
printf("error: incorrect sign \n");
goto zskip;
}
if ((count3==3)&&((oldsgn!=0)||(flag==0)||(C[3]>oldm))) {
printf("error: incorrect sign \n");
goto zskip;
}
}
else
count3=0;
if (((l-oldl)==2)&&(C[3]>oldm)) {
printf("error: incorrect sizes \n");
goto zskip;
}
oldl=l;
oldn=n;
oldsgn=flag;
oldm=C[3];
}
printf("\n");
fprintf(Outfp,"\n");
if (big!=0) {
for (i=0; i<250; i++) {
printf("%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d\n",
v[i*79],v[i*79+1],v[i*79+2],v[i*79+3],v[i*79+4],v[i*79+5],v[i*79+6],v[i*79+7],v[i*79+8],v[i*79+9],
v[i*79+10],v[i*79+11],v[i*79+12],v[i*79+13],v[i*79+14],v[i*79+15],v[i*79+16],v[i*79+17],v[i*79+18],v[i*79+19],v[i*79+20],v[i*79+21],v[i*79+22],v[i*79+23],v[i*79+24],
v[i*79+25],v[i*79+26],v[i*79+27],v[i*79+28],v[i*79+29],v[i*79+30],v[i*79+31],v[i*79+32],v[i*79+33],v[i*79+34],v[i*79+35],v[i*79+36],v[i*79+37],v[i*79+38],v[i*79+39],
v[i*79+40],v[i*79+41],v[i*79+42],v[i*79+43],v[i*79+44],v[i*79+45],v[i*79+46],v[i*79+47],v[i*79+48],v[i*79+49],v[i*79+50],v[i*79+51],v[i*79+52],v[i*79+53],v[i*79+54],
v[i*79+55],v[i*79+56],v[i*79+57],v[i*79+58],v[i*79+59],v[i*79+60],v[i*79+61],v[i*79+62],v[i*79+63],v[i*79+64],v[i*79+65],v[i*79+66],v[i*79+67],v[i*79+68],v[i*79+69],
v[i*79+70],v[i*79+71],v[i*79+72],v[i*79+73],v[i*79+74],v[i*79+75],v[i*79+76],v[i*79+77],v[i*79+78],sv[i]);
fprintf(Outfp,"%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d\n",
v[i*79],v[i*79+1],v[i*79+2],v[i*79+3],v[i*79+4],v[i*79+5],v[i*79+6],v[i*79+7],v[i*79+8],v[i*79+9],
v[i*79+10],v[i*79+11],v[i*79+12],v[i*79+13],v[i*79+14],v[i*79+15],v[i*79+16],v[i*79+17],v[i*79+18],v[i*79+19],v[i*79+20],v[i*79+21],v[i*79+22],v[i*79+23],v[i*79+24],
v[i*79+25],v[i*79+26],v[i*79+27],v[i*79+28],v[i*79+29],v[i*79+30],v[i*79+31],v[i*79+32],v[i*79+33],v[i*79+34],v[i*79+35],v[i*79+36],v[i*79+37],v[i*79+38],v[i*79+39],
v[i*79+40],v[i*79+41],v[i*79+42],v[i*79+43],v[i*79+44],v[i*79+45],v[i*79+46],v[i*79+47],v[i*79+48],v[i*79+49],v[i*79+50],v[i*79+51],v[i*79+52],v[i*79+53],v[i*79+54],
v[i*79+55],v[i*79+56],v[i*79+57],v[i*79+58],v[i*79+59],v[i*79+60],v[i*79+61],v[i*79+62],v[i*79+63],v[i*79+64],v[i*79+65],v[i*79+66],v[i*79+67],v[i*79+68],v[i*79+69],
v[i*79+70],v[i*79+71],v[i*79+72],v[i*79+73],v[i*79+74],v[i*79+75],v[i*79+76],v[i*79+77],v[i*79+78],sv[i]);
}
}
else {
for (i=0; i<125; i++) {
printf(" %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d \n",
v[i*39],v[i*39+1],v[i*39+2],v[i*39+3],v[i*39+4],v[i*39+5],v[i*39+6],v[i*39+7],v[i*39+8],v[i*39+9],
v[i*39+10],v[i*39+11],v[i*39+12],v[i*39+13],v[i*39+14],v[i*39+15],v[i*39+16],v[i*39+17],v[i*39+18],v[i*39+19],v[i*39+20],v[i*39+21],v[i*39+22],v[i*39+23],v[i*39+24],
v[i*39+25],v[i*39+26],v[i*39+27],v[i*39+28],v[i*39+29],v[i*39+30],v[i*39+31],v[i*39+32],v[i*39+33],v[i*39+34],v[i*39+35],v[i*39+36],v[i*39+37],v[i*39+38]);
fprintf(Outfp," %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d \n",
v[i*39],v[i*39+1],v[i*39+2],v[i*39+3],v[i*39+4],v[i*39+5],v[i*39+6],v[i*39+7],v[i*39+8],v[i*39+9],
v[i*39+10],v[i*39+11],v[i*39+12],v[i*39+13],v[i*39+14],v[i*39+15],v[i*39+16],v[i*39+17],v[i*39+18],v[i*39+19],v[i*39+20],v[i*39+21],v[i*39+22],v[i*39+23],v[i*39+24],
v[i*39+25],v[i*39+26],v[i*39+27],v[i*39+28],v[i*39+29],v[i*39+30],v[i*39+31],v[i*39+32],v[i*39+33],v[i*39+34],v[i*39+35],v[i*39+36],v[i*39+37],v[i*39+38]);
}
}
zskip:
fclose(Outfp);
return(0);
}
```