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