907 lines · c
1/*2 * semi.c: test implementations of mathlib seminumerical functions3 *4 * Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.5 * See https://llvm.org/LICENSE.txt for license information.6 * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception7 */8 9#include <stdio.h>10#include "semi.h"11 12static void test_rint(uint32 *in, uint32 *out,13 int isfloor, int isceil) {14 int sign = in[0] & 0x80000000;15 int roundup = (isfloor && sign) || (isceil && !sign);16 uint32 xh, xl, roundword;17 int ex = (in[0] >> 20) & 0x7FF; /* exponent */18 int i;19 20 if ((ex > 0x3ff + 52 - 1) || /* things this big can't be fractional */21 ((in[0] & 0x7FFFFFFF) == 0 && in[1] == 0)) { /* zero */22 /* NaN, Inf, a large integer, or zero: just return the input */23 out[0] = in[0];24 out[1] = in[1];25 return;26 }27 28 /*29 * Special case: ex < 0x3ff, ie our number is in (0,1). Return30 * 1 or 0 according to roundup.31 */32 if (ex < 0x3ff) {33 out[0] = sign | (roundup ? 0x3FF00000 : 0);34 out[1] = 0;35 return;36 }37 38 /*39 * We're not short of time here, so we'll do this the hideously40 * inefficient way. Shift bit by bit so that the units place is41 * somewhere predictable, round, and shift back again.42 */43 xh = in[0];44 xl = in[1];45 roundword = 0;46 for (i = ex; i < 0x3ff + 52; i++) {47 if (roundword & 1)48 roundword |= 2; /* preserve sticky bit */49 roundword = (roundword >> 1) | ((xl & 1) << 31);50 xl = (xl >> 1) | ((xh & 1) << 31);51 xh = xh >> 1;52 }53 if (roundword && roundup) {54 xl++;55 xh += (xl==0);56 }57 for (i = ex; i < 0x3ff + 52; i++) {58 xh = (xh << 1) | ((xl >> 31) & 1);59 xl = (xl & 0x7FFFFFFF) << 1;60 }61 out[0] = xh;62 out[1] = xl;63}64 65char *test_ceil(uint32 *in, uint32 *out) {66 test_rint(in, out, 0, 1);67 return NULL;68}69 70char *test_floor(uint32 *in, uint32 *out) {71 test_rint(in, out, 1, 0);72 return NULL;73}74 75static void test_rintf(uint32 *in, uint32 *out,76 int isfloor, int isceil) {77 int sign = *in & 0x80000000;78 int roundup = (isfloor && sign) || (isceil && !sign);79 uint32 x, roundword;80 int ex = (*in >> 23) & 0xFF; /* exponent */81 int i;82 83 if ((ex > 0x7f + 23 - 1) || /* things this big can't be fractional */84 (*in & 0x7FFFFFFF) == 0) { /* zero */85 /* NaN, Inf, a large integer, or zero: just return the input */86 *out = *in;87 return;88 }89 90 /*91 * Special case: ex < 0x7f, ie our number is in (0,1). Return92 * 1 or 0 according to roundup.93 */94 if (ex < 0x7f) {95 *out = sign | (roundup ? 0x3F800000 : 0);96 return;97 }98 99 /*100 * We're not short of time here, so we'll do this the hideously101 * inefficient way. Shift bit by bit so that the units place is102 * somewhere predictable, round, and shift back again.103 */104 x = *in;105 roundword = 0;106 for (i = ex; i < 0x7F + 23; i++) {107 if (roundword & 1)108 roundword |= 2; /* preserve sticky bit */109 roundword = (roundword >> 1) | ((x & 1) << 31);110 x = x >> 1;111 }112 if (roundword && roundup) {113 x++;114 }115 for (i = ex; i < 0x7F + 23; i++) {116 x = x << 1;117 }118 *out = x;119}120 121char *test_ceilf(uint32 *in, uint32 *out) {122 test_rintf(in, out, 0, 1);123 return NULL;124}125 126char *test_floorf(uint32 *in, uint32 *out) {127 test_rintf(in, out, 1, 0);128 return NULL;129}130 131char *test_fmod(uint32 *a, uint32 *b, uint32 *out) {132 int sign;133 int32 aex, bex;134 uint32 am[2], bm[2];135 136 if (((a[0] & 0x7FFFFFFF) << 1) + !!a[1] > 0xFFE00000 ||137 ((b[0] & 0x7FFFFFFF) << 1) + !!b[1] > 0xFFE00000) {138 /* a or b is NaN: return QNaN, optionally with IVO */139 uint32 an, bn;140 out[0] = 0x7ff80000;141 out[1] = 1;142 an = ((a[0] & 0x7FFFFFFF) << 1) + !!a[1];143 bn = ((b[0] & 0x7FFFFFFF) << 1) + !!b[1];144 if ((an > 0xFFE00000 && an < 0xFFF00000) ||145 (bn > 0xFFE00000 && bn < 0xFFF00000))146 return "i"; /* at least one SNaN: IVO */147 else148 return NULL; /* no SNaNs, but at least 1 QNaN */149 }150 if ((b[0] & 0x7FFFFFFF) == 0 && b[1] == 0) { /* b==0: EDOM */151 out[0] = 0x7ff80000;152 out[1] = 1;153 return "EDOM status=i";154 }155 if ((a[0] & 0x7FF00000) == 0x7FF00000) { /* a==Inf: EDOM */156 out[0] = 0x7ff80000;157 out[1] = 1;158 return "EDOM status=i";159 }160 if ((b[0] & 0x7FF00000) == 0x7FF00000) { /* b==Inf: return a */161 out[0] = a[0];162 out[1] = a[1];163 return NULL;164 }165 if ((a[0] & 0x7FFFFFFF) == 0 && a[1] == 0) { /* a==0: return a */166 out[0] = a[0];167 out[1] = a[1];168 return NULL;169 }170 171 /*172 * OK. That's the special cases cleared out of the way. Now we173 * have finite (though not necessarily normal) a and b.174 */175 sign = a[0] & 0x80000000; /* we discard sign of b */176 test_frexp(a, am, (uint32 *)&aex);177 test_frexp(b, bm, (uint32 *)&bex);178 am[0] &= 0xFFFFF, am[0] |= 0x100000;179 bm[0] &= 0xFFFFF, bm[0] |= 0x100000;180 181 while (aex >= bex) {182 if (am[0] > bm[0] || (am[0] == bm[0] && am[1] >= bm[1])) {183 am[1] -= bm[1];184 am[0] = am[0] - bm[0] - (am[1] > ~bm[1]);185 }186 if (aex > bex) {187 am[0] = (am[0] << 1) | ((am[1] & 0x80000000) >> 31);188 am[1] <<= 1;189 aex--;190 } else191 break;192 }193 194 /*195 * Renormalise final result; this can be cunningly done by196 * passing a denormal to ldexp.197 */198 aex += 0x3fd;199 am[0] |= sign;200 test_ldexp(am, (uint32 *)&aex, out);201 202 return NULL; /* FIXME */203}204 205char *test_fmodf(uint32 *a, uint32 *b, uint32 *out) {206 int sign;207 int32 aex, bex;208 uint32 am, bm;209 210 if ((*a & 0x7FFFFFFF) > 0x7F800000 ||211 (*b & 0x7FFFFFFF) > 0x7F800000) {212 /* a or b is NaN: return QNaN, optionally with IVO */213 uint32 an, bn;214 *out = 0x7fc00001;215 an = *a & 0x7FFFFFFF;216 bn = *b & 0x7FFFFFFF;217 if ((an > 0x7f800000 && an < 0x7fc00000) ||218 (bn > 0x7f800000 && bn < 0x7fc00000))219 return "i"; /* at least one SNaN: IVO */220 else221 return NULL; /* no SNaNs, but at least 1 QNaN */222 }223 if ((*b & 0x7FFFFFFF) == 0) { /* b==0: EDOM */224 *out = 0x7fc00001;225 return "EDOM status=i";226 }227 if ((*a & 0x7F800000) == 0x7F800000) { /* a==Inf: EDOM */228 *out = 0x7fc00001;229 return "EDOM status=i";230 }231 if ((*b & 0x7F800000) == 0x7F800000) { /* b==Inf: return a */232 *out = *a;233 return NULL;234 }235 if ((*a & 0x7FFFFFFF) == 0) { /* a==0: return a */236 *out = *a;237 return NULL;238 }239 240 /*241 * OK. That's the special cases cleared out of the way. Now we242 * have finite (though not necessarily normal) a and b.243 */244 sign = a[0] & 0x80000000; /* we discard sign of b */245 test_frexpf(a, &am, (uint32 *)&aex);246 test_frexpf(b, &bm, (uint32 *)&bex);247 am &= 0x7FFFFF, am |= 0x800000;248 bm &= 0x7FFFFF, bm |= 0x800000;249 250 while (aex >= bex) {251 if (am >= bm) {252 am -= bm;253 }254 if (aex > bex) {255 am <<= 1;256 aex--;257 } else258 break;259 }260 261 /*262 * Renormalise final result; this can be cunningly done by263 * passing a denormal to ldexp.264 */265 aex += 0x7d;266 am |= sign;267 test_ldexpf(&am, (uint32 *)&aex, out);268 269 return NULL; /* FIXME */270}271 272char *test_ldexp(uint32 *x, uint32 *np, uint32 *out) {273 int n = *np;274 int32 n2;275 uint32 y[2];276 int ex = (x[0] >> 20) & 0x7FF; /* exponent */277 int sign = x[0] & 0x80000000;278 279 if (ex == 0x7FF) { /* inf/NaN; just return x */280 out[0] = x[0];281 out[1] = x[1];282 return NULL;283 }284 if ((x[0] & 0x7FFFFFFF) == 0 && x[1] == 0) { /* zero: return x */285 out[0] = x[0];286 out[1] = x[1];287 return NULL;288 }289 290 test_frexp(x, y, (uint32 *)&n2);291 ex = n + n2;292 if (ex > 0x400) { /* overflow */293 out[0] = sign | 0x7FF00000;294 out[1] = 0;295 return "overflow";296 }297 /*298 * Underflow. 2^-1074 is 00000000.00000001; so if ex == -1074299 * then we have something [2^-1075,2^-1074). Under round-to-300 * nearest-even, this whole interval rounds up to 2^-1074,301 * except for the bottom endpoint which rounds to even and is302 * an underflow condition.303 *304 * So, ex < -1074 is definite underflow, and ex == -1074 is305 * underflow iff all mantissa bits are zero.306 */307 if (ex < -1074 || (ex == -1074 && (y[0] & 0xFFFFF) == 0 && y[1] == 0)) {308 out[0] = sign; /* underflow: correctly signed zero */309 out[1] = 0;310 return "underflow";311 }312 313 /*314 * No overflow or underflow; should be nice and simple, unless315 * we have to denormalise and round the result.316 */317 if (ex < -1021) { /* denormalise and round */318 uint32 roundword;319 y[0] &= 0x000FFFFF;320 y[0] |= 0x00100000; /* set leading bit */321 roundword = 0;322 while (ex < -1021) {323 if (roundword & 1)324 roundword |= 2; /* preserve sticky bit */325 roundword = (roundword >> 1) | ((y[1] & 1) << 31);326 y[1] = (y[1] >> 1) | ((y[0] & 1) << 31);327 y[0] = y[0] >> 1;328 ex++;329 }330 if (roundword > 0x80000000 || /* round up */331 (roundword == 0x80000000 && (y[1] & 1))) { /* round up to even */332 y[1]++;333 y[0] += (y[1] == 0);334 }335 out[0] = sign | y[0];336 out[1] = y[1];337 /* Proper ERANGE underflow was handled earlier, but we still338 * expect an IEEE Underflow exception if this partially339 * underflowed result is not exact. */340 if (roundword)341 return "u";342 return NULL; /* underflow was handled earlier */343 } else {344 out[0] = y[0] + (ex << 20);345 out[1] = y[1];346 return NULL;347 }348}349 350char *test_ldexpf(uint32 *x, uint32 *np, uint32 *out) {351 int n = *np;352 int32 n2;353 uint32 y;354 int ex = (*x >> 23) & 0xFF; /* exponent */355 int sign = *x & 0x80000000;356 357 if (ex == 0xFF) { /* inf/NaN; just return x */358 *out = *x;359 return NULL;360 }361 if ((*x & 0x7FFFFFFF) == 0) { /* zero: return x */362 *out = *x;363 return NULL;364 }365 366 test_frexpf(x, &y, (uint32 *)&n2);367 ex = n + n2;368 if (ex > 0x80) { /* overflow */369 *out = sign | 0x7F800000;370 return "overflow";371 }372 /*373 * Underflow. 2^-149 is 00000001; so if ex == -149 then we have374 * something [2^-150,2^-149). Under round-to- nearest-even,375 * this whole interval rounds up to 2^-149, except for the376 * bottom endpoint which rounds to even and is an underflow377 * condition.378 *379 * So, ex < -149 is definite underflow, and ex == -149 is380 * underflow iff all mantissa bits are zero.381 */382 if (ex < -149 || (ex == -149 && (y & 0x7FFFFF) == 0)) {383 *out = sign; /* underflow: correctly signed zero */384 return "underflow";385 }386 387 /*388 * No overflow or underflow; should be nice and simple, unless389 * we have to denormalise and round the result.390 */391 if (ex < -125) { /* denormalise and round */392 uint32 roundword;393 y &= 0x007FFFFF;394 y |= 0x00800000; /* set leading bit */395 roundword = 0;396 while (ex < -125) {397 if (roundword & 1)398 roundword |= 2; /* preserve sticky bit */399 roundword = (roundword >> 1) | ((y & 1) << 31);400 y = y >> 1;401 ex++;402 }403 if (roundword > 0x80000000 || /* round up */404 (roundword == 0x80000000 && (y & 1))) { /* round up to even */405 y++;406 }407 *out = sign | y;408 /* Proper ERANGE underflow was handled earlier, but we still409 * expect an IEEE Underflow exception if this partially410 * underflowed result is not exact. */411 if (roundword)412 return "u";413 return NULL; /* underflow was handled earlier */414 } else {415 *out = y + (ex << 23);416 return NULL;417 }418}419 420char *test_frexp(uint32 *x, uint32 *out, uint32 *nout) {421 int ex = (x[0] >> 20) & 0x7FF; /* exponent */422 if (ex == 0x7FF) { /* inf/NaN; return x/0 */423 out[0] = x[0];424 out[1] = x[1];425 nout[0] = 0;426 return NULL;427 }428 if (ex == 0) { /* denormals/zeros */429 int sign;430 uint32 xh, xl;431 if ((x[0] & 0x7FFFFFFF) == 0 && x[1] == 0) {432 /* zero: return x/0 */433 out[0] = x[0];434 out[1] = x[1];435 nout[0] = 0;436 return NULL;437 }438 sign = x[0] & 0x80000000;439 xh = x[0] & 0x7FFFFFFF;440 xl = x[1];441 ex = 1;442 while (!(xh & 0x100000)) {443 ex--;444 xh = (xh << 1) | ((xl >> 31) & 1);445 xl = (xl & 0x7FFFFFFF) << 1;446 }447 out[0] = sign | 0x3FE00000 | (xh & 0xFFFFF);448 out[1] = xl;449 nout[0] = ex - 0x3FE;450 return NULL;451 }452 out[0] = 0x3FE00000 | (x[0] & 0x800FFFFF);453 out[1] = x[1];454 nout[0] = ex - 0x3FE;455 return NULL; /* ordinary number; no error */456}457 458char *test_frexpf(uint32 *x, uint32 *out, uint32 *nout) {459 int ex = (*x >> 23) & 0xFF; /* exponent */460 if (ex == 0xFF) { /* inf/NaN; return x/0 */461 *out = *x;462 nout[0] = 0;463 return NULL;464 }465 if (ex == 0) { /* denormals/zeros */466 int sign;467 uint32 xv;468 if ((*x & 0x7FFFFFFF) == 0) {469 /* zero: return x/0 */470 *out = *x;471 nout[0] = 0;472 return NULL;473 }474 sign = *x & 0x80000000;475 xv = *x & 0x7FFFFFFF;476 ex = 1;477 while (!(xv & 0x800000)) {478 ex--;479 xv = xv << 1;480 }481 *out = sign | 0x3F000000 | (xv & 0x7FFFFF);482 nout[0] = ex - 0x7E;483 return NULL;484 }485 *out = 0x3F000000 | (*x & 0x807FFFFF);486 nout[0] = ex - 0x7E;487 return NULL; /* ordinary number; no error */488}489 490char *test_modf(uint32 *x, uint32 *fout, uint32 *iout) {491 int ex = (x[0] >> 20) & 0x7FF; /* exponent */492 int sign = x[0] & 0x80000000;493 uint32 fh, fl;494 495 if (((x[0] & 0x7FFFFFFF) | (!!x[1])) > 0x7FF00000) {496 /*497 * NaN input: return the same in _both_ outputs.498 */499 fout[0] = iout[0] = x[0];500 fout[1] = iout[1] = x[1];501 return NULL;502 }503 504 test_rint(x, iout, 0, 0);505 fh = x[0] - iout[0];506 fl = x[1] - iout[1];507 if (!fh && !fl) { /* no fraction part */508 fout[0] = sign;509 fout[1] = 0;510 return NULL;511 }512 if (!(iout[0] & 0x7FFFFFFF) && !iout[1]) { /* no integer part */513 fout[0] = x[0];514 fout[1] = x[1];515 return NULL;516 }517 while (!(fh & 0x100000)) {518 ex--;519 fh = (fh << 1) | ((fl >> 31) & 1);520 fl = (fl & 0x7FFFFFFF) << 1;521 }522 fout[0] = sign | (ex << 20) | (fh & 0xFFFFF);523 fout[1] = fl;524 return NULL;525}526 527char *test_modff(uint32 *x, uint32 *fout, uint32 *iout) {528 int ex = (*x >> 23) & 0xFF; /* exponent */529 int sign = *x & 0x80000000;530 uint32 f;531 532 if ((*x & 0x7FFFFFFF) > 0x7F800000) {533 /*534 * NaN input: return the same in _both_ outputs.535 */536 *fout = *iout = *x;537 return NULL;538 }539 540 test_rintf(x, iout, 0, 0);541 f = *x - *iout;542 if (!f) { /* no fraction part */543 *fout = sign;544 return NULL;545 }546 if (!(*iout & 0x7FFFFFFF)) { /* no integer part */547 *fout = *x;548 return NULL;549 }550 while (!(f & 0x800000)) {551 ex--;552 f = f << 1;553 }554 *fout = sign | (ex << 23) | (f & 0x7FFFFF);555 return NULL;556}557 558char *test_copysign(uint32 *x, uint32 *y, uint32 *out)559{560 int ysign = y[0] & 0x80000000;561 int xhigh = x[0] & 0x7fffffff;562 563 out[0] = ysign | xhigh;564 out[1] = x[1];565 566 /* There can be no error */567 return NULL;568}569 570char *test_copysignf(uint32 *x, uint32 *y, uint32 *out)571{572 int ysign = y[0] & 0x80000000;573 int xhigh = x[0] & 0x7fffffff;574 575 out[0] = ysign | xhigh;576 577 /* There can be no error */578 return NULL;579}580 581char *test_isfinite(uint32 *x, uint32 *out)582{583 int xhigh = x[0];584 /* Being finite means that the exponent is not 0x7ff */585 if ((xhigh & 0x7ff00000) == 0x7ff00000) out[0] = 0;586 else out[0] = 1;587 return NULL;588}589 590char *test_isfinitef(uint32 *x, uint32 *out)591{592 /* Being finite means that the exponent is not 0xff */593 if ((x[0] & 0x7f800000) == 0x7f800000) out[0] = 0;594 else out[0] = 1;595 return NULL;596}597 598char *test_isinff(uint32 *x, uint32 *out)599{600 /* Being infinite means that our bottom 30 bits equate to 0x7f800000 */601 if ((x[0] & 0x7fffffff) == 0x7f800000) out[0] = 1;602 else out[0] = 0;603 return NULL;604}605 606char *test_isinf(uint32 *x, uint32 *out)607{608 int xhigh = x[0];609 int xlow = x[1];610 /* Being infinite means that our fraction is zero and exponent is 0x7ff */611 if (((xhigh & 0x7fffffff) == 0x7ff00000) && (xlow == 0)) out[0] = 1;612 else out[0] = 0;613 return NULL;614}615 616char *test_isnanf(uint32 *x, uint32 *out)617{618 /* Being NaN means that our exponent is 0xff and non-0 fraction */619 int exponent = x[0] & 0x7f800000;620 int fraction = x[0] & 0x007fffff;621 if ((exponent == 0x7f800000) && (fraction != 0)) out[0] = 1;622 else out[0] = 0;623 return NULL;624}625 626char *test_isnan(uint32 *x, uint32 *out)627{628 /* Being NaN means that our exponent is 0x7ff and non-0 fraction */629 int exponent = x[0] & 0x7ff00000;630 int fractionhigh = x[0] & 0x000fffff;631 if ((exponent == 0x7ff00000) && ((fractionhigh != 0) || x[1] != 0))632 out[0] = 1;633 else out[0] = 0;634 return NULL;635}636 637char *test_isnormalf(uint32 *x, uint32 *out)638{639 /* Being normal means exponent is not 0 and is not 0xff */640 int exponent = x[0] & 0x7f800000;641 if (exponent == 0x7f800000) out[0] = 0;642 else if (exponent == 0) out[0] = 0;643 else out[0] = 1;644 return NULL;645}646 647char *test_isnormal(uint32 *x, uint32 *out)648{649 /* Being normal means exponent is not 0 and is not 0x7ff */650 int exponent = x[0] & 0x7ff00000;651 if (exponent == 0x7ff00000) out[0] = 0;652 else if (exponent == 0) out[0] = 0;653 else out[0] = 1;654 return NULL;655}656 657char *test_signbitf(uint32 *x, uint32 *out)658{659 /* Sign bit is bit 31 */660 out[0] = (x[0] >> 31) & 1;661 return NULL;662}663 664char *test_signbit(uint32 *x, uint32 *out)665{666 /* Sign bit is bit 31 */667 out[0] = (x[0] >> 31) & 1;668 return NULL;669}670 671char *test_fpclassify(uint32 *x, uint32 *out)672{673 int exponent = (x[0] & 0x7ff00000) >> 20;674 int fraction = (x[0] & 0x000fffff) | x[1];675 676 if ((exponent == 0x00) && (fraction == 0)) out[0] = 0;677 else if ((exponent == 0x00) && (fraction != 0)) out[0] = 4;678 else if ((exponent == 0x7ff) && (fraction == 0)) out[0] = 3;679 else if ((exponent == 0x7ff) && (fraction != 0)) out[0] = 7;680 else out[0] = 5;681 return NULL;682}683 684char *test_fpclassifyf(uint32 *x, uint32 *out)685{686 int exponent = (x[0] & 0x7f800000) >> 23;687 int fraction = x[0] & 0x007fffff;688 689 if ((exponent == 0x000) && (fraction == 0)) out[0] = 0;690 else if ((exponent == 0x000) && (fraction != 0)) out[0] = 4;691 else if ((exponent == 0xff) && (fraction == 0)) out[0] = 3;692 else if ((exponent == 0xff) && (fraction != 0)) out[0] = 7;693 else out[0] = 5;694 return NULL;695}696 697/*698 * Internal function that compares doubles in x & y and returns -3, -2, -1, 0,699 * 1 if they compare to be signaling, unordered, less than, equal or greater700 * than.701 */702static int fpcmp4(uint32 *x, uint32 *y)703{704 int result = 0;705 706 /*707 * Sort out whether results are ordered or not to begin with708 * NaNs have exponent 0x7ff, and non-zero fraction. Signaling NaNs take709 * higher priority than quiet ones.710 */711 if ((x[0] & 0x7fffffff) >= 0x7ff80000) result = -2;712 else if ((x[0] & 0x7fffffff) > 0x7ff00000) result = -3;713 else if (((x[0] & 0x7fffffff) == 0x7ff00000) && (x[1] != 0)) result = -3;714 if ((y[0] & 0x7fffffff) >= 0x7ff80000 && result != -3) result = -2;715 else if ((y[0] & 0x7fffffff) > 0x7ff00000) result = -3;716 else if (((y[0] & 0x7fffffff) == 0x7ff00000) && (y[1] != 0)) result = -3;717 if (result != 0) return result;718 719 /*720 * The two forms of zero are equal721 */722 if (((x[0] & 0x7fffffff) == 0) && x[1] == 0 &&723 ((y[0] & 0x7fffffff) == 0) && y[1] == 0)724 return 0;725 726 /*727 * If x and y have different signs we can tell that they're not equal728 * If x is +ve we have x > y return 1 - otherwise y is +ve return -1729 */730 if ((x[0] >> 31) != (y[0] >> 31))731 return ((x[0] >> 31) == 0) - ((y[0] >> 31) == 0);732 733 /*734 * Now we have both signs the same, let's do an initial compare of the735 * values.736 *737 * Whoever designed IEEE754's floating point formats is very clever and738 * earns my undying admiration. Once you remove the sign-bit, the739 * floating point numbers can be ordered using the standard <, ==, >740 * operators will treating the fp-numbers as integers with that bit-741 * pattern.742 */743 if ((x[0] & 0x7fffffff) < (y[0] & 0x7fffffff)) result = -1;744 else if ((x[0] & 0x7fffffff) > (y[0] & 0x7fffffff)) result = 1;745 else if (x[1] < y[1]) result = -1;746 else if (x[1] > y[1]) result = 1;747 else result = 0;748 749 /*750 * Now we return the result - is x is positive (and therefore so is y) we751 * return the plain result - otherwise we negate it and return.752 */753 if ((x[0] >> 31) == 0) return result;754 else return -result;755}756 757/*758 * Internal function that compares floats in x & y and returns -3, -2, -1, 0,759 * 1 if they compare to be signaling, unordered, less than, equal or greater760 * than.761 */762static int fpcmp4f(uint32 *x, uint32 *y)763{764 int result = 0;765 766 /*767 * Sort out whether results are ordered or not to begin with768 * NaNs have exponent 0xff, and non-zero fraction - we have to handle all769 * signaling cases over the quiet ones770 */771 if ((x[0] & 0x7fffffff) >= 0x7fc00000) result = -2;772 else if ((x[0] & 0x7fffffff) > 0x7f800000) result = -3;773 if ((y[0] & 0x7fffffff) >= 0x7fc00000 && result != -3) result = -2;774 else if ((y[0] & 0x7fffffff) > 0x7f800000) result = -3;775 if (result != 0) return result;776 777 /*778 * The two forms of zero are equal779 */780 if (((x[0] & 0x7fffffff) == 0) && ((y[0] & 0x7fffffff) == 0))781 return 0;782 783 /*784 * If x and y have different signs we can tell that they're not equal785 * If x is +ve we have x > y return 1 - otherwise y is +ve return -1786 */787 if ((x[0] >> 31) != (y[0] >> 31))788 return ((x[0] >> 31) == 0) - ((y[0] >> 31) == 0);789 790 /*791 * Now we have both signs the same, let's do an initial compare of the792 * values.793 *794 * Whoever designed IEEE754's floating point formats is very clever and795 * earns my undying admiration. Once you remove the sign-bit, the796 * floating point numbers can be ordered using the standard <, ==, >797 * operators will treating the fp-numbers as integers with that bit-798 * pattern.799 */800 if ((x[0] & 0x7fffffff) < (y[0] & 0x7fffffff)) result = -1;801 else if ((x[0] & 0x7fffffff) > (y[0] & 0x7fffffff)) result = 1;802 else result = 0;803 804 /*805 * Now we return the result - is x is positive (and therefore so is y) we806 * return the plain result - otherwise we negate it and return.807 */808 if ((x[0] >> 31) == 0) return result;809 else return -result;810}811 812char *test_isgreater(uint32 *x, uint32 *y, uint32 *out)813{814 int result = fpcmp4(x, y);815 *out = (result == 1);816 return result == -3 ? "i" : NULL;817}818 819char *test_isgreaterequal(uint32 *x, uint32 *y, uint32 *out)820{821 int result = fpcmp4(x, y);822 *out = (result >= 0);823 return result == -3 ? "i" : NULL;824}825 826char *test_isless(uint32 *x, uint32 *y, uint32 *out)827{828 int result = fpcmp4(x, y);829 *out = (result == -1);830 return result == -3 ? "i" : NULL;831}832 833char *test_islessequal(uint32 *x, uint32 *y, uint32 *out)834{835 int result = fpcmp4(x, y);836 *out = (result == -1) || (result == 0);837 return result == -3 ? "i" : NULL;838}839 840char *test_islessgreater(uint32 *x, uint32 *y, uint32 *out)841{842 int result = fpcmp4(x, y);843 *out = (result == -1) || (result == 1);844 return result == -3 ? "i" : NULL;845}846 847char *test_isunordered(uint32 *x, uint32 *y, uint32 *out)848{849 int normal = 0;850 int result = fpcmp4(x, y);851 852 test_isnormal(x, out);853 normal |= *out;854 test_isnormal(y, out);855 normal |= *out;856 *out = (result == -2) || (result == -3);857 return result == -3 ? "i" : NULL;858}859 860char *test_isgreaterf(uint32 *x, uint32 *y, uint32 *out)861{862 int result = fpcmp4f(x, y);863 *out = (result == 1);864 return result == -3 ? "i" : NULL;865}866 867char *test_isgreaterequalf(uint32 *x, uint32 *y, uint32 *out)868{869 int result = fpcmp4f(x, y);870 *out = (result >= 0);871 return result == -3 ? "i" : NULL;872}873 874char *test_islessf(uint32 *x, uint32 *y, uint32 *out)875{876 int result = fpcmp4f(x, y);877 *out = (result == -1);878 return result == -3 ? "i" : NULL;879}880 881char *test_islessequalf(uint32 *x, uint32 *y, uint32 *out)882{883 int result = fpcmp4f(x, y);884 *out = (result == -1) || (result == 0);885 return result == -3 ? "i" : NULL;886}887 888char *test_islessgreaterf(uint32 *x, uint32 *y, uint32 *out)889{890 int result = fpcmp4f(x, y);891 *out = (result == -1) || (result == 1);892 return result == -3 ? "i" : NULL;893}894 895char *test_isunorderedf(uint32 *x, uint32 *y, uint32 *out)896{897 int normal = 0;898 int result = fpcmp4f(x, y);899 900 test_isnormalf(x, out);901 normal |= *out;902 test_isnormalf(y, out);903 normal |= *out;904 *out = (result == -2) || (result == -3);905 return result == -3 ? "i" : NULL;906}907