/******************************************************************************/
/*									      */
/*  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);
	    addn(A, A, m);
	    addn(B, A, 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(U, 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(U, U, m);
      addn(U, U, m);
      addn(U, U, m);
      addn(O, U, m);	 // 9*M
      setn(V, 1, m);
      lshiftn(V, W, oldl-1, m);
      copyn(W, V, m);
      addn(V, 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])
	    goto askip;
	 }
      savek=k;
      goto bskip;
askip:
      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);
      addn(A, A, m);
      addn(B, A, 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);
}