arch/m68k/math-emu/fp_arith.c

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