/*
* semi.c: test implementations of mathlib seminumerical functions
*
* Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
* See https://llvm.org/LICENSE.txt for license information.
* SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
*/
#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;
}