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