/******************************************************************************
*									      *
*  N-WORD ADD								      *
*  04/19/10 (dkc)							      *
*									      *
******************************************************************************/
unsigned int carry(unsigned int a, unsigned int b, unsigned int sum);
void addn(unsigned int *a, unsigned int *b, unsigned int n) {
unsigned int i,s;
unsigned int c[8192];
if (n>8192)
   return;
for (i=n-1; i>0; i--) {
   s=*(a+i)+*(b+i);
   c[i]=carry(*(a+i),*(b+i),s);
   *(b+i)=s;
   }
*b=*a+*b;

for (i=n-2; i>0; i--) {
   s=*(b+i)+c[i+1];
   c[i]+=carry(*(b+i),c[i+1],s);
   *(b+i)=s;
   }
*b=*b+c[1];
return;
}
/*****************************************************************************/
/*									     */
/*  CARRY								     */
/*  01/12/07 (dkc)							     */
/*									     */
/*****************************************************************************/
unsigned int carry(unsigned int a, unsigned int b, unsigned int sum) {
unsigned int c;
if ((a&b)>>31!=0)
   c=1;
else {
   if ((a|b)>>31==0)
      c=0;
   else {
      if (sum>=0x80000000)
	 c=0;
      else
	 c=1;
      }
   }
return c;
}
/******************************************************************************
*									      *
*  32*N-BIT COPY							      *
*  04/02/10 (dkc)							      *
*									      *
******************************************************************************/
void copyn(unsigned int *a, unsigned int *b, unsigned int n) {
unsigned int i;
for (i=0; i<n; i++)
   *(b+i)=*(a+i);
return;
}
/*****************************************************************************/
/*									     */
/*  N-WORD BIT DIVIDE (UNSIGNED, 32-BIT DIVISOR)			     */
/*  02/11/12 (dkc)							     */
/*									     */
/*****************************************************************************/
unsigned int carry(unsigned int a, unsigned int b, unsigned int sum);
unsigned int lmbd(unsigned int mode, unsigned int a);
unsigned int divn_32(unsigned int *an,unsigned int *quotient, unsigned int dn,
	       unsigned int n) {
unsigned int i,j,dshift,ashift,count,flag,shift,c0,c1,temp0;
unsigned int a[1024],d[1024],temp[1024];

if (n>1024)
   return(0);
temp0=0;
for (i=0; i<(n-1); i++) {
   a[i]=an[i];
   d[i]=0;
   temp0=temp0|a[i];
   }
a[n-1]=an[n-1];
if ((temp0==0)&&(a[n-1]<dn)) {
   a[n-1]=0;
   goto zskip;
   }
d[n-1]=dn;
//
// divisor shift count
//
dshift=lmbd(1,dn);
dshift+=32*(n-1);
//
// dividend shift count
//
for (i=0; i<n; i++) {
   if (a[i]!=0) {
      ashift=32*i+lmbd(1,a[i]);
      break;
      }
   }
//
// count
//
shift=dshift-ashift;
count=shift+1;
//
// align dividend and divisor
//
flag=shift/32;
shift=shift-flag*32;
flag=n-1-flag;
if (flag!=(n-1)) {
   d[flag]=d[n-1];
   d[n-1]=0;
   }
if (shift!=0) {
   d[flag-1]=d[flag]>>(32-shift);
   d[flag]=d[flag]<<shift;
   }
//
// negate divisor
//
for (i=0; i<n; i++)
   d[i]=~d[i];
c0=1;
for (i=0; i<n; i++) {
   temp0=d[n-1-i]+c0;
   c0=carry(d[n-1-i],c0,temp0);
   d[n-1-i]=temp0;
   }
//
//  do additions
//
for (i=0; i<count; i++) {
   temp0=a[n-1];
   c0=0;
   for (j=0; j<(n-1); j++) {
      temp[n-1-j]=temp0+d[n-1-j];
      c0+=carry(temp0,d[n-1-j],temp[n-1-j]);
      temp0=a[n-2-j]+c0;
      c1=c0;
      c0=carry(a[n-2-j],c0,temp0);
      }
   temp[0]=a[0]+d[0]+c1;

   if ((temp[0]>>31)==0) {
      for (j=0; j<(n-1); j++) {
	 a[j]=temp[j]<<1;
	 if ((temp[j+1]>>31)!=0)
	    c0=1;
	 else
	    c0=0;
	 a[j]=a[j]+c0;
	 }
      a[n-1]=temp[n-1]<<1;
      a[n-1]=a[n-1]+1;
      }
   else {
      for (j=0; j<(n-1); j++) {
	 a[j]=a[j]<<1;
	 if ((a[j+1]>>31)!=0)
	    c0=1;
	 else
	    c0=0;
	 a[j]=a[j]+c0;
	 }
      a[n-1]=a[n-1]<<1;
      }
   }
//
// shift off extra bits
// (modified to work for divisions with remainders on 08/26/13)
//
for (i=0; i<flag; i++)
   a[i]=0;
if (shift!=31) {
   shift=shift+1;
   shift=32-shift;
   temp0=a[flag]<<shift;
   a[flag]=temp0>>shift;
   }
//
// store quotient
//
zskip:
for (i=0; i<n; i++)
   *(quotient+i)=a[i];
return(1);
}
/*****************************************************************************/
/*									     */
/*  LEFT-MOST BIT DETECTION						     */
/*  01/12/07 (dkc)							     */
/*									     */
/*****************************************************************************/
unsigned int lmbd(unsigned int mode, unsigned int a) {
unsigned int i,mask,count;
if (mode==0)
   a=~a;
if (a==0)
   return(32);
count=0;
mask=0x80000000;
for (i=0; i<32; i++) {
   if ((a&mask)!=0)
      break;
   count+=1;
   mask>>=1;
   }
return(count);
}
/*****************************************************************************/
/*									     */
/*  LEFT-SHIFT 32*N BITS						     */
/*  10/16/09 (dkc)							     */
/*									     */
/*****************************************************************************/
void lshiftn(unsigned int *a, unsigned int *b, unsigned int shift, unsigned
	     int n) {
unsigned int i;
for (i=0; i<n; i++)
   *(b+i)=*(a+i);
while (shift>31) {
   for (i=0; i<n-1; i++)
      *(b+i)=*(b+i+1);
   *(b+n-1)=0;
   shift=shift-32;
   }
if (shift==0)
   return;
for (i=0; i<n-1; i++)
   *(b+i)=(*(b+i)<<shift)|(*(b+i+1)>>(32-shift));
*(b+n-1)=*(b+n-1)<<shift;
return;
}
/******************************************************************************
*									      *
*  NORMALIZE N WORDS							   *
*  04/03/10 (dkc)							      *
*									      *
******************************************************************************/
unsigned int normn(unsigned int *a, unsigned int n) {
unsigned int i,shift,mask;
i=0;
while ((*(a+i)==0)&&(i<n))
   i=i+1;
shift=32*i;
if (i==n)
   return shift;
mask=0x80000000;
while ((*(a+i)&mask)==0) {
   shift=shift+1;
   mask=mask>>1;
   }
return shift;
}
/******************************************************************************
*									      *
*  "OR" N WORDS TOGETHER                                                   *
*  04/02/10 (dkc)							      *
*									      *
******************************************************************************/
unsigned int orn(unsigned int *a, unsigned int n) {
unsigned int i,b;
b=*a;
for (i=1; i<n; i++)
   b=b|*(a+i);
return b;
}
/******************************************************************************
*									      *
*  SET THE LAST OF N WORDS						      *
*  04/02/10 (dkc)							      *
*									      *
******************************************************************************/
void setn(unsigned int *a, unsigned int b, unsigned int n) {
unsigned int i;
for (i=0; i<n-1; i++)
   *(a+i)=0;
*(a+n-1)=b;
return;
}
/*****************************************************************************/
/*									     */
/*  RIGHT-SHIFT 32*N BITS						     */
/*  04/02/10 (dkc)							     */
/*									     */
/*****************************************************************************/
void shiftn(unsigned int *a, unsigned int *b, unsigned int shift, unsigned
	    int n) {
unsigned int i;
for (i=0; i<n; i++)
   *(b+i)=*(a+i);
while (shift>31) {
   for (i=n-1; i>0; i--)
      *(b+i)=*(b+i-1);
   *b=0;
   shift=shift-32;
   }
if (shift==0)
   return;
for (i=n-1; i>0; i--)
   *(b+i)=(*(b+i)>>shift)|(*(b+i-1)<<(32-shift));
*b=*b>>shift;
return;
}
/******************************************************************************
*									      *
*  N-WORD SUBTRACT							      *
*  04/19/10 (dkc)							      *
*									      *
******************************************************************************/
unsigned int carry(unsigned int a, unsigned int b, unsigned int sum);
void subn(unsigned int *a, unsigned int *b, unsigned int n) {
unsigned int i,s,d;
unsigned int c[8192];
if (n>8192)
   return;
for (i=0; i<n; i++) {
   s=*(b+i);
   *(b+i)=~s;
   }

d=1;
for (i=n-1; i>0; i--) {
   s=*(b+i)+d;
   d=carry(*(b+i),d,s);
   *(b+i)=s;
   }
*b=*b+d;

for (i=n-1; i>0; i--) {
   s=*(a+i)+*(b+i);
   c[i]=carry(*(a+i),*(b+i),s);
   *(b+i)=s;
   }
*b=*a+*b;

for (i=n-2; i>0; i--) {
   s=*(b+i)+c[i+1];
   c[i]+=carry(*(b+i),c[i+1],s);
   *(b+i)=s;
   }
*b=*b+c[1];
return;
}