Subversion Repositories planix.SVN

Rev

Rev 2 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
2 - 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 
}