arch/alpha/math-emu/ieee-math.c

/* */
This source file includes following definitions.
sign
cmp128
sll128
srl128
add128
sub128
mul64
div128
normalize
ieee_fpclass
extend_ieee
make_s_ieee
make_t_ieee
round_s_ieee
round_t_ieee
add_kernel_ieee
ieee_CVTST
ieee_CVTTS
ieee_CVTQS
ieee_CVTQT
ieee_CVTTQ
ieee_CMPTEQ
ieee_CMPTLT
ieee_CMPTLE
ieee_CMPTUN
ieee_ADDS
ieee_ADDT
ieee_SUBS
ieee_SUBT
ieee_MULS
ieee_MULT
ieee_DIVS
ieee_DIVT
   1 /*
   2  * ieee-math.c - IEEE floating point emulation code
   3  * Copyright (C) 1989,1990,1991,1995 by
   4  * Digital Equipment Corporation, Maynard, Massachusetts.
   5  *
   6  * Heavily modified for Linux/Alpha.  Changes are Copyright (c) 1995
   7  * by David Mosberger (davidm@azstarnet.com).
   8  *
   9  * This file may be redistributed according to the terms of the
  10  * GNU General Public License.
  11  */
  12 /*
  13  * The original code did not have any comments. I have created many
  14  * comments as I fix the bugs in the code.  My comments are based on
  15  * my observation and interpretation of the code.  If the original
  16  * author would have spend a few minutes to comment the code, we would
  17  * never had a problem of misinterpretation.  -HA
  18  *
  19  * This code could probably be a lot more optimized (especially the
  20  * division routine).  However, my foremost concern was to get the
  21  * IEEE behavior right.  Performance is less critical as these
  22  * functions are used on exceptional numbers only (well, assuming you
  23  * don't turn on the "trap on inexact"...).
  24  */
  25 #include "ieee-math.h"
  26 
  27 #define STICKY_S        0x20000000      /* both in longword 0 of fraction */
  28 #define STICKY_T        1
  29 
  30 /*
  31  * Careful: order matters here!
  32  */
  33 enum {
  34         NaN, QNaN, INFTY, ZERO, DENORM, NORMAL
  35 };
  36 
  37 enum {
  38         SINGLE, DOUBLE
  39 };
  40 
  41 typedef unsigned long fpclass_t;
  42 
  43 #define IEEE_TMAX       0x7fefffffffffffff
  44 #define IEEE_SMAX       0x47efffffe0000000
  45 #define IEEE_SNaN       0xfff00000000f0000
  46 #define IEEE_QNaN       0xfff8000000000000
  47 #define IEEE_PINF       0x7ff0000000000000
  48 #define IEEE_NINF       0xfff0000000000000
  49 
  50 
  51 /*
  52  * The memory format of S floating point numbers differs from the
  53  * register format.  In the following, the bitnumbers above the
  54  * diagram below give the memory format while the numbers below give
  55  * the register format.
  56  *
  57  *        31 30      23 22                              0
  58  *      +-----------------------------------------------+
  59  * S    | s |   exp    |       fraction                 |
  60  *      +-----------------------------------------------+
  61  *        63 62      52 51                            29
  62  *      
  63  * For T floating point numbers, the register and memory formats
  64  * match:
  65  *
  66  *      +-------------------------------------------------------------------+
  67  * T    | s |        exp        |           frac | tion                     |
  68  *      +-------------------------------------------------------------------+
  69  *        63 62               52 51            32 31                       0
  70  */
  71 typedef struct {
  72         unsigned long   f[2];   /* bit 55 in f[0] is the factor of 2^0*/
  73         int             s;      /* 1 bit sign (0 for +, 1 for -) */
  74         int             e;      /* 16 bit signed exponent */
  75 } EXTENDED;
  76 
  77 
  78 /*
  79  * Return the sign of a Q integer, S or T fp number in the register
  80  * format.
  81  */
  82 static inline int
  83 sign (unsigned long a)
     /*  */
  84 {
  85         if ((long) a < 0)
  86                 return -1;
  87         else
  88                 return 1;
  89 }
  90 
  91 
  92 static inline long
  93 cmp128 (const long a[2], const long b[2])
     /*  */
  94 {
  95         if (a[1] < b[1]) return -1;
  96         if (a[1] > b[1]) return  1;
  97         return a[0] - b[0];
  98 }
  99 
 100 
 101 static inline void
 102 sll128 (unsigned long a[2])
     /*  */
 103 {
 104         a[1] = (a[1] << 1) | (a[0] >> 63);
 105         a[0] <<= 1;
 106 }
 107 
 108 
 109 static inline void
 110 srl128 (unsigned long a[2])
     /*  */
 111 {
 112         a[0] = (a[0] >> 1) | (a[1] << 63);
 113         a[1] >>= 1;
 114 }
 115 
 116 
 117 static inline void
 118 add128 (const unsigned long a[2], const unsigned long b[2], unsigned long c[2])
     /*  */
 119 {
 120         unsigned long carry = a[0] > (0xffffffffffffffff - b[0]);
 121 
 122         c[0] = a[0] + b[0];
 123         c[1] = a[1] + b[1] + carry;
 124 }
 125 
 126 
 127 static inline void
 128 sub128 (const unsigned long a[2], const unsigned long b[2], unsigned long c[2])
     /*  */
 129 {
 130         unsigned long borrow = a[0] < b[0];
 131 
 132         c[0] = a[0] - b[0];
 133         c[1] = a[1] - b[1] - borrow;
 134 }
 135 
 136 
 137 static inline void
 138 mul64 (const unsigned long a, const unsigned long b, unsigned long c[2])
     /*  */
 139 {
 140         asm ("mulq  %2,%3,%0\n\t"
 141              "umulh %2,%3,%1"
 142              : "r="(c[0]), "r="(c[1]) : "r"(a), "r"(b));
 143 }
 144 
 145 
 146 static void
 147 div128 (unsigned long a[2], unsigned long b[2], unsigned long c[2])
     /*  */
 148 {
 149         unsigned long mask[2] = {1, 0};
 150 
 151         /*
 152          * Shift b until either the sign bit is set or until it is at
 153          * least as big as the dividend:
 154          */
 155         while (cmp128(b, a) < 0 && sign(b[1]) >= 0) {
 156                 sll128(b);
 157                 sll128(mask);
 158         }
 159         c[0] = c[1] = 0;
 160         do {
 161                 if (cmp128(a, b) >= 0) {
 162                         sub128(a, b, a);
 163                         add128(mask, c, c);
 164                 }
 165                 srl128(mask);
 166                 srl128(b);
 167         } while (mask[0] || mask[1]);
 168 }
 169 
 170 
 171 static void
 172 normalize (EXTENDED *a)
     /*  */
 173 {
 174         if (!a->f[0] && !a->f[1])
 175                 return;         /* zero fraction, unnormalizable... */
 176         /*
 177          * In "extended" format, the "1" in "1.f" is explicit; it is
 178          * in bit 55 of f[0], and the decimal point is understood to
 179          * be between bit 55 and bit 54.  To normalize, shift the
 180          * fraction until we have a "1" in bit 55.
 181          */
 182         if ((a->f[0] & 0xff00000000000000) != 0 || a->f[1] != 0) {
 183                 /*
 184                  * Mantissa is greater than 1.0:
 185                  */
 186                 while ((a->f[0] & 0xff80000000000000) != 0x0080000000000000 ||
 187                        a->f[1] != 0)
 188                 {
 189                         unsigned long sticky;
 190 
 191                         ++a->e;
 192                         sticky = a->f[0] & 1;
 193                         srl128(a->f);
 194                         a->f[0] |= sticky;
 195                 }
 196                 return;
 197         }
 198 
 199         if (!(a->f[0] & 0x0080000000000000)) {
 200                 /*
 201                  * Mantissa is less than 1.0:
 202                  */
 203                 while (!(a->f[0] & 0x0080000000000000)) {
 204                         --a->e;
 205                         a->f[0] <<= 1;
 206                 }
 207                 return;
 208         }
 209 }
 210 
 211 
 212 static inline fpclass_t
 213 ieee_fpclass (unsigned long a)
     /*  */
 214 {
 215         unsigned long exp, fract;
 216 
 217         exp   = (a >> 52) & 0x7ff;      /* 11 bits of exponent */
 218         fract = a & 0x000fffffffffffff; /* 52 bits of fraction */
 219         if (exp == 0) {
 220                 if (fract == 0)
 221                         return ZERO;
 222                 return DENORM;
 223         }
 224         if (exp == 0x7ff) {
 225                 if (fract == 0)
 226                         return INFTY;
 227                 if (((fract >> 51) & 1) != 0)
 228                         return QNaN;
 229                 return NaN;
 230         }
 231         return NORMAL;
 232 }
 233 
 234 
 235 /*
 236  * Translate S/T fp number in register format into extended format.
 237  */
 238 static fpclass_t
 239 extend_ieee (unsigned long a, EXTENDED *b, int prec)
     /*  */
 240 {
 241         fpclass_t result_kind;
 242 
 243         b->s = a >> 63;
 244         b->e = ((a >> 52) & 0x7ff) - 0x3ff;     /* remove bias */
 245         b->f[1] = 0;
 246         /*
 247          * We shift f[1] left three bits so that the higher order bits
 248          * of the fraction will reside in bits 55 through 0 of f[0].
 249          */
 250         b->f[0] = (a & 0x000fffffffffffff) << 3;
 251         result_kind = ieee_fpclass(a);
 252         if (result_kind == NORMAL) {
 253                 /* set implied 1. bit: */
 254                 b->f[0] |= 1UL << 55;
 255         } else if (result_kind == DENORM) {
 256                 if (prec == SINGLE)
 257                         b->e = -126;
 258                 else
 259                         b->e = -1022;
 260         }
 261         return result_kind;
 262 }
 263 
 264 
 265 /*
 266  * INPUT PARAMETERS:
 267  *       a           a number in EXTENDED format to be converted to
 268  *                   s-floating format.
 269  *       f           rounding mode and exception enable bits.
 270  * OUTPUT PARAMETERS:
 271  *       b          will contain the s-floating number that "a" was
 272  *                  converted to (in register format).
 273  */
 274 static unsigned long
 275 make_s_ieee (long f, EXTENDED *a, unsigned long *b)
     /*  */
 276 {
 277         unsigned long res, sticky;
 278 
 279         if (!a->f[0] && !a->f[1]) {
 280                 *b = (unsigned long) a->s << 63;        /* return +/-0 */
 281                 return 0;
 282         }
 283 
 284         normalize(a);
 285         res = 0;
 286 
 287         if (a->e < -0x7e) {
 288                 res = FPCR_INE;
 289                 if (f & IEEE_TRAP_ENABLE_UNF) {
 290                         res |= FPCR_UNF;
 291                         a->e += 0xc0;   /* scale up result by 2^alpha */
 292                 } else {
 293                         /* try making denormalized number: */
 294                         while (a->e < -0x7e) {
 295                                 ++a->e;
 296                                 sticky = a->f[0] & 1;
 297                                 srl128(a->f);
 298                                 if (!a->f[0] && !a->f[0]) {
 299                                         /* underflow: replace with exact 0 */
 300                                         res |= FPCR_UNF;
 301                                         break;
 302                                 }
 303                                 a->f[0] |= sticky;
 304                         }
 305                         a->e = -0x3ff;
 306                 }
 307         }
 308         if (a->e >= 0x80) {
 309                 res = FPCR_OVF | FPCR_INE;
 310                 if (f & IEEE_TRAP_ENABLE_OVF) {
 311                         a->e -= 0xc0;   /* scale down result by 2^alpha */
 312                 } else {
 313                         /*
 314                          * Overflow without trap enabled, substitute
 315                          * result according to rounding mode:
 316                          */
 317                         switch (RM(f)) {
 318                               case ROUND_NEAR:
 319                                 *b = IEEE_PINF;
 320                                 break;
 321 
 322                               case ROUND_CHOP:
 323                                 *b = IEEE_SMAX;
 324                                 break;
 325 
 326                               case ROUND_NINF:
 327                                 if (a->s) {
 328                                         *b = IEEE_PINF;
 329                                 } else {
 330                                         *b = IEEE_SMAX;
 331                                 }
 332                                 break;
 333 
 334                               case ROUND_PINF: 
 335                                 if (a->s) {
 336                                         *b = IEEE_SMAX;
 337                                 } else {
 338                                         *b = IEEE_PINF;
 339                                 }
 340                                 break;
 341                         }
 342                         *b |= ((unsigned long) a->s << 63);
 343                         return res;
 344                 }
 345         }
 346 
 347         *b = (((unsigned long) a->s << 63) |
 348               (((unsigned long) a->e + 0x3ff) << 52) |
 349               ((a->f[0] >> 3) & 0x000fffffe0000000));
 350         return res;
 351 }
 352 
 353 
 354 static unsigned long
 355 make_t_ieee (long f, EXTENDED *a, unsigned long *b)
     /*  */
 356 {
 357         unsigned long res, sticky;
 358 
 359         if (!a->f[0] && !a->f[1]) {
 360                 *b = (unsigned long) a->s << 63;        /* return +/-0 */
 361                 return 0;
 362         }
 363 
 364         normalize(a);
 365         res = 0;
 366         if (a->e < -0x3fe) {
 367                 res = FPCR_INE;
 368                 if (f & IEEE_TRAP_ENABLE_UNF) {
 369                         res |= FPCR_UNF;
 370                         a->e += 0x600;
 371                 } else {
 372                         /* try making denormalized number: */
 373                         while (a->e < -0x3fe) {
 374                                 ++a->e;
 375                                 sticky = a->f[0] & 1;
 376                                 srl128(a->f);
 377                                 if (!a->f[0] && !a->f[0]) {
 378                                         /* underflow: replace with exact 0 */
 379                                         res |= FPCR_UNF;
 380                                         break;
 381                                 }
 382                                 a->f[0] |= sticky;
 383                         }
 384                         a->e = -0x3ff;
 385                 }
 386         }
 387         if (a->e > 0x3ff) {
 388                 res = FPCR_OVF | FPCR_INE;
 389                 if (f & IEEE_TRAP_ENABLE_OVF) {
 390                         a->e -= 0x600;  /* scale down result by 2^alpha */
 391                 } else {
 392                         /*
 393                          * Overflow without trap enabled, substitute
 394                          * result according to rounding mode:
 395                          */
 396                         switch (RM(f)) {
 397                               case ROUND_NEAR:
 398                                 *b = IEEE_PINF;
 399                                 break;
 400 
 401                               case ROUND_CHOP:
 402                                 *b = IEEE_TMAX;
 403                                 break;
 404 
 405                               case ROUND_NINF:
 406                                 if (a->s) {
 407                                         *b = IEEE_PINF;
 408                                 } else {
 409                                         *b = IEEE_TMAX;
 410                                 }
 411                                 break;
 412 
 413                               case ROUND_PINF: 
 414                                 if (a->s) {
 415                                         *b = IEEE_TMAX;
 416                                 } else {
 417                                         *b = IEEE_PINF;
 418                                 }
 419                                 break;
 420                         }
 421                         *b |= ((unsigned long) a->s << 63);
 422                         return res;
 423                 }
 424         }
 425         *b = (((unsigned long) a->s << 63) |
 426               (((unsigned long) a->e + 0x3ff) << 52) |
 427               ((a->f[0] >> 3) & 0x000fffffffffffff));
 428         return res;
 429 }
 430 
 431 
 432 /*
 433  * INPUT PARAMETERS:
 434  *       a          EXTENDED format number to be rounded.
 435  *       rm         integer with value ROUND_NEAR, ROUND_CHOP, etc.
 436  *                  indicates how "a" should be rounded to produce "b".
 437  * OUTPUT PARAMETERS:
 438  *       b          s-floating number produced by rounding "a".
 439  * RETURN VALUE:
 440  *       if no errors occurred, will be zero.  Else will contain flags
 441  *       like FPCR_INE_OP, etc.
 442  */
 443 static unsigned long
 444 round_s_ieee (int f, EXTENDED *a, unsigned long *b)
     /*  */
 445 {
 446         unsigned long diff1, diff2, res = 0;
 447         EXTENDED z1, z2;
 448 
 449         if (!(a->f[0] & 0xffffffff)) {
 450                 return make_s_ieee(f, a, b);    /* no rounding error */
 451         }
 452 
 453         /*
 454          * z1 and z2 are the S-floating numbers with the next smaller/greater
 455          * magnitude than a, respectively.
 456          */
 457         z1.s = z2.s = a->s;
 458         z1.e = z2.e = a->e;
 459         z1.f[0] = z2.f[0] = a->f[0] & 0xffffffff00000000;
 460         z1.f[1] = z2.f[1] = 0;
 461         z2.f[0] += 0x100000000; /* next bigger S float number */
 462 
 463         switch (RM(f)) {
 464               case ROUND_NEAR:
 465                 diff1 = a->f[0] - z1.f[0];
 466                 diff2 = z2.f[0] - a->f[0];
 467                 if (diff1 > diff2)
 468                         res = make_s_ieee(f, &z2, b);
 469                 else if (diff2 > diff1)
 470                         res = make_s_ieee(f, &z1, b);
 471                 else
 472                         /* equal distance: round towards even */
 473                         if (z1.f[0] & 0x100000000)
 474                                 res = make_s_ieee(f, &z2, b);
 475                         else
 476                                 res = make_s_ieee(f, &z1, b);
 477                 break;
 478 
 479               case ROUND_CHOP:
 480                 res = make_s_ieee(f, &z1, b);
 481                 break;
 482 
 483               case ROUND_PINF:
 484                 if (a->s) {
 485                         res = make_s_ieee(f, &z1, b);
 486                 } else {
 487                         res = make_s_ieee(f, &z2, b);
 488                 }
 489                 break;
 490 
 491               case ROUND_NINF:
 492                 if (a->s) {
 493                         res = make_s_ieee(f, &z2, b);
 494                 } else {
 495                         res = make_s_ieee(f, &z1, b);
 496                 }
 497                 break;
 498         }
 499         return FPCR_INE | res;
 500 }
 501 
 502 
 503 static unsigned long
 504 round_t_ieee (int f, EXTENDED *a, unsigned long *b)
     /*  */
 505 {
 506         unsigned long diff1, diff2, res;
 507         EXTENDED z1, z2;
 508 
 509         if (!(a->f[0] & 0x7)) {
 510                 /* no rounding error */
 511                 return make_t_ieee(f, a, b);
 512         }
 513 
 514         z1.s = z2.s = a->s;
 515         z1.e = z2.e = a->e;
 516         z1.f[0] = z2.f[0] = a->f[0] & ~0x7;
 517         z1.f[1] = z2.f[1] = 0;
 518         z2.f[0] += (1 << 3);
 519 
 520         res = 0;
 521         switch (RM(f)) {
 522               case ROUND_NEAR:
 523                 diff1 = a->f[0] - z1.f[0];
 524                 diff2 = z2.f[0] - a->f[0];
 525                 if (diff1 > diff2)
 526                         res = make_t_ieee(f, &z2, b);
 527                 else if (diff2 > diff1)
 528                         res = make_t_ieee(f, &z1, b);
 529                 else
 530                         /* equal distance: round towards even */
 531                         if (z1.f[0] & (1 << 3))
 532                                 res = make_t_ieee(f, &z2, b);
 533                         else
 534                                 res = make_t_ieee(f, &z1, b);
 535                 break;
 536 
 537               case ROUND_CHOP:
 538                 res = make_t_ieee(f, &z1, b);
 539                 break;
 540 
 541               case ROUND_PINF:
 542                 if (a->s) {
 543                         res = make_t_ieee(f, &z1, b);
 544                 } else {
 545                         res = make_t_ieee(f, &z2, b);
 546                 }
 547                 break;
 548 
 549               case ROUND_NINF:
 550                 if (a->s) {
 551                         res = make_t_ieee(f, &z2, b);
 552                 } else {
 553                         res = make_t_ieee(f, &z1, b);
 554                 }
 555                 break;
 556         }
 557         return FPCR_INE | res;
 558 }
 559 
 560 
 561 static fpclass_t
 562 add_kernel_ieee (EXTENDED *op_a, EXTENDED *op_b, EXTENDED *op_c)
     /*  */
 563 {
 564         unsigned long mask, fa, fb, fc;
 565         int diff;
 566 
 567         diff = op_a->e - op_b->e;
 568         fa = op_a->f[0];
 569         fb = op_b->f[0];
 570         if (diff < 0) {
 571                 diff = -diff;
 572                 op_c->e = op_b->e;
 573                 mask = (1UL << diff) - 1;
 574                 fa >>= diff;
 575                 if (op_a->f[0] & mask) {
 576                         fa |= 1;                /* set sticky bit */
 577                 }
 578         } else {
 579                 op_c->e = op_a->e;
 580                 mask = (1UL << diff) - 1;
 581                 fb >>= diff;
 582                 if (op_b->f[0] & mask) {
 583                         fb |= 1;                /* set sticky bit */
 584                 }
 585         }
 586         if (op_a->s)
 587                 fa = -fa;
 588         if (op_b->s)
 589                 fb = -fb;
 590         fc = fa + fb;
 591         op_c->f[1] = 0;
 592         op_c->s = fc >> 63;
 593         if (op_c->s) {
 594                 fc = -fc;
 595         }
 596         op_c->f[0] = fc;
 597         normalize(op_c);
 598         return 0;
 599 }
 600 
 601 
 602 /*
 603  * converts s-floating "a" to t-floating "b".
 604  *
 605  * INPUT PARAMETERS:
 606  *       a           a s-floating number to be converted
 607  *       f           the rounding mode (ROUND_NEAR, etc. )
 608  * OUTPUT PARAMETERS:
 609  *       b           the t-floating number that "a" is converted to.
 610  * RETURN VALUE:
 611  *       error flags - i.e., zero if no errors occurred,
 612  *       FPCR_INV if invalid operation occurred, etc.
 613  */
 614 unsigned long
 615 ieee_CVTST (int f, unsigned long a, unsigned long *b)
     /*  */
 616 {
 617         EXTENDED temp;
 618         fpclass_t a_type;
 619 
 620         a_type = extend_ieee(a, &temp, SINGLE);
 621         if (a_type >= NaN && a_type <= INFTY) {
 622                 *b = a;
 623                 if (a_type == NaN) {
 624                         *b |= (1UL << 51);      /* turn SNaN into QNaN */
 625                         return FPCR_INV;
 626                 }
 627                 return 0;
 628         }
 629         return round_t_ieee(f, &temp, b);
 630 }
 631 
 632 
 633 /*
 634  * converts t-floating "a" to s-floating "b".
 635  *
 636  * INPUT PARAMETERS:
 637  *       a           a t-floating number to be converted
 638  *       f           the rounding mode (ROUND_NEAR, etc. )
 639  * OUTPUT PARAMETERS:
 640  *       b           the s-floating number that "a" is converted to.
 641  * RETURN VALUE:
 642  *       error flags - i.e., zero if no errors occurred,
 643  *       FPCR_INV if invalid operation occurred, etc.
 644  */
 645 unsigned long
 646 ieee_CVTTS (int f, unsigned long a, unsigned long *b)
     /*  */
 647 {
 648         EXTENDED temp;
 649         fpclass_t a_type;
 650 
 651         a_type = extend_ieee(a, &temp, DOUBLE);
 652         if (a_type >= NaN && a_type <= INFTY) {
 653                 *b = a;
 654                 if (a_type == NaN) {
 655                         *b |= (1UL << 51);      /* turn SNaN into QNaN */
 656                         return FPCR_INV;
 657                 }
 658                 return 0;
 659         }
 660         return round_s_ieee(f, &temp, b);
 661 }
 662 
 663 
 664 /*
 665  * converts q-format (64-bit integer) "a" to s-floating "b".
 666  *
 667  * INPUT PARAMETERS:
 668  *       a           an 64-bit integer to be converted.
 669  *       f           the rounding mode (ROUND_NEAR, etc. )
 670  * OUTPUT PARAMETERS:
 671  *       b           the s-floating number "a" is converted to.
 672  * RETURN VALUE:
 673  *       error flags - i.e., zero if no errors occurred,
 674  *       FPCR_INV if invalid operation occurred, etc.
 675  */
 676 unsigned long
 677 ieee_CVTQS (int f, unsigned long a, unsigned long *b)
     /*  */
 678 {
 679         EXTENDED op_b;
 680 
 681         op_b.s    = 0;
 682         op_b.f[0] = a;
 683         op_b.f[1] = 0;
 684         if (sign(a) < 0) {
 685                 op_b.s = 1;
 686                 op_b.f[0] = -a;
 687         }
 688         op_b.e = 55;
 689         normalize(&op_b);
 690         return round_s_ieee(f, &op_b, b);
 691 }
 692 
 693 
 694 /*
 695  * converts 64-bit integer "a" to t-floating "b".
 696  *
 697  * INPUT PARAMETERS:
 698  *       a           a 64-bit integer to be converted.
 699  *       f           the rounding mode (ROUND_NEAR, etc.)
 700  * OUTPUT PARAMETERS:
 701  *       b           the t-floating number "a" is converted to.
 702  * RETURN VALUE:
 703  *       error flags - i.e., zero if no errors occurred,
 704  *       FPCR_INV if invalid operation occurred, etc.
 705  */
 706 unsigned long
 707 ieee_CVTQT (int f, unsigned long a, unsigned long *b)
     /*  */
 708 {
 709         EXTENDED op_b;
 710 
 711         op_b.s    = 0;
 712         op_b.f[0] = a;
 713         op_b.f[1] = 0;
 714         if (sign(a) < 0) {
 715                 op_b.s = 1;
 716                 op_b.f[0] = -a;
 717         }
 718         op_b.e = 55;
 719         normalize(&op_b);
 720         return round_t_ieee(f, &op_b, b);
 721 }
 722 
 723 
 724 /*
 725  * converts t-floating "a" to 64-bit integer (q-format) "b".
 726  *
 727  * INPUT PARAMETERS:
 728  *       a           a t-floating number to be converted.
 729  *       f           the rounding mode (ROUND_NEAR, etc. )
 730  * OUTPUT PARAMETERS:
 731  *       b           the 64-bit integer "a" is converted to.
 732  * RETURN VALUE:
 733  *       error flags - i.e., zero if no errors occurred,
 734  *       FPCR_INV if invalid operation occurred, etc.
 735  */
 736 unsigned long
 737 ieee_CVTTQ (int f, unsigned long a, unsigned long *b)
     /*  */
 738 {
 739         unsigned int midway;
 740         unsigned long ov, uv, res = 0;
 741         fpclass_t a_type;
 742         EXTENDED temp;
 743 
 744         *b = 0;
 745         a_type = extend_ieee(a, &temp, DOUBLE);
 746         if (a_type == NaN || a_type == INFTY)
 747                 return FPCR_INV;
 748         if (a_type == QNaN)
 749                 return 0;
 750 
 751         if (temp.e > 0) {
 752                 ov = 0;
 753                 while (temp.e > 0) {
 754                         --temp.e;
 755                         ov |= temp.f[1] >> 63;
 756                         sll128(temp.f);
 757                 }
 758                 if (ov || (temp.f[1] & 0xffc0000000000000))
 759                         res |= FPCR_IOV | FPCR_INE;
 760         }
 761         if (temp.e < 0) {
 762                 while (temp.e < 0) {
 763                         ++temp.e;
 764                         uv = temp.f[0] & 1;             /* save sticky bit */
 765                         srl128(temp.f);
 766                         temp.f[0] |= uv;
 767                 }
 768         }
 769         *b = ((temp.f[1] << 9) | (temp.f[0] >> 55)) & 0x7fffffffffffffff;
 770         /*
 771          * Notice: the fraction is only 52 bits long.  Thus, rounding
 772          * cannot possibly result in an integer overflow.
 773          */
 774         switch (RM(f)) {
 775               case ROUND_NEAR:
 776                 if (temp.f[0] & 0x0040000000000000) {
 777                         midway = (temp.f[0] & 0x003fffffffffffff) == 0;
 778                         if ((midway && (temp.f[0] & 0x0080000000000000)) ||
 779                             !midway)
 780                                 ++b;
 781                 }
 782                 break;
 783 
 784               case ROUND_PINF:
 785                 if ((temp.f[0] & 0x003fffffffffffff) != 0)
 786                         ++b;
 787                 break;
 788 
 789               case ROUND_NINF:
 790                 if ((temp.f[0] & 0x003fffffffffffff) != 0)
 791                         --b;
 792                 break;
 793 
 794               case ROUND_CHOP:
 795                 /* no action needed */
 796                 break;
 797         }
 798         if ((temp.f[0] & 0x003fffffffffffff) != 0)
 799                 res |= FPCR_INE;
 800 
 801         if (temp.s) {
 802                 *b = -*b;
 803         }
 804         return res;
 805 }
 806 
 807 
 808 unsigned long
 809 ieee_CMPTEQ (unsigned long a, unsigned long b, unsigned long *c)
     /*  */
 810 {
 811         EXTENDED op_a, op_b;
 812         fpclass_t a_type, b_type;
 813 
 814         *c = 0;
 815         a_type = extend_ieee(a, &op_a, DOUBLE);
 816         b_type = extend_ieee(b, &op_b, DOUBLE);
 817         if (a_type == NaN || b_type == NaN)
 818                 return FPCR_INV;
 819         if (a_type == QNaN || b_type == QNaN)
 820                 return 0;
 821 
 822         if ((op_a.e == op_b.e && op_a.s == op_b.s &&
 823              op_a.f[0] == op_b.f[0] && op_a.f[1] == op_b.f[1]) ||
 824             (a_type == ZERO && b_type == ZERO))
 825                 *c = 0x4000000000000000;
 826         return 0;
 827 }
 828 
 829 
 830 unsigned long
 831 ieee_CMPTLT (unsigned long a, unsigned long b, unsigned long *c)
     /*  */
 832 {
 833         fpclass_t a_type, b_type;
 834         EXTENDED op_a, op_b;
 835 
 836         *c = 0;
 837         a_type = extend_ieee(a, &op_a, DOUBLE);
 838         b_type = extend_ieee(b, &op_b, DOUBLE);
 839         if (a_type == NaN || b_type == NaN)
 840                 return FPCR_INV;
 841         if (a_type == QNaN || b_type == QNaN)
 842                 return 0;
 843 
 844         if ((op_a.s == 1 && op_b.s == 0 &&
 845              (a_type != ZERO || b_type != ZERO)) ||
 846             (op_a.s == 1 && op_b.s == 1 &&
 847              (op_a.e > op_b.e || (op_a.e == op_b.e &&
 848                                   cmp128(op_a.f, op_b.f) > 0))) ||
 849             (op_a.s == 0 && op_b.s == 0 &&
 850              (op_a.e < op_b.e || (op_a.e == op_b.e &&
 851                                   cmp128(op_a.f,op_b.f) < 0))))
 852                 *c = 0x4000000000000000;
 853         return 0;
 854 }
 855 
 856 
 857 unsigned long
 858 ieee_CMPTLE (unsigned long a, unsigned long b, unsigned long *c)
     /*  */
 859 {
 860         fpclass_t a_type, b_type;
 861         EXTENDED op_a, op_b;
 862 
 863         *c = 0;
 864         a_type = extend_ieee(a, &op_a, DOUBLE);
 865         b_type = extend_ieee(b, &op_b, DOUBLE);
 866         if (a_type == NaN || b_type == NaN)
 867                 return FPCR_INV;
 868         if (a_type == QNaN || b_type == QNaN)
 869                 return 0;
 870 
 871         if ((a_type == ZERO && b_type == ZERO) ||
 872             (op_a.s == 1 && op_b.s == 0) ||
 873             (op_a.s == 1 && op_b.s == 1 &&
 874              (op_a.e > op_b.e || (op_a.e == op_b.e &&
 875                                   cmp128(op_a.f,op_b.f) >= 0))) ||
 876             (op_a.s == 0 && op_b.s == 0 &&
 877              (op_a.e < op_b.e || (op_a.e == op_b.e &&
 878                                   cmp128(op_a.f,op_b.f) <= 0))))
 879                 *c = 0x4000000000000000;
 880         return 0;
 881 }
 882 
 883 
 884 unsigned long
 885 ieee_CMPTUN (unsigned long a, unsigned long b, unsigned long *c)
     /*  */
 886 {
 887         fpclass_t a_type, b_type;
 888         EXTENDED op_a, op_b;
 889 
 890         *c = 0x4000000000000000;
 891         a_type = extend_ieee(a, &op_a, DOUBLE);
 892         b_type = extend_ieee(b, &op_b, DOUBLE);
 893         if (a_type == NaN || b_type == NaN)
 894                 return FPCR_INV;
 895         if (a_type == QNaN || b_type == QNaN)
 896                 return 0;
 897         *c = 0;
 898         return 0;
 899 }
 900 
 901 
 902 /*
 903  * Add a + b = c, where a, b, and c are ieee s-floating numbers.  "f"
 904  * contains the rounding mode etc.
 905  */
 906 unsigned long
 907 ieee_ADDS (int f, unsigned long a, unsigned long b, unsigned long *c)
     /*  */
 908 {
 909         fpclass_t a_type, b_type;
 910         EXTENDED op_a, op_b, op_c;
 911 
 912         a_type = extend_ieee(a, &op_a, SINGLE);
 913         b_type = extend_ieee(b, &op_b, SINGLE);
 914         if ((a_type >= NaN && a_type <= INFTY) ||
 915             (b_type >= NaN && b_type <= INFTY))
 916         {
 917                 /* propagate NaNs according to arch. ref. handbook: */
 918                 if (b_type == QNaN)
 919                         *c = b;
 920                 else if (b_type == NaN)
 921                         *c = b | (1UL << 51);
 922                 else if (a_type == QNaN)
 923                         *c = a;
 924                 else if (a_type == NaN)
 925                         *c = a | (1UL << 51);
 926 
 927                 if (a_type == NaN || b_type == NaN)
 928                         return FPCR_INV;
 929                 if (a_type == QNaN || b_type == QNaN)
 930                         return 0;
 931 
 932                 if (a_type == INFTY && b_type == INFTY && sign(a) != sign(b)) {
 933                         *c = IEEE_QNaN;
 934                         return FPCR_INV;
 935                 }
 936                 if (a_type == INFTY)
 937                         *c = a;
 938                 else
 939                         *c = b;
 940                 return 0;
 941         }
 942 
 943         add_kernel_ieee(&op_a, &op_b, &op_c);
 944         /* special case for -0 + -0 ==> -0 */
 945         if (a_type == ZERO && b_type == ZERO)
 946                 op_c.s = op_a.s && op_b.s;
 947         return round_s_ieee(f, &op_c, c);
 948 }
 949 
 950 
 951 /*
 952  * Add a + b = c, where a, b, and c are ieee t-floating numbers.  "f"
 953  * contains the rounding mode etc.
 954  */
 955 unsigned long
 956 ieee_ADDT (int f, unsigned long a, unsigned long b, unsigned long *c)
     /*  */
 957 {
 958         fpclass_t a_type, b_type;
 959         EXTENDED op_a, op_b, op_c;
 960 
 961         a_type = extend_ieee(a, &op_a, DOUBLE);
 962         b_type = extend_ieee(b, &op_b, DOUBLE);
 963         if ((a_type >= NaN && a_type <= INFTY) ||
 964             (b_type >= NaN && b_type <= INFTY))
 965         {
 966                 /* propagate NaNs according to arch. ref. handbook: */
 967                 if (b_type == QNaN)
 968                         *c = b;
 969                 else if (b_type == NaN)
 970                         *c = b | (1UL << 51);
 971                 else if (a_type == QNaN)
 972                         *c = a;
 973                 else if (a_type == NaN)
 974                         *c = a | (1UL << 51);
 975 
 976                 if (a_type == NaN || b_type == NaN)
 977                         return FPCR_INV;
 978                 if (a_type == QNaN || b_type == QNaN)
 979                         return 0;
 980 
 981                 if (a_type == INFTY && b_type == INFTY && sign(a) != sign(b)) {
 982                         *c = IEEE_QNaN;
 983                         return FPCR_INV;
 984                 }
 985                 if (a_type == INFTY)
 986                         *c = a;
 987                 else
 988                         *c = b;
 989                 return 0;
 990         }
 991         add_kernel_ieee(&op_a, &op_b, &op_c);
 992         /* special case for -0 + -0 ==> -0 */
 993         if (a_type == ZERO && b_type == ZERO)
 994                 op_c.s = op_a.s && op_b.s;
 995 
 996         return round_t_ieee(f, &op_c, c);
 997 }
 998 
 999 
