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

Python code coverage for Modules/_decimal/libmpdec/sixstep.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 <stdlib.h>
32n/a#include <assert.h>
33n/a#include "bits.h"
34n/a#include "difradix2.h"
35n/a#include "numbertheory.h"
36n/a#include "transpose.h"
37n/a#include "umodarith.h"
38n/a#include "sixstep.h"
39n/a
40n/a
41n/a/* Bignum: Cache efficient Matrix Fourier Transform for arrays of the
42n/a form 2**n (See literature/six-step.txt). */
43n/a
44n/a
45n/a/* forward transform with sign = -1 */
46n/aint
47n/asix_step_fnt(mpd_uint_t *a, mpd_size_t n, int modnum)
48n/a{
49n/a struct fnt_params *tparams;
50n/a mpd_size_t log2n, C, R;
51n/a mpd_uint_t kernel;
52n/a mpd_uint_t umod;
53n/a#ifdef PPRO
54n/a double dmod;
55n/a uint32_t dinvmod[3];
56n/a#endif
57n/a mpd_uint_t *x, w0, w1, wstep;
58n/a mpd_size_t i, k;
59n/a
60n/a
61n/a assert(ispower2(n));
62n/a assert(n >= 16);
63n/a assert(n <= MPD_MAXTRANSFORM_2N);
64n/a
65n/a log2n = mpd_bsr(n);
66n/a C = ((mpd_size_t)1) << (log2n / 2); /* number of columns */
67n/a R = ((mpd_size_t)1) << (log2n - (log2n / 2)); /* number of rows */
68n/a
69n/a
70n/a /* Transpose the matrix. */
71n/a if (!transpose_pow2(a, R, C)) {
72n/a return 0;
73n/a }
74n/a
75n/a /* Length R transform on the rows. */
76n/a if ((tparams = _mpd_init_fnt_params(R, -1, modnum)) == NULL) {
77n/a return 0;
78n/a }
79n/a for (x = a; x < a+n; x += R) {
80n/a fnt_dif2(x, R, tparams);
81n/a }
82n/a
83n/a /* Transpose the matrix. */
84n/a if (!transpose_pow2(a, C, R)) {
85n/a mpd_free(tparams);
86n/a return 0;
87n/a }
88n/a
89n/a /* Multiply each matrix element (addressed by i*C+k) by r**(i*k). */
90n/a SETMODULUS(modnum);
91n/a kernel = _mpd_getkernel(n, -1, modnum);
92n/a for (i = 1; i < R; i++) {
93n/a w0 = 1; /* r**(i*0): initial value for k=0 */
94n/a w1 = POWMOD(kernel, i); /* r**(i*1): initial value for k=1 */
95n/a wstep = MULMOD(w1, w1); /* r**(2*i) */
96n/a for (k = 0; k < C; k += 2) {
97n/a mpd_uint_t x0 = a[i*C+k];
98n/a mpd_uint_t x1 = a[i*C+k+1];
99n/a MULMOD2(&x0, w0, &x1, w1);
100n/a MULMOD2C(&w0, &w1, wstep); /* r**(i*(k+2)) = r**(i*k) * r**(2*i) */
101n/a a[i*C+k] = x0;
102n/a a[i*C+k+1] = x1;
103n/a }
104n/a }
105n/a
106n/a /* Length C transform on the rows. */
107n/a if (C != R) {
108n/a mpd_free(tparams);
109n/a if ((tparams = _mpd_init_fnt_params(C, -1, modnum)) == NULL) {
110n/a return 0;
111n/a }
112n/a }
113n/a for (x = a; x < a+n; x += C) {
114n/a fnt_dif2(x, C, tparams);
115n/a }
116n/a mpd_free(tparams);
117n/a
118n/a#if 0
119n/a /* An unordered transform is sufficient for convolution. */
120n/a /* Transpose the matrix. */
121n/a if (!transpose_pow2(a, R, C)) {
122n/a return 0;
123n/a }
124n/a#endif
125n/a
126n/a return 1;
127n/a}
128n/a
129n/a
130n/a/* reverse transform, sign = 1 */
131n/aint
132n/ainv_six_step_fnt(mpd_uint_t *a, mpd_size_t n, int modnum)
133n/a{
134n/a struct fnt_params *tparams;
135n/a mpd_size_t log2n, C, R;
136n/a mpd_uint_t kernel;
137n/a mpd_uint_t umod;
138n/a#ifdef PPRO
139n/a double dmod;
140n/a uint32_t dinvmod[3];
141n/a#endif
142n/a mpd_uint_t *x, w0, w1, wstep;
143n/a mpd_size_t i, k;
144n/a
145n/a
146n/a assert(ispower2(n));
147n/a assert(n >= 16);
148n/a assert(n <= MPD_MAXTRANSFORM_2N);
149n/a
150n/a log2n = mpd_bsr(n);
151n/a C = ((mpd_size_t)1) << (log2n / 2); /* number of columns */
152n/a R = ((mpd_size_t)1) << (log2n - (log2n / 2)); /* number of rows */
153n/a
154n/a
155n/a#if 0
156n/a /* An unordered transform is sufficient for convolution. */
157n/a /* Transpose the matrix, producing an R*C matrix. */
158n/a if (!transpose_pow2(a, C, R)) {
159n/a return 0;
160n/a }
161n/a#endif
162n/a
163n/a /* Length C transform on the rows. */
164n/a if ((tparams = _mpd_init_fnt_params(C, 1, modnum)) == NULL) {
165n/a return 0;
166n/a }
167n/a for (x = a; x < a+n; x += C) {
168n/a fnt_dif2(x, C, tparams);
169n/a }
170n/a
171n/a /* Multiply each matrix element (addressed by i*C+k) by r**(i*k). */
172n/a SETMODULUS(modnum);
173n/a kernel = _mpd_getkernel(n, 1, modnum);
174n/a for (i = 1; i < R; i++) {
175n/a w0 = 1;
176n/a w1 = POWMOD(kernel, i);
177n/a wstep = MULMOD(w1, w1);
178n/a for (k = 0; k < C; k += 2) {
179n/a mpd_uint_t x0 = a[i*C+k];
180n/a mpd_uint_t x1 = a[i*C+k+1];
181n/a MULMOD2(&x0, w0, &x1, w1);
182n/a MULMOD2C(&w0, &w1, wstep);
183n/a a[i*C+k] = x0;
184n/a a[i*C+k+1] = x1;
185n/a }
186n/a }
187n/a
188n/a /* Transpose the matrix. */
189n/a if (!transpose_pow2(a, R, C)) {
190n/a mpd_free(tparams);
191n/a return 0;
192n/a }
193n/a
194n/a /* Length R transform on the rows. */
195n/a if (R != C) {
196n/a mpd_free(tparams);
197n/a if ((tparams = _mpd_init_fnt_params(R, 1, modnum)) == NULL) {
198n/a return 0;
199n/a }
200n/a }
201n/a for (x = a; x < a+n; x += R) {
202n/a fnt_dif2(x, R, tparams);
203n/a }
204n/a mpd_free(tparams);
205n/a
206n/a /* Transpose the matrix. */
207n/a if (!transpose_pow2(a, C, R)) {
208n/a return 0;
209n/a }
210n/a
211n/a return 1;
212n/a}
213n/a
214n/a