123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112 |
- #include "os.h"
- #include <mp.h>
- #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);
- }
|