aboutsummaryrefslogtreecommitdiff
path: root/math/test/rtest/semi.c
diff options
context:
space:
mode:
Diffstat (limited to 'math/test/rtest/semi.c')
-rw-r--r--math/test/rtest/semi.c905
1 files changed, 905 insertions, 0 deletions
diff --git a/math/test/rtest/semi.c b/math/test/rtest/semi.c
new file mode 100644
index 000000000000..c9f0daf76508
--- /dev/null
+++ b/math/test/rtest/semi.c
@@ -0,0 +1,905 @@
+/*
+ * 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;
+}
generated by cgit v1.3 (git 2.53.0) at 2026-04-20 06:36:37 +0000