/*
 * semi.c: test implementations of mathlib seminumerical functions
 *
 * Copyright (c) 1999-2019, Arm Limited.
 * SPDX-License-Identifier: MIT
 */

#include <stdio.h>
#include "semi.h"

static void test_rint(uint32 *in, uint32 *out,
                       int isfloor, int isceil) {
    int sign = in[0] & 0x80000000;
    int roundup = (isfloor && sign) || (isceil && !sign);
    uint32 xh, xl, roundword;
    int ex = (in[0] >> 20) & 0x7FF;    /* exponent */
    int i;

    if ((ex > 0x3ff + 52 - 1) ||     /* things this big can't be fractional */
        ((in[0] & 0x7FFFFFFF) == 0 && in[1] == 0)) {   /* zero */
        /* NaN, Inf, a large integer, or zero: just return the input */
        out[0] = in[0];
        out[1] = in[1];
        return;
    }

    /*
     * Special case: ex < 0x3ff, ie our number is in (0,1). Return
     * 1 or 0 according to roundup.
     */
    if (ex < 0x3ff) {
        out[0] = sign | (roundup ? 0x3FF00000 : 0);
        out[1] = 0;
        return;
    }

    /*
     * We're not short of time here, so we'll do this the hideously
     * inefficient way. Shift bit by bit so that the units place is
     * somewhere predictable, round, and shift back again.
     */
    xh = in[0];
    xl = in[1];
    roundword = 0;
    for (i = ex; i < 0x3ff + 52; i++) {
        if (roundword & 1)
            roundword |= 2;            /* preserve sticky bit */
        roundword = (roundword >> 1) | ((xl & 1) << 31);
        xl = (xl >> 1) | ((xh & 1) << 31);
        xh = xh >> 1;
    }
    if (roundword && roundup) {
        xl++;
        xh += (xl==0);
    }
    for (i = ex; i < 0x3ff + 52; i++) {
        xh = (xh << 1) | ((xl >> 31) & 1);
        xl = (xl & 0x7FFFFFFF) << 1;
    }
    out[0] = xh;
    out[1] = xl;
}

char *test_ceil(uint32 *in, uint32 *out) {
    test_rint(in, out, 0, 1);
    return NULL;
}

char *test_floor(uint32 *in, uint32 *out) {
    test_rint(in, out, 1, 0);
    return NULL;
}

static void test_rintf(uint32 *in, uint32 *out,
                       int isfloor, int isceil) {
    int sign = *in & 0x80000000;
    int roundup = (isfloor && sign) || (isceil && !sign);
    uint32 x, roundword;
    int ex = (*in >> 23) & 0xFF;       /* exponent */
    int i;

    if ((ex > 0x7f + 23 - 1) ||      /* things this big can't be fractional */
        (*in & 0x7FFFFFFF) == 0) {     /* zero */
        /* NaN, Inf, a large integer, or zero: just return the input */
        *out = *in;
        return;
    }

    /*
     * Special case: ex < 0x7f, ie our number is in (0,1). Return
     * 1 or 0 according to roundup.
     */
    if (ex < 0x7f) {
        *out = sign | (roundup ? 0x3F800000 : 0);
        return;
    }

    /*
     * We're not short of time here, so we'll do this the hideously
     * inefficient way. Shift bit by bit so that the units place is
     * somewhere predictable, round, and shift back again.
     */
    x = *in;
    roundword = 0;
    for (i = ex; i < 0x7F + 23; i++) {
        if (roundword & 1)
            roundword |= 2;            /* preserve sticky bit */
        roundword = (roundword >> 1) | ((x & 1) << 31);
        x = x >> 1;
    }
    if (roundword && roundup) {
        x++;
    }
    for (i = ex; i < 0x7F + 23; i++) {
        x = x << 1;
    }
    *out = x;
}

char *test_ceilf(uint32 *in, uint32 *out) {
    test_rintf(in, out, 0, 1);
    return NULL;
}

