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