Subversion Repositories planix.SVN

Rev

Rev 33 | 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 
{
33 7u83 12 
	int j, s, vn, sign, qsign, rsign;
2 - 13 
	mpdigit qd, *up, *vp, *qp;
14 
	mpint *u, *v, *t;
15 
 
33 7u83 16 
	assert(quotient != remainder);
17 
	assert(divisor->flags & MPnorm);
18 
 
2 - 19 
	// divide bv zero
20 
	if(divisor->top == 0)
21 
		abort();
22 
 
33 7u83 23 
	// division by one or small powers of two
24 
	if(divisor->top == 1 && (divisor->p[0] & divisor->p[0]-1) == 0){
25 
		vlong r = 0;
26 
		if(dividend->top > 0)
27 
			r = (vlong)dividend->sign * (dividend->p[0] & divisor->p[0]-1);
28 
		if(quotient != nil){
29 
			sign = divisor->sign;
30 
			for(s = 0; ((divisor->p[0] >> s) & 1) == 0; s++)
31 
				;
32 
			mpright(dividend, s, quotient);
33 
			if(sign < 0)
34 
				quotient->sign ^= (-mpmagcmp(quotient, mpzero) >> 31) << 1;
35 
		}
36 
		if(remainder != nil){
37 
			remainder->flags |= dividend->flags & MPtimesafe;
38 
			vtomp(r, remainder);
39 
		}
40 
		return;
41 
	}
42 
	assert((dividend->flags & MPtimesafe) == 0);
43 
 
2 - 44 
	// quick check
45 
	if(mpmagcmp(dividend, divisor) < 0){
46 
		if(remainder != nil)
47 
			mpassign(dividend, remainder);
48 
		if(quotient != nil)
49 
			mpassign(mpzero, quotient);
50 
		return;
51 
	}
33 7u83 52 
 
53 
	qsign = divisor->sign * dividend->sign;
54 
	rsign = dividend->sign;
2 - 55 
 
56 
	// D1: shift until divisor, v, has hi bit set (needed to make trial
57 
	//     divisor accurate)
58 
	qd = divisor->p[divisor->top-1];
59 
	for(s = 0; (qd & mpdighi) == 0; s++)
60 
		qd <<= 1;
61 
	u = mpnew((dividend->top+2)*Dbits + s);
62 
	if(s == 0 && divisor != quotient && divisor != remainder) {
63 
		mpassign(dividend, u);
64 
		v = divisor;
65 
	} else {
66 
		mpleft(dividend, s, u);
67 
		v = mpnew(divisor->top*Dbits);
68 
		mpleft(divisor, s, v);
69 
	}
70 
	up = u->p+u->top-1;
71 
	vp = v->p+v->top-1;
72 
	vn = v->top;
73 
 
74 
	// D1a: make sure high digit of dividend is less than high digit of divisor
75 
	if(*up >= *vp){
76 
		*++up = 0;
77 
		u->top++;
78 
	}
79 
 
80 
	// storage for multiplies
81 
	t = mpnew(4*Dbits);
82 
 
83 
	qp = nil;
84 
	if(quotient != nil){
85 
		mpbits(quotient, (u->top - v->top)*Dbits);
86 
		quotient->top = u->top - v->top;
87 
		qp = quotient->p+quotient->top-1;
88 
	}
89 
 
90 
	// D2, D7: loop on length of dividend
91 
	for(j = u->top; j > vn; j--){
92 
 
93 
		// D3: calculate trial divisor
94 
		mpdigdiv(up-1, *vp, &qd);
95 
 
96 
		// D3a: rule out trial divisors 2 greater than real divisor
97 
		if(vn > 1) for(;;){
98 
			memset(t->p, 0, 3*Dbytes);	// mpvecdigmuladd adds to what's there
99 
			mpvecdigmuladd(vp-1, 2, qd, t->p);
100 
			if(mpveccmp(t->p, 3, up-2, 3) > 0)
101 
				qd--;
102 
			else
103 
				break;
104 
		}
105 
 
106 
		// D4: u -= v*qd << j*Dbits
107 
		sign = mpvecdigmulsub(v->p, vn, qd, up-vn);
108 
		if(sign < 0){
109 
 
110 
			// D6: trial divisor was too high, add back borrowed
111 
			//     value and decrease divisor
112 
			mpvecadd(up-vn, vn+1, v->p, vn, up-vn);
113 
			qd--;
114 
		}
115 
 
116 
		// D5: save quotient digit
117 
		if(qp != nil)
118 
			*qp-- = qd;
119 
 
120 
		// push top of u down one
121 
		u->top--;
122 
		*up-- = 0;
123 
	}
124 
	if(qp != nil){
33 7u83 125 
		assert((quotient->flags & MPtimesafe) == 0);
2 - 126 
		mpnorm(quotient);
33 7u83 127 
		if(quotient->top != 0)
128 
			quotient->sign = qsign;
2 - 129 
	}
130 
 
131 
	if(remainder != nil){
33 7u83 132 
		assert((remainder->flags & MPtimesafe) == 0);
2 - 133 
		mpright(u, s, remainder);	// u is the remainder shifted
33 7u83 134 
		if(remainder->top != 0)
135 
			remainder->sign = rsign;
2 - 136 
	}
137 
 
138 
	mpfree(t);
139 
	mpfree(u);
140 
	if(v != divisor)
141 
		mpfree(v);
142 
}