char *test_floorf(uint32 *in, uint32 *out) {
    test_rintf(in, out, 1, 0);
    return NULL;
}

char *test_fmod(uint32 *a, uint32 *b, uint32 *out) {
    int sign;
    int32 aex, bex;
    uint32 am[2], bm[2];

    if (((a[0] & 0x7FFFFFFF) << 1) + !!a[1] > 0xFFE00000 ||
        ((b[0] & 0x7FFFFFFF) << 1) + !!b[1] > 0xFFE00000) {
        /* a or b is NaN: return QNaN, optionally with IVO */
        uint32 an, bn;
        out[0] = 0x7ff80000;
        out[1] = 1;
        an = ((a[0] & 0x7FFFFFFF) << 1) + !!a[1];
        bn = ((b[0] & 0x7FFFFFFF) << 1) + !!b[1];
        if ((an > 0xFFE00000 && an < 0xFFF00000) ||
            (bn > 0xFFE00000 && bn < 0xFFF00000))
            return "i";                /* at least one SNaN: IVO */
        else
            return NULL;               /* no SNaNs, but at least 1 QNaN */
    }
    if ((b[0] & 0x7FFFFFFF) == 0 && b[1] == 0) {   /* b==0: EDOM */
        out[0] = 0x7ff80000;
        out[1] = 1;
        return "EDOM status=i";
    }
    if ((a[0] & 0x7FF00000) == 0x7FF00000) {   /* a==Inf: EDOM */
        out[0] = 0x7ff80000;
        out[1] = 1;
        return "EDOM status=i";
    }
    if ((b[0] & 0x7FF00000) == 0x7FF00000) {   /* b==Inf: return a */
        out[0] = a[0];
        out[1] = a[1];
        return NULL;
    }
    if ((a[0] & 0x7FFFFFFF) == 0 && a[1] == 0) {   /* a==0: return a */
        out[0] = a[0];
        out[1] = a[1];
        return NULL;
    }

    /*
     * OK. That's the special cases cleared out of the way. Now we
     * have finite (though not necessarily normal) a and b.
     */
    sign = a[0] & 0x80000000;          /* we discard sign of b */
    test_frexp(a, am, (uint32 *)&aex);
    test_frexp(b, bm, (uint32 *)&bex);
    am[0] &= 0xFFFFF, am[0] |= 0x100000;
    bm[0] &= 0xFFFFF, bm[0] |= 0x100000;

    while (aex >= bex) {
        if (am[0] > bm[0] || (am[0] == bm[0] && am[1] >= bm[1])) {
            am[1] -= bm[1];
            am[0] = am[0] - bm[0] - (am[1] > ~bm[1]);
        }
        if (aex > bex) {
            am[0] = (am[0] << 1) | ((am[1] & 0x80000000) >> 31);
            am[1] <<= 1;
            aex--;
        } else
            break;
    }

    /*
     * Renormalise final result; this can be cunningly done by
     * passing a denormal to ldexp.
     */
    aex += 0x3fd;
    am[0] |= sign;
    test_ldexp(am, (uint32 *)&aex, out);

    return NULL;                       /* FIXME */
}

