Python code coverage for Modules/_decimal/libmpdec/transpose.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 <string.h>
33	n/a	#include <limits.h>
34	n/a	#include <assert.h>
35	n/a	#include "bits.h"
36	n/a	#include "constants.h"
37	n/a	#include "typearith.h"
38	n/a	#include "transpose.h"
39	n/a
40	n/a
41	n/a	#define BUFSIZE 4096
42	n/a	#define SIDE 128
43	n/a
44	n/a
45	n/a	/* Bignum: The transpose functions are used for very large transforms
46	n/a	in sixstep.c and fourstep.c. */
47	n/a
48	n/a
49	n/a	/* Definition of the matrix transpose */
50	n/a	void
51	n/a	std_trans(mpd_uint_t dest[], mpd_uint_t src[], mpd_size_t rows, mpd_size_t cols)
52	n/a	{
53	n/a	mpd_size_t idest, isrc;
54	n/a	mpd_size_t r, c;
55	n/a
56	n/a	for (r = 0; r < rows; r++) {
57	n/a	isrc = r * cols;
58	n/a	idest = r;
59	n/a	for (c = 0; c < cols; c++) {
60	n/a	dest[idest] = src[isrc];
61	n/a	isrc += 1;
62	n/a	idest += rows;
63	n/a	}
64	n/a	}
65	n/a	}
66	n/a
67	n/a	/*
68	n/a	* Swap half-rows of 2^n * (2*2^n) matrix.
69	n/a	* FORWARD_CYCLE: even/odd permutation of the halfrows.
70	n/a	* BACKWARD_CYCLE: reverse the even/odd permutation.
71	n/a	*/
72	n/a	static int
73	n/a	swap_halfrows_pow2(mpd_uint_t *matrix, mpd_size_t rows, mpd_size_t cols, int dir)
74	n/a	{
75	n/a	mpd_uint_t buf1[BUFSIZE];
76	n/a	mpd_uint_t buf2[BUFSIZE];
77	n/a	mpd_uint_t readbuf, writebuf, *hp;
78	n/a	mpd_size_t *done, dbits;
79	n/a	mpd_size_t b = BUFSIZE, stride;
80	n/a	mpd_size_t hn, hmax; /* halfrow number */
81	n/a	mpd_size_t m, r=0;
82	n/a	mpd_size_t offset;
83	n/a	mpd_size_t next;
84	n/a
85	n/a
86	n/a	assert(cols == mul_size_t(2, rows));
87	n/a
88	n/a	if (dir == FORWARD_CYCLE) {
89	n/a	r = rows;
90	n/a	}
91	n/a	else if (dir == BACKWARD_CYCLE) {
92	n/a	r = 2;
93	n/a	}
94	n/a	else {
95	n/a	abort(); /* GCOV_NOT_REACHED */
96	n/a	}
97	n/a
98	n/a	m = cols - 1;
99	n/a	hmax = rows; /* cycles start at odd halfrows */
100	n/a	dbits = 8 * sizeof *done;
101	n/a	if ((done = mpd_calloc(hmax/(sizeof done) + 1, sizeof done)) == NULL) {
102	n/a	return 0;
103	n/a	}
104	n/a
105	n/a	for (hn = 1; hn <= hmax; hn += 2) {
106	n/a
107	n/a	if (done[hn/dbits] & mpd_bits[hn%dbits]) {
108	n/a	continue;
109	n/a	}
110	n/a
111	n/a	readbuf = buf1; writebuf = buf2;
112	n/a
113	n/a	for (offset = 0; offset < cols/2; offset += b) {
114	n/a
115	n/a	stride = (offset + b < cols/2) ? b : cols/2-offset;
116	n/a
117	n/a	hp = matrix + hn*cols/2;
118	n/a	memcpy(readbuf, hp+offset, stride(sizeof readbuf));
119	n/a	pointerswap(&readbuf, &writebuf);
120	n/a
121	n/a	next = mulmod_size_t(hn, r, m);
122	n/a	hp = matrix + next*cols/2;
123	n/a
124	n/a	while (next != hn) {
125	n/a
126	n/a	memcpy(readbuf, hp+offset, stride(sizeof readbuf));
127	n/a	memcpy(hp+offset, writebuf, stride(sizeof writebuf));
128	n/a	pointerswap(&readbuf, &writebuf);
129	n/a
130	n/a	done[next/dbits] \|= mpd_bits[next%dbits];
131	n/a
132	n/a	next = mulmod_size_t(next, r, m);
133	n/a	hp = matrix + next*cols/2;
134	n/a
135	n/a	}
136	n/a
137	n/a	memcpy(hp+offset, writebuf, stride(sizeof writebuf));
138	n/a
139	n/a	done[hn/dbits] \|= mpd_bits[hn%dbits];
140	n/a	}
141	n/a	}
142	n/a
143	n/a	mpd_free(done);
144	n/a	return 1;
145	n/a	}
146	n/a
147	n/a	/* In-place transpose of a square matrix */
148	n/a	static inline void
149	n/a	squaretrans(mpd_uint_t *buf, mpd_size_t cols)
150	n/a	{
151	n/a	mpd_uint_t tmp;
152	n/a	mpd_size_t idest, isrc;
153	n/a	mpd_size_t r, c;
154	n/a
155	n/a	for (r = 0; r < cols; r++) {
156	n/a	c = r+1;
157	n/a	isrc = r*cols + c;
158	n/a	idest = c*cols + r;
159	n/a	for (c = r+1; c < cols; c++) {
160	n/a	tmp = buf[isrc];
161	n/a	buf[isrc] = buf[idest];
162	n/a	buf[idest] = tmp;
163	n/a	isrc += 1;
164	n/a	idest += cols;
165	n/a	}
166	n/a	}
167	n/a	}
168	n/a
169	n/a	/*
170	n/a	* Transpose 2^n * 2^n matrix. For cache efficiency, the matrix is split into
171	n/a	* square blocks with side length 'SIDE'. First, the blocks are transposed,
172	n/a	* then a square transposition is done on each individual block.
173	n/a	*/
174	n/a	static void
175	n/a	squaretrans_pow2(mpd_uint_t *matrix, mpd_size_t size)
176	n/a	{
177	n/a	mpd_uint_t buf1[SIDE*SIDE];
178	n/a	mpd_uint_t buf2[SIDE*SIDE];
179	n/a	mpd_uint_t to, from;
180	n/a	mpd_size_t b = size;
181	n/a	mpd_size_t r, c;
182	n/a	mpd_size_t i;
183	n/a
184	n/a	while (b > SIDE) b >>= 1;
185	n/a
186	n/a	for (r = 0; r < size; r += b) {
187	n/a
188	n/a	for (c = r; c < size; c += b) {
189	n/a
190	n/a	from = matrix + r*size + c;
191	n/a	to = buf1;
192	n/a	for (i = 0; i < b; i++) {
193	n/a	memcpy(to, from, b(sizeof to));
194	n/a	from += size;
195	n/a	to += b;
196	n/a	}
197	n/a	squaretrans(buf1, b);
198	n/a
199	n/a	if (r == c) {
200	n/a	to = matrix + r*size + c;
201	n/a	from = buf1;
202	n/a	for (i = 0; i < b; i++) {
203	n/a	memcpy(to, from, b(sizeof to));
204	n/a	from += b;
205	n/a	to += size;
206	n/a	}
207	n/a	continue;
208	n/a	}
209	n/a	else {
210	n/a	from = matrix + c*size + r;
211	n/a	to = buf2;
212	n/a	for (i = 0; i < b; i++) {
213	n/a	memcpy(to, from, b(sizeof to));
214	n/a	from += size;
215	n/a	to += b;
216	n/a	}
217	n/a	squaretrans(buf2, b);
218	n/a
219	n/a	to = matrix + c*size + r;
220	n/a	from = buf1;
221	n/a	for (i = 0; i < b; i++) {
222	n/a	memcpy(to, from, b(sizeof to));
223	n/a	from += b;
224	n/a	to += size;
225	n/a	}
226	n/a
227	n/a	to = matrix + r*size + c;
228	n/a	from = buf2;
229	n/a	for (i = 0; i < b; i++) {
230	n/a	memcpy(to, from, b(sizeof to));
231	n/a	from += b;
232	n/a	to += size;
233	n/a	}
234	n/a	}
235	n/a	}
236	n/a	}
237	n/a
238	n/a	}
239	n/a
240	n/a	/*
241	n/a	* In-place transposition of a 2^n x 2^n or a 2^n x (2*2^n)
242	n/a	* or a (2*2^n) x 2^n matrix.
243	n/a	*/
244	n/a	int
245	n/a	transpose_pow2(mpd_uint_t *matrix, mpd_size_t rows, mpd_size_t cols)
246	n/a	{
247	n/a	mpd_size_t size = mul_size_t(rows, cols);
248	n/a
249	n/a	assert(ispower2(rows));
250	n/a	assert(ispower2(cols));
251	n/a
252	n/a	if (cols == rows) {
253	n/a	squaretrans_pow2(matrix, rows);
254	n/a	}
255	n/a	else if (cols == mul_size_t(2, rows)) {
256	n/a	if (!swap_halfrows_pow2(matrix, rows, cols, FORWARD_CYCLE)) {
257	n/a	return 0;
258	n/a	}
259	n/a	squaretrans_pow2(matrix, rows);
260	n/a	squaretrans_pow2(matrix+(size/2), rows);
261	n/a	}
262	n/a	else if (rows == mul_size_t(2, cols)) {
263	n/a	squaretrans_pow2(matrix, cols);
264	n/a	squaretrans_pow2(matrix+(size/2), cols);
265	n/a	if (!swap_halfrows_pow2(matrix, cols, rows, BACKWARD_CYCLE)) {
266	n/a	return 0;
267	n/a	}
268	n/a	}
269	n/a	else {
270	n/a	abort(); /* GCOV_NOT_REACHED */
271	n/a	}
272	n/a
273	n/a	return 1;
274	n/a	}
275	n/a
276	n/a