﻿ n-word arithmetic subroutines
```/******************************************************************************
*									      *
*  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;
}
//
//
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) {
if (mode==0)
a=~a;
if (a==0)
return(32);
count=0;
for (i=0; i<32; i++) {
break;
count+=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) {
i=0;
while ((*(a+i)==0)&&(i<n))
i=i+1;
shift=32*i;
if (i==n)
return shift;
shift=shift+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;
}
```