char *test_fmodf(uint32 *a, uint32 *b, uint32 *out) {
    int sign;
    int32 aex, bex;
    uint32 am, bm;

    if ((*a & 0x7FFFFFFF) > 0x7F800000 ||
        (*b & 0x7FFFFFFF) > 0x7F800000) {
        /* a or b is NaN: return QNaN, optionally with IVO */
        uint32 an, bn;
        *out = 0x7fc00001;
        an = *a & 0x7FFFFFFF;
        bn = *b & 0x7FFFFFFF;
        if ((an > 0x7f800000 && an < 0x7fc00000) ||
            (bn > 0x7f800000 && bn < 0x7fc00000))
            return "i";                /* at least one SNaN: IVO */
        else
            return NULL;               /* no SNaNs, but at least 1 QNaN */
    }
    if ((*b & 0x7FFFFFFF) == 0) {      /* b==0: EDOM */
        *out = 0x7fc00001;
        return "EDOM status=i";
    }
    if ((*a & 0x7F800000) == 0x7F800000) {   /* a==Inf: EDOM */
        *out = 0x7fc00001;
        return "EDOM status=i";
    }
    if ((*b & 0x7F800000) == 0x7F800000) {   /* b==Inf: return a */
        *out = *a;
        return NULL;
    }
    if ((*a & 0x7FFFFFFF) == 0) {      /* a==0: return a */
        *out = *a;
        return NULL;
    }

    /*
     * OK. That's the special cases cleared out of the way. Now we
     * have finite (though not necessarily normal) a and b.
     */
    sign = a[0] & 0x80000000;          /* we discard sign of b */
    test_frexpf(a, &am, (uint32 *)&aex);
    test_frexpf(b, &bm, (uint32 *)&bex);
    am &= 0x7FFFFF, am |= 0x800000;
    bm &= 0x7FFFFF, bm |= 0x800000;

    while (aex >= bex) {
        if (am >= bm) {
            am -= bm;
        }
        if (aex > bex) {
            am <<= 1;
            aex--;
        } else
            break;
    }

    /*
     * Renormalise final result; this can be cunningly done by
     * passing a denormal to ldexp.
     */
    aex += 0x7d;
    am |= sign;
    test_ldexpf(&am, (uint32 *)&aex, out);

    return NULL;                       /* FIXME */
}

char *test_ldexp(uint32 *x, uint32 *np, uint32 *out) {
    int n = *np;
    int32 n2;
    uint32 y[2];
    int ex = (x[0] >> 20) & 0x7FF;     /* exponent */
    int sign = x[0] & 0x80000000;

    if (ex == 0x7FF) {                 /* inf/NaN; just return x */
        out[0] = x[0];
        out[1] = x[1];
        return NULL;
    }
    if ((x[0] & 0x7FFFFFFF) == 0 && x[1] == 0) {   /* zero: return x */
        out[0] = x[0];
        out[1] = x[1];
        return NULL;
    }

    test_frexp(x, y, (uint32 *)&n2);
    ex = n + n2;
    if (ex > 0x400) {                  /* overflow */
        out[0] = sign | 0x7FF00000;
        out[1] = 0;
        return "overflow";
    }
    /*
     * Underflow. 2^-1074 is 00000000.00000001; so if ex == -1074
     * then we have something [2^-1075,2^-1074). Under round-to-
     * nearest-even, this whole interval rounds up to 2^-1074,
     * except for the bottom endpoint which rounds to even and is
     * an underflow condition.
     *
     * So, ex < -1074 is definite underflow, and ex == -1074 is
     * underflow iff all mantissa bits are zero.
     */
    if (ex < -1074 || (ex == -1074 && (y[0] & 0xFFFFF) == 0 && y[1] == 0)) {
        out[0] = sign;                 /* underflow: correctly signed zero */
        out[1] = 0;
        return "underflow";
    }

    /*
     * No overflow or underflow; should be nice and simple, unless
     * we have to denormalise and round the result.
     */
    if (ex < -1021) {                  /* denormalise and round */
        uint32 roundword;
        y[0] &= 0x000FFFFF;
        y[0] |= 0x00100000;            /* set leading bit */
        roundword = 0;
        while (ex < -1021) {
            if (roundword & 1)
                roundword |= 2;        /* preserve sticky bit */
            roundword = (roundword >> 1) | ((y[1] & 1) << 31);
            y[1] = (y[1] >> 1) | ((y[0] & 1) << 31);
            y[0] = y[0] >> 1;
            ex++;
        }
        if (roundword > 0x80000000 ||  /* round up */
            (roundword == 0x80000000 && (y[1] & 1))) {  /* round up to even */
            y[1]++;
            y[0] += (y[1] == 0);
        }
        out[0] = sign | y[0];
        out[1] = y[1];
        /* Proper ERANGE underflow was handled earlier, but we still
         * expect an IEEE Underflow exception if this partially
         * underflowed result is not exact. */
        if (roundword)
            return "u";
        return NULL;                   /* underflow was handled earlier */
    } else {
        out[0] = y[0] + (ex << 20);
        out[1] = y[1];
        return NULL;
    }
}

