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