/*	$OpenBSD: vecast.c,v 1.1 2022/10/22 17:50:28 gkoehler Exp $	*/

/*
 * Copyright (c) 2022 George Koehler <gkoehler@openbsd.org>
 *
 * Permission to use, copy, modify, and distribute this software for any
 * purpose with or without fee is hereby granted, provided that the above
 * copyright notice and this permission notice appear in all copies.
 *
 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
 */

#include <altivec.h>
#include <err.h>
#include <math.h>
#include <stdio.h>

struct double4 {
	double		d[4];
};

union vu {
	vector float	vf;
	vector int	vi;
	vector unsigned	vu;
	float		f[4];
	int		i[4];
	unsigned	u[4];
};

#define AD(a, b, c, d)	(struct double4){a, b, c, d}
#define VF(a, b, c, d)	(vector float)(a, b, c, d)
#define VI(a, b, c, d)	(vector int)(a, b, c, d)
#define VU(a, b, c, d)	(vector unsigned)(a, b, c, d)
#define rsqrt(f)	(1.0 / sqrt(f))

int fail;

void
ck_equal(const char *what, vector float out, vector float answer)
{
	if (vec_any_ne(out, answer)) {
		union vu a, b;

		a.vf = out;
		b.vf = answer;
		warnx("%s: {%a, %a, %a, %a} should be {%a, %a, %a, %a}",
		    what, a.f[0], a.f[1], a.f[2], a.f[3],
		    b.f[0], b.f[1], b.f[2], b.f[3]);
		fail = 1;
	}
}

void
ck_equal_i(const char *what, vector int out, vector int answer)
{
	if (vec_any_ne(out, answer)) {
		union vu a, b;

		a.vi = out;
		b.vi = answer;
		warnx("%s: {%d, %d, %d, %d} should be {%d, %d, %d, %d}",
		    what, a.i[0], a.i[1], a.i[2], a.i[3],
		    b.i[0], b.i[1], b.i[2], b.i[3]);
		fail = 1;
	}
}

void
ck_equal_u(const char *what, vector unsigned out, vector unsigned answer)
{
	if (vec_any_ne(out, answer)) {
		union vu a, b;

		a.vi = out;
		b.vi = answer;
		warnx("%s: {%u, %u, %u, %u} should be {%u, %u, %u, %u}",
		    what, a.u[0], a.u[1], a.u[2], a.u[3],
		    b.u[0], b.u[1], b.u[2], b.u[3]);
		fail = 1;
	}
}

enum error_check {REL_1_IN, ABS_1_IN};

/* Checks that error is at most 1 in err_den. */
void
ck_estimate(const char *what, vector float out, struct double4 answer,
    enum error_check which, double err_den)
{
	union vu u;
	int i, warned = 0;

	u.vf = out;
	for (i = 0; i < 4; i++) {
		double estimate = u.f[i];
		double target = answer.d[i];
		double error;

		switch (which) {
		case REL_1_IN:	/* relative error */
			error = fabs(target / (estimate - target));
			break;
		case ABS_1_IN:	/* absolute error */
			error = fabs(1 / (estimate - target));
			break;
		default:
			errx(1, "invalid check");
		}

		if (error < err_den) {
			if (!warned) {
				warnx("%s: {%a, %a, %a, %a} should be "
				    "near {%a, %a, %a, %a} (1/%g)",
				    what, u.f[0], u.f[1], u.f[2], u.f[3],
				    answer.d[0], answer.d[1], answer.d[2],
				    answer.d[3], err_den);
				warned = 1;
				fail = 1;
			}
			warnx("%a is off %a by 1/%g", estimate, target,
			    error);
		}
	}
}

/*
 * Tries altivec with denormal or subnormal floats.
 * These are single-precision floats f, where
 *   0 < |f| < 2**-126 = 0x1p-126 = 0x10p-130 = 1.17549435E-38F
 */