1000 /*
1001  * Subtract a - b = c, where a, b, and c are ieee s-floating numbers.
1002  * "f" contains the rounding mode etc.
1003  */
1004 unsigned long
1005 ieee_SUBS (int f, unsigned long a, unsigned long b, unsigned long *c)
     /*  */
1006 {
1007         fpclass_t a_type, b_type;
1008         EXTENDED op_a, op_b, op_c;
1009 
1010         a_type = extend_ieee(a, &op_a, SINGLE);
1011         b_type = extend_ieee(b, &op_b, SINGLE);
1012         if ((a_type >= NaN && a_type <= INFTY) ||
1013             (b_type >= NaN && b_type <= INFTY))
1014         {
1015                 /* propagate NaNs according to arch. ref. handbook: */
1016                 if (b_type == QNaN)
1017                         *c = b;
1018                 else if (b_type == NaN)
1019                         *c = b | (1UL << 51);
1020                 else if (a_type == QNaN)
1021                         *c = a;
1022                 else if (a_type == NaN)
1023                         *c = a | (1UL << 51);
1024 
1025                 if (a_type == NaN || b_type == NaN)
1026                         return FPCR_INV;
1027                 if (a_type == QNaN || b_type == QNaN)
1028                         return 0;
1029 
1030                 if (a_type == INFTY && b_type == INFTY && sign(a) == sign(b)) {
1031                         *c = IEEE_QNaN;
1032                         return FPCR_INV;
1033                 }
1034                 if (a_type == INFTY)
1035                         *c = a;
1036                 else
1037                         *c = b ^ (1UL << 63);
1038                 return 0;
1039         }
1040         op_b.s = !op_b.s;
1041         add_kernel_ieee(&op_a, &op_b, &op_c);
1042         /* special case for -0 - +0 ==> -0 */
1043         if (a_type == ZERO && b_type == ZERO)
1044                 op_c.s = op_a.s && op_b.s;
1045 
1046         return round_s_ieee(f, &op_c, c);
1047 }
1048 
1049 
1050 /*
1051  * Subtract a - b = c, where a, b, and c are ieee t-floating numbers.
1052  * "f" contains the rounding mode etc.
1053  */
1054 unsigned long
1055 ieee_SUBT (int f, unsigned long a, unsigned long b, unsigned long *c)
     /*  */
