/*
 * Copyright (c) 2008-2010 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 "mptypes.h"
#include "typearith.h"
#include "transpose.h"
#include "malloc_fail.h"


/*
 * All the following is pasted from transpose.c for the sole reason of
 * running the tests with uint8_t arrays, which allow greater sizes.
 */


#define BUFSIZE 4096
#define SIDE 128


static inline void pointerswap_c(uint8_t **a, uint8_t **b)
{
	uint8_t *tmp;

	tmp = *b;
	*b = *a;
	*a = tmp;
}

/* Definition of the matrix transpose */
static void
std_trans_c(uint8_t dest[], uint8_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_c(uint8_t *matrix, mpd_size_t rows, mpd_size_t cols, int dir)
{
	uint8_t buf1[BUFSIZE];
	uint8_t buf2[BUFSIZE];
	uint8_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 {
		mpd_err_fatal("swap_halfrows_pow2_c: illegal argument for \"dir\"");
	}

	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_c(&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_c(&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];
		}
	}

	__mingw_dfp_get_globals()->mpd_free(done);
	return 1;
}

/* In-place transpose of a square matrix */
static inline void
squaretrans_c(uint8_t *buf, mpd_size_t cols)
{
	uint8_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 tranposition is done on each individual block.
 */
static void
squaretrans_pow2_c(uint8_t *matrix, mpd_size_t size)
{
	uint8_t buf1[SIDE*SIDE];
	uint8_t buf2[SIDE*SIDE];
	uint8_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_c(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_c(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.
 */
static int
transpose_pow2_c(uint8_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_c(matrix, rows);
	}
	else if (cols == mul_size_t(2, rows)) {
		if (!swap_halfrows_pow2_c(matrix, rows, cols, FORWARD_CYCLE)) {
			return 0;
		}
		squaretrans_pow2_c(matrix, rows);
		squaretrans_pow2_c(matrix+(size/2), rows);
	}
	else if (rows == mul_size_t(2, cols)) {
		squaretrans_pow2_c(matrix, cols);
		squaretrans_pow2_c(matrix+(size/2), cols);
		if (!swap_halfrows_pow2_c(matrix, cols, rows, BACKWARD_CYCLE)) {
			return 0;
		}
	}
	else {
		mpd_err_fatal("transpose_pow2_c: illegal matrix size");
	}

	return 1;
}


#include <time.h>
#define UMOD_ARRAY ((mpd_uint_t)1<<28)
#define UCHAR_ARRAY ((mpd_uint_t)1<<29)

/* The faster in-place functions are tested against std_trans(). */
static int
testit_uint(void (* func)(mpd_uint_t *, mpd_size_t, mpd_size_t), mpd_size_t rows, mpd_size_t cols)
{
	mpd_uint_t *a = NULL, *src = NULL, *dest = NULL;
	clock_t start_fast, end_fast, start_std, end_std;
	mpd_size_t msize;
	mpd_size_t i;
	int ret = 1;

	msize = mul_size_t(rows, cols);
	if ((a = mpd_alloc(msize, sizeof *a)) == NULL) {
		goto error;
	}
	if ((src = mpd_alloc(msize, sizeof *src)) == NULL) {
		goto error;
	}
	if ((dest = mpd_alloc(msize, sizeof *dest)) == NULL) {
		goto error;
	}

	for (i = 0; i < msize; i++) {
		a[i] = src[i] = random();
	}

	start_std= clock();
	std_trans(dest, src, rows, cols);
	end_std= clock();

	start_fast = clock();
	func(a, rows, cols);
	end_fast = clock();

	for (i = 0; i < msize; i++) {
		if (a[i] != dest[i]) {
			fprintf(stderr, "FAIL: a[%"PRI_mpd_size_t"] = %"PRI_mpd_uint_t
                                        "\tdest[%"PRI_mpd_size_t"] = %"PRI_mpd_uint_t"\n",
			        i, a[i], i, dest[i]);
			exit(1);
		}
	}

	fprintf(stderr, "size: %10"PRI_mpd_size_t"\tstd_trans: %6.2f sec\t "
                        "in_place_trans: %6.2f sec\n",
		msize,
		(double)(end_std-start_std)/(double)CLOCKS_PER_SEC,
		(double)(end_fast-start_fast)/(double)CLOCKS_PER_SEC);


out:
	if (a) __mingw_dfp_get_globals()->mpd_free(a);
	if (src) __mingw_dfp_get_globals()->mpd_free(src);
	if (dest) __mingw_dfp_get_globals()->mpd_free(dest);
	return ret;

error:
	ret = 0;
	goto out;
}

/* The faster in-place functions are tested against std_trans(). */
static int
testit_uchar(void (* func)(uint8_t *, mpd_size_t, mpd_size_t), mpd_size_t rows, mpd_size_t cols)
{
	uint8_t *a = NULL, *src = NULL, *dest = NULL;
	clock_t start_fast, end_fast, start_std, end_std;
	mpd_size_t msize;
	mpd_size_t i;
	int ret = 1;

	msize = mul_size_t(rows, cols);
	if ((a = mpd_alloc(msize, sizeof *a)) == NULL) {
		goto error;
	}
	if ((src = mpd_alloc(msize, sizeof *src)) == NULL) {
		goto error;
	}
	if ((dest = mpd_alloc(msize, sizeof *dest)) == NULL) {
		goto error;
	}

	for (i = 0; i < msize; i++) {
		a[i] = src[i] = random();
	}

	start_std= clock();
	std_trans_c(dest, src, rows, cols);
	end_std= clock();

	start_fast = clock();
	func(a, rows, cols);
	end_fast = clock();

	for (i = 0; i < msize; i++) {
		if (a[i] != dest[i]) {
			fprintf(stderr, "FAIL: a[%"PRI_mpd_size_t"] = %d\t"
                                        "dest[%"PRI_mpd_size_t"] = %d\n",
			        i, a[i], i, dest[i]);
			exit(1);
		}
	}

	fprintf(stderr, "size: %10"PRI_mpd_size_t"\tstd_trans: %6.2f sec\t "
                        "in_place_trans: %6.2f sec\n",
		msize,
		(double)(end_std-start_std)/(double)CLOCKS_PER_SEC,
		(double)(end_fast-start_fast)/(double)CLOCKS_PER_SEC);


out:
	if (a) __mingw_dfp_get_globals()->mpd_free(a);
	if (src) __mingw_dfp_get_globals()->mpd_free(src);
	if (dest) __mingw_dfp_get_globals()->mpd_free(dest);
	return ret;

error:
	ret = 0;
	goto out;
}


#define transfunc_uint (void (*)(mpd_uint_t *, mpd_size_t, mpd_size_t))
#define transfunc_uchar (void (*)(uint8_t *, mpd_size_t, mpd_size_t))


int main(void)
{
	mpd_size_t rows, cols;


	fprintf(stderr, "Running test_transpose ... \n");

	fprintf(stderr, "\n2^n * 2^n mpd_uint_t matrices:\n\n");
	for (rows = 1; mul_size_t(rows, rows) <= UMOD_ARRAY; rows *= 2) {
		if (!testit_uint(transfunc_uint transpose_pow2, rows, rows)) {
			break;
		}
	}

	fprintf(stderr, "\n2^n * 2*2^n mpd_uint_t matrices:\n\n");
	for (rows = 8, cols = 16; mul_size_t(rows, cols) <= UMOD_ARRAY; rows *= 2, cols*=2) {
		if (!testit_uint(transfunc_uint transpose_pow2, rows, cols)) {
			break;
		}
	}

	fprintf(stderr, "\n2*2^n * 2^n mpd_uint_t matrices:\n\n");
	for (rows = 16, cols = 8; mul_size_t(rows, cols) <= UMOD_ARRAY; rows *= 2, cols*=2) {
		if (!testit_uint(transfunc_uint transpose_pow2, rows, cols)) {
			break;
		}
	}


	fprintf(stderr, "\n2^n * 2^n uint8_t matrices:\n\n");
	for (rows = 1; mul_size_t(rows, rows) <= UCHAR_ARRAY; rows *= 2) {
		if (!testit_uchar(transfunc_uchar transpose_pow2_c, rows, rows)) {
			break;
		}
	}

	fprintf(stderr, "\n2^n * 2*2^n uint8_t matrices:\n\n");
	for (rows = 8, cols = 16; mul_size_t(rows, cols) <= UCHAR_ARRAY; rows *= 2, cols*=2) {
		if (!testit_uchar(transfunc_uchar transpose_pow2_c, rows, cols)) {
			break;
		}
	}

	fprintf(stderr, "\n2*2^n * 2^n uint8_t matrices:\n\n");
	for (rows = 16, cols = 8; mul_size_t(rows, cols) <= UCHAR_ARRAY; rows *= 2, cols*=2) {
		if (!testit_uchar(transfunc_uchar transpose_pow2_c, rows, cols)) {
			break;
		}
	}

	fprintf(stderr, "\ntest_transpose: PASS\n\n");

	return 0;
}