Python code coverage for Modules/_decimal/libmpdec/difradix2.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 <assert.h>
32	n/a	#include "bits.h"
33	n/a	#include "numbertheory.h"
34	n/a	#include "umodarith.h"
35	n/a	#include "difradix2.h"
36	n/a
37	n/a
38	n/a	/* Bignum: The actual transform routine (decimation in frequency). */
39	n/a
40	n/a
41	n/a	/*
42	n/a	* Generate index pairs (x, bitreverse(x)) and carry out the permutation.
43	n/a	* n must be a power of two.
44	n/a	* Algorithm due to Brent/Lehmann, see Joerg Arndt, "Matters Computational",
45	n/a	* Chapter 1.14.4. [http://www.jjj.de/fxt/]
46	n/a	*/
47	n/a	static inline void
48	n/a	bitreverse_permute(mpd_uint_t a[], mpd_size_t n)
49	n/a	{
50	n/a	mpd_size_t x = 0;
51	n/a	mpd_size_t r = 0;
52	n/a	mpd_uint_t t;
53	n/a
54	n/a	do { /* Invariant: r = bitreverse(x) */
55	n/a	if (r > x) {
56	n/a	t = a[x];
57	n/a	a[x] = a[r];
58	n/a	a[r] = t;
59	n/a	}
60	n/a	/* Flip trailing consecutive 1 bits and the first zero bit
61	n/a	* that absorbs a possible carry. */
62	n/a	x += 1;
63	n/a	/* Mirror the operation on r: Flip n_trailing_zeros(x)+1
64	n/a	high bits of r. */
65	n/a	r ^= (n - (n >> (mpd_bsf(x)+1)));
66	n/a	/* The loop invariant is preserved. */
67	n/a	} while (x < n);
68	n/a	}
69	n/a
70	n/a
71	n/a	/* Fast Number Theoretic Transform, decimation in frequency. */
72	n/a	void
73	n/a	fnt_dif2(mpd_uint_t a[], mpd_size_t n, struct fnt_params *tparams)
74	n/a	{
75	n/a	mpd_uint_t *wtable = tparams->wtable;
76	n/a	mpd_uint_t umod;
77	n/a	#ifdef PPRO
78	n/a	double dmod;
79	n/a	uint32_t dinvmod[3];
80	n/a	#endif
81	n/a	mpd_uint_t u0, u1, v0, v1;
82	n/a	mpd_uint_t w, w0, w1, wstep;
83	n/a	mpd_size_t m, mhalf;
84	n/a	mpd_size_t j, r;
85	n/a
86	n/a
87	n/a	assert(ispower2(n));
88	n/a	assert(n >= 4);
89	n/a
90	n/a	SETMODULUS(tparams->modnum);
91	n/a
92	n/a	/* m == n */
93	n/a	mhalf = n / 2;
94	n/a	for (j = 0; j < mhalf; j += 2) {
95	n/a
96	n/a	w0 = wtable[j];
97	n/a	w1 = wtable[j+1];
98	n/a
99	n/a	u0 = a[j];
100	n/a	v0 = a[j+mhalf];
101	n/a
102	n/a	u1 = a[j+1];
103	n/a	v1 = a[j+1+mhalf];
104	n/a
105	n/a	a[j] = addmod(u0, v0, umod);
106	n/a	v0 = submod(u0, v0, umod);
107	n/a
108	n/a	a[j+1] = addmod(u1, v1, umod);
109	n/a	v1 = submod(u1, v1, umod);
110	n/a
111	n/a	MULMOD2(&v0, w0, &v1, w1);
112	n/a
113	n/a	a[j+mhalf] = v0;
114	n/a	a[j+1+mhalf] = v1;
115	n/a
116	n/a	}
117	n/a
118	n/a	wstep = 2;
119	n/a	for (m = n/2; m >= 2; m>>=1, wstep<<=1) {
120	n/a
121	n/a	mhalf = m / 2;
122	n/a
123	n/a	/* j == 0 */
124	n/a	for (r = 0; r < n; r += 2*m) {
125	n/a
126	n/a	u0 = a[r];
127	n/a	v0 = a[r+mhalf];
128	n/a
129	n/a	u1 = a[m+r];
130	n/a	v1 = a[m+r+mhalf];
131	n/a
132	n/a	a[r] = addmod(u0, v0, umod);
133	n/a	v0 = submod(u0, v0, umod);
134	n/a
135	n/a	a[m+r] = addmod(u1, v1, umod);
136	n/a	v1 = submod(u1, v1, umod);
137	n/a
138	n/a	a[r+mhalf] = v0;
139	n/a	a[m+r+mhalf] = v1;
140	n/a	}
141	n/a
142	n/a	for (j = 1; j < mhalf; j++) {
143	n/a
144	n/a	w = wtable[j*wstep];
145	n/a
146	n/a	for (r = 0; r < n; r += 2*m) {
147	n/a
148	n/a	u0 = a[r+j];
149	n/a	v0 = a[r+j+mhalf];
150	n/a
151	n/a	u1 = a[m+r+j];
152	n/a	v1 = a[m+r+j+mhalf];
153	n/a
154	n/a	a[r+j] = addmod(u0, v0, umod);
155	n/a	v0 = submod(u0, v0, umod);
156	n/a
157	n/a	a[m+r+j] = addmod(u1, v1, umod);
158	n/a	v1 = submod(u1, v1, umod);
159	n/a
160	n/a	MULMOD2C(&v0, &v1, w);
161	n/a
162	n/a	a[r+j+mhalf] = v0;
163	n/a	a[m+r+j+mhalf] = v1;
164	n/a	}
165	n/a
166	n/a	}
167	n/a
168	n/a	}
169	n/a
170	n/a	bitreverse_permute(a, n);
171	n/a	}
172	n/a
173	n/a