1056 {
1057         fpclass_t a_type, b_type;
1058         EXTENDED op_a, op_b, op_c;
1059 
1060         a_type = extend_ieee(a, &op_a, DOUBLE);
1061         b_type = extend_ieee(b, &op_b, DOUBLE);
1062         if ((a_type >= NaN && a_type <= INFTY) ||
1063             (b_type >= NaN && b_type <= INFTY))
1064         {
1065                 /* propagate NaNs according to arch. ref. handbook: */
1066                 if (b_type == QNaN)
1067                         *c = b;
1068                 else if (b_type == NaN)
1069                         *c = b | (1UL << 51);
1070                 else if (a_type == QNaN)
1071                         *c = a;
1072                 else if (a_type == NaN)
1073                         *c = a | (1UL << 51);
1074 
1075                 if (a_type == NaN || b_type == NaN)
1076                         return FPCR_INV;
1077                 if (a_type == QNaN || b_type == QNaN)
1078                         return 0;
1079 
1080                 if (a_type == INFTY && b_type == INFTY && sign(a) == sign(b)) {
1081                         *c = IEEE_QNaN;
1082                         return FPCR_INV;
1083                 }
1084                 if (a_type == INFTY)
1085                         *c = a;
1086                 else
1087                         *c = b ^ (1UL << 63);
1088                 return 0;
1089         }
1090         op_b.s = !op_b.s;
1091         add_kernel_ieee(&op_a, &op_b, &op_c);
1092         /* special case for -0 - +0 ==> -0 */
1093         if (a_type == ZERO && b_type == ZERO)
1094                 op_c.s = op_a.s && op_b.s;
1095 
1096         return round_t_ieee(f, &op_c, c);
1097 }
1098 
1099 
1100 /*
1101  * Multiply a x b = c, where a, b, and c are ieee s-floating numbers.
1102  * "f" contains the rounding mode.
1103  */
1104 unsigned long
1105 ieee_MULS (int f, unsigned long a, unsigned long b, unsigned long *c)
     /*  */
