gcc/config/fp-bit.c

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