1 /*
2 * bignum_div.c
3 * BigMath
4 *
5 * Created by Curtis Jones on 2007.08.05.
6 * Copyright 2007 __MyCompanyName__. All rights reserved.
7 *
8 */
9
10 #include "bignum.h"
11
12 // TODO: make sure that divs->oofm == 1
13 //
14 bignum_t *
15 bignum_div_aux (bignum_t *quot, bignum_t *rmdr, bignum_t *divd, bignum_t *divs)
16 {
17 u_int64_t r = 0;
18 u_int64_t n = divd->oofm;
19 u_int64_t j = n - 1;
20 u_int64_t b = divd->base;
21
22 if (divd->sign == BIGNUM_NEGATIVE) {
23 printf("%s().. divd is negative! ", __PRETTY_FUNCTION__); bignum_print(divd);
24 return NULL;
25 }
26
27 if (divs->sign == BIGNUM_NEGATIVE) {
28 printf("%s().. divs is negative! ", __PRETTY_FUNCTION__); bignum_print(divs);
29 return NULL;
30 }
31
32 if (quot != NULL)
33 __bignum_zero(quot);
34
35 do {
36 u_int64_t y = r * b + (u_int64_t)divd->numb[j];
37
38 if (quot != NULL) {
39 __bignum_set(quot, j, y / (u_int64_t)divs->numb[0]);
40
41 if (quot->numb[j] != 0)
42 quot->oofm = BIGNUM_MAX(quot->oofm, (j+1));
43 }
44
45 r = y % (u_int64_t)divs->numb[0];
46 }
47 while (--j >= 0 && j < n);
48
49 if (rmdr != NULL) {
50 __bignum_set(rmdr, 0, r);
51
52 if (rmdr->numb[0] != 0)
53 rmdr->oofm = 1;
54 }
55
56 return quot;
57 }
58
59 // Knuth 4.3.1 - Algorithm D
60 //
61 // a / b = c
62 //
63 // dividend (a) / divisor (b) = quotient (c)
64 //
65 bignum_t *
66 bignum_div (bignum_t *quot, bignum_t *rmdr, bignum_t *divd, bignum_t *divs)
67 {
68 u_int64_t i = 0;
69 u_int64_t m = divd->oofm - divs->oofm;
70 u_int64_t n = divs->oofm;
71 u_int64_t b = divd->base;
72 u_int64_t d = ( b / ((u_int64_t)divs->numb[divs->oofm-1] + 1) );
73 bignum_t *u = NULL, *v = NULL;
74 u_int32_t divd_sign, divs_sign, quot_sign;
75
76 // set divd to positive
77 // set divs to positive
78 //
79 // if divd is negative xor divs is negative
80 // quot is negative
81 //
82 {
83 if (BIGNUM_NEGATIVE == (divd_sign = divd->sign) ^ BIGNUM_NEGATIVE == (divs_sign = divs->sign))
84 quot_sign = BIGNUM_NEGATIVE;
85 else
86 quot_sign = BIGNUM_POSITIVE;
87
88 divd->sign = BIGNUM_POSITIVE;
89 divs->sign = BIGNUM_POSITIVE;
90 }
91
92 if (quot != NULL)
93 __bignum_zero( quot );
94
95 if (rmdr != NULL)
96 __bignum_zero( rmdr );
97
98 // if the dividend is less than the divisor then the
99 // quotient is zero and the remainder is the dividend.
100 //
101 if (bignum_lt(divd, divs)) {
102 if (quot != NULL)
103 bignum_set(quot, 0);
104
105 if (rmdr != NULL)
106 bignum_set_bn(rmdr, divd);
107
108 divd->sign = divd_sign;
109 divs->sign = divs_sign;
110
111 return quot;
112 }
113
114 // if the dividend is equal to the divisor then the
115 // quotient is one and the remainder is zero.
116 //
117 else if (bignum_eq(divd, divs)) {
118 if (quot != NULL)
119 bignum_set(quot, 1);
120
121 if (rmdr != NULL)
122 bignum_set(rmdr, 0);
123
124 divd->sign = divd_sign;
125 divs->sign = divs_sign;
126
127 return quot;
128 }
129
130 // TODO: It is possible that divs->oofm is greater than
131 // one, but that the value is still wholly contained in
132 // the first magnitude. we could test the value of divs
133 // relative to the maximum value that can be held in
134 // a single magnitude....
135 //
136 else if (divs->oofm == 1) {
137 bignum_div_aux(quot, rmdr, divd, divs);
138
139 divd->sign = divd_sign;
140 divs->sign = divs_sign;
141
142 return quot;
143 }
144
145 if (NULL == (u = bignum_alloc())) {
146 printf("%s().. [1] bignum_alloc failed\n", __PRETTY_FUNCTION__);
147 return NULL;
148 }
149
150 if (NULL == (v = bignum_alloc())) {
151 printf("%s().. [2] bignum_alloc failed\n", __PRETTY_FUNCTION__);
152 bignum_free( u );
153 return NULL;
154 }
155
156 // STEP D1 - Normalize
157 //
158 // Set u = divisor * d
159 // Set v = dividend * d
160 //
161 // If u is still of the same magnitude, then add
162 // an additional zero to the end - the high end.
163 // This does not change its value, but the additional
164 // space is required.
165 //
166 {
167 bignum_t * tmpd = bignum_alloc();
168
169 bignum_set(tmpd, d);
170
171 bignum_mul(u, divd, tmpd);
172 bignum_mul(v, divs, tmpd);
173
174 if (u->oofm == divd->oofm)
175 u->oofm++;
176
177 bignum_free( tmpd );
178 }
179
180 // STEP D2 - Initialize j
181 //
182 // Set j to m which is the difference in magnitude between
183 // the dividend and the divisor. We'll be decrementing j
184 // each time through this loop.
185 //
186 {
187 u_int64_t j = m;
188
189 bignum_t * bn_q = bignum_alloc();
190 bignum_t * bn_tmp1 = bignum_alloc();
191 bignum_t * bn_tmp2 = bignum_alloc();
192 bignum_t * bn_tmpu = bignum_alloc();
193
194 // STEP D3 - Part 1 - Calculate q
195 //
196 // q is going to be the j'th digit of the quotient. first
197 // we estimate q and then possibly alter it slightly if
198 // our estimate is off by one or two.
199 //
200 do {
201 u_int64_t x = (u_int64_t)u->numb[j+n] * b + (u_int64_t)u->numb[j+n-1];
202 u_int64_t q = x / (u_int64_t)v->numb[n-1];
203 u_int64_t r = x - (q * (u_int64_t)v->numb[n-1]);
204 u_int64_t o = 0;
205
206 // STEP D3 - Part 2 - Verify q
207 //
208 // It is possible that the estimate of q is off. This
209 // code which I'm not about to explain, alters q as
210 // necessary.
211 //
212 for (i = 0; i < 2; ++i) {
213 if (q < b && q * (u_int64_t)v->numb[n-2] <= b * r + u->numb[j+n-2])
214 break;
215
216 q -= 1;
217 r += (u_int64_t)v->numb[n-1];
218
219 if (r >= b)
220 break;
221 }
222
223 // STEP D4 - Multiply and Subtract
224 //
225 // Now that our q correct, multiply "digits" j through
226 // j + n of u by v. Then replace "digits" j through
227 // j + n of u with the product. Convoluted subtraction.
228 //
229 {
230 int i = 0;
231
232 bignum_set(bn_q, q);
233 bignum_mul(bn_tmp1, v, bn_q);
234 bignum_set(bn_tmp2, 0);
235 bignum_set(bn_tmpu, 0);
236
237 for (i = 0; i <= n; ++i)
238 __bignum_set(bn_tmpu, bn_tmpu->oofm++, u->numb[j+i]);
239
240 if (bignum_lt(bn_tmpu, bn_tmp1)) {
241 bignum_set(bn_q, --q);
242 bignum_mul(bn_tmp1, v, bn_q);
243 }
244
245 if (bignum_lt(bn_tmpu, bn_tmp1))
246 printf("bignum::bignum_div().. warning: q is still too large.\n");
247
248 o = bn_tmpu->oofm;
249 bignum_sub(bn_tmp2, bn_tmpu, bn_tmp1);
250 bignum_set_bn(bn_tmpu, bn_tmp2);
251
252 for (i = 0; i < bn_tmpu->oofm; ++i)
253 __bignum_set(u, i+j, bn_tmpu->numb[i]);
254
255 for (; i <= n; ++i)
256 __bignum_set(u, i+j, 0);
257
258 u->oofm -= o - bn_tmpu->oofm;
259 }
260
261 if (quot != NULL) {
262 __bignum_set(quot, j, q);
263
264 if (quot->numb[j] != 0)
265 quot->oofm = BIGNUM_MAX(quot->oofm, j+1);
266 }
267 }
268 while (--j >= 0 && j <= m);
269
270 bignum_free( bn_q );
271 bignum_free( bn_tmp1 );
272 bignum_free( bn_tmp2 );
273 bignum_free( bn_tmpu );
274 }
275
276 if (rmdr != NULL) {
277 bignum_t *bn_tmpd = bignum_alloc();
278 bignum_set(bn_tmpd, d);
279 bignum_div_aux(rmdr, NULL, u, bn_tmpd);
280 bignum_free( bn_tmpd );
281 }
282
283 bignum_free( u );
284 bignum_free( v );
285
286 divd->sign = divd_sign;
287 divs->sign = divs_sign;
288
289 if (quot != NULL)
290 quot->sign = quot_sign;
291
292 return quot;
293 }
syntax highlighted by Code2HTML, v. 0.9.1