1106 {
1107         fpclass_t a_type, b_type;
1108         EXTENDED op_a, op_b, op_c;
1109 
1110         a_type = extend_ieee(a, &op_a, SINGLE);
1111         b_type = extend_ieee(b, &op_b, SINGLE);
1112         if ((a_type >= NaN && a_type <= INFTY) ||
1113             (b_type >= NaN && b_type <= INFTY))
1114         {
1115                 /* propagate NaNs according to arch. ref. handbook: */
1116                 if (b_type == QNaN)
1117                         *c = b;
1118                 else if (b_type == NaN)
1119                         *c = b | (1UL << 51);
1120                 else if (a_type == QNaN)
1121                         *c = a;
1122                 else if (a_type == NaN)
1123                         *c = a | (1UL << 51);
1124 
1125                 if (a_type == NaN || b_type == NaN)
1126                         return FPCR_INV;
1127                 if (a_type == QNaN || b_type == QNaN)
1128                         return 0;
1129 
1130                 if ((a_type == INFTY && b_type == ZERO) ||
1131                     (b_type == INFTY && a_type == ZERO))
1132                 {
1133                         *c = IEEE_QNaN;         /* return canonical QNaN */
1134                         return FPCR_INV;
1135                 }
1136                 if (a_type == INFTY)
1137                         *c = a ^ ((b >> 63) << 63);
1138                 else if (b_type == INFTY)
1139                         *c = b ^ ((a >> 63) << 63);
1140                 else
1141                         /* either of a and b are +/-0 */
1142                         *c = ((unsigned long) op_a.s ^ op_b.s) << 63;
1143                 return 0;
1144         }
1145         op_c.s = op_a.s ^ op_b.s;
1146         op_c.e = op_a.e + op_b.e;
1147         mul64(op_a.f[0], op_b.f[0], op_c.f);
1148 
1149         normalize(&op_c);
1150         op_c.e -= 55;           /* drop the 55 original bits. */
1151 
1152         return round_s_ieee(f, &op_c, c);
1153 }
1154 
1155 
1156 /*
1157  * Multiply a x b = c, where a, b, and c are ieee t-floating numbers.
1158  * "f" contains the rounding mode.
1159  */
1160 unsigned long
1161 ieee_MULT (int f, unsigned long a, unsigned long b, unsigned long *c)
     /*  */
