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

#	count	content
1	n/a	/*
2	n/a	* Copyright (c) 2008-2016 Stefan Krah. All rights reserved.
3	n/a	*
4	n/a	* Redistribution and use in source and binary forms, with or without
5	n/a	* modification, are permitted provided that the following conditions
6	n/a	* are met:
7	n/a	*
8	n/a	* 1. Redistributions of source code must retain the above copyright
9	n/a	* notice, this list of conditions and the following disclaimer.
10	n/a	*
11	n/a	* 2. Redistributions in binary form must reproduce the above copyright
12	n/a	* notice, this list of conditions and the following disclaimer in the
13	n/a	* documentation and/or other materials provided with the distribution.
14	n/a	*
15	n/a	* THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS "AS IS" AND
16	n/a	* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
17	n/a	* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
18	n/a	* ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
19	n/a	* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
20	n/a	* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
21	n/a	* OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
22	n/a	* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
23	n/a	* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
24	n/a	* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
25	n/a	* SUCH DAMAGE.
26	n/a	*/
27	n/a
28	n/a
29	n/a	#include "mpdecimal.h"
30	n/a	#include <stdio.h>
31	n/a	#include <stdlib.h>
32	n/a	#include <assert.h>
33	n/a	#include "bits.h"
34	n/a	#include "difradix2.h"
35	n/a	#include "numbertheory.h"
36	n/a	#include "transpose.h"
37	n/a	#include "umodarith.h"
38	n/a	#include "sixstep.h"
39	n/a
40	n/a
41	n/a	/* Bignum: Cache efficient Matrix Fourier Transform for arrays of the
42	n/a	form 2*n (See literature/six-step.txt). /
43	n/a
44	n/a
45	n/a	/* forward transform with sign = -1 */
46	n/a	int
47	n/a	six_step_fnt(mpd_uint_t *a, mpd_size_t n, int modnum)
48	n/a	{
49	n/a	struct fnt_params *tparams;
50	n/a	mpd_size_t log2n, C, R;
51	n/a	mpd_uint_t kernel;
52	n/a	mpd_uint_t umod;
53	n/a	#ifdef PPRO
54	n/a	double dmod;
55	n/a	uint32_t dinvmod[3];
56	n/a	#endif
57	n/a	mpd_uint_t *x, w0, w1, wstep;
58	n/a	mpd_size_t i, k;
59	n/a
60	n/a
61	n/a	assert(ispower2(n));
62	n/a	assert(n >= 16);
63	n/a	assert(n <= MPD_MAXTRANSFORM_2N);
64	n/a
65	n/a	log2n = mpd_bsr(n);
66	n/a	C = ((mpd_size_t)1) << (log2n / 2); /* number of columns */
67	n/a	R = ((mpd_size_t)1) << (log2n - (log2n / 2)); /* number of rows */
68	n/a
69	n/a
70	n/a	/* Transpose the matrix. */
71	n/a	if (!transpose_pow2(a, R, C)) {
72	n/a	return 0;
73	n/a	}
74	n/a
75	n/a	/* Length R transform on the rows. */
76	n/a	if ((tparams = _mpd_init_fnt_params(R, -1, modnum)) == NULL) {
77	n/a	return 0;
78	n/a	}
79	n/a	for (x = a; x < a+n; x += R) {
80	n/a	fnt_dif2(x, R, tparams);
81	n/a	}
82	n/a
83	n/a	/* Transpose the matrix. */
84	n/a	if (!transpose_pow2(a, C, R)) {
85	n/a	mpd_free(tparams);
86	n/a	return 0;
87	n/a	}
88	n/a
89	n/a	/* Multiply each matrix element (addressed by iC+k) by r(ik). */
90	n/a	SETMODULUS(modnum);
91	n/a	kernel = _mpd_getkernel(n, -1, modnum);
92	n/a	for (i = 1; i < R; i++) {
93	n/a	w0 = 1; /* r*(i0): initial value for k=0 */
94	n/a	w1 = POWMOD(kernel, i); /* r*(i1): initial value for k=1 */
95	n/a	wstep = MULMOD(w1, w1); /* r*(2i) */
96	n/a	for (k = 0; k < C; k += 2) {
97	n/a	mpd_uint_t x0 = a[i*C+k];
98	n/a	mpd_uint_t x1 = a[i*C+k+1];
99	n/a	MULMOD2(&x0, w0, &x1, w1);
100	n/a	MULMOD2C(&w0, &w1, wstep); /* r*(i(k+2)) = r*(ik) * r*(2i) */
101	n/a	a[i*C+k] = x0;
102	n/a	a[i*C+k+1] = x1;
103	n/a	}
104	n/a	}
105	n/a
106	n/a	/* Length C transform on the rows. */
107	n/a	if (C != R) {
108	n/a	mpd_free(tparams);
109	n/a	if ((tparams = _mpd_init_fnt_params(C, -1, modnum)) == NULL) {
110	n/a	return 0;
111	n/a	}
112	n/a	}
113	n/a	for (x = a; x < a+n; x += C) {
114	n/a	fnt_dif2(x, C, tparams);
115	n/a	}
116	n/a	mpd_free(tparams);
117	n/a
118	n/a	#if 0
119	n/a	/* An unordered transform is sufficient for convolution. */
120	n/a	/* Transpose the matrix. */
121	n/a	if (!transpose_pow2(a, R, C)) {
122	n/a	return 0;
123	n/a	}
124	n/a	#endif
125	n/a
126	n/a	return 1;
127	n/a	}
128	n/a
129	n/a
130	n/a	/* reverse transform, sign = 1 */
131	n/a	int
132	n/a	inv_six_step_fnt(mpd_uint_t *a, mpd_size_t n, int modnum)
133	n/a	{
134	n/a	struct fnt_params *tparams;
135	n/a	mpd_size_t log2n, C, R;
136	n/a	mpd_uint_t kernel;
137	n/a	mpd_uint_t umod;
138	n/a	#ifdef PPRO
139	n/a	double dmod;
140	n/a	uint32_t dinvmod[3];
141	n/a	#endif
142	n/a	mpd_uint_t *x, w0, w1, wstep;
143	n/a	mpd_size_t i, k;
144	n/a
145	n/a
146	n/a	assert(ispower2(n));
147	n/a	assert(n >= 16);
148	n/a	assert(n <= MPD_MAXTRANSFORM_2N);
149	n/a
150	n/a	log2n = mpd_bsr(n);
151	n/a	C = ((mpd_size_t)1) << (log2n / 2); /* number of columns */
152	n/a	R = ((mpd_size_t)1) << (log2n - (log2n / 2)); /* number of rows */
153	n/a
154	n/a
155	n/a	#if 0
156	n/a	/* An unordered transform is sufficient for convolution. */
157	n/a	/* Transpose the matrix, producing an RC matrix. /
158	n/a	if (!transpose_pow2(a, C, R)) {
159	n/a	return 0;
160	n/a	}
161	n/a	#endif
162	n/a
163	n/a	/* Length C transform on the rows. */
164	n/a	if ((tparams = _mpd_init_fnt_params(C, 1, modnum)) == NULL) {
165	n/a	return 0;
166	n/a	}
167	n/a	for (x = a; x < a+n; x += C) {
168	n/a	fnt_dif2(x, C, tparams);
169	n/a	}
170	n/a
171	n/a	/* Multiply each matrix element (addressed by iC+k) by r(ik). */
172	n/a	SETMODULUS(modnum);
173	n/a	kernel = _mpd_getkernel(n, 1, modnum);
174	n/a	for (i = 1; i < R; i++) {
175	n/a	w0 = 1;
176	n/a	w1 = POWMOD(kernel, i);
177	n/a	wstep = MULMOD(w1, w1);
178	n/a	for (k = 0; k < C; k += 2) {
179	n/a	mpd_uint_t x0 = a[i*C+k];
180	n/a	mpd_uint_t x1 = a[i*C+k+1];
181	n/a	MULMOD2(&x0, w0, &x1, w1);
182	n/a	MULMOD2C(&w0, &w1, wstep);
183	n/a	a[i*C+k] = x0;
184	n/a	a[i*C+k+1] = x1;
185	n/a	}
186	n/a	}
187	n/a
188	n/a	/* Transpose the matrix. */
189	n/a	if (!transpose_pow2(a, R, C)) {
190	n/a	mpd_free(tparams);
191	n/a	return 0;
192	n/a	}
193	n/a
194	n/a	/* Length R transform on the rows. */
195	n/a	if (R != C) {
196	n/a	mpd_free(tparams);
197	n/a	if ((tparams = _mpd_init_fnt_params(R, 1, modnum)) == NULL) {
198	n/a	return 0;
199	n/a	}
200	n/a	}
201	n/a	for (x = a; x < a+n; x += R) {
202	n/a	fnt_dif2(x, R, tparams);
203	n/a	}
204	n/a	mpd_free(tparams);
205	n/a
206	n/a	/* Transpose the matrix. */
207	n/a	if (!transpose_pow2(a, C, R)) {
208	n/a	return 0;
209	n/a	}
210	n/a
211	n/a	return 1;
212	n/a	}
213	n/a
214	n/a