arch/m68k/math-emu/fp_arith.c

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