Line data Source code
1 : /*
2 : * Copyright (c) 2008-2012 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 : #include <stdio.h>
31 : #include <stdlib.h>
32 : #include <string.h>
33 : #include <assert.h>
34 : #include "constants.h"
35 : #include "memory.h"
36 : #include "typearith.h"
37 : #include "basearith.h"
38 :
39 :
40 : /*********************************************************************/
41 : /* Calculations in base MPD_RADIX */
42 : /*********************************************************************/
43 :
44 :
45 : /*
46 : * Knuth, TAOCP, Volume 2, 4.3.1:
47 : * w := sum of u (len m) and v (len n)
48 : * n > 0 and m >= n
49 : * The calling function has to handle a possible final carry.
50 : */
51 : mpd_uint_t
52 0 : _mpd_baseadd(mpd_uint_t *w, const mpd_uint_t *u, const mpd_uint_t *v,
53 : mpd_size_t m, mpd_size_t n)
54 : {
55 : mpd_uint_t s;
56 0 : mpd_uint_t carry = 0;
57 : mpd_size_t i;
58 :
59 : assert(n > 0 && m >= n);
60 :
61 : /* add n members of u and v */
62 0 : for (i = 0; i < n; i++) {
63 0 : s = u[i] + (v[i] + carry);
64 0 : carry = (s < u[i]) | (s >= MPD_RADIX);
65 0 : w[i] = carry ? s-MPD_RADIX : s;
66 : }
67 : /* if there is a carry, propagate it */
68 0 : for (; carry && i < m; i++) {
69 0 : s = u[i] + carry;
70 0 : carry = (s == MPD_RADIX);
71 0 : w[i] = carry ? 0 : s;
72 : }
73 : /* copy the rest of u */
74 0 : for (; i < m; i++) {
75 0 : w[i] = u[i];
76 : }
77 :
78 0 : return carry;
79 : }
80 :
81 : /*
82 : * Add the contents of u to w. Carries are propagated further. The caller
83 : * has to make sure that w is big enough.
84 : */
85 : void
86 0 : _mpd_baseaddto(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n)
87 : {
88 : mpd_uint_t s;
89 0 : mpd_uint_t carry = 0;
90 : mpd_size_t i;
91 :
92 0 : if (n == 0) return;
93 :
94 : /* add n members of u to w */
95 0 : for (i = 0; i < n; i++) {
96 0 : s = w[i] + (u[i] + carry);
97 0 : carry = (s < w[i]) | (s >= MPD_RADIX);
98 0 : w[i] = carry ? s-MPD_RADIX : s;
99 : }
100 : /* if there is a carry, propagate it */
101 0 : for (; carry; i++) {
102 0 : s = w[i] + carry;
103 0 : carry = (s == MPD_RADIX);
104 0 : w[i] = carry ? 0 : s;
105 : }
106 : }
107 :
108 : /*
109 : * Add v to w (len m). The calling function has to handle a possible
110 : * final carry. Assumption: m > 0.
111 : */
112 : mpd_uint_t
113 0 : _mpd_shortadd(mpd_uint_t *w, mpd_size_t m, mpd_uint_t v)
114 : {
115 : mpd_uint_t s;
116 : mpd_uint_t carry;
117 : mpd_size_t i;
118 :
119 : assert(m > 0);
120 :
121 : /* add v to w */
122 0 : s = w[0] + v;
123 0 : carry = (s < v) | (s >= MPD_RADIX);
124 0 : w[0] = carry ? s-MPD_RADIX : s;
125 :
126 : /* if there is a carry, propagate it */
127 0 : for (i = 1; carry && i < m; i++) {
128 0 : s = w[i] + carry;
129 0 : carry = (s == MPD_RADIX);
130 0 : w[i] = carry ? 0 : s;
131 : }
132 :
133 0 : return carry;
134 : }
135 :
136 : /* Increment u. The calling function has to handle a possible carry. */
137 : mpd_uint_t
138 0 : _mpd_baseincr(mpd_uint_t *u, mpd_size_t n)
139 : {
140 : mpd_uint_t s;
141 0 : mpd_uint_t carry = 1;
142 : mpd_size_t i;
143 :
144 : assert(n > 0);
145 :
146 : /* if there is a carry, propagate it */
147 0 : for (i = 0; carry && i < n; i++) {
148 0 : s = u[i] + carry;
149 0 : carry = (s == MPD_RADIX);
150 0 : u[i] = carry ? 0 : s;
151 : }
152 :
153 0 : return carry;
154 : }
155 :
156 : /*
157 : * Knuth, TAOCP, Volume 2, 4.3.1:
158 : * w := difference of u (len m) and v (len n).
159 : * number in u >= number in v;
160 : */
161 : void
162 0 : _mpd_basesub(mpd_uint_t *w, const mpd_uint_t *u, const mpd_uint_t *v,
163 : mpd_size_t m, mpd_size_t n)
164 : {
165 : mpd_uint_t d;
166 0 : mpd_uint_t borrow = 0;
167 : mpd_size_t i;
168 :
169 : assert(m > 0 && n > 0);
170 :
171 : /* subtract n members of v from u */
172 0 : for (i = 0; i < n; i++) {
173 0 : d = u[i] - (v[i] + borrow);
174 0 : borrow = (u[i] < d);
175 0 : w[i] = borrow ? d + MPD_RADIX : d;
176 : }
177 : /* if there is a borrow, propagate it */
178 0 : for (; borrow && i < m; i++) {
179 0 : d = u[i] - borrow;
180 0 : borrow = (u[i] == 0);
181 0 : w[i] = borrow ? MPD_RADIX-1 : d;
182 : }
183 : /* copy the rest of u */
184 0 : for (; i < m; i++) {
185 0 : w[i] = u[i];
186 : }
187 0 : }
188 :
189 : /*
190 : * Subtract the contents of u from w. w is larger than u. Borrows are
191 : * propagated further, but eventually w can absorb the final borrow.
192 : */
193 : void
194 0 : _mpd_basesubfrom(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n)
195 : {
196 : mpd_uint_t d;
197 0 : mpd_uint_t borrow = 0;
198 : mpd_size_t i;
199 :
200 0 : if (n == 0) return;
201 :
202 : /* subtract n members of u from w */
203 0 : for (i = 0; i < n; i++) {
204 0 : d = w[i] - (u[i] + borrow);
205 0 : borrow = (w[i] < d);
206 0 : w[i] = borrow ? d + MPD_RADIX : d;
207 : }
208 : /* if there is a borrow, propagate it */
209 0 : for (; borrow; i++) {
210 0 : d = w[i] - borrow;
211 0 : borrow = (w[i] == 0);
212 0 : w[i] = borrow ? MPD_RADIX-1 : d;
213 : }
214 : }
215 :
216 : /* w := product of u (len n) and v (single word) */
217 : void
218 0 : _mpd_shortmul(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n, mpd_uint_t v)
219 : {
220 : mpd_uint_t hi, lo;
221 0 : mpd_uint_t carry = 0;
222 : mpd_size_t i;
223 :
224 : assert(n > 0);
225 :
226 0 : for (i=0; i < n; i++) {
227 :
228 0 : _mpd_mul_words(&hi, &lo, u[i], v);
229 0 : lo = carry + lo;
230 0 : if (lo < carry) hi++;
231 :
232 0 : _mpd_div_words_r(&carry, &w[i], hi, lo);
233 : }
234 0 : w[i] = carry;
235 0 : }
236 :
237 : /*
238 : * Knuth, TAOCP, Volume 2, 4.3.1:
239 : * w := product of u (len m) and v (len n)
240 : * w must be initialized to zero
241 : */
242 : void
243 0 : _mpd_basemul(mpd_uint_t *w, const mpd_uint_t *u, const mpd_uint_t *v,
244 : mpd_size_t m, mpd_size_t n)
245 : {
246 : mpd_uint_t hi, lo;
247 : mpd_uint_t carry;
248 : mpd_size_t i, j;
249 :
250 : assert(m > 0 && n > 0);
251 :
252 0 : for (j=0; j < n; j++) {
253 0 : carry = 0;
254 0 : for (i=0; i < m; i++) {
255 :
256 0 : _mpd_mul_words(&hi, &lo, u[i], v[j]);
257 0 : lo = w[i+j] + lo;
258 0 : if (lo < w[i+j]) hi++;
259 0 : lo = carry + lo;
260 0 : if (lo < carry) hi++;
261 :
262 0 : _mpd_div_words_r(&carry, &w[i+j], hi, lo);
263 : }
264 0 : w[j+m] = carry;
265 : }
266 0 : }
267 :
268 : /*
269 : * Knuth, TAOCP Volume 2, 4.3.1, exercise 16:
270 : * w := quotient of u (len n) divided by a single word v
271 : */
272 : mpd_uint_t
273 0 : _mpd_shortdiv(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n, mpd_uint_t v)
274 : {
275 : mpd_uint_t hi, lo;
276 0 : mpd_uint_t rem = 0;
277 : mpd_size_t i;
278 :
279 : assert(n > 0);
280 :
281 0 : for (i=n-1; i != MPD_SIZE_MAX; i--) {
282 :
283 0 : _mpd_mul_words(&hi, &lo, rem, MPD_RADIX);
284 0 : lo = u[i] + lo;
285 0 : if (lo < u[i]) hi++;
286 :
287 0 : _mpd_div_words(&w[i], &rem, hi, lo, v);
288 : }
289 :
290 0 : return rem;
291 : }
292 :
293 : /*
294 : * Knuth, TAOCP Volume 2, 4.3.1:
295 : * q, r := quotient and remainder of uconst (len nplusm)
296 : * divided by vconst (len n)
297 : * nplusm >= n
298 : *
299 : * If r is not NULL, r will contain the remainder. If r is NULL, the
300 : * return value indicates if there is a remainder: 1 for true, 0 for
301 : * false. A return value of -1 indicates an error.
302 : */
303 : int
304 0 : _mpd_basedivmod(mpd_uint_t *q, mpd_uint_t *r,
305 : const mpd_uint_t *uconst, const mpd_uint_t *vconst,
306 : mpd_size_t nplusm, mpd_size_t n)
307 : {
308 : mpd_uint_t ustatic[MPD_MINALLOC_MAX];
309 : mpd_uint_t vstatic[MPD_MINALLOC_MAX];
310 0 : mpd_uint_t *u = ustatic;
311 0 : mpd_uint_t *v = vstatic;
312 : mpd_uint_t d, qhat, rhat, w2[2];
313 : mpd_uint_t hi, lo, x;
314 : mpd_uint_t carry;
315 : mpd_size_t i, j, m;
316 0 : int retval = 0;
317 :
318 : assert(n > 1 && nplusm >= n);
319 0 : m = sub_size_t(nplusm, n);
320 :
321 : /* D1: normalize */
322 0 : d = MPD_RADIX / (vconst[n-1] + 1);
323 :
324 0 : if (nplusm >= MPD_MINALLOC_MAX) {
325 0 : if ((u = mpd_alloc(nplusm+1, sizeof *u)) == NULL) {
326 0 : return -1;
327 : }
328 : }
329 0 : if (n >= MPD_MINALLOC_MAX) {
330 0 : if ((v = mpd_alloc(n+1, sizeof *v)) == NULL) {
331 0 : mpd_free(u);
332 0 : return -1;
333 : }
334 : }
335 :
336 0 : _mpd_shortmul(u, uconst, nplusm, d);
337 0 : _mpd_shortmul(v, vconst, n, d);
338 :
339 : /* D2: loop */
340 0 : for (j=m; j != MPD_SIZE_MAX; j--) {
341 :
342 : /* D3: calculate qhat and rhat */
343 0 : rhat = _mpd_shortdiv(w2, u+j+n-1, 2, v[n-1]);
344 0 : qhat = w2[1] * MPD_RADIX + w2[0];
345 :
346 : while (1) {
347 0 : if (qhat < MPD_RADIX) {
348 0 : _mpd_singlemul(w2, qhat, v[n-2]);
349 0 : if (w2[1] <= rhat) {
350 0 : if (w2[1] != rhat || w2[0] <= u[j+n-2]) {
351 : break;
352 : }
353 : }
354 : }
355 0 : qhat -= 1;
356 0 : rhat += v[n-1];
357 0 : if (rhat < v[n-1] || rhat >= MPD_RADIX) {
358 : break;
359 : }
360 0 : }
361 : /* D4: multiply and subtract */
362 0 : carry = 0;
363 0 : for (i=0; i <= n; i++) {
364 :
365 0 : _mpd_mul_words(&hi, &lo, qhat, v[i]);
366 :
367 0 : lo = carry + lo;
368 0 : if (lo < carry) hi++;
369 :
370 0 : _mpd_div_words_r(&hi, &lo, hi, lo);
371 :
372 0 : x = u[i+j] - lo;
373 0 : carry = (u[i+j] < x);
374 0 : u[i+j] = carry ? x+MPD_RADIX : x;
375 0 : carry += hi;
376 : }
377 0 : q[j] = qhat;
378 : /* D5: test remainder */
379 0 : if (carry) {
380 0 : q[j] -= 1;
381 : /* D6: add back */
382 0 : (void)_mpd_baseadd(u+j, u+j, v, n+1, n);
383 : }
384 : }
385 :
386 : /* D8: unnormalize */
387 0 : if (r != NULL) {
388 0 : _mpd_shortdiv(r, u, n, d);
389 : /* we are not interested in the return value here */
390 0 : retval = 0;
391 : }
392 : else {
393 0 : retval = !_mpd_isallzero(u, n);
394 : }
395 :
396 :
397 0 : if (u != ustatic) mpd_free(u);
398 0 : if (v != vstatic) mpd_free(v);
399 0 : 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 0 : _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 0 : 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 0 : _mpd_div_word(&q, &r, (mpd_uint_t)shift, MPD_RDIGITS);
442 :
443 0 : if (r != 0) {
444 :
445 0 : ph = mpd_pow10[r];
446 :
447 0 : --m; --n;
448 0 : _mpd_divmod_pow10(&h, &lprev, src[m--], MPD_RDIGITS-r);
449 0 : if (h != 0) { /* r + msdigits > rdigits <==> h != 0 */
450 0 : dest[n--] = h;
451 : }
452 : /* write m-1 shifted words */
453 0 : for (; m != MPD_SIZE_MAX; m--,n--) {
454 0 : _mpd_divmod_pow10(&h, &l, src[m], MPD_RDIGITS-r);
455 0 : dest[n] = ph * lprev + h;
456 0 : lprev = l;
457 : }
458 : /* write least significant word */
459 0 : dest[q] = ph * lprev;
460 : }
461 : else {
462 0 : while (--m != MPD_SIZE_MAX) {
463 0 : dest[m+q] = src[m];
464 : }
465 : }
466 :
467 0 : mpd_uint_zero(dest, q);
468 0 : }
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 0 : _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 0 : 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 0 : _mpd_div_word(&q, &r, (mpd_uint_t)shift, MPD_RDIGITS);
513 :
514 0 : rnd = rest = 0;
515 0 : if (r != 0) {
516 :
517 0 : ph = mpd_pow10[MPD_RDIGITS-r];
518 :
519 0 : _mpd_divmod_pow10(&hprev, &rest, src[q], r);
520 0 : _mpd_divmod_pow10(&rnd, &rest, rest, r-1);
521 :
522 0 : if (rest == 0 && q > 0) {
523 0 : rest = !_mpd_isallzero(src, q);
524 : }
525 : /* write slen-q-1 words */
526 0 : for (j=0,i=q+1; i<slen; i++,j++) {
527 0 : _mpd_divmod_pow10(&h, &l, src[i], r);
528 0 : dest[j] = ph * l + hprev;
529 0 : hprev = h;
530 : }
531 : /* write most significant word */
532 0 : if (hprev != 0) { /* always the case if slen==q-1 */
533 0 : dest[j] = hprev;
534 : }
535 : }
536 : else {
537 0 : if (q > 0) {
538 0 : _mpd_divmod_pow10(&rnd, &rest, src[q-1], MPD_RDIGITS-1);
539 : /* is there any non-zero digit below rnd? */
540 0 : if (rest == 0) rest = !_mpd_isallzero(src, q-1);
541 : }
542 0 : for (j = 0; j < slen-q; j++) {
543 0 : 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 0 : 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 0 : _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 0 : s = w[0] + v;
573 0 : carry = (s < v) | (s >= b);
574 0 : w[0] = carry ? s-b : s;
575 :
576 : /* if there is a carry, propagate it */
577 0 : for (i = 1; carry && i < m; i++) {
578 0 : s = w[i] + carry;
579 0 : carry = (s == b);
580 0 : w[i] = carry ? 0 : s;
581 : }
582 :
583 0 : return carry;
584 : }
585 :
586 : /* w := product of u (len n) and v (single word). Return carry. */
587 : mpd_uint_t
588 0 : _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 0 : mpd_uint_t carry = 0;
592 : mpd_size_t i;
593 :
594 : assert(n > 0);
595 :
596 0 : for (i=0; i < n; i++) {
597 :
598 0 : _mpd_mul_words(&hi, &lo, u[i], v);
599 0 : lo = carry + lo;
600 0 : if (lo < carry) hi++;
601 :
602 0 : _mpd_div_words_r(&carry, &w[i], hi, lo);
603 : }
604 :
605 0 : return carry;
606 : }
607 :
608 : /* w := product of u (len n) and v (single word) */
609 : mpd_uint_t
610 0 : _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 0 : mpd_uint_t carry = 0;
615 : mpd_size_t i;
616 :
617 : assert(n > 0);
618 :
619 0 : for (i=0; i < n; i++) {
620 :
621 0 : _mpd_mul_words(&hi, &lo, u[i], v);
622 0 : lo = carry + lo;
623 0 : if (lo < carry) hi++;
624 :
625 0 : _mpd_div_words(&carry, &w[i], hi, lo, b);
626 : }
627 :
628 0 : 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 0 : _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 0 : mpd_uint_t rem = 0;
641 : mpd_size_t i;
642 :
643 : assert(n > 0);
644 :
645 0 : for (i=n-1; i != MPD_SIZE_MAX; i--) {
646 :
647 0 : _mpd_mul_words(&hi, &lo, rem, b);
648 0 : lo = u[i] + lo;
649 0 : if (lo < u[i]) hi++;
650 :
651 0 : _mpd_div_words(&w[i], &rem, hi, lo, v);
652 : }
653 :
654 0 : return rem;
655 : }
656 :
657 :
658 :
|
