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