int
main(void)
{
	struct double4 dan;
	volatile vector float in1, in2, in3;
	vector float ans;
	vector int ian;
	vector unsigned uan;

	/* in1 + in2 */
	in1 = VF(10, 0x10p-140, 0x20p-130, -0x2000p-134);
	in2 = VF( 4,  0x5p-140, -0x1p-130,  0x1fffp-134);
	ans = VF(14, 0x15p-140, 0x1fp-130,    -0x1p-134);
	ck_equal("vec_add", vec_add(in1, in2), ans);

	/* in1 - in2 */
	in1 = VF(0x4000p-134, 10, 0x10p-140,   0x3p-130);
	in2 = VF(0x3ffep-134,  4,  0x5p-140,  0x40p-130);
	ans = VF(   0x2p-134,  6,  0xbp-140, -0x3dp-130);
	ck_equal("vec_sub", vec_sub(in1, in2), ans);

	/* in1 * in2 + in3 */
	in1 = VF( 0x6p-70,   0x6p-140, 6,   0x6p-100);
	in2 = VF( 0x7p-70,   0x7p50,   7,   0x7p-30);
	in3 = VF(  0,          0,      1, -0x20p-130);
	ans = VF(0x2ap-140, 0x2ap-90, 43,   0xap-130);
	ck_equal("vec_madd", vec_madd(in1, in2, in3), ans);

	/* in3 - in1 * in2 */
	in1 = VF( 0xbp-30,    0xbp-70,   0xbp44,  11);
	in2 = VF( 0x3p-100,   0x3p-70,  -0x3p-138, 3);
	in3 = VF(0x25p-130,     0,         0,     35);
	ans = VF( 0x4p-130, -0x21p-140, 0x21p-94,  2);
	ck_equal("vec_nmsub", vec_nmsub(in1, in2, in3), ans);

	/* 1 / in1 */
	in1 = VF(      3,       0x3p126,       0x3p-126, 0x1p127);
	dan = AD(1.0 / 3, 1.0 / 0x3p126, 1.0 / 0x3p-126, 0x1p-127);
	ck_estimate("vec_re", vec_re(in1), dan, REL_1_IN, 4096);

	/* 1 / sqrt(in1) */
	in1 = VF(1,       2,        0x1p-128,        0x5p-135);
	dan = AD(1, rsqrt(2), rsqrt(0x1p-128), rsqrt(0x5p-135));
	ck_estimate("vec_rsqrt", vec_rsqrte(in1), dan, REL_1_IN, 4096);

	/* log2(in1) */
	in1 = VF(0x1p-130, 0x1p-149, 32, 0x1p-10);
	dan = AD(    -130,     -149,  5,     -10);
	ck_estimate("vec_loge", vec_loge(in1), dan, ABS_1_IN, 32);
	in1 = VF(     0x123p-139,       0xabcp-145,  1, 1);
	dan = AD(log2(0x123p-139), log2(0xabcp-145), 0, 0);
	ck_estimate("vec_loge", vec_loge(in1), dan, ABS_1_IN, 32);

	/* 2**in1 */
	in1 = VF(    -149,     -138,     -127,   10);
	ans = VF(0x1p-149, 0x1p-138, 0x1p-127, 1024);
	ck_equal("vec_expte", vec_expte(in1), ans);
	in1 = VF(    -10,      -145.3,       -136.9,       -127.1);
	dan = AD(0x1p-10, exp2(-145.3), exp2(-136.9), exp2(-127.1));
	ck_estimate("vec_expte", vec_expte(in1), dan, REL_1_IN, 16);

	/* (int)(in1 * 2**exponent) */
	in1 = VF(0x1p-127, 2.34, -0xfedp-140, -19.8);
	ian = VI(  0,      2,         0,      -19);
	ck_equal_i("vec_cts", vec_cts(in1, 0), ian);
	in1 = VF(0x1p-113, -1, -0xabcp-143, 0x1fp-10);
	ian = VI(  0,   -1024,      0,      0x1f);
	ck_equal_i("vec_cts", vec_cts(in1, 10), ian);

	/* (unsigned)(in1 * 2**exponent) */
	in1 = VF(0x1.ap-130, 0x1.ep-140, 24000012, 0);
	uan = VU(  0,          0,      3072001536, 0);
	ck_equal_u("vec_ctu", vec_ctu(in1, 7), uan);

	return fail;
}