1 /* fpu.c --- FPU emulator for stand-alone RX simulator. 2 3 Copyright (C) 2008-2024 Free Software Foundation, Inc. 4 Contributed by Red Hat, Inc. 5 6 This file is part of the GNU simulators. 7 8 This program is free software; you can redistribute it and/or modify 9 it under the terms of the GNU General Public License as published by 10 the Free Software Foundation; either version 3 of the License, or 11 (at your option) any later version. 12 13 This program is distributed in the hope that it will be useful, 14 but WITHOUT ANY WARRANTY; without even the implied warranty of 15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 16 GNU General Public License for more details. 17 18 You should have received a copy of the GNU General Public License 19 along with this program. If not, see <http://www.gnu.org/licenses/>. */ 20 21 /* This must come before any other includes. */ 22 #include "defs.h" 23 24 #include <stdio.h> 25 #include <stdlib.h> 26 27 #include "cpu.h" 28 #include "fpu.h" 29 30 /* FPU encodings are as follows: 31 32 S EXPONENT MANTISSA 33 1 12345678 12345678901234567890123 34 35 0 00000000 00000000000000000000000 +0 36 1 00000000 00000000000000000000000 -0 37 38 X 00000000 00000000000000000000001 Denormals 39 X 00000000 11111111111111111111111 40 41 X 00000001 XXXXXXXXXXXXXXXXXXXXXXX Normals 42 X 11111110 XXXXXXXXXXXXXXXXXXXXXXX 43 44 0 11111111 00000000000000000000000 +Inf 45 1 11111111 00000000000000000000000 -Inf 46 47 X 11111111 0XXXXXXXXXXXXXXXXXXXXXX SNaN (X != 0) 48 X 11111111 1XXXXXXXXXXXXXXXXXXXXXX QNaN (X != 0) 49 50 */ 51 52 #define trace 0 53 #define tprintf if (trace) printf 54 55 /* Some magic numbers. */ 56 #define PLUS_MAX 0x7f7fffffUL 57 #define MINUS_MAX 0xff7fffffUL 58 #define PLUS_INF 0x7f800000UL 59 #define MINUS_INF 0xff800000UL 60 #define PLUS_ZERO 0x00000000UL 61 #define MINUS_ZERO 0x80000000UL 62 63 #define FP_RAISE(e) fp_raise(FPSWBITS_C##e) 64 static void 65 fp_raise (int mask) 66 { 67 regs.r_fpsw |= mask; 68 if (mask != FPSWBITS_CE) 69 { 70 if (regs.r_fpsw & (mask << FPSW_CESH)) 71 regs.r_fpsw |= (mask << FPSW_CFSH); 72 if (regs.r_fpsw & FPSWBITS_FMASK) 73 regs.r_fpsw |= FPSWBITS_FSUM; 74 else 75 regs.r_fpsw &= ~FPSWBITS_FSUM; 76 } 77 } 78 79 /* We classify all numbers as one of these. They correspond to the 80 rows/colums in the exception tables. */ 81 typedef enum { 82 FP_NORMAL, 83 FP_PZERO, 84 FP_NZERO, 85 FP_PINFINITY, 86 FP_NINFINITY, 87 FP_DENORMAL, 88 FP_QNAN, 89 FP_SNAN 90 } FP_Type; 91 92 #if defined DEBUG0 93 static const char *fpt_names[] = { 94 "Normal", "+0", "-0", "+Inf", "-Inf", "Denormal", "QNaN", "SNaN" 95 }; 96 #endif 97 98 #define EXP_BIAS 127 99 #define EXP_ZERO -127 100 #define EXP_INF 128 101 102 #define MANT_BIAS 0x00080000UL 103 104 typedef struct { 105 int exp; 106 unsigned int mant; /* 24 bits */ 107 char type; 108 char sign; 109 fp_t orig_value; 110 } FP_Parts; 111 112 static void 113 fp_explode (fp_t f, FP_Parts *p) 114 { 115 int exp, mant, sign; 116 117 exp = ((f & 0x7f800000UL) >> 23); 118 mant = f & 0x007fffffUL; 119 sign = f & 0x80000000UL; 120 /*printf("explode: %08x %x %2x %6x\n", f, sign, exp, mant);*/ 121 122 p->sign = sign ? -1 : 1; 123 p->exp = exp - EXP_BIAS; 124 p->orig_value = f; 125 p->mant = mant | 0x00800000UL; 126 127 if (p->exp == EXP_ZERO) 128 { 129 if (regs.r_fpsw & FPSWBITS_DN) 130 mant = 0; 131 if (mant) 132 p->type = FP_DENORMAL; 133 else 134 { 135 p->mant = 0; 136 p->type = sign ? FP_NZERO : FP_PZERO; 137 } 138 } 139 else if (p->exp == EXP_INF) 140 { 141 if (mant == 0) 142 p->type = sign ? FP_NINFINITY : FP_PINFINITY; 143 else if (mant & 0x00400000UL) 144 p->type = FP_QNAN; 145 else 146 p->type = FP_SNAN; 147 } 148 else 149 p->type = FP_NORMAL; 150 } 151 152 static fp_t 153 fp_implode (FP_Parts *p) 154 { 155 int exp, mant; 156 157 exp = p->exp + EXP_BIAS; 158 mant = p->mant; 159 /*printf("implode: exp %d mant 0x%x\n", exp, mant);*/ 160 if (p->type == FP_NORMAL) 161 { 162 while (mant 163 && exp > 0 164 && mant < 0x00800000UL) 165 { 166 mant <<= 1; 167 exp --; 168 } 169 while (mant > 0x00ffffffUL) 170 { 171 mant >>= 1; 172 exp ++; 173 } 174 if (exp < 0) 175 { 176 /* underflow */ 177 exp = 0; 178 mant = 0; 179 FP_RAISE (E); 180 } 181 if (exp >= 255) 182 { 183 /* overflow */ 184 exp = 255; 185 mant = 0; 186 FP_RAISE (O); 187 } 188 } 189 mant &= 0x007fffffUL; 190 exp &= 0xff; 191 mant |= exp << 23; 192 if (p->sign < 0) 193 mant |= 0x80000000UL; 194 195 return mant; 196 } 197 198 typedef union { 199 unsigned long long ll; 200 double d; 201 } U_d_ll; 202 203 static int checked_format = 0; 204 205 /* We assume a double format like this: 206 S[1] E[11] M[52] 207 */ 208 209 static double 210 fp_to_double (FP_Parts *p) 211 { 212 U_d_ll u; 213 214 if (!checked_format) 215 { 216 u.d = 1.5; 217 if (u.ll != 0x3ff8000000000000ULL) 218 abort (); 219 u.d = -225; 220 if (u.ll != 0xc06c200000000000ULL) 221 abort (); 222 u.d = 10.1; 223 if (u.ll != 0x4024333333333333ULL) 224 abort (); 225 checked_format = 1; 226 } 227 228 u.ll = 0; 229 if (p->sign < 0) 230 u.ll |= (1ULL << 63); 231 /* Make sure a zero encoding stays a zero. */ 232 if (p->exp != -EXP_BIAS) 233 u.ll |= ((unsigned long long)p->exp + 1023ULL) << 52; 234 u.ll |= (unsigned long long) (p->mant & 0x007fffffUL) << (52 - 23); 235 return u.d; 236 } 237 238 static void 239 double_to_fp (double d, FP_Parts *p) 240 { 241 int exp; 242 U_d_ll u; 243 int sign; 244 245 u.d = d; 246 247 sign = (u.ll & 0x8000000000000000ULL) ? 1 : 0; 248 exp = u.ll >> 52; 249 exp = (exp & 0x7ff); 250 251 if (exp == 0) 252 { 253 /* A generated denormal should show up as an underflow, not 254 here. */ 255 if (sign) 256 fp_explode (MINUS_ZERO, p); 257 else 258 fp_explode (PLUS_ZERO, p); 259 return; 260 } 261 262 exp = exp - 1023; 263 if ((exp + EXP_BIAS) > 254) 264 { 265 FP_RAISE (O); 266 switch (regs.r_fpsw & FPSWBITS_RM) 267 { 268 case FPRM_NEAREST: 269 if (sign) 270 fp_explode (MINUS_INF, p); 271 else 272 fp_explode (PLUS_INF, p); 273 break; 274 case FPRM_ZERO: 275 if (sign) 276 fp_explode (MINUS_MAX, p); 277 else 278 fp_explode (PLUS_MAX, p); 279 break; 280 case FPRM_PINF: 281 if (sign) 282 fp_explode (MINUS_MAX, p); 283 else 284 fp_explode (PLUS_INF, p); 285 break; 286 case FPRM_NINF: 287 if (sign) 288 fp_explode (MINUS_INF, p); 289 else 290 fp_explode (PLUS_MAX, p); 291 break; 292 } 293 return; 294 } 295 if ((exp + EXP_BIAS) < 1) 296 { 297 if (sign) 298 fp_explode (MINUS_ZERO, p); 299 else 300 fp_explode (PLUS_ZERO, p); 301 FP_RAISE (U); 302 } 303 304 p->sign = sign ? -1 : 1; 305 p->exp = exp; 306 p->mant = u.ll >> (52-23) & 0x007fffffUL; 307 p->mant |= 0x00800000UL; 308 p->type = FP_NORMAL; 309 310 if (u.ll & 0x1fffffffULL) 311 { 312 switch (regs.r_fpsw & FPSWBITS_RM) 313 { 314 case FPRM_NEAREST: 315 if (u.ll & 0x10000000ULL) 316 p->mant ++; 317 break; 318 case FPRM_ZERO: 319 break; 320 case FPRM_PINF: 321 if (sign == 1) 322 p->mant ++; 323 break; 324 case FPRM_NINF: 325 if (sign == -1) 326 p->mant ++; 327 break; 328 } 329 FP_RAISE (X); 330 } 331 332 } 333 334 typedef enum { 335 eNR, /* Use the normal result. */ 336 ePZ, eNZ, /* +- zero */ 337 eSZ, /* signed zero - XOR signs of ops together. */ 338 eRZ, /* +- zero depending on rounding mode. */ 339 ePI, eNI, /* +- Infinity */ 340 eSI, /* signed infinity - XOR signs of ops together. */ 341 eQN, eSN, /* Quiet/Signalling NANs */ 342 eIn, /* Invalid. */ 343 eUn, /* Unimplemented. */ 344 eDZ, /* Divide-by-zero. */ 345 eLT, /* less than */ 346 eGT, /* greater than */ 347 eEQ, /* equal to */ 348 } FP_ExceptionCases; 349 350 #if defined DEBUG0 351 static const char *ex_names[] = { 352 "NR", "PZ", "NZ", "SZ", "RZ", "PI", "NI", "SI", "QN", "SN", "IN", "Un", "DZ", "LT", "GT", "EQ" 353 }; 354 #endif 355 356 /* This checks for all exceptional cases (not all FP exceptions) and 357 returns TRUE if it is providing the result in *c. If it returns 358 FALSE, the caller should do the "normal" operation. */ 359 static int 360 check_exceptions (FP_Parts *a, FP_Parts *b, fp_t *c, 361 FP_ExceptionCases ex_tab[5][5], 362 FP_ExceptionCases *case_ret) 363 { 364 FP_ExceptionCases fpec; 365 366 if (a->type == FP_SNAN 367 || b->type == FP_SNAN) 368 fpec = eIn; 369 else if (a->type == FP_QNAN 370 || b->type == FP_QNAN) 371 fpec = eQN; 372 else if (a->type == FP_DENORMAL 373 || b->type == FP_DENORMAL) 374 fpec = eUn; 375 else 376 fpec = ex_tab[(int)(a->type)][(int)(b->type)]; 377 378 /*printf("%s %s -> %s\n", fpt_names[(int)(a->type)], fpt_names[(int)(b->type)], ex_names[(int)(fpec)]);*/ 379 380 if (case_ret) 381 *case_ret = fpec; 382 383 switch (fpec) 384 { 385 case eNR: /* Use the normal result. */ 386 return 0; 387 388 case ePZ: /* + zero */ 389 *c = 0x00000000; 390 return 1; 391 392 case eNZ: /* - zero */ 393 *c = 0x80000000; 394 return 1; 395 396 case eSZ: /* signed zero */ 397 *c = (a->sign == b->sign) ? PLUS_ZERO : MINUS_ZERO; 398 return 1; 399 400 case eRZ: /* +- zero depending on rounding mode. */ 401 if ((regs.r_fpsw & FPSWBITS_RM) == FPRM_NINF) 402 *c = 0x80000000; 403 else 404 *c = 0x00000000; 405 return 1; 406 407 case ePI: /* + Infinity */ 408 *c = 0x7F800000; 409 return 1; 410 411 case eNI: /* - Infinity */ 412 *c = 0xFF800000; 413 return 1; 414 415 case eSI: /* sign Infinity */ 416 *c = (a->sign == b->sign) ? PLUS_INF : MINUS_INF; 417 return 1; 418 419 case eQN: /* Quiet NANs */ 420 if (a->type == FP_QNAN) 421 *c = a->orig_value; 422 else 423 *c = b->orig_value; 424 return 1; 425 426 case eSN: /* Signalling NANs */ 427 if (a->type == FP_SNAN) 428 *c = a->orig_value; 429 else 430 *c = b->orig_value; 431 FP_RAISE (V); 432 return 1; 433 434 case eIn: /* Invalid. */ 435 FP_RAISE (V); 436 if (a->type == FP_SNAN) 437 *c = a->orig_value | 0x00400000; 438 else if (b->type == FP_SNAN) 439 *c = b->orig_value | 0x00400000; 440 else 441 *c = 0x7fc00000; 442 return 1; 443 444 case eUn: /* Unimplemented. */ 445 FP_RAISE (E); 446 return 1; 447 448 case eDZ: /* Division-by-zero. */ 449 *c = (a->sign == b->sign) ? PLUS_INF : MINUS_INF; 450 FP_RAISE (Z); 451 return 1; 452 453 default: 454 return 0; 455 } 456 } 457 458 #define CHECK_EXCEPTIONS(FPPa, FPPb, fpc, ex_tab) \ 459 if (check_exceptions (&FPPa, &FPPb, &fpc, ex_tab, 0)) \ 460 return fpc; 461 462 /* For each operation, we have two tables of how nonnormal cases are 463 handled. The DN=0 case is first, followed by the DN=1 case, with 464 each table using the following layout: */ 465 466 static FP_ExceptionCases ex_add_tab[5][5] = { 467 /* N +0 -0 +In -In */ 468 { eNR, eNR, eNR, ePI, eNI }, /* Normal */ 469 { eNR, ePZ, eRZ, ePI, eNI }, /* +0 */ 470 { eNR, eRZ, eNZ, ePI, eNI }, /* -0 */ 471 { ePI, ePI, ePI, ePI, eIn }, /* +Inf */ 472 { eNI, eNI, eNI, eIn, eNI }, /* -Inf */ 473 }; 474 475 fp_t 476 rxfp_add (fp_t fa, fp_t fb) 477 { 478 FP_Parts a, b, c; 479 fp_t rv; 480 double da, db; 481 482 fp_explode (fa, &a); 483 fp_explode (fb, &b); 484 CHECK_EXCEPTIONS (a, b, rv, ex_add_tab); 485 486 da = fp_to_double (&a); 487 db = fp_to_double (&b); 488 tprintf("%g + %g = %g\n", da, db, da+db); 489 490 double_to_fp (da+db, &c); 491 rv = fp_implode (&c); 492 return rv; 493 } 494 495 static FP_ExceptionCases ex_sub_tab[5][5] = { 496 /* N +0 -0 +In -In */ 497 { eNR, eNR, eNR, eNI, ePI }, /* Normal */ 498 { eNR, eRZ, ePZ, eNI, ePI }, /* +0 */ 499 { eNR, eNZ, eRZ, eNI, ePI }, /* -0 */ 500 { ePI, ePI, ePI, eIn, ePI }, /* +Inf */ 501 { eNI, eNI, eNI, eNI, eIn }, /* -Inf */ 502 }; 503 504 fp_t 505 rxfp_sub (fp_t fa, fp_t fb) 506 { 507 FP_Parts a, b, c; 508 fp_t rv; 509 double da, db; 510 511 fp_explode (fa, &a); 512 fp_explode (fb, &b); 513 CHECK_EXCEPTIONS (a, b, rv, ex_sub_tab); 514 515 da = fp_to_double (&a); 516 db = fp_to_double (&b); 517 tprintf("%g - %g = %g\n", da, db, da-db); 518 519 double_to_fp (da-db, &c); 520 rv = fp_implode (&c); 521 522 return rv; 523 } 524 525 static FP_ExceptionCases ex_mul_tab[5][5] = { 526 /* N +0 -0 +In -In */ 527 { eNR, eNR, eNR, eSI, eSI }, /* Normal */ 528 { eNR, ePZ, eNZ, eIn, eIn }, /* +0 */ 529 { eNR, eNZ, ePZ, eIn, eIn }, /* -0 */ 530 { eSI, eIn, eIn, ePI, eNI }, /* +Inf */ 531 { eSI, eIn, eIn, eNI, ePI }, /* -Inf */ 532 }; 533 534 fp_t 535 rxfp_mul (fp_t fa, fp_t fb) 536 { 537 FP_Parts a, b, c; 538 fp_t rv; 539 double da, db; 540 541 fp_explode (fa, &a); 542 fp_explode (fb, &b); 543 CHECK_EXCEPTIONS (a, b, rv, ex_mul_tab); 544 545 da = fp_to_double (&a); 546 db = fp_to_double (&b); 547 tprintf("%g x %g = %g\n", da, db, da*db); 548 549 double_to_fp (da*db, &c); 550 rv = fp_implode (&c); 551 552 return rv; 553 } 554 555 static FP_ExceptionCases ex_div_tab[5][5] = { 556 /* N +0 -0 +In -In */ 557 { eNR, eDZ, eDZ, eSZ, eSZ }, /* Normal */ 558 { eSZ, eIn, eIn, ePZ, eNZ }, /* +0 */ 559 { eSZ, eIn, eIn, eNZ, ePZ }, /* -0 */ 560 { eSI, ePI, eNI, eIn, eIn }, /* +Inf */ 561 { eSI, eNI, ePI, eIn, eIn }, /* -Inf */ 562 }; 563 564 fp_t 565 rxfp_div (fp_t fa, fp_t fb) 566 { 567 FP_Parts a, b, c; 568 fp_t rv; 569 double da, db; 570 571 fp_explode (fa, &a); 572 fp_explode (fb, &b); 573 CHECK_EXCEPTIONS (a, b, rv, ex_div_tab); 574 575 da = fp_to_double (&a); 576 db = fp_to_double (&b); 577 tprintf("%g / %g = %g\n", da, db, da/db); 578 579 double_to_fp (da/db, &c); 580 rv = fp_implode (&c); 581 582 return rv; 583 } 584 585 static FP_ExceptionCases ex_cmp_tab[5][5] = { 586 /* N +0 -0 +In -In */ 587 { eNR, eNR, eNR, eLT, eGT }, /* Normal */ 588 { eNR, eEQ, eEQ, eLT, eGT }, /* +0 */ 589 { eNR, eEQ, eEQ, eLT, eGT }, /* -0 */ 590 { eGT, eGT, eGT, eEQ, eGT }, /* +Inf */ 591 { eLT, eLT, eLT, eLT, eEQ }, /* -Inf */ 592 }; 593 594 void 595 rxfp_cmp (fp_t fa, fp_t fb) 596 { 597 FP_Parts a, b; 598 fp_t c; 599 FP_ExceptionCases reason; 600 int flags = 0; 601 double da, db; 602 603 fp_explode (fa, &a); 604 fp_explode (fb, &b); 605 606 if (check_exceptions (&a, &b, &c, ex_cmp_tab, &reason)) 607 { 608 if (reason == eQN) 609 { 610 /* Special case - incomparable. */ 611 set_flags (FLAGBIT_Z | FLAGBIT_S | FLAGBIT_O, FLAGBIT_O); 612 return; 613 } 614 return; 615 } 616 617 switch (reason) 618 { 619 case eEQ: 620 flags = FLAGBIT_Z; 621 break; 622 case eLT: 623 flags = FLAGBIT_S; 624 break; 625 case eGT: 626 flags = 0; 627 break; 628 case eNR: 629 da = fp_to_double (&a); 630 db = fp_to_double (&b); 631 tprintf("fcmp: %g cmp %g\n", da, db); 632 if (da < db) 633 flags = FLAGBIT_S; 634 else if (da == db) 635 flags = FLAGBIT_Z; 636 else 637 flags = 0; 638 break; 639 default: 640 abort(); 641 } 642 643 set_flags (FLAGBIT_Z | FLAGBIT_S | FLAGBIT_O, flags); 644 } 645 646 long 647 rxfp_ftoi (fp_t fa, int round_mode) 648 { 649 FP_Parts a; 650 fp_t rv; 651 int sign; 652 int whole_bits, frac_bits; 653 654 fp_explode (fa, &a); 655 sign = fa & 0x80000000UL; 656 657 switch (a.type) 658 { 659 case FP_NORMAL: 660 break; 661 case FP_PZERO: 662 case FP_NZERO: 663 return 0; 664 case FP_PINFINITY: 665 FP_RAISE (V); 666 return 0x7fffffffL; 667 case FP_NINFINITY: 668 FP_RAISE (V); 669 return 0x80000000L; 670 case FP_DENORMAL: 671 FP_RAISE (E); 672 return 0; 673 case FP_QNAN: 674 case FP_SNAN: 675 FP_RAISE (V); 676 return sign ? 0x80000000U : 0x7fffffff; 677 } 678 679 if (a.exp >= 31) 680 { 681 FP_RAISE (V); 682 return sign ? 0x80000000U : 0x7fffffff; 683 } 684 685 a.exp -= 23; 686 687 if (a.exp <= -25) 688 { 689 /* Less than 0.49999 */ 690 frac_bits = a.mant; 691 whole_bits = 0; 692 } 693 else if (a.exp < 0) 694 { 695 frac_bits = a.mant << (32 + a.exp); 696 whole_bits = a.mant >> (-a.exp); 697 } 698 else 699 { 700 frac_bits = 0; 701 whole_bits = a.mant << a.exp; 702 } 703 704 if (frac_bits) 705 { 706 switch (round_mode & 3) 707 { 708 case FPRM_NEAREST: 709 if (frac_bits & 0x80000000UL) 710 whole_bits ++; 711 break; 712 case FPRM_ZERO: 713 break; 714 case FPRM_PINF: 715 if (!sign) 716 whole_bits ++; 717 break; 718 case FPRM_NINF: 719 if (sign) 720 whole_bits ++; 721 break; 722 } 723 } 724 725 rv = sign ? -whole_bits : whole_bits; 726 727 return rv; 728 } 729 730 fp_t 731 rxfp_itof (long fa, int round_mode) 732 { 733 fp_t rv; 734 int sign = 0; 735 unsigned int frac_bits; 736 volatile unsigned int whole_bits; 737 FP_Parts a = {0}; 738 739 if (fa == 0) 740 return PLUS_ZERO; 741 742 if (fa < 0) 743 { 744 fa = -fa; 745 sign = 1; 746 a.sign = -1; 747 } 748 else 749 a.sign = 1; 750 751 whole_bits = fa; 752 a.exp = 31; 753 754 while (! (whole_bits & 0x80000000UL)) 755 { 756 a.exp --; 757 whole_bits <<= 1; 758 } 759 frac_bits = whole_bits & 0xff; 760 whole_bits = whole_bits >> 8; 761 762 if (frac_bits) 763 { 764 /* We must round */ 765 switch (round_mode & 3) 766 { 767 case FPRM_NEAREST: 768 if (frac_bits & 0x80) 769 whole_bits ++; 770 break; 771 case FPRM_ZERO: 772 break; 773 case FPRM_PINF: 774 if (!sign) 775 whole_bits ++; 776 break; 777 case FPRM_NINF: 778 if (sign) 779 whole_bits ++; 780 break; 781 } 782 } 783 784 a.mant = whole_bits; 785 if (whole_bits & 0xff000000UL) 786 { 787 a.mant >>= 1; 788 a.exp ++; 789 } 790 791 rv = fp_implode (&a); 792 return rv; 793 } 794 795