char *test_ldexpf(uint32 *x, uint32 *np, uint32 *out) {
    int n = *np;
    int32 n2;
    uint32 y;
    int ex = (*x >> 23) & 0xFF;     /* exponent */
    int sign = *x & 0x80000000;

    if (ex == 0xFF) {                 /* inf/NaN; just return x */
        *out = *x;
        return NULL;
    }
    if ((*x & 0x7FFFFFFF) == 0) {      /* zero: return x */
        *out = *x;
        return NULL;
    }

    test_frexpf(x, &y, (uint32 *)&n2);
    ex = n + n2;
    if (ex > 0x80) {                  /* overflow */
        *out = sign | 0x7F800000;
        return "overflow";
    }
    /*
     * Underflow. 2^-149 is 00000001; so if ex == -149 then we have
     * something [2^-150,2^-149). Under round-to- nearest-even,
     * this whole interval rounds up to 2^-149, except for the
     * bottom endpoint which rounds to even and is an underflow
     * condition.
     *
     * So, ex < -149 is definite underflow, and ex == -149 is
     * underflow iff all mantissa bits are zero.
     */
    if (ex < -149 || (ex == -149 && (y & 0x7FFFFF) == 0)) {
        *out = sign;                 /* underflow: correctly signed zero */
        return "underflow";
    }

    /*
     * No overflow or underflow; should be nice and simple, unless
     * we have to denormalise and round the result.
     */
    if (ex < -125) {                  /* denormalise and round */
        uint32 roundword;
        y &= 0x007FFFFF;
        y |= 0x00800000;               /* set leading bit */
        roundword = 0;
        while (ex < -125) {
            if (roundword & 1)
                roundword |= 2;        /* preserve sticky bit */
            roundword = (roundword >> 1) | ((y & 1) << 31);
            y = y >> 1;
            ex++;
        }
        if (roundword > 0x80000000 ||  /* round up */
            (roundword == 0x80000000 && (y & 1))) {  /* round up to even */
            y++;
        }
        *out = sign | y;
        /* Proper ERANGE underflow was handled earlier, but we still
         * expect an IEEE Underflow exception if this partially
         * underflowed result is not exact. */
        if (roundword)
            return "u";
        return NULL;                   /* underflow was handled earlier */
    } else {
        *out = y + (ex << 23);
        return NULL;
    }
}

char *test_frexp(uint32 *x, uint32 *out, uint32 *nout) {
    int ex = (x[0] >> 20) & 0x7FF;     /* exponent */
    if (ex == 0x7FF) {                 /* inf/NaN; return x/0 */
        out[0] = x[0];
        out[1] = x[1];
        nout[0] = 0;
        return NULL;
    }
    if (ex == 0) {                     /* denormals/zeros */
        int sign;
        uint32 xh, xl;
        if ((x[0] & 0x7FFFFFFF) == 0 && x[1] == 0) {
            /* zero: return x/0 */
            out[0] = x[0];
            out[1] = x[1];
            nout[0] = 0;
            return NULL;
        }
        sign = x[0] & 0x80000000;
        xh = x[0] & 0x7FFFFFFF;
        xl = x[1];
        ex = 1;
        while (!(xh & 0x100000)) {
            ex--;
            xh = (xh << 1) | ((xl >> 31) & 1);
            xl = (xl & 0x7FFFFFFF) << 1;
        }
        out[0] = sign | 0x3FE00000 | (xh & 0xFFFFF);
        out[1] = xl;
        nout[0] = ex - 0x3FE;
        return NULL;
    }
    out[0] = 0x3FE00000 | (x[0] & 0x800FFFFF);
    out[1] = x[1];
    nout[0] = ex - 0x3FE;
    return NULL;                       /* ordinary number; no error */
}