1162 {
1163         fpclass_t a_type, b_type;
1164         EXTENDED op_a, op_b, op_c;
1165 
1166         *c = IEEE_QNaN;
1167         a_type = extend_ieee(a, &op_a, DOUBLE);
1168         b_type = extend_ieee(b, &op_b, DOUBLE);
1169         if ((a_type >= NaN && a_type <= ZERO) ||
1170             (b_type >= NaN && b_type <= ZERO))
1171         {
1172                 /* propagate NaNs according to arch. ref. handbook: */
1173                 if (b_type == QNaN)
1174                         *c = b;
1175                 else if (b_type == NaN)
1176                         *c = b | (1UL << 51);
1177                 else if (a_type == QNaN)
1178                         *c = a;
1179                 else if (a_type == NaN)
1180                         *c = a | (1UL << 51);
1181 
1182                 if (a_type == NaN || b_type == NaN)
1183                         return FPCR_INV;
1184                 if (a_type == QNaN || b_type == QNaN)
1185                         return 0;
1186 
1187                 if ((a_type == INFTY && b_type == ZERO) ||
1188                     (b_type == INFTY && a_type == ZERO))
1189                 {
1190                         *c = IEEE_QNaN;         /* return canonical QNaN */
1191                         return FPCR_INV;
1192                 }
1193                 if (a_type == INFTY)
1194                         *c = a ^ ((b >> 63) << 63);
1195                 else if (b_type == INFTY)
1196                         *c = b ^ ((a >> 63) << 63);
1197                 else
1198                         /* either of a and b are +/-0 */
1199                         *c = ((unsigned long) op_a.s ^ op_b.s) << 63;
1200                 return 0;
1201         }
1202         op_c.s = op_a.s ^ op_b.s;
1203         op_c.e = op_a.e + op_b.e;
1204         mul64(op_a.f[0], op_b.f[0], op_c.f);
1205 
1206         normalize(&op_c);
1207         op_c.e -= 55;   /* drop the 55 original bits. */
1208 
1209         return round_t_ieee(f, &op_c, c);
1210 }
1211 
1212 
1213 /*
1214  * Divide a / b = c, where a, b, and c are ieee s-floating numbers.
1215  * "f" contains the rounding mode etc.
1216  */
1217 unsigned long
1218 ieee_DIVS (int f, unsigned long a, unsigned long b, unsigned long *c)
     /*  */
