1 /* 2 3 fp_arith.c: floating-point math routines for the Linux-m68k 4 floating point emulator. 5 6 Copyright (c) 1998-1999 David Huggins-Daines. 7 8 Somewhat based on the AlphaLinux floating point emulator, by David 9 Mosberger-Tang. 10 11 You may copy, modify, and redistribute this file under the terms of 12 the GNU General Public License, version 2, or any later version, at 13 your convenience. 14 */ 15 16 #include "fp_emu.h" 17 #include "multi_arith.h" 18 #include "fp_arith.h" 19 20 const struct fp_ext fp_QNaN = 21 { 22 .exp = 0x7fff, 23 .mant = { .m64 = ~0 } 24 }; 25 26 const struct fp_ext fp_Inf = 27 { 28 .exp = 0x7fff, 29 }; 30 31 /* let's start with the easy ones */ 32 33 struct fp_ext * 34 fp_fabs(struct fp_ext *dest, struct fp_ext *src) 35 { 36 dprint(PINSTR, "fabs\n"); 37 38 fp_monadic_check(dest, src); 39 40 dest->sign = 0; 41 42 return dest; 43 } 44 45 struct fp_ext * 46 fp_fneg(struct fp_ext *dest, struct fp_ext *src) 47 { 48 dprint(PINSTR, "fneg\n"); 49 50 fp_monadic_check(dest, src); 51 52 dest->sign = !dest->sign; 53 54 return dest; 55 } 56 57 /* Now, the slightly harder ones */ 58 59 /* fp_fadd: Implements the kernel of the FADD, FSADD, FDADD, FSUB, 60 FDSUB, and FCMP instructions. */ 61 62 struct fp_ext * 63 fp_fadd(struct fp_ext *dest, struct fp_ext *src) 64 { 65 int diff; 66 67 dprint(PINSTR, "fadd\n"); 68 69 fp_dyadic_check(dest, src); 70 71 if (IS_INF(dest)) { 72 /* infinity - infinity == NaN */ 73 if (IS_INF(src) && (src->sign != dest->sign)) 74 fp_set_nan(dest); 75 return dest; 76 } 77 if (IS_INF(src)) { 78 fp_copy_ext(dest, src); 79 return dest; 80 } 81 82 if (IS_ZERO(dest)) { 83 if (IS_ZERO(src)) { 84 if (src->sign != dest->sign) { 85 if (FPDATA->rnd == FPCR_ROUND_RM) 86 dest->sign = 1; 87 else 88 dest->sign = 0; 89 } 90 } else 91 fp_copy_ext(dest, src); 92 return dest; 93 } 94 95 dest->lowmant = src->lowmant = 0; 96 97 if ((diff = dest->exp - src->exp) > 0) 98 fp_denormalize(src, diff); 99 else if ((diff = -diff) > 0) 100 fp_denormalize(dest, diff); 101 102 if (dest->sign == src->sign) { 103 if (fp_addmant(dest, src)) 104 if (!fp_addcarry(dest)) 105 return dest; 106 } else { 107 if (dest->mant.m64 < src->mant.m64) { 108 fp_submant(dest, src, dest); 109 dest->sign = !dest->sign; 110 } else 111 fp_submant(dest, dest, src); 112 } 113 114 return dest; 115 } 116 117 /* fp_fsub: Implements the kernel of the FSUB, FSSUB, and FDSUB 118 instructions. 119 120 Remember that the arguments are in assembler-syntax order! */ 121 122 struct fp_ext * 123 fp_fsub(struct fp_ext *dest, struct fp_ext *src) 124 { 125 dprint(PINSTR, "fsub "); 126 127 src->sign = !src->sign; 128 return fp_fadd(dest, src); 129 } 130 131 132 struct fp_ext * 133 fp_fcmp(struct fp_ext *dest, struct fp_ext *src) 134 { 135 dprint(PINSTR, "fcmp "); 136 137 FPDATA->temp[1] = *dest; 138 src->sign = !src->sign; 139 return fp_fadd(&FPDATA->temp[1], src); 140 } 141 142 struct fp_ext * 143 fp_ftst(struct fp_ext *dest, struct fp_ext *src) 144 { 145 dprint(PINSTR, "ftst\n"); 146 147 (void)dest; 148 149 return src; 150 } 151 152 struct fp_ext * 153 fp_fmul(struct fp_ext *dest, struct fp_ext *src) 154 { 155 union fp_mant128 temp; 156 int exp; 157 158 dprint(PINSTR, "fmul\n"); 159 160 fp_dyadic_check(dest, src); 161 162 /* calculate the correct sign now, as it's necessary for infinities */ 163 dest->sign = src->sign ^ dest->sign; 164 165 /* Handle infinities */ 166 if (IS_INF(dest)) { 167 if (IS_ZERO(src)) 168 fp_set_nan(dest); 169 return dest; 170 } 171 if (IS_INF(src)) { 172 if (IS_ZERO(dest)) 173 fp_set_nan(dest); 174 else 175 fp_copy_ext(dest, src); 176 return dest; 177 } 178 179 /* Of course, as we all know, zero * anything = zero. You may 180 not have known that it might be a positive or negative 181 zero... */ 182 if (IS_ZERO(dest) || IS_ZERO(src)) { 183 dest->exp = 0; 184 dest->mant.m64 = 0; 185 dest->lowmant = 0; 186 187 return dest; 188 } 189 190 exp = dest->exp + src->exp - 0x3ffe; 191 192 /* shift up the mantissa for denormalized numbers, 193 so that the highest bit is set, this makes the 194 shift of the result below easier */ 195 if ((long)dest->mant.m32[0] >= 0) 196 exp -= fp_overnormalize(dest); 197 if ((long)src->mant.m32[0] >= 0) 198 exp -= fp_overnormalize(src); 199 200 /* now, do a 64-bit multiply with expansion */ 201 fp_multiplymant(&temp, dest, src); 202 203 /* normalize it back to 64 bits and stuff it back into the 204 destination struct */ 205 if ((long)temp.m32[0] > 0) { 206 exp--; 207 fp_putmant128(dest, &temp, 1); 208 } else 209 fp_putmant128(dest, &temp, 0); 210 211 if (exp >= 0x7fff) { 212 fp_set_ovrflw(dest); 213 return dest; 214 } 215 dest->exp = exp; 216 if (exp < 0) { 217 fp_set_sr(FPSR_EXC_UNFL); 218 fp_denormalize(dest, -exp); 219 } 220 221 return dest; 222 } 223 224 /* fp_fdiv: Implements the "kernel" of the FDIV, FSDIV, FDDIV and 225 FSGLDIV instructions. 226 227 Note that the order of the operands is counter-intuitive: instead 228 of src / dest, the result is actually dest / src. */ 229 230 struct fp_ext * 231 fp_fdiv(struct fp_ext *dest, struct fp_ext *src) 232 { 233 union fp_mant128 temp; 234 int exp; 235 236 dprint(PINSTR, "fdiv\n"); 237 238 fp_dyadic_check(dest, src); 239 240 /* calculate the correct sign now, as it's necessary for infinities */ 241 dest->sign = src->sign ^ dest->sign; 242 243 /* Handle infinities */ 244 if (IS_INF(dest)) { 245 /* infinity / infinity = NaN (quiet, as always) */ 246 if (IS_INF(src)) 247 fp_set_nan(dest); 248 /* infinity / anything else = infinity (with approprate sign) */ 249 return dest; 250 } 251 if (IS_INF(src)) { 252 /* anything / infinity = zero (with appropriate sign) */ 253 dest->exp = 0; 254 dest->mant.m64 = 0; 255 dest->lowmant = 0; 256 257 return dest; 258 } 259 260 /* zeroes */ 261 if (IS_ZERO(dest)) { 262 /* zero / zero = NaN */ 263 if (IS_ZERO(src)) 264 fp_set_nan(dest); 265 /* zero / anything else = zero */ 266 return dest; 267 } 268 if (IS_ZERO(src)) { 269 /* anything / zero = infinity (with appropriate sign) */ 270 fp_set_sr(FPSR_EXC_DZ); 271 dest->exp = 0x7fff; 272 dest->mant.m64 = 0; 273 274 return dest; 275 } 276 277 exp = dest->exp - src->exp + 0x3fff; 278 279 /* shift up the mantissa for denormalized numbers, 280 so that the highest bit is set, this makes lots 281 of things below easier */ 282 if ((long)dest->mant.m32[0] >= 0) 283 exp -= fp_overnormalize(dest); 284 if ((long)src->mant.m32[0] >= 0) 285 exp -= fp_overnormalize(src); 286 287 /* now, do the 64-bit divide */ 288 fp_dividemant(&temp, dest, src); 289 290 /* normalize it back to 64 bits and stuff it back into the 291 destination struct */ 292 if (!temp.m32[0]) { 293 exp--; 294 fp_putmant128(dest, &temp, 32); 295 } else 296 fp_putmant128(dest, &temp, 31); 297 298 if (exp >= 0x7fff) { 299 fp_set_ovrflw(dest); 300 return dest; 301 } 302 dest->exp = exp; 303 if (exp < 0) { 304 fp_set_sr(FPSR_EXC_UNFL); 305 fp_denormalize(dest, -exp); 306 } 307 308 return dest; 309 } 310 311 struct fp_ext * 312 fp_fsglmul(struct fp_ext *dest, struct fp_ext *src) 313 { 314 int exp; 315 316 dprint(PINSTR, "fsglmul\n"); 317 318 fp_dyadic_check(dest, src); 319 320 /* calculate the correct sign now, as it's necessary for infinities */ 321 dest->sign = src->sign ^ dest->sign; 322 323 /* Handle infinities */ 324 if (IS_INF(dest)) { 325 if (IS_ZERO(src)) 326 fp_set_nan(dest); 327 return dest; 328 } 329 if (IS_INF(src)) { 330 if (IS_ZERO(dest)) 331 fp_set_nan(dest); 332 else 333 fp_copy_ext(dest, src); 334 return dest; 335 } 336 337 /* Of course, as we all know, zero * anything = zero. You may 338 not have known that it might be a positive or negative 339 zero... */ 340 if (IS_ZERO(dest) || IS_ZERO(src)) { 341 dest->exp = 0; 342 dest->mant.m64 = 0; 343 dest->lowmant = 0; 344 345 return dest; 346 } 347 348 exp = dest->exp + src->exp - 0x3ffe; 349 350 /* do a 32-bit multiply */ 351 fp_mul64(dest->mant.m32[0], dest->mant.m32[1], 352 dest->mant.m32[0] & 0xffffff00, 353 src->mant.m32[0] & 0xffffff00); 354 355 if (exp >= 0x7fff) { 356 fp_set_ovrflw(dest); 357 return dest; 358 } 359 dest->exp = exp; 360 if (exp < 0) { 361 fp_set_sr(FPSR_EXC_UNFL); 362 fp_denormalize(dest, -exp); 363 } 364 365 return dest; 366 } 367 368 struct fp_ext * 369 fp_fsgldiv(struct fp_ext *dest, struct fp_ext *src) 370 { 371 int exp; 372 unsigned long quot, rem; 373 374 dprint(PINSTR, "fsgldiv\n"); 375 376 fp_dyadic_check(dest, src); 377 378 /* calculate the correct sign now, as it's necessary for infinities */ 379 dest->sign = src->sign ^ dest->sign; 380 381 /* Handle infinities */ 382 if (IS_INF(dest)) { 383 /* infinity / infinity = NaN (quiet, as always) */ 384 if (IS_INF(src)) 385 fp_set_nan(dest); 386 /* infinity / anything else = infinity (with approprate sign) */ 387 return dest; 388 } 389 if (IS_INF(src)) { 390 /* anything / infinity = zero (with appropriate sign) */ 391 dest->exp = 0; 392 dest->mant.m64 = 0; 393 dest->lowmant = 0; 394 395 return dest; 396 } 397 398 /* zeroes */ 399 if (IS_ZERO(dest)) { 400 /* zero / zero = NaN */ 401 if (IS_ZERO(src)) 402 fp_set_nan(dest); 403 /* zero / anything else = zero */ 404 return dest; 405 } 406 if (IS_ZERO(src)) { 407 /* anything / zero = infinity (with appropriate sign) */ 408 fp_set_sr(FPSR_EXC_DZ); 409 dest->exp = 0x7fff; 410 dest->mant.m64 = 0; 411 412 return dest; 413 } 414 415 exp = dest->exp - src->exp + 0x3fff; 416 417 dest->mant.m32[0] &= 0xffffff00; 418 src->mant.m32[0] &= 0xffffff00; 419 420 /* do the 32-bit divide */ 421 if (dest->mant.m32[0] >= src->mant.m32[0]) { 422 fp_sub64(dest->mant, src->mant); 423 fp_div64(quot, rem, dest->mant.m32[0], 0, src->mant.m32[0]); 424 dest->mant.m32[0] = 0x80000000 | (quot >> 1); 425 dest->mant.m32[1] = (quot & 1) | rem; /* only for rounding */ 426 } else { 427 fp_div64(quot, rem, dest->mant.m32[0], 0, src->mant.m32[0]); 428 dest->mant.m32[0] = quot; 429 dest->mant.m32[1] = rem; /* only for rounding */ 430 exp--; 431 } 432 433 if (exp >= 0x7fff) { 434 fp_set_ovrflw(dest); 435 return dest; 436 } 437 dest->exp = exp; 438 if (exp < 0) { 439 fp_set_sr(FPSR_EXC_UNFL); 440 fp_denormalize(dest, -exp); 441 } 442 443 return dest; 444 } 445 446 /* fp_roundint: Internal rounding function for use by several of these 447 emulated instructions. 448 449 This one rounds off the fractional part using the rounding mode 450 specified. */ 451 452 static void fp_roundint(struct fp_ext *dest, int mode) 453 { 454 union fp_mant64 oldmant; 455 unsigned long mask; 456 457 if (!fp_normalize_ext(dest)) 458 return; 459 460 /* infinities and zeroes */ 461 if (IS_INF(dest) || IS_ZERO(dest)) 462 return; 463 464 /* first truncate the lower bits */ 465 oldmant = dest->mant; 466 switch (dest->exp) { 467 case 0 ... 0x3ffe: 468 dest->mant.m64 = 0; 469 break; 470 case 0x3fff ... 0x401e: 471 dest->mant.m32[0] &= 0xffffffffU << (0x401e - dest->exp); 472 dest->mant.m32[1] = 0; 473 if (oldmant.m64 == dest->mant.m64) 474 return; 475 break; 476 case 0x401f ... 0x403e: 477 dest->mant.m32[1] &= 0xffffffffU << (0x403e - dest->exp); 478 if (oldmant.m32[1] == dest->mant.m32[1]) 479 return; 480 break; 481 default: 482 return; 483 } 484 fp_set_sr(FPSR_EXC_INEX2); 485 486 /* We might want to normalize upwards here... however, since 487 we know that this is only called on the output of fp_fdiv, 488 or with the input to fp_fint or fp_fintrz, and the inputs 489 to all these functions are either normal or denormalized 490 (no subnormals allowed!), there's really no need. 491 492 In the case of fp_fdiv, observe that 0x80000000 / 0xffff = 493 0xffff8000, and the same holds for 128-bit / 64-bit. (i.e. the 494 smallest possible normal dividend and the largest possible normal 495 divisor will still produce a normal quotient, therefore, (normal 496 << 64) / normal is normal in all cases) */ 497 498 switch (mode) { 499 case FPCR_ROUND_RN: 500 switch (dest->exp) { 501 case 0 ... 0x3ffd: 502 return; 503 case 0x3ffe: 504 /* As noted above, the input is always normal, so the 505 guard bit (bit 63) is always set. therefore, the 506 only case in which we will NOT round to 1.0 is when 507 the input is exactly 0.5. */ 508 if (oldmant.m64 == (1ULL << 63)) 509 return; 510 break; 511 case 0x3fff ... 0x401d: 512 mask = 1 << (0x401d - dest->exp); 513 if (!(oldmant.m32[0] & mask)) 514 return; 515 if (oldmant.m32[0] & (mask << 1)) 516 break; 517 if (!(oldmant.m32[0] << (dest->exp - 0x3ffd)) && 518 !oldmant.m32[1]) 519 return; 520 break; 521 case 0x401e: 522 if (oldmant.m32[1] & 0x80000000) 523 return; 524 if (oldmant.m32[0] & 1) 525 break; 526 if (!(oldmant.m32[1] << 1)) 527 return; 528 break; 529 case 0x401f ... 0x403d: 530 mask = 1 << (0x403d - dest->exp); 531 if (!(oldmant.m32[1] & mask)) 532 return; 533 if (oldmant.m32[1] & (mask << 1)) 534 break; 535 if (!(oldmant.m32[1] << (dest->exp - 0x401d))) 536 return; 537 break; 538 default: 539 return; 540 } 541 break; 542 case FPCR_ROUND_RZ: 543 return; 544 default: 545 if (dest->sign ^ (mode - FPCR_ROUND_RM)) 546 break; 547 return; 548 } 549 550 switch (dest->exp) { 551 case 0 ... 0x3ffe: 552 dest->exp = 0x3fff; 553 dest->mant.m64 = 1ULL << 63; 554 break; 555 case 0x3fff ... 0x401e: 556 mask = 1 << (0x401e - dest->exp); 557 if (dest->mant.m32[0] += mask) 558 break; 559 dest->mant.m32[0] = 0x80000000; 560 dest->exp++; 561 break; 562 case 0x401f ... 0x403e: 563 mask = 1 << (0x403e - dest->exp); 564 if (dest->mant.m32[1] += mask) 565 break; 566 if (dest->mant.m32[0] += 1) 567 break; 568 dest->mant.m32[0] = 0x80000000; 569 dest->exp++; 570 break; 571 } 572 } 573 574 /* modrem_kernel: Implementation of the FREM and FMOD instructions 575 (which are exactly the same, except for the rounding used on the 576 intermediate value) */ 577 578 static struct fp_ext * 579 modrem_kernel(struct fp_ext *dest, struct fp_ext *src, int mode) 580 { 581 struct fp_ext tmp; 582 583 fp_dyadic_check(dest, src); 584 585 /* Infinities and zeros */ 586 if (IS_INF(dest) || IS_ZERO(src)) { 587 fp_set_nan(dest); 588 return dest; 589 } 590 if (IS_ZERO(dest) || IS_INF(src)) 591 return dest; 592 593 /* FIXME: there is almost certainly a smarter way to do this */ 594 fp_copy_ext(&tmp, dest); 595 fp_fdiv(&tmp, src); /* NOTE: src might be modified */ 596 fp_roundint(&tmp, mode); 597 fp_fmul(&tmp, src); 598 fp_fsub(dest, &tmp); 599 600 /* set the quotient byte */ 601 fp_set_quotient((dest->mant.m64 & 0x7f) | (dest->sign << 7)); 602 return dest; 603 } 604 605 /* fp_fmod: Implements the kernel of the FMOD instruction. 606 607 Again, the argument order is backwards. The result, as defined in 608 the Motorola manuals, is: 609 610 fmod(src,dest) = (dest - (src * floor(dest / src))) */ 611 612 struct fp_ext * 613 fp_fmod(struct fp_ext *dest, struct fp_ext *src) 614 { 615 dprint(PINSTR, "fmod\n"); 616 return modrem_kernel(dest, src, FPCR_ROUND_RZ); 617 } 618 619 /* fp_frem: Implements the kernel of the FREM instruction. 620 621 frem(src,dest) = (dest - (src * round(dest / src))) 622 */ 623 624 struct fp_ext * 625 fp_frem(struct fp_ext *dest, struct fp_ext *src) 626 { 627 dprint(PINSTR, "frem\n"); 628 return modrem_kernel(dest, src, FPCR_ROUND_RN); 629 } 630 631 struct fp_ext * 632 fp_fint(struct fp_ext *dest, struct fp_ext *src) 633 { 634 dprint(PINSTR, "fint\n"); 635 636 fp_copy_ext(dest, src); 637 638 fp_roundint(dest, FPDATA->rnd); 639 640 return dest; 641 } 642 643 struct fp_ext * 644 fp_fintrz(struct fp_ext *dest, struct fp_ext *src) 645 { 646 dprint(PINSTR, "fintrz\n"); 647 648 fp_copy_ext(dest, src); 649 650 fp_roundint(dest, FPCR_ROUND_RZ); 651 652 return dest; 653 } 654 655 struct fp_ext * 656 fp_fscale(struct fp_ext *dest, struct fp_ext *src) 657 { 658 int scale, oldround; 659 660 dprint(PINSTR, "fscale\n"); 661 662 fp_dyadic_check(dest, src); 663 664 /* Infinities */ 665 if (IS_INF(src)) { 666 fp_set_nan(dest); 667 return dest; 668 } 669 if (IS_INF(dest)) 670 return dest; 671 672 /* zeroes */ 673 if (IS_ZERO(src) || IS_ZERO(dest)) 674 return dest; 675 676 /* Source exponent out of range */ 677 if (src->exp >= 0x400c) { 678 fp_set_ovrflw(dest); 679 return dest; 680 } 681 682 /* src must be rounded with round to zero. */ 683 oldround = FPDATA->rnd; 684 FPDATA->rnd = FPCR_ROUND_RZ; 685 scale = fp_conv_ext2long(src); 686 FPDATA->rnd = oldround; 687 688 /* new exponent */ 689 scale += dest->exp; 690 691 if (scale >= 0x7fff) { 692 fp_set_ovrflw(dest); 693 } else if (scale <= 0) { 694 fp_set_sr(FPSR_EXC_UNFL); 695 fp_denormalize(dest, -scale); 696 } else 697 dest->exp = scale; 698 699 return dest; 700 } 701 702