char *test_frexpf(uint32 *x, uint32 *out, uint32 *nout) {
    int ex = (*x >> 23) & 0xFF;        /* exponent */
    if (ex == 0xFF) {                  /* inf/NaN; return x/0 */
        *out = *x;
        nout[0] = 0;
        return NULL;
    }
    if (ex == 0) {                     /* denormals/zeros */
        int sign;
        uint32 xv;
        if ((*x & 0x7FFFFFFF) == 0) {
            /* zero: return x/0 */
            *out = *x;
            nout[0] = 0;
            return NULL;
        }
        sign = *x & 0x80000000;
        xv = *x & 0x7FFFFFFF;
        ex = 1;
        while (!(xv & 0x800000)) {
            ex--;
            xv = xv << 1;
        }
        *out = sign | 0x3F000000 | (xv & 0x7FFFFF);
        nout[0] = ex - 0x7E;
        return NULL;
    }
    *out = 0x3F000000 | (*x & 0x807FFFFF);
    nout[0] = ex - 0x7E;
    return NULL;                       /* ordinary number; no error */
}

char *test_modf(uint32 *x, uint32 *fout, uint32 *iout) {
    int ex = (x[0] >> 20) & 0x7FF;     /* exponent */
    int sign = x[0] & 0x80000000;
    uint32 fh, fl;

    if (((x[0] & 0x7FFFFFFF) | (!!x[1])) > 0x7FF00000) {
        /*
         * NaN input: return the same in _both_ outputs.
         */
        fout[0] = iout[0] = x[0];
        fout[1] = iout[1] = x[1];
        return NULL;
    }

    test_rint(x, iout, 0, 0);
    fh = x[0] - iout[0];
    fl = x[1] - iout[1];
    if (!fh && !fl) {                  /* no fraction part */
        fout[0] = sign;
        fout[1] = 0;
        return NULL;
    }
    if (!(iout[0] & 0x7FFFFFFF) && !iout[1]) {   /* no integer part */
        fout[0] = x[0];
        fout[1] = x[1];
        return NULL;
    }
    while (!(fh & 0x100000)) {
        ex--;
        fh = (fh << 1) | ((fl >> 31) & 1);
        fl = (fl & 0x7FFFFFFF) << 1;
    }
    fout[0] = sign | (ex << 20) | (fh & 0xFFFFF);
    fout[1] = fl;
    return NULL;
}

char *test_modff(uint32 *x, uint32 *fout, uint32 *iout) {
    int ex = (*x >> 23) & 0xFF;        /* exponent */
    int sign = *x & 0x80000000;
    uint32 f;

    if ((*x & 0x7FFFFFFF) > 0x7F800000) {
        /*
         * NaN input: return the same in _both_ outputs.
         */
        *fout = *iout = *x;
        return NULL;
    }

    test_rintf(x, iout, 0, 0);
    f = *x - *iout;
    if (!f) {                          /* no fraction part */
        *fout = sign;
        return NULL;
    }
    if (!(*iout & 0x7FFFFFFF)) {       /* no integer part */
        *fout = *x;
        return NULL;
    }
    while (!(f & 0x800000)) {
        ex--;
        f = f << 1;
    }
    *fout = sign | (ex << 23) | (f & 0x7FFFFF);
    return NULL;
}

char *test_copysign(uint32 *x, uint32 *y, uint32 *out)
{
    int ysign = y[0] & 0x80000000;
    int xhigh = x[0] & 0x7fffffff;

    out[0] = ysign | xhigh;
    out[1] = x[1];

    /* There can be no error */
    return NULL;
}

char *test_copysignf(uint32 *x, uint32 *y, uint32 *out)
{
    int ysign = y[0] & 0x80000000;
    int xhigh = x[0] & 0x7fffffff;

    out[0] = ysign | xhigh;

    /* There can be no error */
    return NULL;
}

