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