1 | /* |
2 | * Copyright (c) 2008-2020 Stefan Krah. All rights reserved. |
3 | * |
4 | * Redistribution and use in source and binary forms, with or without |
5 | * modification, are permitted provided that the following conditions |
6 | * are met: |
7 | * |
8 | * 1. Redistributions of source code must retain the above copyright |
9 | * notice, this list of conditions and the following disclaimer. |
10 | * |
11 | * 2. Redistributions in binary form must reproduce the above copyright |
12 | * notice, this list of conditions and the following disclaimer in the |
13 | * documentation and/or other materials provided with the distribution. |
14 | * |
15 | * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS "AS IS" AND |
16 | * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE |
17 | * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE |
18 | * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE |
19 | * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL |
20 | * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS |
21 | * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) |
22 | * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT |
23 | * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY |
24 | * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF |
25 | * SUCH DAMAGE. |
26 | */ |
27 | |
28 | |
29 | #include "mpdecimal.h" |
30 | |
31 | #include <assert.h> |
32 | #include <stdio.h> |
33 | |
34 | #include "basearith.h" |
35 | #include "constants.h" |
36 | #include "typearith.h" |
37 | |
38 | |
39 | /*********************************************************************/ |
40 | /* Calculations in base MPD_RADIX */ |
41 | /*********************************************************************/ |
42 | |
43 | |
44 | /* |
45 | * Knuth, TAOCP, Volume 2, 4.3.1: |
46 | * w := sum of u (len m) and v (len n) |
47 | * n > 0 and m >= n |
48 | * The calling function has to handle a possible final carry. |
49 | */ |
50 | mpd_uint_t |
51 | _mpd_baseadd(mpd_uint_t *w, const mpd_uint_t *u, const mpd_uint_t *v, |
52 | mpd_size_t m, mpd_size_t n) |
53 | { |
54 | mpd_uint_t s; |
55 | mpd_uint_t carry = 0; |
56 | mpd_size_t i; |
57 | |
58 | assert(n > 0 && m >= n); |
59 | |
60 | /* add n members of u and v */ |
61 | for (i = 0; i < n; i++) { |
62 | s = u[i] + (v[i] + carry); |
63 | carry = (s < u[i]) | (s >= MPD_RADIX); |
64 | w[i] = carry ? s-MPD_RADIX : s; |
65 | } |
66 | /* if there is a carry, propagate it */ |
67 | for (; carry && i < m; i++) { |
68 | s = u[i] + carry; |
69 | carry = (s == MPD_RADIX); |
70 | w[i] = carry ? 0 : s; |
71 | } |
72 | /* copy the rest of u */ |
73 | for (; i < m; i++) { |
74 | w[i] = u[i]; |
75 | } |
76 | |
77 | return carry; |
78 | } |
79 | |
80 | /* |
81 | * Add the contents of u to w. Carries are propagated further. The caller |
82 | * has to make sure that w is big enough. |
83 | */ |
84 | void |
85 | _mpd_baseaddto(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n) |
86 | { |
87 | mpd_uint_t s; |
88 | mpd_uint_t carry = 0; |
89 | mpd_size_t i; |
90 | |
91 | if (n == 0) return; |
92 | |
93 | /* add n members of u to w */ |
94 | for (i = 0; i < n; i++) { |
95 | s = w[i] + (u[i] + carry); |
96 | carry = (s < w[i]) | (s >= MPD_RADIX); |
97 | w[i] = carry ? s-MPD_RADIX : s; |
98 | } |
99 | /* if there is a carry, propagate it */ |
100 | for (; carry; i++) { |
101 | s = w[i] + carry; |
102 | carry = (s == MPD_RADIX); |
103 | w[i] = carry ? 0 : s; |
104 | } |
105 | } |
106 | |
107 | /* |
108 | * Add v to w (len m). The calling function has to handle a possible |
109 | * final carry. Assumption: m > 0. |
110 | */ |
111 | mpd_uint_t |
112 | _mpd_shortadd(mpd_uint_t *w, mpd_size_t m, mpd_uint_t v) |
113 | { |
114 | mpd_uint_t s; |
115 | mpd_uint_t carry; |
116 | mpd_size_t i; |
117 | |
118 | assert(m > 0); |
119 | |
120 | /* add v to w */ |
121 | s = w[0] + v; |
122 | carry = (s < v) | (s >= MPD_RADIX); |
123 | w[0] = carry ? s-MPD_RADIX : s; |
124 | |
125 | /* if there is a carry, propagate it */ |
126 | for (i = 1; carry && i < m; i++) { |
127 | s = w[i] + carry; |
128 | carry = (s == MPD_RADIX); |
129 | w[i] = carry ? 0 : s; |
130 | } |
131 | |
132 | return carry; |
133 | } |
134 | |
135 | /* Increment u. The calling function has to handle a possible carry. */ |
136 | mpd_uint_t |
137 | _mpd_baseincr(mpd_uint_t *u, mpd_size_t n) |
138 | { |
139 | mpd_uint_t s; |
140 | mpd_uint_t carry = 1; |
141 | mpd_size_t i; |
142 | |
143 | assert(n > 0); |
144 | |
145 | /* if there is a carry, propagate it */ |
146 | for (i = 0; carry && i < n; i++) { |
147 | s = u[i] + carry; |
148 | carry = (s == MPD_RADIX); |
149 | u[i] = carry ? 0 : s; |
150 | } |
151 | |
152 | return carry; |
153 | } |
154 | |
155 | /* |
156 | * Knuth, TAOCP, Volume 2, 4.3.1: |
157 | * w := difference of u (len m) and v (len n). |
158 | * number in u >= number in v; |
159 | */ |
160 | void |
161 | _mpd_basesub(mpd_uint_t *w, const mpd_uint_t *u, const mpd_uint_t *v, |
162 | mpd_size_t m, mpd_size_t n) |
163 | { |
164 | mpd_uint_t d; |
165 | mpd_uint_t borrow = 0; |
166 | mpd_size_t i; |
167 | |
168 | assert(m > 0 && n > 0); |
169 | |
170 | /* subtract n members of v from u */ |
171 | for (i = 0; i < n; i++) { |
172 | d = u[i] - (v[i] + borrow); |
173 | borrow = (u[i] < d); |
174 | w[i] = borrow ? d + MPD_RADIX : d; |
175 | } |
176 | /* if there is a borrow, propagate it */ |
177 | for (; borrow && i < m; i++) { |
178 | d = u[i] - borrow; |
179 | borrow = (u[i] == 0); |
180 | w[i] = borrow ? MPD_RADIX-1 : d; |
181 | } |
182 | /* copy the rest of u */ |
183 | for (; i < m; i++) { |
184 | w[i] = u[i]; |
185 | } |
186 | } |
187 | |
188 | /* |
189 | * Subtract the contents of u from w. w is larger than u. Borrows are |
190 | * propagated further, but eventually w can absorb the final borrow. |
191 | */ |
192 | void |
193 | _mpd_basesubfrom(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n) |
194 | { |
195 | mpd_uint_t d; |
196 | mpd_uint_t borrow = 0; |
197 | mpd_size_t i; |
198 | |
199 | if (n == 0) return; |
200 | |
201 | /* subtract n members of u from w */ |
202 | for (i = 0; i < n; i++) { |
203 | d = w[i] - (u[i] + borrow); |
204 | borrow = (w[i] < d); |
205 | w[i] = borrow ? d + MPD_RADIX : d; |
206 | } |
207 | /* if there is a borrow, propagate it */ |
208 | for (; borrow; i++) { |
209 | d = w[i] - borrow; |
210 | borrow = (w[i] == 0); |
211 | w[i] = borrow ? MPD_RADIX-1 : d; |
212 | } |
213 | } |
214 | |
215 | /* w := product of u (len n) and v (single word) */ |
216 | void |
217 | _mpd_shortmul(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n, mpd_uint_t v) |
218 | { |
219 | mpd_uint_t hi, lo; |
220 | mpd_uint_t carry = 0; |
221 | mpd_size_t i; |
222 | |
223 | assert(n > 0); |
224 | |
225 | for (i=0; i < n; i++) { |
226 | |
227 | _mpd_mul_words(&hi, &lo, u[i], v); |
228 | lo = carry + lo; |
229 | if (lo < carry) hi++; |
230 | |
231 | _mpd_div_words_r(&carry, &w[i], hi, lo); |
232 | } |
233 | w[i] = carry; |
234 | } |
235 | |
236 | /* |
237 | * Knuth, TAOCP, Volume 2, 4.3.1: |
238 | * w := product of u (len m) and v (len n) |
239 | * w must be initialized to zero |
240 | */ |
241 | void |
242 | _mpd_basemul(mpd_uint_t *w, const mpd_uint_t *u, const mpd_uint_t *v, |
243 | mpd_size_t m, mpd_size_t n) |
244 | { |
245 | mpd_uint_t hi, lo; |
246 | mpd_uint_t carry; |
247 | mpd_size_t i, j; |
248 | |
249 | assert(m > 0 && n > 0); |
250 | |
251 | for (j=0; j < n; j++) { |
252 | carry = 0; |
253 | for (i=0; i < m; i++) { |
254 | |
255 | _mpd_mul_words(&hi, &lo, u[i], v[j]); |
256 | lo = w[i+j] + lo; |
257 | if (lo < w[i+j]) hi++; |
258 | lo = carry + lo; |
259 | if (lo < carry) hi++; |
260 | |
261 | _mpd_div_words_r(&carry, &w[i+j], hi, lo); |
262 | } |
263 | w[j+m] = carry; |
264 | } |
265 | } |
266 | |
267 | /* |
268 | * Knuth, TAOCP Volume 2, 4.3.1, exercise 16: |
269 | * w := quotient of u (len n) divided by a single word v |
270 | */ |
271 | mpd_uint_t |
272 | _mpd_shortdiv(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n, mpd_uint_t v) |
273 | { |
274 | mpd_uint_t hi, lo; |
275 | mpd_uint_t rem = 0; |
276 | mpd_size_t i; |
277 | |
278 | assert(n > 0); |
279 | |
280 | for (i=n-1; i != MPD_SIZE_MAX; i--) { |
281 | |
282 | _mpd_mul_words(&hi, &lo, rem, MPD_RADIX); |
283 | lo = u[i] + lo; |
284 | if (lo < u[i]) hi++; |
285 | |
286 | _mpd_div_words(&w[i], &rem, hi, lo, v); |
287 | } |
288 | |
289 | return rem; |
290 | } |
291 | |
292 | /* |
293 | * Knuth, TAOCP Volume 2, 4.3.1: |
294 | * q, r := quotient and remainder of uconst (len nplusm) |
295 | * divided by vconst (len n) |
296 | * nplusm >= n |
297 | * |
298 | * If r is not NULL, r will contain the remainder. If r is NULL, the |
299 | * return value indicates if there is a remainder: 1 for true, 0 for |
300 | * false. A return value of -1 indicates an error. |
301 | */ |
302 | int |
303 | _mpd_basedivmod(mpd_uint_t *q, mpd_uint_t *r, |
304 | const mpd_uint_t *uconst, const mpd_uint_t *vconst, |
305 | mpd_size_t nplusm, mpd_size_t n) |
306 | { |
307 | mpd_uint_t ustatic[MPD_MINALLOC_MAX]; |
308 | mpd_uint_t vstatic[MPD_MINALLOC_MAX]; |
309 | mpd_uint_t *u = ustatic; |
310 | mpd_uint_t *v = vstatic; |
311 | mpd_uint_t d, qhat, rhat, w2[2]; |
312 | mpd_uint_t hi, lo, x; |
313 | mpd_uint_t carry; |
314 | mpd_size_t i, j, m; |
315 | int retval = 0; |
316 | |
317 | assert(n > 1 && nplusm >= n); |
318 | m = sub_size_t(nplusm, n); |
319 | |
320 | /* D1: normalize */ |
321 | d = MPD_RADIX / (vconst[n-1] + 1); |
322 | |
323 | if (nplusm >= MPD_MINALLOC_MAX) { |
324 | if ((u = mpd_alloc(nplusm+1, sizeof *u)) == NULL) { |
325 | return -1; |
326 | } |
327 | } |
328 | if (n >= MPD_MINALLOC_MAX) { |
329 | if ((v = mpd_alloc(n+1, sizeof *v)) == NULL) { |
330 | mpd_free(u); |
331 | return -1; |
332 | } |
333 | } |
334 | |
335 | _mpd_shortmul(u, uconst, nplusm, d); |
336 | _mpd_shortmul(v, vconst, n, d); |
337 | |
338 | /* D2: loop */ |
339 | for (j=m; j != MPD_SIZE_MAX; j--) { |
340 | assert(2 <= j+n && j+n <= nplusm); /* annotation for scan-build */ |
341 | |
342 | /* D3: calculate qhat and rhat */ |
343 | rhat = _mpd_shortdiv(w2, u+j+n-1, 2, v[n-1]); |
344 | qhat = w2[1] * MPD_RADIX + w2[0]; |
345 | |
346 | while (1) { |
347 | if (qhat < MPD_RADIX) { |
348 | _mpd_singlemul(w2, qhat, v[n-2]); |
349 | if (w2[1] <= rhat) { |
350 | if (w2[1] != rhat || w2[0] <= u[j+n-2]) { |
351 | break; |
352 | } |
353 | } |
354 | } |
355 | qhat -= 1; |
356 | rhat += v[n-1]; |
357 | if (rhat < v[n-1] || rhat >= MPD_RADIX) { |
358 | break; |
359 | } |
360 | } |
361 | /* D4: multiply and subtract */ |
362 | carry = 0; |
363 | for (i=0; i <= n; i++) { |
364 | |
365 | _mpd_mul_words(&hi, &lo, qhat, v[i]); |
366 | |
367 | lo = carry + lo; |
368 | if (lo < carry) hi++; |
369 | |
370 | _mpd_div_words_r(&hi, &lo, hi, lo); |
371 | |
372 | x = u[i+j] - lo; |
373 | carry = (u[i+j] < x); |
374 | u[i+j] = carry ? x+MPD_RADIX : x; |
375 | carry += hi; |
376 | } |
377 | q[j] = qhat; |
378 | /* D5: test remainder */ |
379 | if (carry) { |
380 | q[j] -= 1; |
381 | /* D6: add back */ |
382 | (void)_mpd_baseadd(u+j, u+j, v, n+1, n); |
383 | } |
384 | } |
385 | |
386 | /* D8: unnormalize */ |
387 | if (r != NULL) { |
388 | _mpd_shortdiv(r, u, n, d); |
389 | /* we are not interested in the return value here */ |
390 | retval = 0; |
391 | } |
392 | else { |
393 | retval = !_mpd_isallzero(u, n); |
394 | } |
395 | |
396 | |
397 | if (u != ustatic) mpd_free(u); |
398 | if (v != vstatic) mpd_free(v); |
399 | return retval; |
400 | } |
401 | |
402 | /* |
403 | * Left shift of src by 'shift' digits; src may equal dest. |
404 | * |
405 | * dest := area of n mpd_uint_t with space for srcdigits+shift digits. |
406 | * src := coefficient with length m. |
407 | * |
408 | * The case splits in the function are non-obvious. The following |
409 | * equations might help: |
410 | * |
411 | * Let msdigits denote the number of digits in the most significant |
412 | * word of src. Then 1 <= msdigits <= rdigits. |
413 | * |
414 | * 1) shift = q * rdigits + r |
415 | * 2) srcdigits = qsrc * rdigits + msdigits |
416 | * 3) destdigits = shift + srcdigits |
417 | * = q * rdigits + r + qsrc * rdigits + msdigits |
418 | * = q * rdigits + (qsrc * rdigits + (r + msdigits)) |
419 | * |
420 | * The result has q zero words, followed by the coefficient that |
421 | * is left-shifted by r. The case r == 0 is trivial. For r > 0, it |
422 | * is important to keep in mind that we always read m source words, |
423 | * but write m+1 destination words if r + msdigits > rdigits, m words |
424 | * otherwise. |
425 | */ |
426 | void |
427 | _mpd_baseshiftl(mpd_uint_t *dest, mpd_uint_t *src, mpd_size_t n, mpd_size_t m, |
428 | mpd_size_t shift) |
429 | { |
430 | #if defined(__GNUC__) && !defined(__INTEL_COMPILER) && !defined(__clang__) |
431 | /* spurious uninitialized warnings */ |
432 | mpd_uint_t l=l, lprev=lprev, h=h; |
433 | #else |
434 | mpd_uint_t l, lprev, h; |
435 | #endif |
436 | mpd_uint_t q, r; |
437 | mpd_uint_t ph; |
438 | |
439 | assert(m > 0 && n >= m); |
440 | |
441 | _mpd_div_word(&q, &r, (mpd_uint_t)shift, MPD_RDIGITS); |
442 | |
443 | if (r != 0) { |
444 | |
445 | ph = mpd_pow10[r]; |
446 | |
447 | --m; --n; |
448 | _mpd_divmod_pow10(&h, &lprev, src[m--], MPD_RDIGITS-r); |
449 | if (h != 0) { /* r + msdigits > rdigits <==> h != 0 */ |
450 | dest[n--] = h; |
451 | } |
452 | /* write m-1 shifted words */ |
453 | for (; m != MPD_SIZE_MAX; m--,n--) { |
454 | _mpd_divmod_pow10(&h, &l, src[m], MPD_RDIGITS-r); |
455 | dest[n] = ph * lprev + h; |
456 | lprev = l; |
457 | } |
458 | /* write least significant word */ |
459 | dest[q] = ph * lprev; |
460 | } |
461 | else { |
462 | while (--m != MPD_SIZE_MAX) { |
463 | dest[m+q] = src[m]; |
464 | } |
465 | } |
466 | |
467 | mpd_uint_zero(dest, q); |
468 | } |
469 | |
470 | /* |
471 | * Right shift of src by 'shift' digits; src may equal dest. |
472 | * Assumption: srcdigits-shift > 0. |
473 | * |
474 | * dest := area with space for srcdigits-shift digits. |
475 | * src := coefficient with length 'slen'. |
476 | * |
477 | * The case splits in the function rely on the following equations: |
478 | * |
479 | * Let msdigits denote the number of digits in the most significant |
480 | * word of src. Then 1 <= msdigits <= rdigits. |
481 | * |
482 | * 1) shift = q * rdigits + r |
483 | * 2) srcdigits = qsrc * rdigits + msdigits |
484 | * 3) destdigits = srcdigits - shift |
485 | * = qsrc * rdigits + msdigits - (q * rdigits + r) |
486 | * = (qsrc - q) * rdigits + msdigits - r |
487 | * |
488 | * Since destdigits > 0 and 1 <= msdigits <= rdigits: |
489 | * |
490 | * 4) qsrc >= q |
491 | * 5) qsrc == q ==> msdigits > r |
492 | * |
493 | * The result has slen-q words if msdigits > r, slen-q-1 words otherwise. |
494 | */ |
495 | mpd_uint_t |
496 | _mpd_baseshiftr(mpd_uint_t *dest, mpd_uint_t *src, mpd_size_t slen, |
497 | mpd_size_t shift) |
498 | { |
499 | #if defined(__GNUC__) && !defined(__INTEL_COMPILER) && !defined(__clang__) |
500 | /* spurious uninitialized warnings */ |
501 | mpd_uint_t l=l, h=h, hprev=hprev; /* low, high, previous high */ |
502 | #else |
503 | mpd_uint_t l, h, hprev; /* low, high, previous high */ |
504 | #endif |
505 | mpd_uint_t rnd, rest; /* rounding digit, rest */ |
506 | mpd_uint_t q, r; |
507 | mpd_size_t i, j; |
508 | mpd_uint_t ph; |
509 | |
510 | assert(slen > 0); |
511 | |
512 | _mpd_div_word(&q, &r, (mpd_uint_t)shift, MPD_RDIGITS); |
513 | |
514 | rnd = rest = 0; |
515 | if (r != 0) { |
516 | |
517 | ph = mpd_pow10[MPD_RDIGITS-r]; |
518 | |
519 | _mpd_divmod_pow10(&hprev, &rest, src[q], r); |
520 | _mpd_divmod_pow10(&rnd, &rest, rest, r-1); |
521 | |
522 | if (rest == 0 && q > 0) { |
523 | rest = !_mpd_isallzero(src, q); |
524 | } |
525 | /* write slen-q-1 words */ |
526 | for (j=0,i=q+1; i<slen; i++,j++) { |
527 | _mpd_divmod_pow10(&h, &l, src[i], r); |
528 | dest[j] = ph * l + hprev; |
529 | hprev = h; |
530 | } |
531 | /* write most significant word */ |
532 | if (hprev != 0) { /* always the case if slen==q-1 */ |
533 | dest[j] = hprev; |
534 | } |
535 | } |
536 | else { |
537 | if (q > 0) { |
538 | _mpd_divmod_pow10(&rnd, &rest, src[q-1], MPD_RDIGITS-1); |
539 | /* is there any non-zero digit below rnd? */ |
540 | if (rest == 0) rest = !_mpd_isallzero(src, q-1); |
541 | } |
542 | for (j = 0; j < slen-q; j++) { |
543 | dest[j] = src[q+j]; |
544 | } |
545 | } |
546 | |
547 | /* 0-4 ==> rnd+rest < 0.5 */ |
548 | /* 5 ==> rnd+rest == 0.5 */ |
549 | /* 6-9 ==> rnd+rest > 0.5 */ |
550 | return (rnd == 0 || rnd == 5) ? rnd + !!rest : rnd; |
551 | } |
552 | |
553 | |
554 | /*********************************************************************/ |
555 | /* Calculations in base b */ |
556 | /*********************************************************************/ |
557 | |
558 | /* |
559 | * Add v to w (len m). The calling function has to handle a possible |
560 | * final carry. Assumption: m > 0. |
561 | */ |
562 | mpd_uint_t |
563 | _mpd_shortadd_b(mpd_uint_t *w, mpd_size_t m, mpd_uint_t v, mpd_uint_t b) |
564 | { |
565 | mpd_uint_t s; |
566 | mpd_uint_t carry; |
567 | mpd_size_t i; |
568 | |
569 | assert(m > 0); |
570 | |
571 | /* add v to w */ |
572 | s = w[0] + v; |
573 | carry = (s < v) | (s >= b); |
574 | w[0] = carry ? s-b : s; |
575 | |
576 | /* if there is a carry, propagate it */ |
577 | for (i = 1; carry && i < m; i++) { |
578 | s = w[i] + carry; |
579 | carry = (s == b); |
580 | w[i] = carry ? 0 : s; |
581 | } |
582 | |
583 | return carry; |
584 | } |
585 | |
586 | /* w := product of u (len n) and v (single word). Return carry. */ |
587 | mpd_uint_t |
588 | _mpd_shortmul_c(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n, mpd_uint_t v) |
589 | { |
590 | mpd_uint_t hi, lo; |
591 | mpd_uint_t carry = 0; |
592 | mpd_size_t i; |
593 | |
594 | assert(n > 0); |
595 | |
596 | for (i=0; i < n; i++) { |
597 | |
598 | _mpd_mul_words(&hi, &lo, u[i], v); |
599 | lo = carry + lo; |
600 | if (lo < carry) hi++; |
601 | |
602 | _mpd_div_words_r(&carry, &w[i], hi, lo); |
603 | } |
604 | |
605 | return carry; |
606 | } |
607 | |
608 | /* w := product of u (len n) and v (single word) */ |
609 | mpd_uint_t |
610 | _mpd_shortmul_b(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n, |
611 | mpd_uint_t v, mpd_uint_t b) |
612 | { |
613 | mpd_uint_t hi, lo; |
614 | mpd_uint_t carry = 0; |
615 | mpd_size_t i; |
616 | |
617 | assert(n > 0); |
618 | |
619 | for (i=0; i < n; i++) { |
620 | |
621 | _mpd_mul_words(&hi, &lo, u[i], v); |
622 | lo = carry + lo; |
623 | if (lo < carry) hi++; |
624 | |
625 | _mpd_div_words(&carry, &w[i], hi, lo, b); |
626 | } |
627 | |
628 | return carry; |
629 | } |
630 | |
631 | /* |
632 | * Knuth, TAOCP Volume 2, 4.3.1, exercise 16: |
633 | * w := quotient of u (len n) divided by a single word v |
634 | */ |
635 | mpd_uint_t |
636 | _mpd_shortdiv_b(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n, |
637 | mpd_uint_t v, mpd_uint_t b) |
638 | { |
639 | mpd_uint_t hi, lo; |
640 | mpd_uint_t rem = 0; |
641 | mpd_size_t i; |
642 | |
643 | assert(n > 0); |
644 | |
645 | for (i=n-1; i != MPD_SIZE_MAX; i--) { |
646 | |
647 | _mpd_mul_words(&hi, &lo, rem, b); |
648 | lo = u[i] + lo; |
649 | if (lo < u[i]) hi++; |
650 | |
651 | _mpd_div_words(&w[i], &rem, hi, lo, v); |
652 | } |
653 | |
654 | return rem; |
655 | } |
656 | |