virtual_mul.c

/*-
 * Copyright (c) 2017 Hans Petter Selasky
 *
 * 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 <stdint.h>
#include <string.h>

#include <sys/queue.h>

#include "virtual_int.h"

#ifndef VOSS_X3_LOG2_COMBA
#define	VOSS_X3_LOG2_COMBA 5
#endif

#if (VOSS_X3_LOG2_COMBA < 2)
#error "VOSS_X3_LOG2_COMBA must be greater than 1"
#endif

struct voss_x3_input_double {
	double	a;
	double	b;
} __aligned(16);

/*
 * <input size> = "stride"
 * <output size> = 2 * "stride"
 */
static void
voss_x3_multiply_sub_double(struct voss_x3_input_double *input, double *ptr_low, double *ptr_high,
    const size_t stride, const uint8_t toggle)
{
	size_t x;
	size_t y;

	if (stride >= (1UL << VOSS_X3_LOG2_COMBA)) {
		const size_t strideh = stride >> 1;

		if (toggle) {

			/* inverse step */
			for (x = 0; x != strideh; x++) {
				double a, b, c, d;

				a = ptr_low[x];
				b = ptr_low[x + strideh];
				c = ptr_high[x];
				d = ptr_high[x + strideh];

				ptr_low[x + strideh] = a + b;
				ptr_high[x] = a + b + c + d;
			}

			voss_x3_multiply_sub_double(input, ptr_low, ptr_low + strideh, strideh, 1);

			for (x = 0; x != strideh; x++)
				ptr_low[x + strideh] = -ptr_low[x + strideh];

			voss_x3_multiply_sub_double(input + strideh, ptr_low + strideh, ptr_high + strideh, strideh, 1);

			/* forward step */
			for (x = 0; x != strideh; x++) {
				double a, b, c, d;

				a = ptr_low[x];
				b = ptr_low[x + strideh];
				c = ptr_high[x];
				d = ptr_high[x + strideh];

				ptr_low[x + strideh] = -a - b;
				ptr_high[x] = c + b - d;

				input[x + strideh].a += input[x].a;
				input[x + strideh].b += input[x].b;
			}

			voss_x3_multiply_sub_double(input + strideh, ptr_low + strideh, ptr_high, strideh, 0);
		} else {
			voss_x3_multiply_sub_double(input + strideh, ptr_low + strideh, ptr_high, strideh, 1);

			/* inverse step */
			for (x = 0; x != strideh; x++) {
				double a, b, c, d;

				a = ptr_low[x];
				b = ptr_low[x + strideh];
				c = ptr_high[x];
				d = ptr_high[x + strideh];

				ptr_low[x + strideh] = -a - b;
				ptr_high[x] = a + b + c + d;

				input[x + strideh].a -= input[x].a;
				input[x + strideh].b -= input[x].b;
			}

			voss_x3_multiply_sub_double(input + strideh, ptr_low + strideh, ptr_high + strideh, strideh, 0);

			for (x = 0; x != strideh; x++)
				ptr_low[x + strideh] = -ptr_low[x + strideh];

			voss_x3_multiply_sub_double(input, ptr_low, ptr_low + strideh, strideh, 0);

			/* forward step */
			for (x = 0; x != strideh; x++) {
				double a, b, c, d;

				a = ptr_low[x];
				b = ptr_low[x + strideh];
				c = ptr_high[x];
				d = ptr_high[x + strideh];

				ptr_low[x + strideh] = b - a;
				ptr_high[x] = c - b - d;
			}
		}
	} else {
		for (x = 0; x != stride; x++) {
			double value = input[x].a;

			for (y = 0; y != (stride - x); y++) {
				ptr_low[x + y] += input[y].b * value;
			}

			for (; y != stride; y++) {
				ptr_high[x + y - stride] += input[y].b * value;
			}
		}
	}
}

/*
 * <input size> = "max"
 * <output size> = 2 * "max"
 */
void
voss_x3_multiply_double(const int64_t *va, const double *vb, double *pc, const size_t max)
{
	struct voss_x3_input_double input[max];
	size_t x;

	/* check for non-power of two */
	if (max & (max - 1))
		return;

	/* setup input vector */
	for (x = 0; x != max; x++) {
		input[x].a = va[x];
		input[x].b = vb[x];
	}

	/* do multiplication */
	voss_x3_multiply_sub_double(input, pc, pc + max, max, 1);
}