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 }