1 b3f61791 2004-03-21 devnull #include "os.h"
  2 b3f61791 2004-03-21 devnull #include <mp.h>
  3 b3f61791 2004-03-21 devnull #include "dat.h"
  4 b3f61791 2004-03-21 devnull 
  5 cbeb0b26 2006-04-01 devnull /* division ala knuth, seminumerical algorithms, pp 237-238 */
  6 cbeb0b26 2006-04-01 devnull /* the numbers are stored backwards to what knuth expects so j */
  7 cbeb0b26 2006-04-01 devnull /* counts down rather than up. */
  8 b3f61791 2004-03-21 devnull 
  9 b3f61791 2004-03-21 devnull void
 10 b3f61791 2004-03-21 devnull mpdiv(mpint *dividend, mpint *divisor, mpint *quotient, mpint *remainder)
 11 b3f61791 2004-03-21 devnull {
 12 b3f61791 2004-03-21 devnull 	int j, s, vn, sign;
 13 b3f61791 2004-03-21 devnull 	mpdigit qd, *up, *vp, *qp;
 14 b3f61791 2004-03-21 devnull 	mpint *u, *v, *t;
 15 b3f61791 2004-03-21 devnull 
 16 cbeb0b26 2006-04-01 devnull 	/* divide bv zero */
 17 b3f61791 2004-03-21 devnull 	if(divisor->top == 0)
 18 b3f61791 2004-03-21 devnull 		abort();
 19 b3f61791 2004-03-21 devnull 
 20 cbeb0b26 2006-04-01 devnull 	/* quick check */
 21 b3f61791 2004-03-21 devnull 	if(mpmagcmp(dividend, divisor) < 0){
 22 b3f61791 2004-03-21 devnull 		if(remainder != nil)
 23 b3f61791 2004-03-21 devnull 			mpassign(dividend, remainder);
 24 b3f61791 2004-03-21 devnull 		if(quotient != nil)
 25 b3f61791 2004-03-21 devnull 			mpassign(mpzero, quotient);
 26 b3f61791 2004-03-21 devnull 		return;
 27 b3f61791 2004-03-21 devnull 	}
 28 b3f61791 2004-03-21 devnull 
 29 cbeb0b26 2006-04-01 devnull 	/* D1: shift until divisor, v, has hi bit set (needed to make trial */
 30 cbeb0b26 2006-04-01 devnull 	/*     divisor accurate) */
 31 b3f61791 2004-03-21 devnull 	qd = divisor->p[divisor->top-1];
 32 b3f61791 2004-03-21 devnull 	for(s = 0; (qd & mpdighi) == 0; s++)
 33 b3f61791 2004-03-21 devnull 		qd <<= 1;
 34 b3f61791 2004-03-21 devnull 	u = mpnew((dividend->top+2)*Dbits + s);
 35 b3f61791 2004-03-21 devnull 	if(s == 0 && divisor != quotient && divisor != remainder) {
 36 b3f61791 2004-03-21 devnull 		mpassign(dividend, u);
 37 b3f61791 2004-03-21 devnull 		v = divisor;
 38 b3f61791 2004-03-21 devnull 	} else {
 39 b3f61791 2004-03-21 devnull 		mpleft(dividend, s, u);
 40 b3f61791 2004-03-21 devnull 		v = mpnew(divisor->top*Dbits);
 41 b3f61791 2004-03-21 devnull 		mpleft(divisor, s, v);
 42 b3f61791 2004-03-21 devnull 	}
 43 b3f61791 2004-03-21 devnull 	up = u->p+u->top-1;
 44 b3f61791 2004-03-21 devnull 	vp = v->p+v->top-1;
 45 b3f61791 2004-03-21 devnull 	vn = v->top;
 46 b3f61791 2004-03-21 devnull 
 47 cbeb0b26 2006-04-01 devnull 	/* D1a: make sure high digit of dividend is less than high digit of divisor */
 48 b3f61791 2004-03-21 devnull 	if(*up >= *vp){
 49 b3f61791 2004-03-21 devnull 		*++up = 0;
 50 b3f61791 2004-03-21 devnull 		u->top++;
 51 b3f61791 2004-03-21 devnull 	}
 52 b3f61791 2004-03-21 devnull 
 53 cbeb0b26 2006-04-01 devnull 	/* storage for multiplies */
 54 b3f61791 2004-03-21 devnull 	t = mpnew(4*Dbits);
 55 b3f61791 2004-03-21 devnull 
 56 b3f61791 2004-03-21 devnull 	qp = nil;
 57 b3f61791 2004-03-21 devnull 	if(quotient != nil){
 58 b3f61791 2004-03-21 devnull 		mpbits(quotient, (u->top - v->top)*Dbits);
 59 b3f61791 2004-03-21 devnull 		quotient->top = u->top - v->top;
 60 b3f61791 2004-03-21 devnull 		qp = quotient->p+quotient->top-1;
 61 b3f61791 2004-03-21 devnull 	}
 62 b3f61791 2004-03-21 devnull 
 63 cbeb0b26 2006-04-01 devnull 	/* D2, D7: loop on length of dividend */
 64 b3f61791 2004-03-21 devnull 	for(j = u->top; j > vn; j--){
 65 b3f61791 2004-03-21 devnull 
 66 cbeb0b26 2006-04-01 devnull 		/* D3: calculate trial divisor */
 67 b3f61791 2004-03-21 devnull 		mpdigdiv(up-1, *vp, &qd);
 68 b3f61791 2004-03-21 devnull 
 69 cbeb0b26 2006-04-01 devnull 		/* D3a: rule out trial divisors 2 greater than real divisor */
 70 b3f61791 2004-03-21 devnull 		if(vn > 1) for(;;){
 71 cbeb0b26 2006-04-01 devnull 			memset(t->p, 0, 3*Dbytes);	/* mpvecdigmuladd adds to what's there */
 72 b3f61791 2004-03-21 devnull 			mpvecdigmuladd(vp-1, 2, qd, t->p);
 73 b3f61791 2004-03-21 devnull 			if(mpveccmp(t->p, 3, up-2, 3) > 0)
 74 b3f61791 2004-03-21 devnull 				qd--;
 75 b3f61791 2004-03-21 devnull 			else
 76 b3f61791 2004-03-21 devnull 				break;
 77 b3f61791 2004-03-21 devnull 		}
 78 b3f61791 2004-03-21 devnull 
 79 cbeb0b26 2006-04-01 devnull 		/* D4: u -= v*qd << j*Dbits */
 80 b3f61791 2004-03-21 devnull 		sign = mpvecdigmulsub(v->p, vn, qd, up-vn);
 81 b3f61791 2004-03-21 devnull 		if(sign < 0){
 82 b3f61791 2004-03-21 devnull 
 83 cbeb0b26 2006-04-01 devnull 			/* D6: trial divisor was too high, add back borrowed */
 84 cbeb0b26 2006-04-01 devnull 			/*     value and decrease divisor */
 85 b3f61791 2004-03-21 devnull 			mpvecadd(up-vn, vn+1, v->p, vn, up-vn);
 86 b3f61791 2004-03-21 devnull 			qd--;
 87 b3f61791 2004-03-21 devnull 		}
 88 b3f61791 2004-03-21 devnull 
 89 cbeb0b26 2006-04-01 devnull 		/* D5: save quotient digit */
 90 b3f61791 2004-03-21 devnull 		if(qp != nil)
 91 b3f61791 2004-03-21 devnull 			*qp-- = qd;
 92 b3f61791 2004-03-21 devnull 
 93 cbeb0b26 2006-04-01 devnull 		/* push top of u down one */
 94 b3f61791 2004-03-21 devnull 		u->top--;
 95 b3f61791 2004-03-21 devnull 		*up-- = 0;
 96 b3f61791 2004-03-21 devnull 	}
 97 b3f61791 2004-03-21 devnull 	if(qp != nil){
 98 b3f61791 2004-03-21 devnull 		mpnorm(quotient);
 99 b3f61791 2004-03-21 devnull 		if(dividend->sign != divisor->sign)
100 b3f61791 2004-03-21 devnull 			quotient->sign = -1;
101 b3f61791 2004-03-21 devnull 	}
102 b3f61791 2004-03-21 devnull 
103 b3f61791 2004-03-21 devnull 	if(remainder != nil){
104 cbeb0b26 2006-04-01 devnull 		mpright(u, s, remainder);	/* u is the remainder shifted */
105 b3f61791 2004-03-21 devnull 		remainder->sign = dividend->sign;
106 b3f61791 2004-03-21 devnull 	}
107 b3f61791 2004-03-21 devnull 
108 b3f61791 2004-03-21 devnull 	mpfree(t);
109 b3f61791 2004-03-21 devnull 	mpfree(u);
110 b3f61791 2004-03-21 devnull 	if(v != divisor)
111 b3f61791 2004-03-21 devnull 		mpfree(v);
112 b3f61791 2004-03-21 devnull }