char *test_isfinite(uint32 *x, uint32 *out)
{
    int xhigh = x[0];
    /* Being finite means that the exponent is not 0x7ff */
    if ((xhigh & 0x7ff00000) == 0x7ff00000) out[0] = 0;
    else out[0] = 1;
    return NULL;
}

char *test_isfinitef(uint32 *x, uint32 *out)
{
    /* Being finite means that the exponent is not 0xff */
    if ((x[0] & 0x7f800000) == 0x7f800000) out[0] = 0;
    else out[0] = 1;
    return NULL;
}

char *test_isinff(uint32 *x, uint32 *out)
{
    /* Being infinite means that our bottom 30 bits equate to 0x7f800000 */
    if ((x[0] & 0x7fffffff) == 0x7f800000) out[0] = 1;
    else out[0] = 0;
    return NULL;
}

char *test_isinf(uint32 *x, uint32 *out)
{
    int xhigh = x[0];
    int xlow = x[1];
    /* Being infinite means that our fraction is zero and exponent is 0x7ff */
    if (((xhigh & 0x7fffffff) == 0x7ff00000) && (xlow == 0)) out[0] = 1;
    else out[0] = 0;
    return NULL;
}

char *test_isnanf(uint32 *x, uint32 *out)
{
    /* Being NaN means that our exponent is 0xff and non-0 fraction */
    int exponent = x[0] & 0x7f800000;
    int fraction = x[0] & 0x007fffff;
    if ((exponent == 0x7f800000) && (fraction != 0)) out[0] = 1;
    else out[0] = 0;
    return NULL;
}

char *test_isnan(uint32 *x, uint32 *out)
{
    /* Being NaN means that our exponent is 0x7ff and non-0 fraction */
    int exponent = x[0] & 0x7ff00000;
    int fractionhigh = x[0] & 0x000fffff;
    if ((exponent == 0x7ff00000) && ((fractionhigh != 0) || x[1] != 0))
        out[0] = 1;
    else out[0] = 0;
    return NULL;
}

char *test_isnormalf(uint32 *x, uint32 *out)
{
    /* Being normal means exponent is not 0 and is not 0xff */
    int exponent = x[0] & 0x7f800000;
    if (exponent == 0x7f800000) out[0] = 0;
    else if (exponent == 0) out[0] = 0;
    else out[0] = 1;
    return NULL;
}

char *test_isnormal(uint32 *x, uint32 *out)
{
    /* Being normal means exponent is not 0 and is not 0x7ff */
    int exponent = x[0] & 0x7ff00000;
    if (exponent == 0x7ff00000) out[0] = 0;
    else if (exponent == 0) out[0] = 0;
    else out[0] = 1;
    return NULL;
}

char *test_signbitf(uint32 *x, uint32 *out)
{
    /* Sign bit is bit 31 */
    out[0] = (x[0] >> 31) & 1;
    return NULL;
}

char *test_signbit(uint32 *x, uint32 *out)
{
    /* Sign bit is bit 31 */
    out[0] = (x[0] >> 31) & 1;
    return NULL;
}

char *test_fpclassify(uint32 *x, uint32 *out)
{
    int exponent = (x[0] & 0x7ff00000) >> 20;
    int fraction = (x[0] & 0x000fffff) | x[1];

    if ((exponent == 0x00) && (fraction == 0)) out[0] = 0;
    else if ((exponent == 0x00) && (fraction != 0)) out[0] = 4;
    else if ((exponent == 0x7ff) && (fraction == 0)) out[0] = 3;
    else if ((exponent == 0x7ff) && (fraction != 0)) out[0] = 7;
    else out[0] = 5;
    return NULL;
}

