/*****************************************************************************/
/*									     */
/*  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);
}