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