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