char *test_fpclassifyf(uint32 *x, uint32 *out)
{
    int exponent = (x[0] & 0x7f800000) >> 23;
    int fraction = x[0] & 0x007fffff;

    if ((exponent == 0x000) && (fraction == 0)) out[0] = 0;
    else if ((exponent == 0x000) && (fraction != 0)) out[0] = 4;
    else if ((exponent == 0xff) && (fraction == 0)) out[0] = 3;
    else if ((exponent == 0xff) && (fraction != 0)) out[0] = 7;
    else out[0] = 5;
    return NULL;
}

/*
 * Internal function that compares doubles in x & y and returns -3, -2, -1, 0,
 * 1 if they compare to be signaling, unordered, less than, equal or greater
 * than.
 */
static int fpcmp4(uint32 *x, uint32 *y)
{
    int result = 0;

    /*
     * Sort out whether results are ordered or not to begin with
     * NaNs have exponent 0x7ff, and non-zero fraction. Signaling NaNs take
     * higher priority than quiet ones.
     */
    if ((x[0] & 0x7fffffff) >= 0x7ff80000) result = -2;
    else if ((x[0] & 0x7fffffff) > 0x7ff00000) result = -3;
    else if (((x[0] & 0x7fffffff) == 0x7ff00000) && (x[1] != 0)) result = -3;
    if ((y[0] & 0x7fffffff) >= 0x7ff80000 && result != -3) result = -2;
    else if ((y[0] & 0x7fffffff) > 0x7ff00000) result = -3;
    else if (((y[0] & 0x7fffffff) == 0x7ff00000) && (y[1] != 0)) result = -3;
    if (result != 0) return result;

    /*
     * The two forms of zero are equal
     */
    if (((x[0] & 0x7fffffff) == 0) && x[1] == 0 &&
        ((y[0] & 0x7fffffff) == 0) && y[1] == 0)
        return 0;

    /*
     * If x and y have different signs we can tell that they're not equal
     * If x is +ve we have x > y return 1 - otherwise y is +ve return -1
     */
    if ((x[0] >> 31) != (y[0] >> 31))
        return ((x[0] >> 31) == 0) - ((y[0] >> 31) == 0);

    /*
     * Now we have both signs the same, let's do an initial compare of the
     * values.
     *
     * Whoever designed IEEE754's floating point formats is very clever and
     * earns my undying admiration.  Once you remove the sign-bit, the
     * floating point numbers can be ordered using the standard <, ==, >
     * operators will treating the fp-numbers as integers with that bit-
     * pattern.
     */
    if ((x[0] & 0x7fffffff) < (y[0] & 0x7fffffff)) result = -1;
    else if ((x[0] & 0x7fffffff) > (y[0] & 0x7fffffff)) result = 1;
    else if (x[1] < y[1]) result = -1;
    else if (x[1] > y[1]) result = 1;
    else result = 0;

    /*
     * Now we return the result - is x is positive (and therefore so is y) we
     * return the plain result - otherwise we negate it and return.
     */
    if ((x[0] >> 31) == 0) return result;
    else return -result;
}

/*
 * Internal function that compares floats in x & y and returns -3, -2, -1, 0,
 * 1 if they compare to be signaling, unordered, less than, equal or greater
 * than.
 */
