mpdiv.c (2480B)
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 }