1219 {
1220         fpclass_t a_type, b_type;
1221         EXTENDED op_a, op_b, op_c;
1222 
1223         a_type = extend_ieee(a, &op_a, SINGLE);
1224         b_type = extend_ieee(b, &op_b, SINGLE);
1225         if ((a_type >= NaN && a_type <= ZERO) ||
1226             (b_type >= NaN && b_type <= ZERO))
1227         {
1228                 unsigned long res;
1229 
1230                 /* propagate NaNs according to arch. ref. handbook: */
1231                 if (b_type == QNaN)
1232                         *c = b;
1233                 else if (b_type == NaN)
1234                         *c = b | (1UL << 51);
1235                 else if (a_type == QNaN)
1236                         *c = a;
1237                 else if (a_type == NaN)
1238                         *c = a | (1UL << 51);
1239 
1240                 if (a_type == NaN || b_type == NaN)
1241                         return FPCR_INV;
1242                 if (a_type == QNaN || b_type == QNaN)
1243                         return 0;
1244 
1245                 res = 0;
1246                 *c = IEEE_PINF;
1247                 if (a_type == INFTY) {
1248                         if (b_type == INFTY) {
1249                                 *c = IEEE_QNaN;
1250                                 return FPCR_INV;
1251                         }
1252                 } else if (b_type == ZERO) {
1253                         if (a_type == ZERO) {
1254                                 *c = IEEE_QNaN;
1255                                 return FPCR_INV;
1256                         }
1257                         res = FPCR_DZE;
1258                 } else
1259                         /* a_type == ZERO || b_type == INFTY */
1260                         *c = 0;
1261                 *c |= (unsigned long) (op_a.s ^ op_b.s) << 63;
1262                 return res;
1263         }
1264         op_c.s = op_a.s ^ op_b.s;
1265         op_c.e = op_a.e - op_b.e;
1266 
1267         op_a.f[1] = op_a.f[0];
1268         op_a.f[0] = 0;
1269         div128(op_a.f, op_b.f, op_c.f);
1270         if (a_type != ZERO)
1271                 /* force a sticky bit because DIVs never hit exact .5: */
1272                 op_c.f[0] |= STICKY_S;
1273         normalize(&op_c);
1274         op_c.e -= 9;            /* remove excess exp from original shift */
1275         return round_s_ieee(f, &op_c, c);
1276 }
1277 
1278 
1279 /*
1280  * Divide a/b = c, where a, b, and c are ieee t-floating numbers.  "f"
1281  * contains the rounding mode etc.
1282  */
1283 unsigned long
1284 ieee_DIVT (int f, unsigned long a, unsigned long b, unsigned long *c)
     /*  */