static int fpcmp4f(uint32 *x, uint32 *y)
{
    int result = 0;

    /*
     * Sort out whether results are ordered or not to begin with
     * NaNs have exponent 0xff, and non-zero fraction - we have to handle all
     * signaling cases over the quiet ones
     */
    if ((x[0] & 0x7fffffff) >= 0x7fc00000) result = -2;
    else if ((x[0] & 0x7fffffff) > 0x7f800000) result = -3;
    if ((y[0] & 0x7fffffff) >= 0x7fc00000 && result != -3) result = -2;
    else if ((y[0] & 0x7fffffff) > 0x7f800000) result = -3;
    if (result != 0) return result;

    /*
     * The two forms of zero are equal
     */
    if (((x[0] & 0x7fffffff) == 0) && ((y[0] & 0x7fffffff) == 0))
        return 0;

    /*
     * If x and y have different signs we can tell that they're not equal
     * If x is +ve we have x > y return 1 - otherwise y is +ve return -1
     */
    if ((x[0] >> 31) != (y[0] >> 31))
        return ((x[0] >> 31) == 0) - ((y[0] >> 31) == 0);

    /*
     * Now we have both signs the same, let's do an initial compare of the
     * values.
     *
     * Whoever designed IEEE754's floating point formats is very clever and
     * earns my undying admiration.  Once you remove the sign-bit, the
     * floating point numbers can be ordered using the standard <, ==, >
     * operators will treating the fp-numbers as integers with that bit-
     * pattern.
     */
    if ((x[0] & 0x7fffffff) < (y[0] & 0x7fffffff)) result = -1;
    else if ((x[0] & 0x7fffffff) > (y[0] & 0x7fffffff)) result = 1;
    else result = 0;

    /*
     * Now we return the result - is x is positive (and therefore so is y) we
     * return the plain result - otherwise we negate it and return.
     */
    if ((x[0] >> 31) == 0) return result;
    else return -result;
}

char *test_isgreater(uint32 *x, uint32 *y, uint32 *out)
{
    int result = fpcmp4(x, y);
    *out = (result == 1);
    return result == -3 ? "i" : NULL;
}

char *test_isgreaterequal(uint32 *x, uint32 *y, uint32 *out)
{
    int result = fpcmp4(x, y);
    *out = (result >= 0);
    return result == -3 ? "i" : NULL;
}

char *test_isless(uint32 *x, uint32 *y, uint32 *out)
{
    int result = fpcmp4(x, y);
    *out = (result == -1);
    return result == -3 ? "i" : NULL;
}

char *test_islessequal(uint32 *x, uint32 *y, uint32 *out)
{
    int result = fpcmp4(x, y);
    *out = (result == -1) || (result == 0);
    return result == -3 ? "i" : NULL;
}

char *test_islessgreater(uint32 *x, uint32 *y, uint32 *out)
{
    int result = fpcmp4(x, y);
    *out = (result == -1) || (result == 1);
    return result == -3 ? "i" : NULL;
}

char *test_isunordered(uint32 *x, uint32 *y, uint32 *out)
{
    int normal = 0;
    int result = fpcmp4(x, y);

    test_isnormal(x, out);
    normal |= *out;
    test_isnormal(y, out);
    normal |= *out;
    *out = (result == -2) || (result == -3);
    return result == -3 ? "i" : NULL;
}

char *test_isgreaterf(uint32 *x, uint32 *y, uint32 *out)
{
    int result = fpcmp4f(x, y);
    *out = (result == 1);
    return result == -3 ? "i" : NULL;
}

char *test_isgreaterequalf(uint32 *x, uint32 *y, uint32 *out)
{
    int result = fpcmp4f(x, y);
    *out = (result >= 0);
    return result == -3 ? "i" : NULL;
}

char *test_islessf(uint32 *x, uint32 *y, uint32 *out)
{
    int result = fpcmp4f(x, y);
    *out = (result == -1);
    return result == -3 ? "i" : NULL;
}

char *test_islessequalf(uint32 *x, uint32 *y, uint32 *out)
{
    int result = fpcmp4f(x, y);
    *out = (result == -1) || (result == 0);
    return result == -3 ? "i" : NULL;
}

char *test_islessgreaterf(uint32 *x, uint32 *y, uint32 *out)
{
    int result = fpcmp4f(x, y);
    *out = (result == -1) || (result == 1);
    return result == -3 ? "i" : NULL;
}

char *test_isunorderedf(uint32 *x, uint32 *y, uint32 *out)
{
    int normal = 0;
    int result = fpcmp4f(x, y);

    test_isnormalf(x, out);
    normal |= *out;
    test_isnormalf(y, out);
    normal |= *out;
    *out = (result == -2) || (result == -3);
    return result == -3 ? "i" : NULL;
}