blob: 90ab4ba5a26318e0c012e45ce303959c97dca843 [file] [log] [blame]
 #include "os.h" #include #include "dat.h" /* division ala knuth, seminumerical algorithms, pp 237-238 */ /* the numbers are stored backwards to what knuth expects so j */ /* counts down rather than up. */ void mpdiv(mpint *dividend, mpint *divisor, mpint *quotient, mpint *remainder) { int j, s, vn, sign; mpdigit qd, *up, *vp, *qp; mpint *u, *v, *t; /* divide bv zero */ if(divisor->top == 0) abort(); /* quick check */ if(mpmagcmp(dividend, divisor) < 0){ if(remainder != nil) mpassign(dividend, remainder); if(quotient != nil) mpassign(mpzero, quotient); return; } /* D1: shift until divisor, v, has hi bit set (needed to make trial */ /* divisor accurate) */ qd = divisor->p[divisor->top-1]; for(s = 0; (qd & mpdighi) == 0; s++) qd <<= 1; u = mpnew((dividend->top+2)*Dbits + s); if(s == 0 && divisor != quotient && divisor != remainder) { mpassign(dividend, u); v = divisor; } else { mpleft(dividend, s, u); v = mpnew(divisor->top*Dbits); mpleft(divisor, s, v); } up = u->p+u->top-1; vp = v->p+v->top-1; vn = v->top; /* D1a: make sure high digit of dividend is less than high digit of divisor */ if(*up >= *vp){ *++up = 0; u->top++; } /* storage for multiplies */ t = mpnew(4*Dbits); qp = nil; if(quotient != nil){ mpbits(quotient, (u->top - v->top)*Dbits); quotient->top = u->top - v->top; qp = quotient->p+quotient->top-1; } /* D2, D7: loop on length of dividend */ for(j = u->top; j > vn; j--){ /* D3: calculate trial divisor */ mpdigdiv(up-1, *vp, &qd); /* D3a: rule out trial divisors 2 greater than real divisor */ if(vn > 1) for(;;){ memset(t->p, 0, 3*Dbytes); /* mpvecdigmuladd adds to what's there */ mpvecdigmuladd(vp-1, 2, qd, t->p); if(mpveccmp(t->p, 3, up-2, 3) > 0) qd--; else break; } /* D4: u -= v*qd << j*Dbits */ sign = mpvecdigmulsub(v->p, vn, qd, up-vn); if(sign < 0){ /* D6: trial divisor was too high, add back borrowed */ /* value and decrease divisor */ mpvecadd(up-vn, vn+1, v->p, vn, up-vn); qd--; } /* D5: save quotient digit */ if(qp != nil) *qp-- = qd; /* push top of u down one */ u->top--; *up-- = 0; } if(qp != nil){ mpnorm(quotient); if(dividend->sign != divisor->sign) quotient->sign = -1; } if(remainder != nil){ mpright(u, s, remainder); /* u is the remainder shifted */ remainder->sign = dividend->sign; } mpfree(t); mpfree(u); if(v != divisor) mpfree(v); }