| /* |
| * Copyright (c) 2008-2016 Stefan Krah. All rights reserved. |
| * |
| * Redistribution and use in source and binary forms, with or without |
| * modification, are permitted provided that the following conditions |
| * are met: |
| * |
| * 1. Redistributions of source code must retain the above copyright |
| * notice, this list of conditions and the following disclaimer. |
| * |
| * 2. Redistributions in binary form must reproduce the above copyright |
| * notice, this list of conditions and the following disclaimer in the |
| * documentation and/or other materials provided with the distribution. |
| * |
| * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS "AS IS" AND |
| * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE |
| * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE |
| * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE |
| * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL |
| * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS |
| * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) |
| * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT |
| * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY |
| * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF |
| * SUCH DAMAGE. |
| */ |
| |
| |
| #include mpdecimal_header |
| #include <stdio.h> |
| #include <stdlib.h> |
| #include <string.h> |
| #include <assert.h> |
| #include "constants.h" |
| #include "memory.h" |
| #include "typearith.h" |
| #include "basearith.h" |
| |
| |
| /*********************************************************************/ |
| /* Calculations in base MPD_RADIX */ |
| /*********************************************************************/ |
| |
| |
| /* |
| * Knuth, TAOCP, Volume 2, 4.3.1: |
| * w := sum of u (len m) and v (len n) |
| * n > 0 and m >= n |
| * The calling function has to handle a possible final carry. |
| */ |
| mpd_uint_t |
| _mpd_baseadd(mpd_uint_t *w, const mpd_uint_t *u, const mpd_uint_t *v, |
| mpd_size_t m, mpd_size_t n) |
| { |
| mpd_uint_t s; |
| mpd_uint_t carry = 0; |
| mpd_size_t i; |
| |
| assert(n > 0 && m >= n); |
| |
| /* add n members of u and v */ |
| for (i = 0; i < n; i++) { |
| s = u[i] + (v[i] + carry); |
| carry = (s < u[i]) | (s >= MPD_RADIX); |
| w[i] = carry ? s-MPD_RADIX : s; |
| } |
| /* if there is a carry, propagate it */ |
| for (; carry && i < m; i++) { |
| s = u[i] + carry; |
| carry = (s == MPD_RADIX); |
| w[i] = carry ? 0 : s; |
| } |
| /* copy the rest of u */ |
| for (; i < m; i++) { |
| w[i] = u[i]; |
| } |
| |
| return carry; |
| } |
| |
| /* |
| * Add the contents of u to w. Carries are propagated further. The caller |
| * has to make sure that w is big enough. |
| */ |
| void |
| _mpd_baseaddto(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n) |
| { |
| mpd_uint_t s; |
| mpd_uint_t carry = 0; |
| mpd_size_t i; |
| |
| if (n == 0) return; |
| |
| /* add n members of u to w */ |
| for (i = 0; i < n; i++) { |
| s = w[i] + (u[i] + carry); |
| carry = (s < w[i]) | (s >= MPD_RADIX); |
| w[i] = carry ? s-MPD_RADIX : s; |
| } |
| /* if there is a carry, propagate it */ |
| for (; carry; i++) { |
| s = w[i] + carry; |
| carry = (s == MPD_RADIX); |
| w[i] = carry ? 0 : s; |
| } |
| } |
| |
| /* |
| * Add v to w (len m). The calling function has to handle a possible |
| * final carry. Assumption: m > 0. |
| */ |
| mpd_uint_t |
| _mpd_shortadd(mpd_uint_t *w, mpd_size_t m, mpd_uint_t v) |
| { |
| mpd_uint_t s; |
| mpd_uint_t carry; |
| mpd_size_t i; |
| |
| assert(m > 0); |
| |
| /* add v to w */ |
| s = w[0] + v; |
| carry = (s < v) | (s >= MPD_RADIX); |
| w[0] = carry ? s-MPD_RADIX : s; |
| |
| /* if there is a carry, propagate it */ |
| for (i = 1; carry && i < m; i++) { |
| s = w[i] + carry; |
| carry = (s == MPD_RADIX); |
| w[i] = carry ? 0 : s; |
| } |
| |
| return carry; |
| } |
| |
| /* Increment u. The calling function has to handle a possible carry. */ |
| mpd_uint_t |
| _mpd_baseincr(mpd_uint_t *u, mpd_size_t n) |
| { |
| mpd_uint_t s; |
| mpd_uint_t carry = 1; |
| mpd_size_t i; |
| |
| assert(n > 0); |
| |
| /* if there is a carry, propagate it */ |
| for (i = 0; carry && i < n; i++) { |
| s = u[i] + carry; |
| carry = (s == MPD_RADIX); |
| u[i] = carry ? 0 : s; |
| } |
| |
| return carry; |
| } |
| |
| /* |
| * Knuth, TAOCP, Volume 2, 4.3.1: |
| * w := difference of u (len m) and v (len n). |
| * number in u >= number in v; |
| */ |
| void |
| _mpd_basesub(mpd_uint_t *w, const mpd_uint_t *u, const mpd_uint_t *v, |
| mpd_size_t m, mpd_size_t n) |
| { |
| mpd_uint_t d; |
| mpd_uint_t borrow = 0; |
| mpd_size_t i; |
| |
| assert(m > 0 && n > 0); |
| |
| /* subtract n members of v from u */ |
| for (i = 0; i < n; i++) { |
| d = u[i] - (v[i] + borrow); |
| borrow = (u[i] < d); |
| w[i] = borrow ? d + MPD_RADIX : d; |
| } |
| /* if there is a borrow, propagate it */ |
| for (; borrow && i < m; i++) { |
| d = u[i] - borrow; |
| borrow = (u[i] == 0); |
| w[i] = borrow ? MPD_RADIX-1 : d; |
| } |
| /* copy the rest of u */ |
| for (; i < m; i++) { |
| w[i] = u[i]; |
| } |
| } |
| |
| /* |
| * Subtract the contents of u from w. w is larger than u. Borrows are |
| * propagated further, but eventually w can absorb the final borrow. |
| */ |
| void |
| _mpd_basesubfrom(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n) |
| { |
| mpd_uint_t d; |
| mpd_uint_t borrow = 0; |
| mpd_size_t i; |
| |
| if (n == 0) return; |
| |
| /* subtract n members of u from w */ |
| for (i = 0; i < n; i++) { |
| d = w[i] - (u[i] + borrow); |
| borrow = (w[i] < d); |
| w[i] = borrow ? d + MPD_RADIX : d; |
| } |
| /* if there is a borrow, propagate it */ |
| for (; borrow; i++) { |
| d = w[i] - borrow; |
| borrow = (w[i] == 0); |
| w[i] = borrow ? MPD_RADIX-1 : d; |
| } |
| } |
| |
| /* w := product of u (len n) and v (single word) */ |
| void |
| _mpd_shortmul(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n, mpd_uint_t v) |
| { |
| mpd_uint_t hi, lo; |
| mpd_uint_t carry = 0; |
| mpd_size_t i; |
| |
| assert(n > 0); |
| |
| for (i=0; i < n; i++) { |
| |
| _mpd_mul_words(&hi, &lo, u[i], v); |
| lo = carry + lo; |
| if (lo < carry) hi++; |
| |
| _mpd_div_words_r(&carry, &w[i], hi, lo); |
| } |
| w[i] = carry; |
| } |
| |
| /* |
| * Knuth, TAOCP, Volume 2, 4.3.1: |
| * w := product of u (len m) and v (len n) |
| * w must be initialized to zero |
| */ |
| void |
| _mpd_basemul(mpd_uint_t *w, const mpd_uint_t *u, const mpd_uint_t *v, |
| mpd_size_t m, mpd_size_t n) |
| { |
| mpd_uint_t hi, lo; |
| mpd_uint_t carry; |
| mpd_size_t i, j; |
| |
| assert(m > 0 && n > 0); |
| |
| for (j=0; j < n; j++) { |
| carry = 0; |
| for (i=0; i < m; i++) { |
| |
| _mpd_mul_words(&hi, &lo, u[i], v[j]); |
| lo = w[i+j] + lo; |
| if (lo < w[i+j]) hi++; |
| lo = carry + lo; |
| if (lo < carry) hi++; |
| |
| _mpd_div_words_r(&carry, &w[i+j], hi, lo); |
| } |
| w[j+m] = carry; |
| } |
| } |
| |
| /* |
| * Knuth, TAOCP Volume 2, 4.3.1, exercise 16: |
| * w := quotient of u (len n) divided by a single word v |
| */ |
| mpd_uint_t |
| _mpd_shortdiv(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n, mpd_uint_t v) |
| { |
| mpd_uint_t hi, lo; |
| mpd_uint_t rem = 0; |
| mpd_size_t i; |
| |
| assert(n > 0); |
| |
| for (i=n-1; i != MPD_SIZE_MAX; i--) { |
| |
| _mpd_mul_words(&hi, &lo, rem, MPD_RADIX); |
| lo = u[i] + lo; |
| if (lo < u[i]) hi++; |
| |
| _mpd_div_words(&w[i], &rem, hi, lo, v); |
| } |
| |
| return rem; |
| } |
| |
| /* |
| * Knuth, TAOCP Volume 2, 4.3.1: |
| * q, r := quotient and remainder of uconst (len nplusm) |
| * divided by vconst (len n) |
| * nplusm >= n |
| * |
| * If r is not NULL, r will contain the remainder. If r is NULL, the |
| * return value indicates if there is a remainder: 1 for true, 0 for |
| * false. A return value of -1 indicates an error. |
| */ |
| int |
| _mpd_basedivmod(mpd_uint_t *q, mpd_uint_t *r, |
| const mpd_uint_t *uconst, const mpd_uint_t *vconst, |
| mpd_size_t nplusm, mpd_size_t n) |
| { |
| mpd_uint_t ustatic[MPD_MINALLOC_MAX]; |
| mpd_uint_t vstatic[MPD_MINALLOC_MAX]; |
| mpd_uint_t *u = ustatic; |
| mpd_uint_t *v = vstatic; |
| mpd_uint_t d, qhat, rhat, w2[2]; |
| mpd_uint_t hi, lo, x; |
| mpd_uint_t carry; |
| mpd_size_t i, j, m; |
| int retval = 0; |
| |
| assert(n > 1 && nplusm >= n); |
| m = sub_size_t(nplusm, n); |
| |
| /* D1: normalize */ |
| d = MPD_RADIX / (vconst[n-1] + 1); |
| |
| if (nplusm >= MPD_MINALLOC_MAX) { |
| if ((u = mpd_alloc(nplusm+1, sizeof *u)) == NULL) { |
| return -1; |
| } |
| } |
| if (n >= MPD_MINALLOC_MAX) { |
| if ((v = mpd_alloc(n+1, sizeof *v)) == NULL) { |
| mpd_free(u); |
| return -1; |
| } |
| } |
| |
| _mpd_shortmul(u, uconst, nplusm, d); |
| _mpd_shortmul(v, vconst, n, d); |
| |
| /* D2: loop */ |
| for (j=m; j != MPD_SIZE_MAX; j--) { |
| |
| /* D3: calculate qhat and rhat */ |
| rhat = _mpd_shortdiv(w2, u+j+n-1, 2, v[n-1]); |
| qhat = w2[1] * MPD_RADIX + w2[0]; |
| |
| while (1) { |
| if (qhat < MPD_RADIX) { |
| _mpd_singlemul(w2, qhat, v[n-2]); |
| if (w2[1] <= rhat) { |
| if (w2[1] != rhat || w2[0] <= u[j+n-2]) { |
| break; |
| } |
| } |
| } |
| qhat -= 1; |
| rhat += v[n-1]; |
| if (rhat < v[n-1] || rhat >= MPD_RADIX) { |
| break; |
| } |
| } |
| /* D4: multiply and subtract */ |
| carry = 0; |
| for (i=0; i <= n; i++) { |
| |
| _mpd_mul_words(&hi, &lo, qhat, v[i]); |
| |
| lo = carry + lo; |
| if (lo < carry) hi++; |
| |
| _mpd_div_words_r(&hi, &lo, hi, lo); |
| |
| x = u[i+j] - lo; |
| carry = (u[i+j] < x); |
| u[i+j] = carry ? x+MPD_RADIX : x; |
| carry += hi; |
| } |
| q[j] = qhat; |
| /* D5: test remainder */ |
| if (carry) { |
| q[j] -= 1; |
| /* D6: add back */ |
| (void)_mpd_baseadd(u+j, u+j, v, n+1, n); |
| } |
| } |
| |
| /* D8: unnormalize */ |
| if (r != NULL) { |
| _mpd_shortdiv(r, u, n, d); |
| /* we are not interested in the return value here */ |
| retval = 0; |
| } |
| else { |
| retval = !_mpd_isallzero(u, n); |
| } |
| |
| |
| if (u != ustatic) mpd_free(u); |
| if (v != vstatic) mpd_free(v); |
| return retval; |
| } |
| |
| /* |
| * Left shift of src by 'shift' digits; src may equal dest. |
| * |
| * dest := area of n mpd_uint_t with space for srcdigits+shift digits. |
| * src := coefficient with length m. |
| * |
| * The case splits in the function are non-obvious. The following |
| * equations might help: |
| * |
| * Let msdigits denote the number of digits in the most significant |
| * word of src. Then 1 <= msdigits <= rdigits. |
| * |
| * 1) shift = q * rdigits + r |
| * 2) srcdigits = qsrc * rdigits + msdigits |
| * 3) destdigits = shift + srcdigits |
| * = q * rdigits + r + qsrc * rdigits + msdigits |
| * = q * rdigits + (qsrc * rdigits + (r + msdigits)) |
| * |
| * The result has q zero words, followed by the coefficient that |
| * is left-shifted by r. The case r == 0 is trivial. For r > 0, it |
| * is important to keep in mind that we always read m source words, |
| * but write m+1 destination words if r + msdigits > rdigits, m words |
| * otherwise. |
| */ |
| void |
| _mpd_baseshiftl(mpd_uint_t *dest, mpd_uint_t *src, mpd_size_t n, mpd_size_t m, |
| mpd_size_t shift) |
| { |
| #if defined(__GNUC__) && !defined(__INTEL_COMPILER) && !defined(__clang__) |
| /* spurious uninitialized warnings */ |
| mpd_uint_t l=l, lprev=lprev, h=h; |
| #else |
| mpd_uint_t l, lprev, h; |
| #endif |
| mpd_uint_t q, r; |
| mpd_uint_t ph; |
| |
| assert(m > 0 && n >= m); |
| |
| _mpd_div_word(&q, &r, (mpd_uint_t)shift, MPD_RDIGITS); |
| |
| if (r != 0) { |
| |
| ph = mpd_pow10[r]; |
| |
| --m; --n; |
| _mpd_divmod_pow10(&h, &lprev, src[m--], MPD_RDIGITS-r); |
| if (h != 0) { /* r + msdigits > rdigits <==> h != 0 */ |
| dest[n--] = h; |
| } |
| /* write m-1 shifted words */ |
| for (; m != MPD_SIZE_MAX; m--,n--) { |
| _mpd_divmod_pow10(&h, &l, src[m], MPD_RDIGITS-r); |
| dest[n] = ph * lprev + h; |
| lprev = l; |
| } |
| /* write least significant word */ |
| dest[q] = ph * lprev; |
| } |
| else { |
| while (--m != MPD_SIZE_MAX) { |
| dest[m+q] = src[m]; |
| } |
| } |
| |
| mpd_uint_zero(dest, q); |
| } |
| |
| /* |
| * Right shift of src by 'shift' digits; src may equal dest. |
| * Assumption: srcdigits-shift > 0. |
| * |
| * dest := area with space for srcdigits-shift digits. |
| * src := coefficient with length 'slen'. |
| * |
| * The case splits in the function rely on the following equations: |
| * |
| * Let msdigits denote the number of digits in the most significant |
| * word of src. Then 1 <= msdigits <= rdigits. |
| * |
| * 1) shift = q * rdigits + r |
| * 2) srcdigits = qsrc * rdigits + msdigits |
| * 3) destdigits = srcdigits - shift |
| * = qsrc * rdigits + msdigits - (q * rdigits + r) |
| * = (qsrc - q) * rdigits + msdigits - r |
| * |
| * Since destdigits > 0 and 1 <= msdigits <= rdigits: |
| * |
| * 4) qsrc >= q |
| * 5) qsrc == q ==> msdigits > r |
| * |
| * The result has slen-q words if msdigits > r, slen-q-1 words otherwise. |
| */ |
| mpd_uint_t |
| _mpd_baseshiftr(mpd_uint_t *dest, mpd_uint_t *src, mpd_size_t slen, |
| mpd_size_t shift) |
| { |
| #if defined(__GNUC__) && !defined(__INTEL_COMPILER) && !defined(__clang__) |
| /* spurious uninitialized warnings */ |
| mpd_uint_t l=l, h=h, hprev=hprev; /* low, high, previous high */ |
| #else |
| mpd_uint_t l, h, hprev; /* low, high, previous high */ |
| #endif |
| mpd_uint_t rnd, rest; /* rounding digit, rest */ |
| mpd_uint_t q, r; |
| mpd_size_t i, j; |
| mpd_uint_t ph; |
| |
| assert(slen > 0); |
| |
| _mpd_div_word(&q, &r, (mpd_uint_t)shift, MPD_RDIGITS); |
| |
| rnd = rest = 0; |
| if (r != 0) { |
| |
| ph = mpd_pow10[MPD_RDIGITS-r]; |
| |
| _mpd_divmod_pow10(&hprev, &rest, src[q], r); |
| _mpd_divmod_pow10(&rnd, &rest, rest, r-1); |
| |
| if (rest == 0 && q > 0) { |
| rest = !_mpd_isallzero(src, q); |
| } |
| /* write slen-q-1 words */ |
| for (j=0,i=q+1; i<slen; i++,j++) { |
| _mpd_divmod_pow10(&h, &l, src[i], r); |
| dest[j] = ph * l + hprev; |
| hprev = h; |
| } |
| /* write most significant word */ |
| if (hprev != 0) { /* always the case if slen==q-1 */ |
| dest[j] = hprev; |
| } |
| } |
| else { |
| if (q > 0) { |
| _mpd_divmod_pow10(&rnd, &rest, src[q-1], MPD_RDIGITS-1); |
| /* is there any non-zero digit below rnd? */ |
| if (rest == 0) rest = !_mpd_isallzero(src, q-1); |
| } |
| for (j = 0; j < slen-q; j++) { |
| dest[j] = src[q+j]; |
| } |
| } |
| |
| /* 0-4 ==> rnd+rest < 0.5 */ |
| /* 5 ==> rnd+rest == 0.5 */ |
| /* 6-9 ==> rnd+rest > 0.5 */ |
| return (rnd == 0 || rnd == 5) ? rnd + !!rest : rnd; |
| } |
| |
| |
| /*********************************************************************/ |
| /* Calculations in base b */ |
| /*********************************************************************/ |
| |
| /* |
| * Add v to w (len m). The calling function has to handle a possible |
| * final carry. Assumption: m > 0. |
| */ |
| mpd_uint_t |
| _mpd_shortadd_b(mpd_uint_t *w, mpd_size_t m, mpd_uint_t v, mpd_uint_t b) |
| { |
| mpd_uint_t s; |
| mpd_uint_t carry; |
| mpd_size_t i; |
| |
| assert(m > 0); |
| |
| /* add v to w */ |
| s = w[0] + v; |
| carry = (s < v) | (s >= b); |
| w[0] = carry ? s-b : s; |
| |
| /* if there is a carry, propagate it */ |
| for (i = 1; carry && i < m; i++) { |
| s = w[i] + carry; |
| carry = (s == b); |
| w[i] = carry ? 0 : s; |
| } |
| |
| return carry; |
| } |
| |
| /* w := product of u (len n) and v (single word). Return carry. */ |
| mpd_uint_t |
| _mpd_shortmul_c(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n, mpd_uint_t v) |
| { |
| mpd_uint_t hi, lo; |
| mpd_uint_t carry = 0; |
| mpd_size_t i; |
| |
| assert(n > 0); |
| |
| for (i=0; i < n; i++) { |
| |
| _mpd_mul_words(&hi, &lo, u[i], v); |
| lo = carry + lo; |
| if (lo < carry) hi++; |
| |
| _mpd_div_words_r(&carry, &w[i], hi, lo); |
| } |
| |
| return carry; |
| } |
| |
| /* w := product of u (len n) and v (single word) */ |
| mpd_uint_t |
| _mpd_shortmul_b(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n, |
| mpd_uint_t v, mpd_uint_t b) |
| { |
| mpd_uint_t hi, lo; |
| mpd_uint_t carry = 0; |
| mpd_size_t i; |
| |
| assert(n > 0); |
| |
| for (i=0; i < n; i++) { |
| |
| _mpd_mul_words(&hi, &lo, u[i], v); |
| lo = carry + lo; |
| if (lo < carry) hi++; |
| |
| _mpd_div_words(&carry, &w[i], hi, lo, b); |
| } |
| |
| return carry; |
| } |
| |
| /* |
| * Knuth, TAOCP Volume 2, 4.3.1, exercise 16: |
| * w := quotient of u (len n) divided by a single word v |
| */ |
| mpd_uint_t |
| _mpd_shortdiv_b(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n, |
| mpd_uint_t v, mpd_uint_t b) |
| { |
| mpd_uint_t hi, lo; |
| mpd_uint_t rem = 0; |
| mpd_size_t i; |
| |
| assert(n > 0); |
| |
| for (i=n-1; i != MPD_SIZE_MAX; i--) { |
| |
| _mpd_mul_words(&hi, &lo, rem, b); |
| lo = u[i] + lo; |
| if (lo < u[i]) hi++; |
| |
| _mpd_div_words(&w[i], &rem, hi, lo, v); |
| } |
| |
| return rem; |
| } |
| |
| |
| |
