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 sysfatal("mpdiv: divide by zero"); 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 } 113