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