1285 {
1286         fpclass_t a_type, b_type;
1287         EXTENDED op_a, op_b, op_c;
1288 
1289         *c = IEEE_QNaN;
1290         a_type = extend_ieee(a, &op_a, DOUBLE);
1291         b_type = extend_ieee(b, &op_b, DOUBLE);
1292         if ((a_type >= NaN && a_type <= ZERO) ||
1293             (b_type >= NaN && b_type <= ZERO))
1294         {
1295                 unsigned long res;
1296 
1297                 /* propagate NaNs according to arch. ref. handbook: */
1298                 if (b_type == QNaN)
1299                         *c = b;
1300                 else if (b_type == NaN)
1301                         *c = b | (1UL << 51);
1302                 else if (a_type == QNaN)
1303                         *c = a;
1304                 else if (a_type == NaN)
1305                         *c = a | (1UL << 51);
1306 
1307                 if (a_type == NaN || b_type == NaN)
1308                         return FPCR_INV;
1309                 if (a_type == QNaN || b_type == QNaN)
1310                         return 0;
1311 
1312                 res = 0;
1313                 *c = IEEE_PINF;
1314                 if (a_type == INFTY) {
1315                         if (b_type == INFTY) {
1316                                 *c = IEEE_QNaN;
1317                                 return FPCR_INV;
1318                         }
1319                 } else if (b_type == ZERO) {
1320                         if (a_type == ZERO) {
1321                                 *c = IEEE_QNaN;
1322                                 return FPCR_INV;
1323                         }
1324                         res = FPCR_DZE;
1325                 } else
1326                         /* a_type == ZERO || b_type == INFTY */
1327                         *c = 0;
1328                 *c |= (unsigned long) (op_a.s ^ op_b.s) << 63;
1329                 return res;
1330         }
1331         op_c.s = op_a.s ^ op_b.s;
1332         op_c.e = op_a.e - op_b.e;
1333 
1334         op_a.f[1] = op_a.f[0];
1335         op_a.f[0] = 0;
1336         div128(op_a.f, op_b.f, op_c.f);
1337         if (a_type != ZERO)
1338                 /* force a sticky bit because DIVs never hit exact .5 */
1339                 op_c.f[0] |= STICKY_T;
1340         normalize(&op_c);
1341         op_c.e -= 9;            /* remove excess exp from original shift */
1342         return round_t_ieee(f, &op_c, c);
1343 }
/* */
root/arch/alpha/math-emu/ieee-math.c

DEFINITIONS
