diff options
Diffstat (limited to 'Modules/_decimal/libmpdec/crt.c')
| -rw-r--r-- | Modules/_decimal/libmpdec/crt.c | 179 |
1 files changed, 0 insertions, 179 deletions
diff --git a/Modules/_decimal/libmpdec/crt.c b/Modules/_decimal/libmpdec/crt.c deleted file mode 100644 index 4a1e80a..0000000 --- a/Modules/_decimal/libmpdec/crt.c +++ /dev/null @@ -1,179 +0,0 @@ -/* - * 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.h" -#include <stdio.h> -#include <assert.h> -#include "numbertheory.h" -#include "umodarith.h" -#include "crt.h" - - -/* Bignum: Chinese Remainder Theorem, extends the maximum transform length. */ - - -/* Multiply P1P2 by v, store result in w. */ -static inline void -_crt_mulP1P2_3(mpd_uint_t w[3], mpd_uint_t v) -{ - mpd_uint_t hi1, hi2, lo; - - _mpd_mul_words(&hi1, &lo, LH_P1P2, v); - w[0] = lo; - - _mpd_mul_words(&hi2, &lo, UH_P1P2, v); - lo = hi1 + lo; - if (lo < hi1) hi2++; - - w[1] = lo; - w[2] = hi2; -} - -/* Add 3 words from v to w. The result is known to fit in w. */ -static inline void -_crt_add3(mpd_uint_t w[3], mpd_uint_t v[3]) -{ - mpd_uint_t carry; - mpd_uint_t s; - - s = w[0] + v[0]; - carry = (s < w[0]); - w[0] = s; - - s = w[1] + (v[1] + carry); - carry = (s < w[1]); - w[1] = s; - - w[2] = w[2] + (v[2] + carry); -} - -/* Divide 3 words in u by v, store result in w, return remainder. */ -static inline mpd_uint_t -_crt_div3(mpd_uint_t *w, const mpd_uint_t *u, mpd_uint_t v) -{ - mpd_uint_t r1 = u[2]; - mpd_uint_t r2; - - if (r1 < v) { - w[2] = 0; - } - else { - _mpd_div_word(&w[2], &r1, u[2], v); /* GCOV_NOT_REACHED */ - } - - _mpd_div_words(&w[1], &r2, r1, u[1], v); - _mpd_div_words(&w[0], &r1, r2, u[0], v); - - return r1; -} - - -/* - * Chinese Remainder Theorem: - * Algorithm from Joerg Arndt, "Matters Computational", - * Chapter 37.4.1 [http://www.jjj.de/fxt/] - * - * See also Knuth, TAOCP, Volume 2, 4.3.2, exercise 7. - */ - -/* - * CRT with carry: x1, x2, x3 contain numbers modulo p1, p2, p3. For each - * triple of members of the arrays, find the unique z modulo p1*p2*p3, with - * zmax = p1*p2*p3 - 1. - * - * In each iteration of the loop, split z into result[i] = z % MPD_RADIX - * and carry = z / MPD_RADIX. Let N be the size of carry[] and cmax the - * maximum carry. - * - * Limits for the 32-bit build: - * - * N = 2**96 - * cmax = 7711435591312380274 - * - * Limits for the 64 bit build: - * - * N = 2**192 - * cmax = 627710135393475385904124401220046371710 - * - * The following statements hold for both versions: - * - * 1) cmax + zmax < N, so the addition does not overflow. - * - * 2) (cmax + zmax) / MPD_RADIX == cmax. - * - * 3) If c <= cmax, then c_next = (c + zmax) / MPD_RADIX <= cmax. - */ -void -crt3(mpd_uint_t *x1, mpd_uint_t *x2, mpd_uint_t *x3, mpd_size_t rsize) -{ - mpd_uint_t p1 = mpd_moduli[P1]; - mpd_uint_t umod; -#ifdef PPRO - double dmod; - uint32_t dinvmod[3]; -#endif - mpd_uint_t a1, a2, a3; - mpd_uint_t s; - mpd_uint_t z[3], t[3]; - mpd_uint_t carry[3] = {0,0,0}; - mpd_uint_t hi, lo; - mpd_size_t i; - - for (i = 0; i < rsize; i++) { - - a1 = x1[i]; - a2 = x2[i]; - a3 = x3[i]; - - SETMODULUS(P2); - s = ext_submod(a2, a1, umod); - s = MULMOD(s, INV_P1_MOD_P2); - - _mpd_mul_words(&hi, &lo, s, p1); - lo = lo + a1; - if (lo < a1) hi++; - - SETMODULUS(P3); - s = dw_submod(a3, hi, lo, umod); - s = MULMOD(s, INV_P1P2_MOD_P3); - - z[0] = lo; - z[1] = hi; - z[2] = 0; - - _crt_mulP1P2_3(t, s); - _crt_add3(z, t); - _crt_add3(carry, z); - - x1[i] = _crt_div3(carry, carry, MPD_RADIX); - } - - assert(carry[0] == 0 && carry[1] == 0 && carry[2] == 0); -} - - |