1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
|
#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);
}
|