1 /* 2 * Copyright (c) 2008-2016 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 "typearith.h" 36 #include "basearith.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 341 /* D3: calculate qhat and rhat */ 342 rhat = _mpd_shortdiv(w2, u+j+n-1, 2, v[n-1]); 343 qhat = w2[1] * MPD_RADIX + w2[0]; 344 345 while (1) { 346 if (qhat < MPD_RADIX) { 347 _mpd_singlemul(w2, qhat, v[n-2]); 348 if (w2[1] <= rhat) { 349 if (w2[1] != rhat || w2[0] <= u[j+n-2]) { 350 break; 351 } 352 } 353 } 354 qhat -= 1; 355 rhat += v[n-1]; 356 if (rhat < v[n-1] || rhat >= MPD_RADIX) { 357 break; 358 } 359 } 360 /* D4: multiply and subtract */ 361 carry = 0; 362 for (i=0; i <= n; i++) { 363 364 _mpd_mul_words(&hi, &lo, qhat, v[i]); 365 366 lo = carry + lo; 367 if (lo < carry) hi++; 368 369 _mpd_div_words_r(&hi, &lo, hi, lo); 370 371 x = u[i+j] - lo; 372 carry = (u[i+j] < x); 373 u[i+j] = carry ? x+MPD_RADIX : x; 374 carry += hi; 375 } 376 q[j] = qhat; 377 /* D5: test remainder */ 378 if (carry) { 379 q[j] -= 1; 380 /* D6: add back */ 381 (void)_mpd_baseadd(u+j, u+j, v, n+1, n); 382 } 383 } 384 385 /* D8: unnormalize */ 386 if (r != NULL) { 387 _mpd_shortdiv(r, u, n, d); 388 /* we are not interested in the return value here */ 389 retval = 0; 390 } 391 else { 392 retval = !_mpd_isallzero(u, n); 393 } 394 395 396 if (u != ustatic) mpd_free(u); 397 if (v != vstatic) mpd_free(v); 398 return retval; 399 } 400 401 /* 402 * Left shift of src by 'shift' digits; src may equal dest. 403 * 404 * dest := area of n mpd_uint_t with space for srcdigits+shift digits. 405 * src := coefficient with length m. 406 * 407 * The case splits in the function are non-obvious. The following 408 * equations might help: 409 * 410 * Let msdigits denote the number of digits in the most significant 411 * word of src. Then 1 <= msdigits <= rdigits. 412 * 413 * 1) shift = q * rdigits + r 414 * 2) srcdigits = qsrc * rdigits + msdigits 415 * 3) destdigits = shift + srcdigits 416 * = q * rdigits + r + qsrc * rdigits + msdigits 417 * = q * rdigits + (qsrc * rdigits + (r + msdigits)) 418 * 419 * The result has q zero words, followed by the coefficient that 420 * is left-shifted by r. The case r == 0 is trivial. For r > 0, it 421 * is important to keep in mind that we always read m source words, 422 * but write m+1 destination words if r + msdigits > rdigits, m words 423 * otherwise. 424 */ 425 void 426 _mpd_baseshiftl(mpd_uint_t *dest, mpd_uint_t *src, mpd_size_t n, mpd_size_t m, 427 mpd_size_t shift) 428 { 429 #if defined(__GNUC__) && !defined(__INTEL_COMPILER) && !defined(__clang__) 430 /* spurious uninitialized warnings */ 431 mpd_uint_t l=l, lprev=lprev, h=h; 432 #else 433 mpd_uint_t l, lprev, h; 434 #endif 435 mpd_uint_t q, r; 436 mpd_uint_t ph; 437 438 assert(m > 0 && n >= m); 439 440 _mpd_div_word(&q, &r, (mpd_uint_t)shift, MPD_RDIGITS); 441 442 if (r != 0) { 443 444 ph = mpd_pow10[r]; 445 446 --m; --n; 447 _mpd_divmod_pow10(&h, &lprev, src[m--], MPD_RDIGITS-r); 448 if (h != 0) { /* r + msdigits > rdigits <==> h != 0 */ 449 dest[n--] = h; 450 } 451 /* write m-1 shifted words */ 452 for (; m != MPD_SIZE_MAX; m--,n--) { 453 _mpd_divmod_pow10(&h, &l, src[m], MPD_RDIGITS-r); 454 dest[n] = ph * lprev + h; 455 lprev = l; 456 } 457 /* write least significant word */ 458 dest[q] = ph * lprev; 459 } 460 else { 461 while (--m != MPD_SIZE_MAX) { 462 dest[m+q] = src[m]; 463 } 464 } 465 466 mpd_uint_zero(dest, q); 467 } 468 469 /* 470 * Right shift of src by 'shift' digits; src may equal dest. 471 * Assumption: srcdigits-shift > 0. 472 * 473 * dest := area with space for srcdigits-shift digits. 474 * src := coefficient with length 'slen'. 475 * 476 * The case splits in the function rely on the following equations: 477 * 478 * Let msdigits denote the number of digits in the most significant 479 * word of src. Then 1 <= msdigits <= rdigits. 480 * 481 * 1) shift = q * rdigits + r 482 * 2) srcdigits = qsrc * rdigits + msdigits 483 * 3) destdigits = srcdigits - shift 484 * = qsrc * rdigits + msdigits - (q * rdigits + r) 485 * = (qsrc - q) * rdigits + msdigits - r 486 * 487 * Since destdigits > 0 and 1 <= msdigits <= rdigits: 488 * 489 * 4) qsrc >= q 490 * 5) qsrc == q ==> msdigits > r 491 * 492 * The result has slen-q words if msdigits > r, slen-q-1 words otherwise. 493 */ 494 mpd_uint_t 495 _mpd_baseshiftr(mpd_uint_t *dest, mpd_uint_t *src, mpd_size_t slen, 496 mpd_size_t shift) 497 { 498 #if defined(__GNUC__) && !defined(__INTEL_COMPILER) && !defined(__clang__) 499 /* spurious uninitialized warnings */ 500 mpd_uint_t l=l, h=h, hprev=hprev; /* low, high, previous high */ 501 #else 502 mpd_uint_t l, h, hprev; /* low, high, previous high */ 503 #endif 504 mpd_uint_t rnd, rest; /* rounding digit, rest */ 505 mpd_uint_t q, r; 506 mpd_size_t i, j; 507 mpd_uint_t ph; 508 509 assert(slen > 0); 510 511 _mpd_div_word(&q, &r, (mpd_uint_t)shift, MPD_RDIGITS); 512 513 rnd = rest = 0; 514 if (r != 0) { 515 516 ph = mpd_pow10[MPD_RDIGITS-r]; 517 518 _mpd_divmod_pow10(&hprev, &rest, src[q], r); 519 _mpd_divmod_pow10(&rnd, &rest, rest, r-1); 520 521 if (rest == 0 && q > 0) { 522 rest = !_mpd_isallzero(src, q); 523 } 524 /* write slen-q-1 words */ 525 for (j=0,i=q+1; i<slen; i++,j++) { 526 _mpd_divmod_pow10(&h, &l, src[i], r); 527 dest[j] = ph * l + hprev; 528 hprev = h; 529 } 530 /* write most significant word */ 531 if (hprev != 0) { /* always the case if slen==q-1 */ 532 dest[j] = hprev; 533 } 534 } 535 else { 536 if (q > 0) { 537 _mpd_divmod_pow10(&rnd, &rest, src[q-1], MPD_RDIGITS-1); 538 /* is there any non-zero digit below rnd? */ 539 if (rest == 0) rest = !_mpd_isallzero(src, q-1); 540 } 541 for (j = 0; j < slen-q; j++) { 542 dest[j] = src[q+j]; 543 } 544 } 545 546 /* 0-4 ==> rnd+rest < 0.5 */ 547 /* 5 ==> rnd+rest == 0.5 */ 548 /* 6-9 ==> rnd+rest > 0.5 */ 549 return (rnd == 0 || rnd == 5) ? rnd + !!rest : rnd; 550 } 551 552 553 /*********************************************************************/ 554 /* Calculations in base b */ 555 /*********************************************************************/ 556 557 /* 558 * Add v to w (len m). The calling function has to handle a possible 559 * final carry. Assumption: m > 0. 560 */ 561 mpd_uint_t 562 _mpd_shortadd_b(mpd_uint_t *w, mpd_size_t m, mpd_uint_t v, mpd_uint_t b) 563 { 564 mpd_uint_t s; 565 mpd_uint_t carry; 566 mpd_size_t i; 567 568 assert(m > 0); 569 570 /* add v to w */ 571 s = w[0] + v; 572 carry = (s < v) | (s >= b); 573 w[0] = carry ? s-b : s; 574 575 /* if there is a carry, propagate it */ 576 for (i = 1; carry && i < m; i++) { 577 s = w[i] + carry; 578 carry = (s == b); 579 w[i] = carry ? 0 : s; 580 } 581 582 return carry; 583 } 584 585 /* w := product of u (len n) and v (single word). Return carry. */ 586 mpd_uint_t 587 _mpd_shortmul_c(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n, mpd_uint_t v) 588 { 589 mpd_uint_t hi, lo; 590 mpd_uint_t carry = 0; 591 mpd_size_t i; 592 593 assert(n > 0); 594 595 for (i=0; i < n; i++) { 596 597 _mpd_mul_words(&hi, &lo, u[i], v); 598 lo = carry + lo; 599 if (lo < carry) hi++; 600 601 _mpd_div_words_r(&carry, &w[i], hi, lo); 602 } 603 604 return carry; 605 } 606 607 /* w := product of u (len n) and v (single word) */ 608 mpd_uint_t 609 _mpd_shortmul_b(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n, 610 mpd_uint_t v, mpd_uint_t b) 611 { 612 mpd_uint_t hi, lo; 613 mpd_uint_t carry = 0; 614 mpd_size_t i; 615 616 assert(n > 0); 617 618 for (i=0; i < n; i++) { 619 620 _mpd_mul_words(&hi, &lo, u[i], v); 621 lo = carry + lo; 622 if (lo < carry) hi++; 623 624 _mpd_div_words(&carry, &w[i], hi, lo, b); 625 } 626 627 return carry; 628 } 629 630 /* 631 * Knuth, TAOCP Volume 2, 4.3.1, exercise 16: 632 * w := quotient of u (len n) divided by a single word v 633 */ 634 mpd_uint_t 635 _mpd_shortdiv_b(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n, 636 mpd_uint_t v, mpd_uint_t b) 637 { 638 mpd_uint_t hi, lo; 639 mpd_uint_t rem = 0; 640 mpd_size_t i; 641 642 assert(n > 0); 643 644 for (i=n-1; i != MPD_SIZE_MAX; i--) { 645 646 _mpd_mul_words(&hi, &lo, rem, b); 647 lo = u[i] + lo; 648 if (lo < u[i]) hi++; 649 650 _mpd_div_words(&w[i], &rem, hi, lo, v); 651 } 652 653 return rem; 654 } 655 656 657 658