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