gcc/config/fp-bit.c

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