/*      $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;
}