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