gcc/config/fp-bit.c

   1 /* This is a software floating point library which can be used
   2    for targets without hardware floating point.
   3    Copyright (C) 1994, 1995, 1996, 1997, 1998, 2000, 2001, 2002, 2003,
   4    2004, 2005, 2008, 2009 Free Software Foundation, Inc.
   5
   6 This file is part of GCC.
   7
   8 GCC is free software; you can redistribute it and/or modify it under
   9 the terms of the GNU General Public License as published by the Free
  10 Software Foundation; either version 3, or (at your option) any later
  11 version.
  12
  13 GCC is distributed in the hope that it will be useful, but WITHOUT ANY
  14 WARRANTY; without even the implied warranty of MERCHANTABILITY or
  15 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  16 for more details.
  17
  18 Under Section 7 of GPL version 3, you are granted additional
  19 permissions described in the GCC Runtime Library Exception, version
  20 3.1, as published by the Free Software Foundation.
  21
  22 You should have received a copy of the GNU General Public License and
  23 a copy of the GCC Runtime Library Exception along with this program;
  24 see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
  25 <http://www.gnu.org/licenses/>.  */
  26
  27 /* This implements IEEE 754 format arithmetic, but does not provide a
  28    mechanism for setting the rounding mode, or for generating or handling
  29    exceptions.
  30
  31    The original code by Steve Chamberlain, hacked by Mark Eichin and Jim
  32    Wilson, all of Cygnus Support.  */
  33
  34 /* The intended way to use this file is to make two copies, add `#define FLOAT'
  35    to one copy, then compile both copies and add them to libgcc.a.  */
  36
  37 #include "tconfig.h"
  38 #include "coretypes.h"
  39 #include "tm.h"
  40 #include "config/fp-bit.h"
  41
  42 /* The following macros can be defined to change the behavior of this file:
  43    FLOAT: Implement a `float', aka SFmode, fp library.  If this is not
  44      defined, then this file implements a `double', aka DFmode, fp library.
  45    FLOAT_ONLY: Used with FLOAT, to implement a `float' only library, i.e.
  46      don't include float->double conversion which requires the double library.
  47      This is useful only for machines which can't support doubles, e.g. some
  48      8-bit processors.
  49    CMPtype: Specify the type that floating point compares should return.
  50      This defaults to SItype, aka int.
  51    US_SOFTWARE_GOFAST: This makes all entry points use the same names as the
  52      US Software goFast library.
  53    _DEBUG_BITFLOAT: This makes debugging the code a little easier, by adding
  54      two integers to the FLO_union_type.
  55    NO_DENORMALS: Disable handling of denormals.
  56    NO_NANS: Disable nan and infinity handling
  57    SMALL_MACHINE: Useful when operations on QIs and HIs are faster
  58      than on an SI */
  59
  60 /* We don't currently support extended floats (long doubles) on machines
  61    without hardware to deal with them.
  62
  63    These stubs are just to keep the linker from complaining about unresolved
  64    references which can be pulled in from libio & libstdc++, even if the
  65    user isn't using long doubles.  However, they may generate an unresolved
  66    external to abort if abort is not used by the function, and the stubs
  67    are referenced from within libc, since libgcc goes before and after the
  68    system library.  */
  69
  70 #ifdef DECLARE_LIBRARY_RENAMES
  71   DECLARE_LIBRARY_RENAMES
  72 #endif
  73
  74 #ifdef EXTENDED_FLOAT_STUBS
  75 extern void abort (void);
  76 void __extendsfxf2 (void) { abort(); }
  77 void __extenddfxf2 (void) { abort(); }
  78 void __truncxfdf2 (void) { abort(); }
  79 void __truncxfsf2 (void) { abort(); }
  80 void __fixxfsi (void) { abort(); }
  81 void __floatsixf (void) { abort(); }
  82 void __addxf3 (void) { abort(); }
  83 void __subxf3 (void) { abort(); }
  84 void __mulxf3 (void) { abort(); }
  85 void __divxf3 (void) { abort(); }
  86 void __negxf2 (void) { abort(); }
  87 void __eqxf2 (void) { abort(); }
  88 void __nexf2 (void) { abort(); }
  89 void __gtxf2 (void) { abort(); }
  90 void __gexf2 (void) { abort(); }
  91 void __lexf2 (void) { abort(); }
  92 void __ltxf2 (void) { abort(); }
  93
  94 void __extendsftf2 (void) { abort(); }
  95 void __extenddftf2 (void) { abort(); }
  96 void __trunctfdf2 (void) { abort(); }
  97 void __trunctfsf2 (void) { abort(); }
  98 void __fixtfsi (void) { abort(); }
  99 void __floatsitf (void) { abort(); }
 100 void __addtf3 (void) { abort(); }
 101 void __subtf3 (void) { abort(); }
 102 void __multf3 (void) { abort(); }
 103 void __divtf3 (void) { abort(); }
 104 void __negtf2 (void) { abort(); }
 105 void __eqtf2 (void) { abort(); }
 106 void __netf2 (void) { abort(); }
 107 void __gttf2 (void) { abort(); }
 108 void __getf2 (void) { abort(); }
 109 void __letf2 (void) { abort(); }
 110 void __lttf2 (void) { abort(); }
 111 #else   /* !EXTENDED_FLOAT_STUBS, rest of file */
 112
 113 /* IEEE "special" number predicates */
 114
 115 #ifdef NO_NANS
 116
 117 #define nan() 0
 118 #define isnan(x) 0
 119 #define isinf(x) 0
 120 #else
 121
 122 #if   defined L_thenan_sf
 123 const fp_number_type __thenan_sf = { CLASS_SNAN, 0, 0, {(fractype) 0} };
 124 #elif defined L_thenan_df
 125 const fp_number_type __thenan_df = { CLASS_SNAN, 0, 0, {(fractype) 0} };
 126 #elif defined L_thenan_tf
 127 const fp_number_type __thenan_tf = { CLASS_SNAN, 0, 0, {(fractype) 0} };
 128 #elif defined TFLOAT
 129 extern const fp_number_type __thenan_tf;
 130 #elif defined FLOAT
 131 extern const fp_number_type __thenan_sf;
 132 #else
 133 extern const fp_number_type __thenan_df;
 134 #endif
 135
 136 INLINE
 137 static const fp_number_type *
 138 makenan (void)
 139 {
 140 #ifdef TFLOAT
 141   return & __thenan_tf;
 142 #elif defined FLOAT
 143   return & __thenan_sf;
 144 #else
 145   return & __thenan_df;
 146 #endif
 147 }
 148
 149 INLINE
 150 static int
 151 isnan (const fp_number_type *x)
 152 {
 153   return __builtin_expect (x->class == CLASS_SNAN || x->class == CLASS_QNAN,
 154                            0);
 155 }
 156
 157 INLINE
 158 static int
 159 isinf (const fp_number_type *  x)
 160 {
 161   return __builtin_expect (x->class == CLASS_INFINITY, 0);
 162 }
 163
 164 #endif /* NO_NANS */
 165
 166 INLINE
 167 static int
 168 iszero (const fp_number_type *  x)
 169 {
 170   return x->class == CLASS_ZERO;
 171 }
 172
 173 INLINE
 174 static void
 175 flip_sign ( fp_number_type *  x)
 176 {
 177   x->sign = !x->sign;
 178 }
 179
 180 /* Count leading zeroes in N.  */
 181 INLINE
 182 static int
 183 clzusi (USItype n)
 184 {
 185   extern int __clzsi2 (USItype);
 186   if (sizeof (USItype) == sizeof (unsigned int))
 187     return __builtin_clz (n);
 188   else if (sizeof (USItype) == sizeof (unsigned long))
 189     return __builtin_clzl (n);
 190   else if (sizeof (USItype) == sizeof (unsigned long long))
 191     return __builtin_clzll (n);
 192   else
 193     return __clzsi2 (n);
 194 }
 195
 196 extern FLO_type pack_d (const fp_number_type * );
 197
 198 #if defined(L_pack_df) || defined(L_pack_sf) || defined(L_pack_tf)
 199 FLO_type
 200 pack_d (const fp_number_type *src)
 201 {
 202   FLO_union_type dst;
 203   fractype fraction = src->fraction.ll; /* wasn't unsigned before? */
 204   int sign = src->sign;
 205   int exp = 0;
 206
 207   if (LARGEST_EXPONENT_IS_NORMAL (FRAC_NBITS) && (isnan (src) || isinf (src)))
 208     {
 209       /* We can't represent these values accurately.  By using the
 210          largest possible magnitude, we guarantee that the conversion
 211          of infinity is at least as big as any finite number.  */
 212       exp = EXPMAX;
 213       fraction = ((fractype) 1 << FRACBITS) - 1;
 214     }
 215   else if (isnan (src))
 216     {
 217       exp = EXPMAX;
 218       if (src->class == CLASS_QNAN || 1)
 219         {
 220 #ifdef QUIET_NAN_NEGATED
 221           fraction |= QUIET_NAN - 1;
 222 #else
 223           fraction |= QUIET_NAN;
 224 #endif
 225         }
 226     }
 227   else if (isinf (src))
 228     {
 229       exp = EXPMAX;
 230       fraction = 0;
 231     }
 232   else if (iszero (src))
 233     {
 234       exp = 0;
 235       fraction = 0;
 236     }
 237   else if (fraction == 0)
 238     {
 239       exp = 0;
 240     }
 241   else
 242     {
 243       if (__builtin_expect (src->normal_exp < NORMAL_EXPMIN, 0))
 244         {
 245 #ifdef NO_DENORMALS
 246           /* Go straight to a zero representation if denormals are not
 247              supported.  The denormal handling would be harmless but
 248              isn't unnecessary.  */
 249           exp = 0;
 250           fraction = 0;
 251 #else /* NO_DENORMALS */
 252           /* This number's exponent is too low to fit into the bits
 253              available in the number, so we'll store 0 in the exponent and
 254              shift the fraction to the right to make up for it.  */
 255
 256           int shift = NORMAL_EXPMIN - src->normal_exp;
 257
 258           exp = 0;
 259
 260           if (shift > FRAC_NBITS - NGARDS)
 261             {
 262               /* No point shifting, since it's more that 64 out.  */
 263               fraction = 0;
 264             }
 265           else
 266             {
 267               int lowbit = (fraction & (((fractype)1 << shift) - 1)) ? 1 : 0;
 268               fraction = (fraction >> shift) | lowbit;
 269             }
 270           if ((fraction & GARDMASK) == GARDMSB)
 271             {
 272               if ((fraction & (1 << NGARDS)))
 273                 fraction += GARDROUND + 1;
 274             }
 275           else
 276             {
 277               /* Add to the guards to round up.  */
 278               fraction += GARDROUND;
 279             }
 280           /* Perhaps the rounding means we now need to change the
 281              exponent, because the fraction is no longer denormal.  */
 282           if (fraction >= IMPLICIT_1)
 283             {
 284               exp += 1;
 285             }
 286           fraction >>= NGARDS;
 287 #endif /* NO_DENORMALS */
 288         }
 289       else if (!LARGEST_EXPONENT_IS_NORMAL (FRAC_NBITS)
 290                && __builtin_expect (src->normal_exp > EXPBIAS, 0))
 291         {
 292           exp = EXPMAX;
 293           fraction = 0;
 294         }
 295       else
 296         {
 297           exp = src->normal_exp + EXPBIAS;
 298           if (!ROUND_TOWARDS_ZERO)
 299             {
 300               /* IF the gard bits are the all zero, but the first, then we're
 301                  half way between two numbers, choose the one which makes the
 302                  lsb of the answer 0.  */
 303               if ((fraction & GARDMASK) == GARDMSB)
 304                 {
 305                   if (fraction & (1 << NGARDS))
 306                     fraction += GARDROUND + 1;
 307                 }
 308               else
 309                 {
 310                   /* Add a one to the guards to round up */
 311                   fraction += GARDROUND;
 312                 }
 313               if (fraction >= IMPLICIT_2)
 314                 {
 315                   fraction >>= 1;
 316                   exp += 1;
 317                 }
 318             }
 319           fraction >>= NGARDS;
 320
 321           if (LARGEST_EXPONENT_IS_NORMAL (FRAC_NBITS) && exp > EXPMAX)
 322             {
 323               /* Saturate on overflow.  */
 324               exp = EXPMAX;
 325               fraction = ((fractype) 1 << FRACBITS) - 1;
 326             }
 327         }
 328     }
 329
 330   /* We previously used bitfields to store the number, but this doesn't
 331      handle little/big endian systems conveniently, so use shifts and
 332      masks */
 333 #ifdef FLOAT_BIT_ORDER_MISMATCH
 334   dst.bits.fraction = fraction;
 335   dst.bits.exp = exp;
 336   dst.bits.sign = sign;
 337 #else
 338 # if defined TFLOAT && defined HALFFRACBITS
 339  {
 340    halffractype high, low, unity;
 341    int lowsign, lowexp;
 342
 343    unity = (halffractype) 1 << HALFFRACBITS;
 344
 345    /* Set HIGH to the high double's significand, masking out the implicit 1.
 346       Set LOW to the low double's full significand.  */
 347    high = (fraction >> (FRACBITS - HALFFRACBITS)) & (unity - 1);
 348    low = fraction & (unity * 2 - 1);
 349
 350    /* Get the initial sign and exponent of the low double.  */
 351    lowexp = exp - HALFFRACBITS - 1;
 352    lowsign = sign;
 353
 354    /* HIGH should be rounded like a normal double, making |LOW| <=
 355       0.5 ULP of HIGH.  Assume round-to-nearest.  */
 356    if (exp < EXPMAX)
 357      if (low > unity || (low == unity && (high & 1) == 1))
 358        {
 359          /* Round HIGH up and adjust LOW to match.  */
 360          high++;
 361          if (high == unity)
 362            {
 363              /* May make it infinite, but that's OK.  */
 364              high = 0;
 365              exp++;
 366            }
 367          low = unity * 2 - low;
 368          lowsign ^= 1;
 369        }
 370
 371    high |= (halffractype) exp << HALFFRACBITS;
 372    high |= (halffractype) sign << (HALFFRACBITS + EXPBITS);
 373
 374    if (exp == EXPMAX || exp == 0 || low == 0)
 375      low = 0;
 376    else
 377      {
 378        while (lowexp > 0 && low < unity)
 379          {
 380            low <<= 1;
 381            lowexp--;
 382          }
 383
 384        if (lowexp <= 0)
 385          {
 386            halffractype roundmsb, round;
 387            int shift;
 388
 389            shift = 1 - lowexp;
 390            roundmsb = (1 << (shift - 1));
 391            round = low & ((roundmsb << 1) - 1);
 392
 393            low >>= shift;
 394            lowexp = 0;
 395
 396            if (round > roundmsb || (round == roundmsb && (low & 1) == 1))
 397              {
 398                low++;
 399                if (low == unity)
 400                  /* LOW rounds up to the smallest normal number.  */
 401                  lowexp++;
 402              }
 403          }
 404
 405        low &= unity - 1;
 406        low |= (halffractype) lowexp << HALFFRACBITS;
 407        low |= (halffractype) lowsign << (HALFFRACBITS + EXPBITS);
 408      }
 409    dst.value_raw = ((fractype) high << HALFSHIFT) | low;
 410  }
 411 # else
 412   dst.value_raw = fraction & ((((fractype)1) << FRACBITS) - (fractype)1);
 413   dst.value_raw |= ((fractype) (exp & ((1 << EXPBITS) - 1))) << FRACBITS;
 414   dst.value_raw |= ((fractype) (sign & 1)) << (FRACBITS | EXPBITS);
 415 # endif
 416 #endif
 417
 418 #if defined(FLOAT_WORD_ORDER_MISMATCH) && !defined(FLOAT)
 419 #ifdef TFLOAT
 420   {
 421     qrtrfractype tmp1 = dst.words[0];
 422     qrtrfractype tmp2 = dst.words[1];
 423     dst.words[0] = dst.words[3];
 424     dst.words[1] = dst.words[2];
 425     dst.words[2] = tmp2;
 426     dst.words[3] = tmp1;
 427   }
 428 #else
 429   {
 430     halffractype tmp = dst.words[0];
 431     dst.words[0] = dst.words[1];
 432     dst.words[1] = tmp;
 433   }
 434 #endif
 435 #endif
 436
 437   return dst.value;
 438 }
 439 #endif
 440
 441 #if defined(L_unpack_df) || defined(L_unpack_sf) || defined(L_unpack_tf)
 442 void
 443 unpack_d (FLO_union_type * src, fp_number_type * dst)
 444 {
 445   /* We previously used bitfields to store the number, but this doesn't
 446      handle little/big endian systems conveniently, so use shifts and
 447      masks */
 448   fractype fraction;
 449   int exp;
 450   int sign;
 451
 452 #if defined(FLOAT_WORD_ORDER_MISMATCH) && !defined(FLOAT)
 453   FLO_union_type swapped;
 454
 455 #ifdef TFLOAT
 456   swapped.words[0] = src->words[3];
 457   swapped.words[1] = src->words[2];
 458   swapped.words[2] = src->words[1];
 459   swapped.words[3] = src->words[0];
 460 #else
 461   swapped.words[0] = src->words[1];
 462   swapped.words[1] = src->words[0];
 463 #endif
 464   src = &swapped;
 465 #endif
 466
 467 #ifdef FLOAT_BIT_ORDER_MISMATCH
 468   fraction = src->bits.fraction;
 469   exp = src->bits.exp;
 470   sign = src->bits.sign;
 471 #else
 472 # if defined TFLOAT && defined HALFFRACBITS
 473  {
 474    halffractype high, low;
 475
 476    high = src->value_raw >> HALFSHIFT;
 477    low = src->value_raw & (((fractype)1 << HALFSHIFT) - 1);
 478
 479    fraction = high & ((((fractype)1) << HALFFRACBITS) - 1);
 480    fraction <<= FRACBITS - HALFFRACBITS;
 481    exp = ((int)(high >> HALFFRACBITS)) & ((1 << EXPBITS) - 1);
 482    sign = ((int)(high >> (((HALFFRACBITS + EXPBITS))))) & 1;
 483
 484    if (exp != EXPMAX && exp != 0 && low != 0)
 485      {
 486        int lowexp = ((int)(low >> HALFFRACBITS)) & ((1 << EXPBITS) - 1);
 487        int lowsign = ((int)(low >> (((HALFFRACBITS + EXPBITS))))) & 1;
 488        int shift;
 489        fractype xlow;
 490
 491        xlow = low & ((((fractype)1) << HALFFRACBITS) - 1);
 492        if (lowexp)
 493          xlow |= (((halffractype)1) << HALFFRACBITS);
 494        else
 495          lowexp = 1;
 496        shift = (FRACBITS - HALFFRACBITS) - (exp - lowexp);
 497        if (shift > 0)
 498          xlow <<= shift;
 499        else if (shift < 0)
 500          xlow >>= -shift;
 501        if (sign == lowsign)
 502          fraction += xlow;
 503        else if (fraction >= xlow)
 504          fraction -= xlow;
 505        else
 506          {
 507            /* The high part is a power of two but the full number is lower.
 508               This code will leave the implicit 1 in FRACTION, but we'd
 509               have added that below anyway.  */
 510            fraction = (((fractype) 1 << FRACBITS) - xlow) << 1;
 511            exp--;
 512          }
 513      }
 514  }
 515 # else
 516   fraction = src->value_raw & ((((fractype)1) << FRACBITS) - 1);
 517   exp = ((int)(src->value_raw >> FRACBITS)) & ((1 << EXPBITS) - 1);
 518   sign = ((int)(src->value_raw >> (FRACBITS + EXPBITS))) & 1;
 519 # endif
 520 #endif
 521
 522   dst->sign = sign;
 523   if (exp == 0)
 524     {
 525       /* Hmm.  Looks like 0 */
 526       if (fraction == 0
 527 #ifdef NO_DENORMALS
 528           || 1
 529 #endif
 530           )
 531         {
 532           /* tastes like zero */
 533           dst->class = CLASS_ZERO;
 534         }
 535       else
 536         {
 537           /* Zero exponent with nonzero fraction - it's denormalized,
 538              so there isn't a leading implicit one - we'll shift it so
 539              it gets one.  */
 540           dst->normal_exp = exp - EXPBIAS + 1;
 541           fraction <<= NGARDS;
 542
 543           dst->class = CLASS_NUMBER;
 544 #if 1
 545           while (fraction < IMPLICIT_1)
 546             {
 547               fraction <<= 1;
 548               dst->normal_exp--;
 549             }
 550 #endif
 551           dst->fraction.ll = fraction;
 552         }
 553     }
 554   else if (!LARGEST_EXPONENT_IS_NORMAL (FRAC_NBITS)
 555            && __builtin_expect (exp == EXPMAX, 0))
 556     {
 557       /* Huge exponent*/
 558       if (fraction == 0)
 559         {
 560           /* Attached to a zero fraction - means infinity */
 561           dst->class = CLASS_INFINITY;
 562         }
 563       else
 564         {
 565           /* Nonzero fraction, means nan */
 566 #ifdef QUIET_NAN_NEGATED
 567           if ((fraction & QUIET_NAN) == 0)
 568 #else
 569           if (fraction & QUIET_NAN)
 570 #endif
 571             {
 572               dst->class = CLASS_QNAN;
 573             }
 574           else
 575             {
 576               dst->class = CLASS_SNAN;
 577             }
 578           /* Keep the fraction part as the nan number */
 579           dst->fraction.ll = fraction;
 580         }
 581     }
 582   else
 583     {
 584       /* Nothing strange about this number */
 585       dst->normal_exp = exp - EXPBIAS;
 586       dst->class = CLASS_NUMBER;
 587       dst->fraction.ll = (fraction << NGARDS) | IMPLICIT_1;
 588     }
 589 }
 590 #endif /* L_unpack_df || L_unpack_sf */
 591
 592 #if defined(L_addsub_sf) || defined(L_addsub_df) || defined(L_addsub_tf)
 593 static const fp_number_type *
 594 _fpadd_parts (fp_number_type * a,
 595               fp_number_type * b,
 596               fp_number_type * tmp)
 597 {
 598   intfrac tfraction;
 599
 600   /* Put commonly used fields in local variables.  */
 601   int a_normal_exp;
 602   int b_normal_exp;
 603   fractype a_fraction;
 604   fractype b_fraction;
 605
 606   if (isnan (a))
 607     {
 608       return a;
 609     }
 610   if (isnan (b))
 611     {
 612       return b;
 613     }
 614   if (isinf (a))
 615     {
 616       /* Adding infinities with opposite signs yields a NaN.  */
 617       if (isinf (b) && a->sign != b->sign)
 618         return makenan ();
 619       return a;
 620     }
 621   if (isinf (b))
 622     {
 623       return b;
 624     }
 625   if (iszero (b))
 626     {
 627       if (iszero (a))
 628         {
 629           *tmp = *a;
 630           tmp->sign = a->sign & b->sign;
 631           return tmp;
 632         }
 633       return a;
 634     }
 635   if (iszero (a))
 636     {
 637       return b;
 638     }
 639
 640   /* Got two numbers. shift the smaller and increment the exponent till
 641      they're the same */
 642   {
 643     int diff;
 644     int sdiff;
 645
 646     a_normal_exp = a->normal_exp;
 647     b_normal_exp = b->normal_exp;
 648     a_fraction = a->fraction.ll;
 649     b_fraction = b->fraction.ll;
 650
 651     diff = a_normal_exp - b_normal_exp;
 652     sdiff = diff;
 653
 654     if (diff < 0)
 655       diff = -diff;
 656     if (diff < FRAC_NBITS)
 657       {
 658         if (sdiff > 0)
 659           {
 660             b_normal_exp += diff;
 661             LSHIFT (b_fraction, diff);
 662           }
 663         else if (sdiff < 0)
 664           {
 665             a_normal_exp += diff;
 666             LSHIFT (a_fraction, diff);
 667           }
 668       }
 669     else
 670       {
 671         /* Somethings's up.. choose the biggest */
 672         if (a_normal_exp > b_normal_exp)
 673           {
 674             b_normal_exp = a_normal_exp;
 675             b_fraction = 0;
 676           }
 677         else
 678           {
 679             a_normal_exp = b_normal_exp;
 680             a_fraction = 0;
 681           }
 682       }
 683   }
 684
 685   if (a->sign != b->sign)
 686     {
 687       if (a->sign)
 688         {
 689           tfraction = -a_fraction + b_fraction;
 690         }
 691       else
 692         {
 693           tfraction = a_fraction - b_fraction;
 694         }
 695       if (tfraction >= 0)
 696         {
 697           tmp->sign = 0;
 698           tmp->normal_exp = a_normal_exp;
 699           tmp->fraction.ll = tfraction;
 700         }
 701       else
 702         {
 703           tmp->sign = 1;
 704           tmp->normal_exp = a_normal_exp;
 705           tmp->fraction.ll = -tfraction;
 706         }
 707       /* and renormalize it */
 708
 709       while (tmp->fraction.ll < IMPLICIT_1 && tmp->fraction.ll)
 710         {
 711           tmp->fraction.ll <<= 1;
 712           tmp->normal_exp--;
 713         }
 714     }
 715   else
 716     {
 717       tmp->sign = a->sign;
 718       tmp->normal_exp = a_normal_exp;
 719       tmp->fraction.ll = a_fraction + b_fraction;
 720     }
 721   tmp->class = CLASS_NUMBER;
 722   /* Now the fraction is added, we have to shift down to renormalize the
 723      number */
 724
 725   if (tmp->fraction.ll >= IMPLICIT_2)
 726     {
 727       LSHIFT (tmp->fraction.ll, 1);
 728       tmp->normal_exp++;
 729     }
 730   return tmp;
 731 }
 732
 733 FLO_type
 734 add (FLO_type arg_a, FLO_type arg_b)
 735 {
 736   fp_number_type a;
 737   fp_number_type b;
 738   fp_number_type tmp;
 739   const fp_number_type *res;
 740   FLO_union_type au, bu;
 741
 742   au.value = arg_a;
 743   bu.value = arg_b;
 744
 745   unpack_d (&au, &a);
 746   unpack_d (&bu, &b);
 747
 748   res = _fpadd_parts (&a, &b, &tmp);
 749
 750   return pack_d (res);
 751 }
 752
 753 FLO_type
 754 sub (FLO_type arg_a, FLO_type arg_b)
 755 {
 756   fp_number_type a;
 757   fp_number_type b;
 758   fp_number_type tmp;
 759   const fp_number_type *res;
 760   FLO_union_type au, bu;
 761
 762   au.value = arg_a;
 763   bu.value = arg_b;
 764
 765   unpack_d (&au, &a);
 766   unpack_d (&bu, &b);
 767
 768   b.sign ^= 1;
 769
 770   res = _fpadd_parts (&a, &b, &tmp);
 771
 772   return pack_d (res);
 773 }
 774 #endif /* L_addsub_sf || L_addsub_df */
 775
 776 #if defined(L_mul_sf) || defined(L_mul_df) || defined(L_mul_tf)
 777 static inline __attribute__ ((__always_inline__)) const fp_number_type *
 778 _fpmul_parts ( fp_number_type *  a,
 779                fp_number_type *  b,
 780                fp_number_type * tmp)
 781 {
 782   fractype low = 0;
 783   fractype high = 0;
 784
 785   if (isnan (a))
 786     {
 787       a->sign = a->sign != b->sign;
 788       return a;
 789     }
 790   if (isnan (b))
 791     {
 792       b->sign = a->sign != b->sign;
 793       return b;
 794     }
 795   if (isinf (a))
 796     {
 797       if (iszero (b))
 798         return makenan ();
 799       a->sign = a->sign != b->sign;
 800       return a;
 801     }
 802   if (isinf (b))
 803     {
 804       if (iszero (a))
 805         {
 806           return makenan ();
 807         }
 808       b->sign = a->sign != b->sign;
 809       return b;
 810     }
 811   if (iszero (a))
 812     {
 813       a->sign = a->sign != b->sign;
 814       return a;
 815     }
 816   if (iszero (b))
 817     {
 818       b->sign = a->sign != b->sign;
 819       return b;
 820     }
 821
 822   /* Calculate the mantissa by multiplying both numbers to get a
 823      twice-as-wide number.  */
 824   {
 825 #if defined(NO_DI_MODE) || defined(TFLOAT)
 826     {
 827       fractype x = a->fraction.ll;
 828       fractype ylow = b->fraction.ll;
 829       fractype yhigh = 0;
 830       int bit;
 831
 832       /* ??? This does multiplies one bit at a time.  Optimize.  */
 833       for (bit = 0; bit < FRAC_NBITS; bit++)
 834         {
 835           int carry;
 836
 837           if (x & 1)
 838             {
 839               carry = (low += ylow) < ylow;
 840               high += yhigh + carry;
 841             }
 842           yhigh <<= 1;
 843           if (ylow & FRACHIGH)
 844             {
 845               yhigh |= 1;
 846             }
 847           ylow <<= 1;
 848           x >>= 1;
 849         }
 850     }
 851 #elif defined(FLOAT)
 852     /* Multiplying two USIs to get a UDI, we're safe.  */
 853     {
 854       UDItype answer = (UDItype)a->fraction.ll * (UDItype)b->fraction.ll;
 855
 856       high = answer >> BITS_PER_SI;
 857       low = answer;
 858     }
 859 #else
 860     /* fractype is DImode, but we need the result to be twice as wide.
 861        Assuming a widening multiply from DImode to TImode is not
 862        available, build one by hand.  */
 863     {
 864       USItype nl = a->fraction.ll;
 865       USItype nh = a->fraction.ll >> BITS_PER_SI;
 866       USItype ml = b->fraction.ll;
 867       USItype mh = b->fraction.ll >> BITS_PER_SI;
 868       UDItype pp_ll = (UDItype) ml * nl;
 869       UDItype pp_hl = (UDItype) mh * nl;
 870       UDItype pp_lh = (UDItype) ml * nh;
 871       UDItype pp_hh = (UDItype) mh * nh;
 872       UDItype res2 = 0;
 873       UDItype res0 = 0;
 874       UDItype ps_hh__ = pp_hl + pp_lh;
 875       if (ps_hh__ < pp_hl)
 876         res2 += (UDItype)1 << BITS_PER_SI;
 877       pp_hl = (UDItype)(USItype)ps_hh__ << BITS_PER_SI;
 878       res0 = pp_ll + pp_hl;
 879       if (res0 < pp_ll)
 880         res2++;
 881       res2 += (ps_hh__ >> BITS_PER_SI) + pp_hh;
 882       high = res2;
 883       low = res0;
 884     }
 885 #endif
 886   }
 887
 888   tmp->normal_exp = a->normal_exp + b->normal_exp
 889     + FRAC_NBITS - (FRACBITS + NGARDS);
 890   tmp->sign = a->sign != b->sign;
 891   while (high >= IMPLICIT_2)
 892     {
 893       tmp->normal_exp++;
 894       if (high & 1)
 895         {
 896           low >>= 1;
 897           low |= FRACHIGH;
 898         }
 899       high >>= 1;
 900     }
 901   while (high < IMPLICIT_1)
 902     {
 903       tmp->normal_exp--;
 904
 905       high <<= 1;
 906       if (low & FRACHIGH)
 907         high |= 1;
 908       low <<= 1;
 909     }
 910
 911   if (!ROUND_TOWARDS_ZERO && (high & GARDMASK) == GARDMSB)
 912     {
 913       if (high & (1 << NGARDS))
 914         {
 915           /* Because we're half way, we would round to even by adding
 916              GARDROUND + 1, except that's also done in the packing
 917              function, and rounding twice will lose precision and cause
 918              the result to be too far off.  Example: 32-bit floats with
 919              bit patterns 0xfff * 0x3f800400 ~= 0xfff (less than 0.5ulp
 920              off), not 0x1000 (more than 0.5ulp off).  */
 921         }
 922       else if (low)
 923         {
 924           /* We're a further than half way by a small amount corresponding
 925              to the bits set in "low".  Knowing that, we round here and
 926              not in pack_d, because there we don't have "low" available
 927              anymore.  */
 928           high += GARDROUND + 1;
 929
 930           /* Avoid further rounding in pack_d.  */
 931           high &= ~(fractype) GARDMASK;
 932         }
 933     }
 934   tmp->fraction.ll = high;
 935   tmp->class = CLASS_NUMBER;
 936   return tmp;
 937 }
 938
 939 FLO_type
 940 multiply (FLO_type arg_a, FLO_type arg_b)
 941 {
 942   fp_number_type a;
 943   fp_number_type b;
 944   fp_number_type tmp;
 945   const fp_number_type *res;
 946   FLO_union_type au, bu;
 947
 948   au.value = arg_a;
 949   bu.value = arg_b;
 950
 951   unpack_d (&au, &a);
 952   unpack_d (&bu, &b);
 953
 954   res = _fpmul_parts (&a, &b, &tmp);
 955
 956   return pack_d (res);
 957 }
 958 #endif /* L_mul_sf || L_mul_df || L_mul_tf */
 959
 960 #if defined(L_div_sf) || defined(L_div_df) || defined(L_div_tf)
 961 static inline __attribute__ ((__always_inline__)) const fp_number_type *
 962 _fpdiv_parts (fp_number_type * a,
 963               fp_number_type * b)
 964 {
 965   fractype bit;
 966   fractype numerator;
 967   fractype denominator;
 968   fractype quotient;
 969
 970   if (isnan (a))
 971     {
 972       return a;
 973     }
 974   if (isnan (b))
 975     {
 976       return b;
 977     }
 978
 979   a->sign = a->sign ^ b->sign;
 980
 981   if (isinf (a) || iszero (a))
 982     {
 983       if (a->class == b->class)
 984         return makenan ();
 985       return a;
 986     }
 987
 988   if (isinf (b))
 989     {
 990       a->fraction.ll = 0;
 991       a->normal_exp = 0;
 992       return a;
 993     }
 994   if (iszero (b))
 995     {
 996       a->class = CLASS_INFINITY;
 997       return a;
 998     }
 999
1000   /* Calculate the mantissa by multiplying both 64bit numbers to get a
1001      128 bit number */
1002   {
1003     /* quotient =
1004        ( numerator / denominator) * 2^(numerator exponent -  denominator exponent)
1005      */
1006
1007     a->normal_exp = a->normal_exp - b->normal_exp;
1008     numerator = a->fraction.ll;
1009     denominator = b->fraction.ll;
1010
1011     if (numerator < denominator)
1012       {
1013         /* Fraction will be less than 1.0 */
1014         numerator *= 2;
1015         a->normal_exp--;
1016       }
1017     bit = IMPLICIT_1;
1018     quotient = 0;
1019     /* ??? Does divide one bit at a time.  Optimize.  */
1020     while (bit)
1021       {
1022         if (numerator >= denominator)
1023           {
1024             quotient |= bit;
1025             numerator -= denominator;
1026           }
1027         bit >>= 1;
1028         numerator *= 2;
1029       }
1030
1031     if (!ROUND_TOWARDS_ZERO && (quotient & GARDMASK) == GARDMSB)
1032       {
1033         if (quotient & (1 << NGARDS))
1034           {
1035             /* Because we're half way, we would round to even by adding
1036                GARDROUND + 1, except that's also done in the packing
1037                function, and rounding twice will lose precision and cause
1038                the result to be too far off.  */
1039           }
1040         else if (numerator)
1041           {
1042             /* We're a further than half way by the small amount
1043                corresponding to the bits set in "numerator".  Knowing
1044                that, we round here and not in pack_d, because there we
1045                don't have "numerator" available anymore.  */
1046             quotient += GARDROUND + 1;
1047
1048             /* Avoid further rounding in pack_d.  */
1049             quotient &= ~(fractype) GARDMASK;
1050           }
1051       }
1052
1053     a->fraction.ll = quotient;
1054     return (a);
1055   }
1056 }
1057
1058 FLO_type
1059 divide (FLO_type arg_a, FLO_type arg_b)
1060 {
1061   fp_number_type a;
1062   fp_number_type b;
1063   const fp_number_type *res;
1064   FLO_union_type au, bu;
1065
1066   au.value = arg_a;
1067   bu.value = arg_b;
1068
1069   unpack_d (&au, &a);
1070   unpack_d (&bu, &b);
1071
1072   res = _fpdiv_parts (&a, &b);
1073
1074   return pack_d (res);
1075 }
1076 #endif /* L_div_sf || L_div_df */
1077
1078 #if defined(L_fpcmp_parts_sf) || defined(L_fpcmp_parts_df) \
1079     || defined(L_fpcmp_parts_tf)
1080 /* according to the demo, fpcmp returns a comparison with 0... thus
1081    a<b -> -1
1082    a==b -> 0
1083    a>b -> +1
1084  */
1085
1086 int
1087 __fpcmp_parts (fp_number_type * a, fp_number_type * b)
1088 {
1089 #if 0
1090   /* either nan -> unordered. Must be checked outside of this routine.  */
1091   if (isnan (a) && isnan (b))
1092     {
1093       return 1;                 /* still unordered! */
1094     }
1095 #endif
1096
1097   if (isnan (a) || isnan (b))
1098     {
1099       return 1;                 /* how to indicate unordered compare? */
1100     }
1101   if (isinf (a) && isinf (b))
1102     {
1103       /* +inf > -inf, but +inf != +inf */
1104       /* b    \a| +inf(0)| -inf(1)
1105        ______\+--------+--------
1106        +inf(0)| a==b(0)| a<b(-1)
1107        -------+--------+--------
1108        -inf(1)| a>b(1) | a==b(0)
1109        -------+--------+--------
1110        So since unordered must be nonzero, just line up the columns...
1111        */
1112       return b->sign - a->sign;
1113     }
1114   /* but not both...  */
1115   if (isinf (a))
1116     {
1117       return a->sign ? -1 : 1;
1118     }
1119   if (isinf (b))
1120     {
1121       return b->sign ? 1 : -1;
1122     }
1123   if (iszero (a) && iszero (b))
1124     {
1125       return 0;
1126     }
1127   if (iszero (a))
1128     {
1129       return b->sign ? 1 : -1;
1130     }
1131   if (iszero (b))
1132     {
1133       return a->sign ? -1 : 1;
1134     }
1135   /* now both are "normal".  */
1136   if (a->sign != b->sign)
1137     {
1138       /* opposite signs */
1139       return a->sign ? -1 : 1;
1140     }
1141   /* same sign; exponents? */
1142   if (a->normal_exp > b->normal_exp)
1143     {
1144       return a->sign ? -1 : 1;
1145     }
1146   if (a->normal_exp < b->normal_exp)
1147     {
1148       return a->sign ? 1 : -1;
1149     }
1150   /* same exponents; check size.  */
1151   if (a->fraction.ll > b->fraction.ll)
1152     {
1153       return a->sign ? -1 : 1;
1154     }
1155   if (a->fraction.ll < b->fraction.ll)
1156     {
1157       return a->sign ? 1 : -1;
1158     }
1159   /* after all that, they're equal.  */
1160   return 0;
1161 }
1162 #endif
1163
1164 #if defined(L_compare_sf) || defined(L_compare_df) || defined(L_compoare_tf)
1165 CMPtype
1166 compare (FLO_type arg_a, FLO_type arg_b)
1167 {
1168   fp_number_type a;
1169   fp_number_type b;
1170   FLO_union_type au, bu;
1171
1172   au.value = arg_a;
1173   bu.value = arg_b;
1174
1175   unpack_d (&au, &a);
1176   unpack_d (&bu, &b);
1177
1178   return __fpcmp_parts (&a, &b);
1179 }
1180 #endif /* L_compare_sf || L_compare_df */
1181
1182 #ifndef US_SOFTWARE_GOFAST
1183
1184 /* These should be optimized for their specific tasks someday.  */
1185
1186 #if defined(L_eq_sf) || defined(L_eq_df) || defined(L_eq_tf)
1187 CMPtype
1188 _eq_f2 (FLO_type arg_a, FLO_type arg_b)
1189 {
1190   fp_number_type a;
1191   fp_number_type b;
1192   FLO_union_type au, bu;
1193
1194   au.value = arg_a;
1195   bu.value = arg_b;
1196
1197   unpack_d (&au, &a);
1198   unpack_d (&bu, &b);
1199
1200   if (isnan (&a) || isnan (&b))
1201     return 1;                   /* false, truth == 0 */
1202
1203   return __fpcmp_parts (&a, &b) ;
1204 }
1205 #endif /* L_eq_sf || L_eq_df */
1206
1207 #if defined(L_ne_sf) || defined(L_ne_df) || defined(L_ne_tf)
1208 CMPtype
1209 _ne_f2 (FLO_type arg_a, FLO_type arg_b)
1210 {
1211   fp_number_type a;
1212   fp_number_type b;
1213   FLO_union_type au, bu;
1214
1215   au.value = arg_a;
1216   bu.value = arg_b;
1217
1218   unpack_d (&au, &a);
1219   unpack_d (&bu, &b);
1220
1221   if (isnan (&a) || isnan (&b))
1222     return 1;                   /* true, truth != 0 */
1223
1224   return  __fpcmp_parts (&a, &b) ;
1225 }
1226 #endif /* L_ne_sf || L_ne_df */
1227
1228 #if defined(L_gt_sf) || defined(L_gt_df) || defined(L_gt_tf)
1229 CMPtype
1230 _gt_f2 (FLO_type arg_a, FLO_type arg_b)
1231 {
1232   fp_number_type a;
1233   fp_number_type b;
1234   FLO_union_type au, bu;
1235
1236   au.value = arg_a;
1237   bu.value = arg_b;
1238
1239   unpack_d (&au, &a);
1240   unpack_d (&bu, &b);
1241
1242   if (isnan (&a) || isnan (&b))
1243     return -1;                  /* false, truth > 0 */
1244
1245   return __fpcmp_parts (&a, &b);
1246 }
1247 #endif /* L_gt_sf || L_gt_df */
1248
1249 #if defined(L_ge_sf) || defined(L_ge_df) || defined(L_ge_tf)
1250 CMPtype
1251 _ge_f2 (FLO_type arg_a, FLO_type arg_b)
1252 {
1253   fp_number_type a;
1254   fp_number_type b;
1255   FLO_union_type au, bu;
1256
1257   au.value = arg_a;
1258   bu.value = arg_b;
1259
1260   unpack_d (&au, &a);
1261   unpack_d (&bu, &b);
1262
1263   if (isnan (&a) || isnan (&b))
1264     return -1;                  /* false, truth >= 0 */
1265   return __fpcmp_parts (&a, &b) ;
1266 }
1267 #endif /* L_ge_sf || L_ge_df */
1268
1269 #if defined(L_lt_sf) || defined(L_lt_df) || defined(L_lt_tf)
1270 CMPtype
1271 _lt_f2 (FLO_type arg_a, FLO_type arg_b)
1272 {
1273   fp_number_type a;
1274   fp_number_type b;
1275   FLO_union_type au, bu;
1276
1277   au.value = arg_a;
1278   bu.value = arg_b;
1279
1280   unpack_d (&au, &a);
1281   unpack_d (&bu, &b);
1282
1283   if (isnan (&a) || isnan (&b))
1284     return 1;                   /* false, truth < 0 */
1285
1286   return __fpcmp_parts (&a, &b);
1287 }
1288 #endif /* L_lt_sf || L_lt_df */
1289
1290 #if defined(L_le_sf) || defined(L_le_df) || defined(L_le_tf)
1291 CMPtype
1292 _le_f2 (FLO_type arg_a, FLO_type arg_b)
1293 {
1294   fp_number_type a;
1295   fp_number_type b;
1296   FLO_union_type au, bu;
1297
1298   au.value = arg_a;
1299   bu.value = arg_b;
1300
1301   unpack_d (&au, &a);
1302   unpack_d (&bu, &b);
1303
1304   if (isnan (&a) || isnan (&b))
1305     return 1;                   /* false, truth <= 0 */
1306
1307   return __fpcmp_parts (&a, &b) ;
1308 }
1309 #endif /* L_le_sf || L_le_df */
1310
1311 #endif /* ! US_SOFTWARE_GOFAST */
1312
1313 #if defined(L_unord_sf) || defined(L_unord_df) || defined(L_unord_tf)
1314 CMPtype
1315 _unord_f2 (FLO_type arg_a, FLO_type arg_b)
1316 {
1317   fp_number_type a;
1318   fp_number_type b;
1319   FLO_union_type au, bu;
1320
1321   au.value = arg_a;
1322   bu.value = arg_b;
1323
1324   unpack_d (&au, &a);
1325   unpack_d (&bu, &b);
1326
1327   return (isnan (&a) || isnan (&b));
1328 }
1329 #endif /* L_unord_sf || L_unord_df */
1330
1331 #if defined(L_si_to_sf) || defined(L_si_to_df) || defined(L_si_to_tf)
1332 FLO_type
1333 si_to_float (SItype arg_a)
1334 {
1335   fp_number_type in;
1336
1337   in.class = CLASS_NUMBER;
1338   in.sign = arg_a < 0;
1339   if (!arg_a)
1340     {
1341       in.class = CLASS_ZERO;
1342     }
1343   else
1344     {
1345       USItype uarg;
1346       int shift;
1347       in.normal_exp = FRACBITS + NGARDS;
1348       if (in.sign)
1349         {
1350           /* Special case for minint, since there is no +ve integer
1351              representation for it */
1352           if (arg_a == (- MAX_SI_INT - 1))
1353             {
1354               return (FLO_type)(- MAX_SI_INT - 1);
1355             }
1356           uarg = (-arg_a);
1357         }
1358       else
1359         uarg = arg_a;
1360
1361       in.fraction.ll = uarg;
1362       shift = clzusi (uarg) - (BITS_PER_SI - 1 - FRACBITS - NGARDS);
1363       if (shift > 0)
1364         {
1365           in.fraction.ll <<= shift;
1366           in.normal_exp -= shift;
1367         }
1368     }
1369   return pack_d (&in);
1370 }
1371 #endif /* L_si_to_sf || L_si_to_df */
1372
1373 #if defined(L_usi_to_sf) || defined(L_usi_to_df) || defined(L_usi_to_tf)
1374 FLO_type
1375 usi_to_float (USItype arg_a)
1376 {
1377   fp_number_type in;
1378
1379   in.sign = 0;
1380   if (!arg_a)
1381     {
1382       in.class = CLASS_ZERO;
1383     }
1384   else
1385     {
1386       int shift;
1387       in.class = CLASS_NUMBER;
1388       in.normal_exp = FRACBITS + NGARDS;
1389       in.fraction.ll = arg_a;
1390
1391       shift = clzusi (arg_a) - (BITS_PER_SI - 1 - FRACBITS - NGARDS);
1392       if (shift < 0)
1393         {
1394           fractype guard = in.fraction.ll & (((fractype)1 << -shift) - 1);
1395           in.fraction.ll >>= -shift;
1396           in.fraction.ll |= (guard != 0);
1397           in.normal_exp -= shift;
1398         }
1399       else if (shift > 0)
1400         {
1401           in.fraction.ll <<= shift;
1402           in.normal_exp -= shift;
1403         }
1404     }
1405   return pack_d (&in);
1406 }
1407 #endif
1408
1409 #if defined(L_sf_to_si) || defined(L_df_to_si) || defined(L_tf_to_si)
1410 SItype
1411 float_to_si (FLO_type arg_a)
1412 {
1413   fp_number_type a;
1414   SItype tmp;
1415   FLO_union_type au;
1416
1417   au.value = arg_a;
1418   unpack_d (&au, &a);
1419
1420   if (iszero (&a))
1421     return 0;
1422   if (isnan (&a))
1423     return 0;
1424   /* get reasonable MAX_SI_INT...  */
1425   if (isinf (&a))
1426     return a.sign ? (-MAX_SI_INT)-1 : MAX_SI_INT;
1427   /* it is a number, but a small one */
1428   if (a.normal_exp < 0)
1429     return 0;
1430   if (a.normal_exp > BITS_PER_SI - 2)
1431     return a.sign ? (-MAX_SI_INT)-1 : MAX_SI_INT;
1432   tmp = a.fraction.ll >> ((FRACBITS + NGARDS) - a.normal_exp);
1433   return a.sign ? (-tmp) : (tmp);
1434 }
1435 #endif /* L_sf_to_si || L_df_to_si */
1436
1437 #if defined(L_sf_to_usi) || defined(L_df_to_usi) || defined(L_tf_to_usi)
1438 #if defined US_SOFTWARE_GOFAST || defined(L_tf_to_usi)
1439 /* While libgcc2.c defines its own __fixunssfsi and __fixunsdfsi routines,
1440    we also define them for GOFAST because the ones in libgcc2.c have the
1441    wrong names and I'd rather define these here and keep GOFAST CYG-LOC's
1442    out of libgcc2.c.  We can't define these here if not GOFAST because then
1443    there'd be duplicate copies.  */
1444
1445 USItype
1446 float_to_usi (FLO_type arg_a)
1447 {
1448   fp_number_type a;
1449   FLO_union_type au;
1450
1451   au.value = arg_a;
1452   unpack_d (&au, &a);
1453
1454   if (iszero (&a))
1455     return 0;
1456   if (isnan (&a))
1457     return 0;
1458   /* it is a negative number */
1459   if (a.sign)
1460     return 0;
1461   /* get reasonable MAX_USI_INT...  */
1462   if (isinf (&a))
1463     return MAX_USI_INT;
1464   /* it is a number, but a small one */
1465   if (a.normal_exp < 0)
1466     return 0;
1467   if (a.normal_exp > BITS_PER_SI - 1)
1468     return MAX_USI_INT;
1469   else if (a.normal_exp > (FRACBITS + NGARDS))
1470     return a.fraction.ll << (a.normal_exp - (FRACBITS + NGARDS));
1471   else
1472     return a.fraction.ll >> ((FRACBITS + NGARDS) - a.normal_exp);
1473 }
1474 #endif /* US_SOFTWARE_GOFAST */
1475 #endif /* L_sf_to_usi || L_df_to_usi */
1476
1477 #if defined(L_negate_sf) || defined(L_negate_df) || defined(L_negate_tf)
1478 FLO_type
1479 negate (FLO_type arg_a)
1480 {
1481   fp_number_type a;
1482   FLO_union_type au;
1483
1484   au.value = arg_a;
1485   unpack_d (&au, &a);
1486
1487   flip_sign (&a);
1488   return pack_d (&a);
1489 }
1490 #endif /* L_negate_sf || L_negate_df */
1491
1492 #ifdef FLOAT
1493
1494 #if defined(L_make_sf)
1495 SFtype
1496 __make_fp(fp_class_type class,
1497              unsigned int sign,
1498              int exp,
1499              USItype frac)
1500 {
1501   fp_number_type in;
1502
1503   in.class = class;
1504   in.sign = sign;
1505   in.normal_exp = exp;
1506   in.fraction.ll = frac;
1507   return pack_d (&in);
1508 }
1509 #endif /* L_make_sf */
1510
1511 #ifndef FLOAT_ONLY
1512
1513 /* This enables one to build an fp library that supports float but not double.
1514    Otherwise, we would get an undefined reference to __make_dp.
1515    This is needed for some 8-bit ports that can't handle well values that
1516    are 8-bytes in size, so we just don't support double for them at all.  */
1517
1518 #if defined(L_sf_to_df)
1519 DFtype
1520 sf_to_df (SFtype arg_a)
1521 {
1522   fp_number_type in;
1523   FLO_union_type au;
1524
1525   au.value = arg_a;
1526   unpack_d (&au, &in);
1527
1528   return __make_dp (in.class, in.sign, in.normal_exp,
1529                     ((UDItype) in.fraction.ll) << F_D_BITOFF);
1530 }
1531 #endif /* L_sf_to_df */
1532
1533 #if defined(L_sf_to_tf) && defined(TMODES)
1534 TFtype
1535 sf_to_tf (SFtype arg_a)
1536 {
1537   fp_number_type in;
1538   FLO_union_type au;
1539
1540   au.value = arg_a;
1541   unpack_d (&au, &in);
1542
1543   return __make_tp (in.class, in.sign, in.normal_exp,
1544                     ((UTItype) in.fraction.ll) << F_T_BITOFF);
1545 }
1546 #endif /* L_sf_to_df */
1547
1548 #endif /* ! FLOAT_ONLY */
1549 #endif /* FLOAT */
1550
1551 #ifndef FLOAT
1552
1553 extern SFtype __make_fp (fp_class_type, unsigned int, int, USItype);
1554
1555 #if defined(L_make_df)
1556 DFtype
1557 __make_dp (fp_class_type class, unsigned int sign, int exp, UDItype frac)
1558 {
1559   fp_number_type in;
1560
1561   in.class = class;
1562   in.sign = sign;
1563   in.normal_exp = exp;
1564   in.fraction.ll = frac;
1565   return pack_d (&in);
1566 }
1567 #endif /* L_make_df */
1568
1569 #if defined(L_df_to_sf)
1570 SFtype
1571 df_to_sf (DFtype arg_a)
1572 {
1573   fp_number_type in;
1574   USItype sffrac;
1575   FLO_union_type au;
1576
1577   au.value = arg_a;
1578   unpack_d (&au, &in);
1579
1580   sffrac = in.fraction.ll >> F_D_BITOFF;
1581
1582   /* We set the lowest guard bit in SFFRAC if we discarded any non
1583      zero bits.  */
1584   if ((in.fraction.ll & (((USItype) 1 << F_D_BITOFF) - 1)) != 0)
1585     sffrac |= 1;
1586
1587   return __make_fp (in.class, in.sign, in.normal_exp, sffrac);
1588 }
1589 #endif /* L_df_to_sf */
1590
1591 #if defined(L_df_to_tf) && defined(TMODES) \
1592     && !defined(FLOAT) && !defined(TFLOAT)
1593 TFtype
1594 df_to_tf (DFtype arg_a)
1595 {
1596   fp_number_type in;
1597   FLO_union_type au;
1598
1599   au.value = arg_a;
1600   unpack_d (&au, &in);
1601
1602   return __make_tp (in.class, in.sign, in.normal_exp,
1603                     ((UTItype) in.fraction.ll) << D_T_BITOFF);
1604 }
1605 #endif /* L_sf_to_df */
1606
1607 #ifdef TFLOAT
1608 #if defined(L_make_tf)
1609 TFtype
1610 __make_tp(fp_class_type class,
1611              unsigned int sign,
1612              int exp,
1613              UTItype frac)
1614 {
1615   fp_number_type in;
1616
1617   in.class = class;
1618   in.sign = sign;
1619   in.normal_exp = exp;
1620   in.fraction.ll = frac;
1621   return pack_d (&in);
1622 }
1623 #endif /* L_make_tf */
1624
1625 #if defined(L_tf_to_df)
1626 DFtype
1627 tf_to_df (TFtype arg_a)
1628 {
1629   fp_number_type in;
1630   UDItype sffrac;
1631   FLO_union_type au;
1632
1633   au.value = arg_a;
1634   unpack_d (&au, &in);
1635
1636   sffrac = in.fraction.ll >> D_T_BITOFF;
1637
1638   /* We set the lowest guard bit in SFFRAC if we discarded any non
1639      zero bits.  */
1640   if ((in.fraction.ll & (((UTItype) 1 << D_T_BITOFF) - 1)) != 0)
1641     sffrac |= 1;
1642
1643   return __make_dp (in.class, in.sign, in.normal_exp, sffrac);
1644 }
1645 #endif /* L_tf_to_df */
1646
1647 #if defined(L_tf_to_sf)
1648 SFtype
1649 tf_to_sf (TFtype arg_a)
1650 {
1651   fp_number_type in;
1652   USItype sffrac;
1653   FLO_union_type au;
1654
1655   au.value = arg_a;
1656   unpack_d (&au, &in);
1657
1658   sffrac = in.fraction.ll >> F_T_BITOFF;
1659
1660   /* We set the lowest guard bit in SFFRAC if we discarded any non
1661      zero bits.  */
1662   if ((in.fraction.ll & (((UTItype) 1 << F_T_BITOFF) - 1)) != 0)
1663     sffrac |= 1;
1664
1665   return __make_fp (in.class, in.sign, in.normal_exp, sffrac);
1666 }
1667 #endif /* L_tf_to_sf */
1668 #endif /* TFLOAT */
1669
1670 #endif /* ! FLOAT */
1671 #endif /* !EXTENDED_FLOAT_STUBS */