summaryrefslogtreecommitdiffstats
path: root/Modules/_decimal/libmpdec/transpose.c
diff options
context:
space:
mode:
Diffstat (limited to 'Modules/_decimal/libmpdec/transpose.c')
-rw-r--r--Modules/_decimal/libmpdec/transpose.c276
1 files changed, 0 insertions, 276 deletions
diff --git a/Modules/_decimal/libmpdec/transpose.c b/Modules/_decimal/libmpdec/transpose.c
deleted file mode 100644
index 55d6d89..0000000
--- a/Modules/_decimal/libmpdec/transpose.c
+++ /dev/null
@@ -1,276 +0,0 @@
-/*
- * Copyright (c) 2008-2016 Stefan Krah. All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- *
- * 1. Redistributions of source code must retain the above copyright
- * notice, this list of conditions and the following disclaimer.
- *
- * 2. Redistributions in binary form must reproduce the above copyright
- * notice, this list of conditions and the following disclaimer in the
- * documentation and/or other materials provided with the distribution.
- *
- * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS "AS IS" AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- */
-
-
-#include "mpdecimal.h"
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include <limits.h>
-#include <assert.h>
-#include "bits.h"
-#include "constants.h"
-#include "typearith.h"
-#include "transpose.h"
-
-
-#define BUFSIZE 4096
-#define SIDE 128
-
-
-/* Bignum: The transpose functions are used for very large transforms
- in sixstep.c and fourstep.c. */
-
-
-/* Definition of the matrix transpose */
-void
-std_trans(mpd_uint_t dest[], mpd_uint_t src[], mpd_size_t rows, mpd_size_t cols)
-{
- mpd_size_t idest, isrc;
- mpd_size_t r, c;
-
- for (r = 0; r < rows; r++) {
- isrc = r * cols;
- idest = r;
- for (c = 0; c < cols; c++) {
- dest[idest] = src[isrc];
- isrc += 1;
- idest += rows;
- }
- }
-}
-
-/*
- * Swap half-rows of 2^n * (2*2^n) matrix.
- * FORWARD_CYCLE: even/odd permutation of the halfrows.
- * BACKWARD_CYCLE: reverse the even/odd permutation.
- */
-static int
-swap_halfrows_pow2(mpd_uint_t *matrix, mpd_size_t rows, mpd_size_t cols, int dir)
-{
- mpd_uint_t buf1[BUFSIZE];
- mpd_uint_t buf2[BUFSIZE];
- mpd_uint_t *readbuf, *writebuf, *hp;
- mpd_size_t *done, dbits;
- mpd_size_t b = BUFSIZE, stride;
- mpd_size_t hn, hmax; /* halfrow number */
- mpd_size_t m, r=0;
- mpd_size_t offset;
- mpd_size_t next;
-
-
- assert(cols == mul_size_t(2, rows));
-
- if (dir == FORWARD_CYCLE) {
- r = rows;
- }
- else if (dir == BACKWARD_CYCLE) {
- r = 2;
- }
- else {
- abort(); /* GCOV_NOT_REACHED */
- }
-
- m = cols - 1;
- hmax = rows; /* cycles start at odd halfrows */
- dbits = 8 * sizeof *done;
- if ((done = mpd_calloc(hmax/(sizeof *done) + 1, sizeof *done)) == NULL) {
- return 0;
- }
-
- for (hn = 1; hn <= hmax; hn += 2) {
-
- if (done[hn/dbits] & mpd_bits[hn%dbits]) {
- continue;
- }
-
- readbuf = buf1; writebuf = buf2;
-
- for (offset = 0; offset < cols/2; offset += b) {
-
- stride = (offset + b < cols/2) ? b : cols/2-offset;
-
- hp = matrix + hn*cols/2;
- memcpy(readbuf, hp+offset, stride*(sizeof *readbuf));
- pointerswap(&readbuf, &writebuf);
-
- next = mulmod_size_t(hn, r, m);
- hp = matrix + next*cols/2;
-
- while (next != hn) {
-
- memcpy(readbuf, hp+offset, stride*(sizeof *readbuf));
- memcpy(hp+offset, writebuf, stride*(sizeof *writebuf));
- pointerswap(&readbuf, &writebuf);
-
- done[next/dbits] |= mpd_bits[next%dbits];
-
- next = mulmod_size_t(next, r, m);
- hp = matrix + next*cols/2;
-
- }
-
- memcpy(hp+offset, writebuf, stride*(sizeof *writebuf));
-
- done[hn/dbits] |= mpd_bits[hn%dbits];
- }
- }
-
- mpd_free(done);
- return 1;
-}
-
-/* In-place transpose of a square matrix */
-static inline void
-squaretrans(mpd_uint_t *buf, mpd_size_t cols)
-{
- mpd_uint_t tmp;
- mpd_size_t idest, isrc;
- mpd_size_t r, c;
-
- for (r = 0; r < cols; r++) {
- c = r+1;
- isrc = r*cols + c;
- idest = c*cols + r;
- for (c = r+1; c < cols; c++) {
- tmp = buf[isrc];
- buf[isrc] = buf[idest];
- buf[idest] = tmp;
- isrc += 1;
- idest += cols;
- }
- }
-}
-
-/*
- * Transpose 2^n * 2^n matrix. For cache efficiency, the matrix is split into
- * square blocks with side length 'SIDE'. First, the blocks are transposed,
- * then a square transposition is done on each individual block.
- */
-static void
-squaretrans_pow2(mpd_uint_t *matrix, mpd_size_t size)
-{
- mpd_uint_t buf1[SIDE*SIDE];
- mpd_uint_t buf2[SIDE*SIDE];
- mpd_uint_t *to, *from;
- mpd_size_t b = size;
- mpd_size_t r, c;
- mpd_size_t i;
-
- while (b > SIDE) b >>= 1;
-
- for (r = 0; r < size; r += b) {
-
- for (c = r; c < size; c += b) {
-
- from = matrix + r*size + c;
- to = buf1;
- for (i = 0; i < b; i++) {
- memcpy(to, from, b*(sizeof *to));
- from += size;
- to += b;
- }
- squaretrans(buf1, b);
-
- if (r == c) {
- to = matrix + r*size + c;
- from = buf1;
- for (i = 0; i < b; i++) {
- memcpy(to, from, b*(sizeof *to));
- from += b;
- to += size;
- }
- continue;
- }
- else {
- from = matrix + c*size + r;
- to = buf2;
- for (i = 0; i < b; i++) {
- memcpy(to, from, b*(sizeof *to));
- from += size;
- to += b;
- }
- squaretrans(buf2, b);
-
- to = matrix + c*size + r;
- from = buf1;
- for (i = 0; i < b; i++) {
- memcpy(to, from, b*(sizeof *to));
- from += b;
- to += size;
- }
-
- to = matrix + r*size + c;
- from = buf2;
- for (i = 0; i < b; i++) {
- memcpy(to, from, b*(sizeof *to));
- from += b;
- to += size;
- }
- }
- }
- }
-
-}
-
-/*
- * In-place transposition of a 2^n x 2^n or a 2^n x (2*2^n)
- * or a (2*2^n) x 2^n matrix.
- */
-int
-transpose_pow2(mpd_uint_t *matrix, mpd_size_t rows, mpd_size_t cols)
-{
- mpd_size_t size = mul_size_t(rows, cols);
-
- assert(ispower2(rows));
- assert(ispower2(cols));
-
- if (cols == rows) {
- squaretrans_pow2(matrix, rows);
- }
- else if (cols == mul_size_t(2, rows)) {
- if (!swap_halfrows_pow2(matrix, rows, cols, FORWARD_CYCLE)) {
- return 0;
- }
- squaretrans_pow2(matrix, rows);
- squaretrans_pow2(matrix+(size/2), rows);
- }
- else if (rows == mul_size_t(2, cols)) {
- squaretrans_pow2(matrix, cols);
- squaretrans_pow2(matrix+(size/2), cols);
- if (!swap_halfrows_pow2(matrix, cols, rows, BACKWARD_CYCLE)) {
- return 0;
- }
- }
- else {
- abort(); /* GCOV_NOT_REACHED */
- }
-
- return 1;
-}
-
-