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