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

Python code coverage for Modules/_decimal/libmpdec/fourstep.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 <assert.h>
31n/a#include "numbertheory.h"
32n/a#include "sixstep.h"
33n/a#include "transpose.h"
34n/a#include "umodarith.h"
35n/a#include "fourstep.h"
36n/a
37n/a
38n/a/* Bignum: Cache efficient Matrix Fourier Transform for arrays of the
39n/a form 3 * 2**n (See literature/matrix-transform.txt). */
40n/a
41n/a
42n/a#ifndef PPRO
43n/astatic inline void
44n/astd_size3_ntt(mpd_uint_t *x1, mpd_uint_t *x2, mpd_uint_t *x3,
45n/a mpd_uint_t w3table[3], mpd_uint_t umod)
46n/a{
47n/a mpd_uint_t r1, r2;
48n/a mpd_uint_t w;
49n/a mpd_uint_t s, tmp;
50n/a
51n/a
52n/a /* k = 0 -> w = 1 */
53n/a s = *x1;
54n/a s = addmod(s, *x2, umod);
55n/a s = addmod(s, *x3, umod);
56n/a
57n/a r1 = s;
58n/a
59n/a /* k = 1 */
60n/a s = *x1;
61n/a
62n/a w = w3table[1];
63n/a tmp = MULMOD(*x2, w);
64n/a s = addmod(s, tmp, umod);
65n/a
66n/a w = w3table[2];
67n/a tmp = MULMOD(*x3, w);
68n/a s = addmod(s, tmp, umod);
69n/a
70n/a r2 = s;
71n/a
72n/a /* k = 2 */
73n/a s = *x1;
74n/a
75n/a w = w3table[2];
76n/a tmp = MULMOD(*x2, w);
77n/a s = addmod(s, tmp, umod);
78n/a
79n/a w = w3table[1];
80n/a tmp = MULMOD(*x3, w);
81n/a s = addmod(s, tmp, umod);
82n/a
83n/a *x3 = s;
84n/a *x2 = r2;
85n/a *x1 = r1;
86n/a}
87n/a#else /* PPRO */
88n/astatic inline void
89n/appro_size3_ntt(mpd_uint_t *x1, mpd_uint_t *x2, mpd_uint_t *x3, mpd_uint_t w3table[3],
90n/a mpd_uint_t umod, double *dmod, uint32_t dinvmod[3])
91n/a{
92n/a mpd_uint_t r1, r2;
93n/a mpd_uint_t w;
94n/a mpd_uint_t s, tmp;
95n/a
96n/a
97n/a /* k = 0 -> w = 1 */
98n/a s = *x1;
99n/a s = addmod(s, *x2, umod);
100n/a s = addmod(s, *x3, umod);
101n/a
102n/a r1 = s;
103n/a
104n/a /* k = 1 */
105n/a s = *x1;
106n/a
107n/a w = w3table[1];
108n/a tmp = ppro_mulmod(*x2, w, dmod, dinvmod);
109n/a s = addmod(s, tmp, umod);
110n/a
111n/a w = w3table[2];
112n/a tmp = ppro_mulmod(*x3, w, dmod, dinvmod);
113n/a s = addmod(s, tmp, umod);
114n/a
115n/a r2 = s;
116n/a
117n/a /* k = 2 */
118n/a s = *x1;
119n/a
120n/a w = w3table[2];
121n/a tmp = ppro_mulmod(*x2, w, dmod, dinvmod);
122n/a s = addmod(s, tmp, umod);
123n/a
124n/a w = w3table[1];
125n/a tmp = ppro_mulmod(*x3, w, dmod, dinvmod);
126n/a s = addmod(s, tmp, umod);
127n/a
128n/a *x3 = s;
129n/a *x2 = r2;
130n/a *x1 = r1;
131n/a}
132n/a#endif
133n/a
134n/a
135n/a/* forward transform, sign = -1; transform length = 3 * 2**n */
136n/aint
137n/afour_step_fnt(mpd_uint_t *a, mpd_size_t n, int modnum)
138n/a{
139n/a mpd_size_t R = 3; /* number of rows */
140n/a mpd_size_t C = n / 3; /* number of columns */
141n/a mpd_uint_t w3table[3];
142n/a mpd_uint_t kernel, w0, w1, wstep;
143n/a mpd_uint_t *s, *p0, *p1, *p2;
144n/a mpd_uint_t umod;
145n/a#ifdef PPRO
146n/a double dmod;
147n/a uint32_t dinvmod[3];
148n/a#endif
149n/a mpd_size_t i, k;
150n/a
151n/a
152n/a assert(n >= 48);
153n/a assert(n <= 3*MPD_MAXTRANSFORM_2N);
154n/a
155n/a
156n/a /* Length R transform on the columns. */
157n/a SETMODULUS(modnum);
158n/a _mpd_init_w3table(w3table, -1, modnum);
159n/a for (p0=a, p1=p0+C, p2=p0+2*C; p0<a+C; p0++,p1++,p2++) {
160n/a
161n/a SIZE3_NTT(p0, p1, p2, w3table);
162n/a }
163n/a
164n/a /* Multiply each matrix element (addressed by i*C+k) by r**(i*k). */
165n/a kernel = _mpd_getkernel(n, -1, modnum);
166n/a for (i = 1; i < R; i++) {
167n/a w0 = 1; /* r**(i*0): initial value for k=0 */
168n/a w1 = POWMOD(kernel, i); /* r**(i*1): initial value for k=1 */
169n/a wstep = MULMOD(w1, w1); /* r**(2*i) */
170n/a for (k = 0; k < C-1; k += 2) {
171n/a mpd_uint_t x0 = a[i*C+k];
172n/a mpd_uint_t x1 = a[i*C+k+1];
173n/a MULMOD2(&x0, w0, &x1, w1);
174n/a MULMOD2C(&w0, &w1, wstep); /* r**(i*(k+2)) = r**(i*k) * r**(2*i) */
175n/a a[i*C+k] = x0;
176n/a a[i*C+k+1] = x1;
177n/a }
178n/a }
179n/a
180n/a /* Length C transform on the rows. */
181n/a for (s = a; s < a+n; s += C) {
182n/a if (!six_step_fnt(s, C, modnum)) {
183n/a return 0;
184n/a }
185n/a }
186n/a
187n/a#if 0
188n/a /* An unordered transform is sufficient for convolution. */
189n/a /* Transpose the matrix. */
190n/a transpose_3xpow2(a, R, C);
191n/a#endif
192n/a
193n/a return 1;
194n/a}
195n/a
196n/a/* backward transform, sign = 1; transform length = 3 * 2**n */
197n/aint
198n/ainv_four_step_fnt(mpd_uint_t *a, mpd_size_t n, int modnum)
199n/a{
200n/a mpd_size_t R = 3; /* number of rows */
201n/a mpd_size_t C = n / 3; /* number of columns */
202n/a mpd_uint_t w3table[3];
203n/a mpd_uint_t kernel, w0, w1, wstep;
204n/a mpd_uint_t *s, *p0, *p1, *p2;
205n/a mpd_uint_t umod;
206n/a#ifdef PPRO
207n/a double dmod;
208n/a uint32_t dinvmod[3];
209n/a#endif
210n/a mpd_size_t i, k;
211n/a
212n/a
213n/a assert(n >= 48);
214n/a assert(n <= 3*MPD_MAXTRANSFORM_2N);
215n/a
216n/a
217n/a#if 0
218n/a /* An unordered transform is sufficient for convolution. */
219n/a /* Transpose the matrix, producing an R*C matrix. */
220n/a transpose_3xpow2(a, C, R);
221n/a#endif
222n/a
223n/a /* Length C transform on the rows. */
224n/a for (s = a; s < a+n; s += C) {
225n/a if (!inv_six_step_fnt(s, C, modnum)) {
226n/a return 0;
227n/a }
228n/a }
229n/a
230n/a /* Multiply each matrix element (addressed by i*C+k) by r**(i*k). */
231n/a SETMODULUS(modnum);
232n/a kernel = _mpd_getkernel(n, 1, modnum);
233n/a for (i = 1; i < R; i++) {
234n/a w0 = 1;
235n/a w1 = POWMOD(kernel, i);
236n/a wstep = MULMOD(w1, w1);
237n/a for (k = 0; k < C; k += 2) {
238n/a mpd_uint_t x0 = a[i*C+k];
239n/a mpd_uint_t x1 = a[i*C+k+1];
240n/a MULMOD2(&x0, w0, &x1, w1);
241n/a MULMOD2C(&w0, &w1, wstep);
242n/a a[i*C+k] = x0;
243n/a a[i*C+k+1] = x1;
244n/a }
245n/a }
246n/a
247n/a /* Length R transform on the columns. */
248n/a _mpd_init_w3table(w3table, 1, modnum);
249n/a for (p0=a, p1=p0+C, p2=p0+2*C; p0<a+C; p0++,p1++,p2++) {
250n/a
251n/a SIZE3_NTT(p0, p1, p2, w3table);
252n/a }
253n/a
254n/a return 1;
255n/a}
256n/a
257n/a