ยปCore Development>Code coverage>Modules/_decimal/libmpdec/crt.c

Python code coverage for Modules/_decimal/libmpdec/crt.c

#countcontent
1n/a/*
2n/a * Copyright (c) 2008-2016 Stefan Krah. All rights reserved.
3n/a *
4n/a * Redistribution and use in source and binary forms, with or without
5n/a * modification, are permitted provided that the following conditions
6n/a * are met:
7n/a *
8n/a * 1. Redistributions of source code must retain the above copyright
9n/a * notice, this list of conditions and the following disclaimer.
10n/a *
11n/a * 2. Redistributions in binary form must reproduce the above copyright
12n/a * notice, this list of conditions and the following disclaimer in the
13n/a * documentation and/or other materials provided with the distribution.
14n/a *
15n/a * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS "AS IS" AND
16n/a * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
17n/a * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
18n/a * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
19n/a * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
20n/a * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
21n/a * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
22n/a * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
23n/a * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
24n/a * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
25n/a * SUCH DAMAGE.
26n/a */
27n/a
28n/a
29n/a#include "mpdecimal.h"
30n/a#include <stdio.h>
31n/a#include <assert.h>
32n/a#include "numbertheory.h"
33n/a#include "umodarith.h"
34n/a#include "crt.h"
35n/a
36n/a
37n/a/* Bignum: Chinese Remainder Theorem, extends the maximum transform length. */
38n/a
39n/a
40n/a/* Multiply P1P2 by v, store result in w. */
41n/astatic inline void
42n/a_crt_mulP1P2_3(mpd_uint_t w[3], mpd_uint_t v)
43n/a{
44n/a mpd_uint_t hi1, hi2, lo;
45n/a
46n/a _mpd_mul_words(&hi1, &lo, LH_P1P2, v);
47n/a w[0] = lo;
48n/a
49n/a _mpd_mul_words(&hi2, &lo, UH_P1P2, v);
50n/a lo = hi1 + lo;
51n/a if (lo < hi1) hi2++;
52n/a
53n/a w[1] = lo;
54n/a w[2] = hi2;
55n/a}
56n/a
57n/a/* Add 3 words from v to w. The result is known to fit in w. */
58n/astatic inline void
59n/a_crt_add3(mpd_uint_t w[3], mpd_uint_t v[3])
60n/a{
61n/a mpd_uint_t carry;
62n/a mpd_uint_t s;
63n/a
64n/a s = w[0] + v[0];
65n/a carry = (s < w[0]);
66n/a w[0] = s;
67n/a
68n/a s = w[1] + (v[1] + carry);
69n/a carry = (s < w[1]);
70n/a w[1] = s;
71n/a
72n/a w[2] = w[2] + (v[2] + carry);
73n/a}
74n/a
75n/a/* Divide 3 words in u by v, store result in w, return remainder. */
76n/astatic inline mpd_uint_t
77n/a_crt_div3(mpd_uint_t *w, const mpd_uint_t *u, mpd_uint_t v)
78n/a{
79n/a mpd_uint_t r1 = u[2];
80n/a mpd_uint_t r2;
81n/a
82n/a if (r1 < v) {
83n/a w[2] = 0;
84n/a }
85n/a else {
86n/a _mpd_div_word(&w[2], &r1, u[2], v); /* GCOV_NOT_REACHED */
87n/a }
88n/a
89n/a _mpd_div_words(&w[1], &r2, r1, u[1], v);
90n/a _mpd_div_words(&w[0], &r1, r2, u[0], v);
91n/a
92n/a return r1;
93n/a}
94n/a
95n/a
96n/a/*
97n/a * Chinese Remainder Theorem:
98n/a * Algorithm from Joerg Arndt, "Matters Computational",
99n/a * Chapter 37.4.1 [http://www.jjj.de/fxt/]
100n/a *
101n/a * See also Knuth, TAOCP, Volume 2, 4.3.2, exercise 7.
102n/a */
103n/a
104n/a/*
105n/a * CRT with carry: x1, x2, x3 contain numbers modulo p1, p2, p3. For each
106n/a * triple of members of the arrays, find the unique z modulo p1*p2*p3, with
107n/a * zmax = p1*p2*p3 - 1.
108n/a *
109n/a * In each iteration of the loop, split z into result[i] = z % MPD_RADIX
110n/a * and carry = z / MPD_RADIX. Let N be the size of carry[] and cmax the
111n/a * maximum carry.
112n/a *
113n/a * Limits for the 32-bit build:
114n/a *
115n/a * N = 2**96
116n/a * cmax = 7711435591312380274
117n/a *
118n/a * Limits for the 64 bit build:
119n/a *
120n/a * N = 2**192
121n/a * cmax = 627710135393475385904124401220046371710
122n/a *
123n/a * The following statements hold for both versions:
124n/a *
125n/a * 1) cmax + zmax < N, so the addition does not overflow.
126n/a *
127n/a * 2) (cmax + zmax) / MPD_RADIX == cmax.
128n/a *
129n/a * 3) If c <= cmax, then c_next = (c + zmax) / MPD_RADIX <= cmax.
130n/a */
131n/avoid
132n/acrt3(mpd_uint_t *x1, mpd_uint_t *x2, mpd_uint_t *x3, mpd_size_t rsize)
133n/a{
134n/a mpd_uint_t p1 = mpd_moduli[P1];
135n/a mpd_uint_t umod;
136n/a#ifdef PPRO
137n/a double dmod;
138n/a uint32_t dinvmod[3];
139n/a#endif
140n/a mpd_uint_t a1, a2, a3;
141n/a mpd_uint_t s;
142n/a mpd_uint_t z[3], t[3];
143n/a mpd_uint_t carry[3] = {0,0,0};
144n/a mpd_uint_t hi, lo;
145n/a mpd_size_t i;
146n/a
147n/a for (i = 0; i < rsize; i++) {
148n/a
149n/a a1 = x1[i];
150n/a a2 = x2[i];
151n/a a3 = x3[i];
152n/a
153n/a SETMODULUS(P2);
154n/a s = ext_submod(a2, a1, umod);
155n/a s = MULMOD(s, INV_P1_MOD_P2);
156n/a
157n/a _mpd_mul_words(&hi, &lo, s, p1);
158n/a lo = lo + a1;
159n/a if (lo < a1) hi++;
160n/a
161n/a SETMODULUS(P3);
162n/a s = dw_submod(a3, hi, lo, umod);
163n/a s = MULMOD(s, INV_P1P2_MOD_P3);
164n/a
165n/a z[0] = lo;
166n/a z[1] = hi;
167n/a z[2] = 0;
168n/a
169n/a _crt_mulP1P2_3(t, s);
170n/a _crt_add3(z, t);
171n/a _crt_add3(carry, z);
172n/a
173n/a x1[i] = _crt_div3(carry, carry, MPD_RADIX);
174n/a }
175n/a
176n/a assert(carry[0] == 0 && carry[1] == 0 && carry[2] == 0);
177n/a}
178n/a
179n/a