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
   5 /* Copyright (C) 1994 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, 675 Mass Ave, Cambridge, MA 02139, 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 /* The following macros can be defined to change the behaviour of this file:
  47    FLOAT: Implement a `float', aka SFmode, fp library.  If this is not
  48      defined, then this file implements a `double', aka DFmode, fp library.
  49    FLOAT_ONLY: Used with FLOAT, to implement a `float' only library, i.e.
  50      don't include float->double conversion which requires the double library.
  51      This is useful only for machines which can't support doubles, e.g. some
  52      8-bit processors.
  53    CMPtype: Specify the type that floating point compares should return.
  54      This defaults to SItype, aka int.
  55    US_SOFTWARE_GOFAST: This makes all entry points use the same names as the
  56      US Software goFast library.  If this is not defined, the entry points use
  57      the same names as libgcc1.c.
  58    _DEBUG_BITFLOAT: This makes debugging the code a little easier, by adding
  59      two integers to the FLO_union_type.
  60    NO_NANS: Disable nan and infinity handling
  61    SMALL_MACHINE: Useful when operations on QIs and HIs are faster
  62      than on an SI */
  63
  64 typedef SFtype __attribute__ ((mode (SF)));
  65 typedef DFtype __attribute__ ((mode (DF)));
  66
  67 typedef int HItype __attribute__ ((mode (HI)));
  68 typedef int SItype __attribute__ ((mode (SI)));
  69 typedef int DItype __attribute__ ((mode (DI)));
  70
  71 /* The type of the result of a fp compare */
  72 #ifndef CMPtype
  73 #define CMPtype SItype
  74 #endif
  75
  76 typedef unsigned int UHItype __attribute__ ((mode (HI)));
  77 typedef unsigned int USItype __attribute__ ((mode (SI)));
  78 typedef unsigned int UDItype __attribute__ ((mode (DI)));
  79
  80 #define MAX_SI_INT   ((SItype) ((unsigned) (~0)>>1))
  81 #define MAX_USI_INT  ((USItype) ~0)
  82
  83
  84 #ifdef FLOAT_ONLY
  85 #define NO_DI_MODE
  86 #endif
  87
  88 #ifdef FLOAT
  89 #       define NGARDS    7L
  90 #       define GARDROUND 0x3f
  91 #       define GARDMASK  0x7f
  92 #       define GARDMSB   0x40
  93 #       define EXPBITS 8
  94 #       define EXPBIAS 127
  95 #       define FRACBITS 23
  96 #       define EXPMAX (0xff)
  97 #       define QUIET_NAN 0x100000L
  98 #       define FRAC_NBITS 32
  99 #       define FRACHIGH  0x80000000L
 100 #       define FRACHIGH2 0xc0000000L
 101         typedef USItype fractype;
 102         typedef UHItype halffractype;
 103         typedef SFtype FLO_type;
 104         typedef SItype intfrac;
 105
 106 #else
 107 #       define PREFIXFPDP dp
 108 #       define PREFIXSFDF df
 109 #       define NGARDS 8L
 110 #       define GARDROUND 0x7f
 111 #       define GARDMASK  0xff
 112 #       define GARDMSB   0x80
 113 #       define EXPBITS 11
 114 #       define EXPBIAS 1023
 115 #       define FRACBITS 52
 116 #       define EXPMAX (0x7ff)
 117 #       define QUIET_NAN 0x8000000000000LL
 118 #       define FRAC_NBITS 64
 119 #       define FRACHIGH  0x8000000000000000LL
 120 #       define FRACHIGH2 0xc000000000000000LL
 121         typedef UDItype fractype;
 122         typedef USItype halffractype;
 123         typedef DFtype FLO_type;
 124         typedef DItype intfrac;
 125 #endif
 126
 127 #ifdef US_SOFTWARE_GOFAST
 128 #       ifdef FLOAT
 129 #               define add              fpadd
 130 #               define sub              fpsub
 131 #               define multiply         fpmul
 132 #               define divide           fpdiv
 133 #               define compare          fpcmp
 134 #               define si_to_float      sitofp
 135 #               define float_to_si      fptosi
 136 #               define float_to_usi     fptoui
 137 #               define negate           __negsf2
 138 #               define sf_to_df         fptodp
 139 #               define dptofp           dptofp
 140 #else
 141 #               define add              dpadd
 142 #               define sub              dpsub
 143 #               define multiply         dpmul
 144 #               define divide           dpdiv
 145 #               define compare          dpcmp
 146 #               define si_to_float      litodp
 147 #               define float_to_si      dptoli
 148 #               define float_to_usi     dptoul
 149 #               define negate           __negdf2
 150 #               define df_to_sf         dptofp
 151 #endif
 152 #else
 153 #       ifdef FLOAT
 154 #               define add              __addsf3
 155 #               define sub              __subsf3
 156 #               define multiply         __mulsf3
 157 #               define divide           __divsf3
 158 #               define compare          __cmpsf2
 159 #               define _eq_f2           __eqsf2
 160 #               define _ne_f2           __nesf2
 161 #               define _gt_f2           __gtsf2
 162 #               define _ge_f2           __gesf2
 163 #               define _lt_f2           __ltsf2
 164 #               define _le_f2           __lesf2
 165 #               define si_to_float      __floatsisf
 166 #               define float_to_si      __fixsfsi
 167 #               define float_to_usi     __fixunssfsi
 168 #               define negate           __negsf2
 169 #               define sf_to_df         __extendsfdf2
 170 #else
 171 #               define add              __adddf3
 172 #               define sub              __subdf3
 173 #               define multiply         __muldf3
 174 #               define divide           __divdf3
 175 #               define compare          __cmpdf2
 176 #               define _eq_f2           __eqdf2
 177 #               define _ne_f2           __nedf2
 178 #               define _gt_f2           __gtdf2
 179 #               define _ge_f2           __gedf2
 180 #               define _lt_f2           __ltdf2
 181 #               define _le_f2           __ledf2
 182 #               define si_to_float      __floatsidf
 183 #               define float_to_si      __fixdfsi
 184 #               define float_to_usi     __fixunsdfsi
 185 #               define negate           __negdf2
 186 #               define df_to_sf         __truncdfsf2
 187 #       endif
 188 #endif
 189
 190
 191 #define INLINE __inline__
 192
 193 /* Preserve the sticky-bit when shifting fractions to the right.  */
 194 #define LSHIFT(a) { a = (a & 1) | (a >> 1); }
 195
 196 /* numeric parameters */
 197 /* F_D_BITOFF is the number of bits offset between the MSB of the mantissa
 198    of a float and of a double. Assumes there are only two float types.
 199    (double::FRAC_BITS+double::NGARGS-(float::FRAC_BITS-float::NGARDS))
 200  */
 201 #define F_D_BITOFF (52+8-(23+7))
 202
 203
 204 #define NORMAL_EXPMIN (-(EXPBIAS)+1)
 205 #define IMPLICIT_1 (1LL<<(FRACBITS+NGARDS))
 206 #define IMPLICIT_2 (1LL<<(FRACBITS+1+NGARDS))
 207
 208 /* common types */
 209
 210 typedef enum
 211 {
 212   CLASS_SNAN,
 213   CLASS_QNAN,
 214   CLASS_ZERO,
 215   CLASS_NUMBER,
 216   CLASS_INFINITY
 217 } fp_class_type;
 218
 219 typedef struct
 220 {
 221 #ifdef SMALL_MACHINE
 222   char class;
 223   unsigned char sign;
 224   short normal_exp;
 225 #else
 226   fp_class_type class;
 227   unsigned int sign;
 228   int normal_exp;
 229 #endif
 230
 231   union
 232     {
 233       fractype ll;
 234       halffractype l[2];
 235     } fraction;
 236 } fp_number_type;
 237
 238 typedef union
 239 {
 240   FLO_type value;
 241 #ifdef _DEBUG_BITFLOAT
 242   int l[2];
 243 #endif
 244   struct
 245     {
 246 #ifndef FLOAT_BIT_ORDER_MISMATCH
 247       unsigned int sign:1 __attribute__ ((packed));
 248       unsigned int exp:EXPBITS __attribute__ ((packed));
 249       fractype fraction:FRACBITS __attribute__ ((packed));
 250 #else
 251       fractype fraction:FRACBITS __attribute__ ((packed));
 252       unsigned int exp:EXPBITS __attribute__ ((packed));
 253       unsigned int sign:1 __attribute__ ((packed));
 254 #endif
 255     }
 256   bits;
 257 }
 258 FLO_union_type;
 259
 260
 261 /* end of header */
 262
 263 /* IEEE "special" number predicates */
 264
 265 #ifdef NO_NANS
 266
 267 #define nan() 0
 268 #define isnan(x) 0
 269 #define isinf(x) 0
 270 #else
 271
 272 INLINE
 273 static fp_number_type *
 274 nan ()
 275 {
 276   static fp_number_type thenan;
 277
 278   return &thenan;
 279 }
 280
 281 INLINE
 282 static int
 283 isnan ( fp_number_type *  x)
 284 {
 285   return x->class == CLASS_SNAN || x->class == CLASS_QNAN;
 286 }
 287
 288 INLINE
 289 static int
 290 isinf ( fp_number_type *  x)
 291 {
 292   return x->class == CLASS_INFINITY;
 293 }
 294
 295 #endif
 296
 297 INLINE
 298 static int
 299 iszero ( fp_number_type *  x)
 300 {
 301   return x->class == CLASS_ZERO;
 302 }
 303
 304 INLINE
 305 static void
 306 flip_sign ( fp_number_type *  x)
 307 {
 308   x->sign = !x->sign;
 309 }
 310
 311 static FLO_type
 312 pack_d ( fp_number_type *  src)
 313 {
 314   FLO_union_type dst;
 315   fractype fraction = src->fraction.ll; /* wasn't unsigned before? */
 316
 317   dst.bits.sign = src->sign;
 318
 319   if (isnan (src))
 320     {
 321       dst.bits.exp = EXPMAX;
 322       dst.bits.fraction = src->fraction.ll;
 323       if (src->class == CLASS_QNAN || 1)
 324         {
 325           dst.bits.fraction |= QUIET_NAN;
 326         }
 327     }
 328   else if (isinf (src))
 329     {
 330       dst.bits.exp = EXPMAX;
 331       dst.bits.fraction = 0;
 332     }
 333   else if (iszero (src))
 334     {
 335       dst.bits.exp = 0;
 336       dst.bits.fraction = 0;
 337     }
 338   else if (fraction == 0)
 339     {
 340       dst.value = 0;
 341     }
 342   else
 343     {
 344       if (src->normal_exp < NORMAL_EXPMIN)
 345         {
 346           /* This number's exponent is too low to fit into the bits
 347              available in the number, so we'll store 0 in the exponent and
 348              shift the fraction to the right to make up for it.  */
 349
 350           int shift = NORMAL_EXPMIN - src->normal_exp;
 351
 352           dst.bits.exp = 0;
 353
 354           if (shift > FRAC_NBITS - NGARDS)
 355             {
 356               /* No point shifting, since it's more that 64 out.  */
 357               fraction = 0;
 358             }
 359           else
 360             {
 361               /* Shift by the value */
 362               fraction >>= shift;
 363             }
 364           fraction >>= NGARDS;
 365           dst.bits.fraction = fraction;
 366         }
 367       else if (src->normal_exp > EXPBIAS)
 368         {
 369           dst.bits.exp = EXPMAX;
 370           dst.bits.fraction = 0;
 371         }
 372       else
 373         {
 374           dst.bits.exp = src->normal_exp + EXPBIAS;
 375           /* IF the gard bits are the all zero, but the first, then we're
 376              half way between two numbers, choose the one which makes the
 377              lsb of the answer 0.  */
 378           if ((fraction & GARDMASK) == GARDMSB)
 379             {
 380               if (fraction & (1 << NGARDS))
 381                 fraction += GARDROUND + 1;
 382             }
 383           else
 384             {
 385               /* Add a one to the guards to round up */
 386               fraction += GARDROUND;
 387             }
 388           if (fraction >= IMPLICIT_2)
 389             {
 390               fraction >>= 1;
 391               dst.bits.exp += 1;
 392             }
 393           fraction >>= NGARDS;
 394           dst.bits.fraction = fraction;
 395         }
 396     }
 397   return dst.value;
 398 }
 399
 400 static void
 401 unpack_d (FLO_union_type * src, fp_number_type * dst)
 402 {
 403   fractype fraction = src->bits.fraction;
 404
 405   dst->sign = src->bits.sign;
 406   if (src->bits.exp == 0)
 407     {
 408       /* Hmm.  Looks like 0 */
 409       if (fraction == 0)
 410         {
 411           /* tastes like zero */
 412           dst->class = CLASS_ZERO;
 413         }
 414       else
 415         {
 416           /* Zero exponent with non zero fraction - it's denormalized,
 417              so there isn't a leading implicit one - we'll shift it so
 418              it gets one.  */
 419           dst->normal_exp = src->bits.exp - EXPBIAS + 1;
 420           fraction <<= NGARDS;
 421
 422           dst->class = CLASS_NUMBER;
 423 #if 1
 424           while (fraction < IMPLICIT_1)
 425             {
 426               fraction <<= 1;
 427               dst->normal_exp--;
 428             }
 429 #endif
 430           dst->fraction.ll = fraction;
 431         }
 432     }
 433   else if (src->bits.exp == EXPMAX)
 434     {
 435       /* Huge exponent*/
 436       if (fraction == 0)
 437         {
 438           /* Attatched to a zero fraction - means infinity */
 439           dst->class = CLASS_INFINITY;
 440         }
 441       else
 442         {
 443           /* Non zero fraction, means nan */
 444           if (dst->sign)
 445             {
 446               dst->class = CLASS_SNAN;
 447             }
 448           else
 449             {
 450               dst->class = CLASS_QNAN;
 451             }
 452           /* Keep the fraction part as the nan number */
 453           dst->fraction.ll = fraction;
 454         }
 455     }
 456   else
 457     {
 458       /* Nothing strange about this number */
 459       dst->normal_exp = src->bits.exp - EXPBIAS;
 460       dst->class = CLASS_NUMBER;
 461       dst->fraction.ll = (fraction << NGARDS) | IMPLICIT_1;
 462     }
 463 }
 464
 465 static fp_number_type *
 466 _fpadd_parts (fp_number_type * a,
 467               fp_number_type * b,
 468               fp_number_type * tmp)
 469 {
 470   intfrac tfraction;
 471
 472   /* Put commonly used fields in local variables.  */
 473   int a_normal_exp;
 474   int b_normal_exp;
 475   fractype a_fraction;
 476   fractype b_fraction;
 477
 478   if (isnan (a))
 479     {
 480       return a;
 481     }
 482   if (isnan (b))
 483     {
 484       return b;
 485     }
 486   if (isinf (a))
 487     {
 488       /* Adding infinities with opposite signs yields a NaN.  */
 489       if (isinf (b) && a->sign != b->sign)
 490         return nan ();
 491       return a;
 492     }
 493   if (isinf (b))
 494     {
 495       return b;
 496     }
 497   if (iszero (b))
 498     {
 499       return a;
 500     }
 501   if (iszero (a))
 502     {
 503       return b;
 504     }
 505
 506   /* Got two numbers. shift the smaller and increment the exponent till
 507      they're the same */
 508   {
 509     int diff;
 510
 511     a_normal_exp = a->normal_exp;
 512     b_normal_exp = b->normal_exp;
 513     a_fraction = a->fraction.ll;
 514     b_fraction = b->fraction.ll;
 515
 516     diff = a_normal_exp - b_normal_exp;
 517
 518     if (diff < 0)
 519       diff = -diff;
 520     if (diff < FRAC_NBITS)
 521       {
 522         /* ??? This does shifts one bit at a time.  Optimize.  */
 523         while (a_normal_exp > b_normal_exp)
 524           {
 525             b_normal_exp++;
 526             LSHIFT (b_fraction);
 527           }
 528         while (b_normal_exp > a_normal_exp)
 529           {
 530             a_normal_exp++;
 531             LSHIFT (a_fraction);
 532           }
 533       }
 534     else
 535       {
 536         /* Somethings's up.. choose the biggest */
 537         if (a_normal_exp > b_normal_exp)
 538           {
 539             b_normal_exp = a_normal_exp;
 540             b_fraction = 0;
 541           }
 542         else
 543           {
 544             a_normal_exp = b_normal_exp;
 545             a_fraction = 0;
 546           }
 547       }
 548   }
 549
 550   if (a->sign != b->sign)
 551     {
 552       if (a->sign)
 553         {
 554           tfraction = -a_fraction + b_fraction;
 555         }
 556       else
 557         {
 558           tfraction = a_fraction - b_fraction;
 559         }
 560       if (tfraction > 0)
 561         {
 562           tmp->sign = 0;
 563           tmp->normal_exp = a_normal_exp;
 564           tmp->fraction.ll = tfraction;
 565         }
 566       else
 567         {
 568           tmp->sign = 1;
 569           tmp->normal_exp = a_normal_exp;
 570           tmp->fraction.ll = -tfraction;
 571         }
 572       /* and renomalize it */
 573
 574       while (tmp->fraction.ll < IMPLICIT_1 && tmp->fraction.ll)
 575         {
 576           tmp->fraction.ll <<= 1;
 577           tmp->normal_exp--;
 578         }
 579     }
 580   else
 581     {
 582       tmp->sign = a->sign;
 583       tmp->normal_exp = a_normal_exp;
 584       tmp->fraction.ll = a_fraction + b_fraction;
 585     }
 586   tmp->class = CLASS_NUMBER;
 587   /* Now the fraction is added, we have to shift down to renormalize the
 588      number */
 589
 590   if (tmp->fraction.ll >= IMPLICIT_2)
 591     {
 592       LSHIFT (tmp->fraction.ll);
 593       tmp->normal_exp++;
 594     }
 595   return tmp;
 596
 597 }
 598
 599 FLO_type
 600 add (FLO_type arg_a, FLO_type arg_b)
 601 {
 602   fp_number_type a;
 603   fp_number_type b;
 604   fp_number_type tmp;
 605   fp_number_type *res;
 606
 607   unpack_d ((FLO_union_type *) & arg_a, &a);
 608   unpack_d ((FLO_union_type *) & arg_b, &b);
 609
 610   res = _fpadd_parts (&a, &b, &tmp);
 611
 612   return pack_d (res);
 613 }
 614
 615 FLO_type
 616 sub (FLO_type arg_a, FLO_type arg_b)
 617 {
 618   fp_number_type a;
 619   fp_number_type b;
 620   fp_number_type tmp;
 621   fp_number_type *res;
 622
 623   unpack_d ((FLO_union_type *) & arg_a, &a);
 624   unpack_d ((FLO_union_type *) & arg_b, &b);
 625
 626   b.sign ^= 1;
 627
 628   res = _fpadd_parts (&a, &b, &tmp);
 629
 630   return pack_d (res);
 631 }
 632
 633 static fp_number_type *
 634 _fpmul_parts ( fp_number_type *  a,
 635                fp_number_type *  b,
 636                fp_number_type * tmp)
 637 {
 638   fractype low = 0;
 639   fractype high = 0;
 640
 641   if (isnan (a))
 642     {
 643       a->sign = a->sign != b->sign;
 644       return a;
 645     }
 646   if (isnan (b))
 647     {
 648       b->sign = a->sign != b->sign;
 649       return b;
 650     }
 651   if (isinf (a))
 652     {
 653       if (iszero (b))
 654         return nan ();
 655       a->sign = a->sign != b->sign;
 656       return a;
 657     }
 658   if (isinf (b))
 659     {
 660       if (iszero (a))
 661         {
 662           return nan ();
 663         }
 664       b->sign = a->sign != b->sign;
 665       return b;
 666     }
 667   if (iszero (a))
 668     {
 669       a->sign = a->sign != b->sign;
 670       return a;
 671     }
 672   if (iszero (b))
 673     {
 674       b->sign = a->sign != b->sign;
 675       return b;
 676     }
 677
 678   /* Calculate the mantissa by multiplying both 64bit numbers to get a
 679      128 bit number */
 680   {
 681     fractype x = a->fraction.ll;
 682     fractype ylow = b->fraction.ll;
 683     fractype yhigh = 0;
 684     int bit;
 685
 686 #if defined(NO_DI_MODE)
 687     {
 688       /* ??? This does multiplies one bit at a time.  Optimize.  */
 689       for (bit = 0; bit < FRAC_NBITS; bit++)
 690         {
 691           int carry;
 692
 693           if (x & 1)
 694             {
 695               carry = (low += ylow) < ylow;
 696               high += yhigh + carry;
 697             }
 698           yhigh <<= 1;
 699           if (ylow & FRACHIGH)
 700             {
 701               yhigh |= 1;
 702             }
 703           ylow <<= 1;
 704           x >>= 1;
 705         }
 706     }
 707 #elif defined(FLOAT)
 708     {
 709       /* Multiplying two 32 bit numbers to get a 64 bit number  on
 710         a machine with DI, so we're safe */
 711
 712       DItype answer = (DItype)(a->fraction.ll) * (DItype)(b->fraction.ll);
 713
 714       high = answer >> 32;
 715       low = answer;
 716     }
 717 #else
 718     /* Doing a 64*64 to 128 */
 719     {
 720       UDItype nl = a->fraction.ll & 0xffffffff;
 721       UDItype nh = a->fraction.ll >> 32;
 722       UDItype ml = b->fraction.ll & 0xffffffff;
 723       UDItype mh = b->fraction.ll >>32;
 724       UDItype pp_ll = ml * nl;
 725       UDItype pp_hl = mh * nl;
 726       UDItype pp_lh = ml * nh;
 727       UDItype pp_hh = mh * nh;
 728       UDItype res2 = 0;
 729       UDItype res0 = 0;
 730       UDItype ps_hh__ = pp_hl + pp_lh;
 731       if (ps_hh__ < pp_hl)
 732         res2 += 0x100000000LL;
 733       pp_hl = (ps_hh__ << 32) & 0xffffffff00000000LL;
 734       res0 = pp_ll + pp_hl;
 735       if (res0 < pp_ll)
 736         res2++;
 737       res2 += ((ps_hh__ >> 32) & 0xffffffffL) + pp_hh;
 738       high = res2;
 739       low = res0;
 740     }
 741 #endif
 742   }
 743
 744   tmp->normal_exp = a->normal_exp + b->normal_exp;
 745   tmp->sign = a->sign != b->sign;
 746 #ifdef FLOAT
 747   tmp->normal_exp += 2;         /* ??????????????? */
 748 #else
 749   tmp->normal_exp += 4;         /* ??????????????? */
 750 #endif
 751   while (high >= IMPLICIT_2)
 752     {
 753       tmp->normal_exp++;
 754       if (high & 1)
 755         {
 756           low >>= 1;
 757           low |= FRACHIGH;
 758         }
 759       high >>= 1;
 760     }
 761   while (high < IMPLICIT_1)
 762     {
 763       tmp->normal_exp--;
 764
 765       high <<= 1;
 766       if (low & FRACHIGH)
 767         high |= 1;
 768       low <<= 1;
 769     }
 770   /* rounding is tricky. if we only round if it won't make us round later. */
 771 #if 0
 772   if (low & FRACHIGH2)
 773     {
 774       if (((high & GARDMASK) != GARDMSB)
 775           && (((high + 1) & GARDMASK) == GARDMSB))
 776         {
 777           /* don't round, it gets done again later. */
 778         }
 779       else
 780         {
 781           high++;
 782         }
 783     }
 784 #endif
 785   if ((high & GARDMASK) == GARDMSB)
 786     {
 787       if (high & (1 << NGARDS))
 788         {
 789           /* half way, so round to even */
 790           high += GARDROUND + 1;
 791         }
 792       else if (low)
 793         {
 794           /* but we really weren't half way */
 795           high += GARDROUND + 1;
 796         }
 797     }
 798   tmp->fraction.ll = high;
 799   tmp->class = CLASS_NUMBER;
 800   return tmp;
 801 }
 802
 803 FLO_type
 804 multiply (FLO_type arg_a, FLO_type arg_b)
 805 {
 806   fp_number_type a;
 807   fp_number_type b;
 808   fp_number_type tmp;
 809   fp_number_type *res;
 810
 811   unpack_d ((FLO_union_type *) & arg_a, &a);
 812   unpack_d ((FLO_union_type *) & arg_b, &b);
 813
 814   res = _fpmul_parts (&a, &b, &tmp);
 815
 816   return pack_d (res);
 817 }
 818
 819 static fp_number_type *
 820 _fpdiv_parts (fp_number_type * a,
 821               fp_number_type * b,
 822               fp_number_type * tmp)
 823 {
 824   fractype low = 0;
 825   fractype high = 0;
 826   fractype r0, r1, y0, y1, bit;
 827   fractype q;
 828   fractype numerator;
 829   fractype denominator;
 830   fractype quotient;
 831   fractype remainder;
 832
 833   if (isnan (a))
 834     {
 835       return a;
 836     }
 837   if (isnan (b))
 838     {
 839       return b;
 840     }
 841   if (isinf (a) || iszero (a))
 842     {
 843       if (a->class == b->class)
 844         return nan ();
 845       return a;
 846     }
 847   a->sign = a->sign ^ b->sign;
 848
 849   if (isinf (b))
 850     {
 851       a->fraction.ll = 0;
 852       a->normal_exp = 0;
 853       return a;
 854     }
 855   if (iszero (b))
 856     {
 857       a->class = CLASS_INFINITY;
 858       return b;
 859     }
 860
 861   /* Calculate the mantissa by multiplying both 64bit numbers to get a
 862      128 bit number */
 863   {
 864     int carry;
 865     intfrac d0, d1;             /* weren't unsigned before ??? */
 866
 867     /* quotient =
 868        ( numerator / denominator) * 2^(numerator exponent -  denominator exponent)
 869      */
 870
 871     a->normal_exp = a->normal_exp - b->normal_exp;
 872     numerator = a->fraction.ll;
 873     denominator = b->fraction.ll;
 874
 875     if (numerator < denominator)
 876       {
 877         /* Fraction will be less than 1.0 */
 878         numerator *= 2;
 879         a->normal_exp--;
 880       }
 881     bit = IMPLICIT_1;
 882     quotient = 0;
 883     /* ??? Does divide one bit at a time.  Optimize.  */
 884     while (bit)
 885       {
 886         if (numerator >= denominator)
 887           {
 888             quotient |= bit;
 889             numerator -= denominator;
 890           }
 891         bit >>= 1;
 892         numerator *= 2;
 893       }
 894
 895     if ((quotient & GARDMASK) == GARDMSB)
 896       {
 897         if (quotient & (1 << NGARDS))
 898           {
 899             /* half way, so round to even */
 900             quotient += GARDROUND + 1;
 901           }
 902         else if (numerator)
 903           {
 904             /* but we really weren't half way, more bits exist */
 905             quotient += GARDROUND + 1;
 906           }
 907       }
 908
 909     a->fraction.ll = quotient;
 910     return (a);
 911   }
 912 }
 913
 914 FLO_type
 915 divide (FLO_type arg_a, FLO_type arg_b)
 916 {
 917   fp_number_type a;
 918   fp_number_type b;
 919   fp_number_type tmp;
 920   fp_number_type *res;
 921
 922   unpack_d ((FLO_union_type *) & arg_a, &a);
 923   unpack_d ((FLO_union_type *) & arg_b, &b);
 924
 925   res = _fpdiv_parts (&a, &b, &tmp);
 926
 927   return pack_d (res);
 928 }
 929
 930 /* according to the demo, fpcmp returns a comparison with 0... thus
 931    a<b -> -1
 932    a==b -> 0
 933    a>b -> +1
 934  */
 935
 936 static int
 937 _fpcmp_parts (fp_number_type * a, fp_number_type * b)
 938 {
 939 #if 0
 940   /* either nan -> unordered. Must be checked outside of this routine. */
 941   if (isnan (a) && isnan (b))
 942     {
 943       return 1;                 /* still unordered! */
 944     }
 945 #endif
 946
 947   if (isnan (a) || isnan (b))
 948     {
 949       return 1;                 /* how to indicate unordered compare? */
 950     }
 951   if (isinf (a) && isinf (b))
 952     {
 953       /* +inf > -inf, but +inf != +inf */
 954       /* b    \a| +inf(0)| -inf(1)
 955        ______\+--------+--------
 956        +inf(0)| a==b(0)| a<b(-1)
 957        -------+--------+--------
 958        -inf(1)| a>b(1) | a==b(0)
 959        -------+--------+--------
 960        So since unordered must be non zero, just line up the columns...
 961        */
 962       return b->sign - a->sign;
 963     }
 964   /* but not both... */
 965   if (isinf (a))
 966     {
 967       return a->sign ? -1 : 1;
 968     }
 969   if (isinf (b))
 970     {
 971       return b->sign ? 1 : -1;
 972     }
 973   if (iszero (a) && iszero (b))
 974     {
 975       return 0;
 976     }
 977   if (iszero (a))
 978     {
 979       return b->sign ? 1 : -1;
 980     }
 981   if (iszero (b))
 982     {
 983       return a->sign ? -1 : 1;
 984     }
 985   /* now both are "normal". */
 986   if (a->sign != b->sign)
 987     {
 988       /* opposite signs */
 989       return a->sign ? -1 : 1;
 990     }
 991   /* same sign; exponents? */
 992   if (a->normal_exp > b->normal_exp)
 993     {
 994       return a->sign ? -1 : 1;
 995     }
 996   if (a->normal_exp < b->normal_exp)
 997     {
 998       return a->sign ? 1 : -1;
 999     }
1000   /* same exponents; check size. */
1001   if (a->fraction.ll > b->fraction.ll)
1002     {
1003       return a->sign ? -1 : 1;
1004     }
1005   if (a->fraction.ll < b->fraction.ll)
1006     {
1007       return a->sign ? 1 : -1;
1008     }
1009   /* after all that, they're equal. */
1010   return 0;
1011 }
1012
1013 CMPtype
1014 compare (FLO_type arg_a, FLO_type arg_b)
1015 {
1016   fp_number_type a;
1017   fp_number_type b;
1018
1019   unpack_d ((FLO_union_type *) & arg_a, &a);
1020   unpack_d ((FLO_union_type *) & arg_b, &b);
1021
1022   return _fpcmp_parts (&a, &b);
1023 }
1024
1025 #ifndef US_SOFTWARE_GOFAST
1026
1027 /* These should be optimized for their specific tasks someday.  */
1028
1029 CMPtype
1030 _eq_f2 (FLO_type arg_a, FLO_type arg_b)
1031 {
1032   fp_number_type a;
1033   fp_number_type b;
1034
1035   unpack_d ((FLO_union_type *) & arg_a, &a);
1036   unpack_d ((FLO_union_type *) & arg_b, &b);
1037
1038   if (isnan (&a) || isnan (&b))
1039     return 1;                   /* false, truth == 0 */
1040
1041   return _fpcmp_parts (&a, &b) ;
1042 }
1043
1044 CMPtype
1045 _ne_f2 (FLO_type arg_a, FLO_type arg_b)
1046 {
1047   fp_number_type a;
1048   fp_number_type b;
1049
1050   unpack_d ((FLO_union_type *) & arg_a, &a);
1051   unpack_d ((FLO_union_type *) & arg_b, &b);
1052
1053   if (isnan (&a) || isnan (&b))
1054     return 1;                   /* true, truth != 0 */
1055
1056   return  _fpcmp_parts (&a, &b) ;
1057 }
1058
1059 CMPtype
1060 _gt_f2 (FLO_type arg_a, FLO_type arg_b)
1061 {
1062   fp_number_type a;
1063   fp_number_type b;
1064
1065   unpack_d ((FLO_union_type *) & arg_a, &a);
1066   unpack_d ((FLO_union_type *) & arg_b, &b);
1067
1068   if (isnan (&a) || isnan (&b))
1069     return -1;                  /* false, truth > 0 */
1070
1071   return _fpcmp_parts (&a, &b);
1072 }
1073
1074 CMPtype
1075 _ge_f2 (FLO_type arg_a, FLO_type arg_b)
1076 {
1077   fp_number_type a;
1078   fp_number_type b;
1079
1080   unpack_d ((FLO_union_type *) & arg_a, &a);
1081   unpack_d ((FLO_union_type *) & arg_b, &b);
1082
1083   if (isnan (&a) || isnan (&b))
1084     return -1;                  /* false, truth >= 0 */
1085   return _fpcmp_parts (&a, &b) ;
1086 }
1087
1088 CMPtype
1089 _lt_f2 (FLO_type arg_a, FLO_type arg_b)
1090 {
1091   fp_number_type a;
1092   fp_number_type b;
1093
1094   unpack_d ((FLO_union_type *) & arg_a, &a);
1095   unpack_d ((FLO_union_type *) & arg_b, &b);
1096
1097   if (isnan (&a) || isnan (&b))
1098     return 1;                   /* false, truth < 0 */
1099
1100   return _fpcmp_parts (&a, &b);
1101 }
1102
1103 CMPtype
1104 _le_f2 (FLO_type arg_a, FLO_type arg_b)
1105 {
1106   fp_number_type a;
1107   fp_number_type b;
1108
1109   unpack_d ((FLO_union_type *) & arg_a, &a);
1110   unpack_d ((FLO_union_type *) & arg_b, &b);
1111
1112   if (isnan (&a) || isnan (&b))
1113     return 1;                   /* false, truth <= 0 */
1114
1115   return _fpcmp_parts (&a, &b) ;
1116 }
1117
1118 #endif /* ! US_SOFTWARE_GOFAST */
1119
1120 FLO_type
1121 si_to_float (SItype arg_a)
1122 {
1123   fp_number_type in;
1124
1125   in.class = CLASS_NUMBER;
1126   in.sign = arg_a < 0;
1127   if (!arg_a)
1128     {
1129       in.class = CLASS_ZERO;
1130     }
1131   else
1132     {
1133       in.normal_exp = FRACBITS + NGARDS;
1134       if (in.sign)
1135         {
1136           /* Special case for minint, since there is no +ve integer
1137              representation for it */
1138           if (arg_a == 0x80000000)
1139             {
1140               return -2147483648.0;
1141             }
1142           in.fraction.ll = (-arg_a);
1143         }
1144       else
1145         in.fraction.ll = arg_a;
1146
1147       while (in.fraction.ll < (1LL << (FRACBITS + NGARDS)))
1148         {
1149           in.fraction.ll <<= 1;
1150           in.normal_exp -= 1;
1151         }
1152     }
1153   return pack_d (&in);
1154 }
1155
1156 SItype
1157 float_to_si (FLO_type arg_a)
1158 {
1159   fp_number_type a;
1160   SItype tmp;
1161
1162   unpack_d ((FLO_union_type *) & arg_a, &a);
1163   if (iszero (&a))
1164     return 0;
1165   if (isnan (&a))
1166     return 0;
1167   /* get reasonable MAX_SI_INT... */
1168   if (isinf (&a))
1169     return a.sign ? MAX_SI_INT : (-MAX_SI_INT)-1;
1170   /* it is a number, but a small one */
1171   if (a.normal_exp < 0)
1172     return 0;
1173   if (a.normal_exp > 30)
1174     return a.sign ? (-MAX_SI_INT)-1 : MAX_SI_INT;
1175   tmp = a.fraction.ll >> ((FRACBITS + NGARDS) - a.normal_exp);
1176   return a.sign ? (-tmp) : (tmp);
1177 }
1178
1179 #ifdef US_SOFTWARE_GOFAST
1180 /* While libgcc2.c defines its own __fixunssfsi and __fixunsdfsi routines,
1181    we also define them for GOFAST because the ones in libgcc2.c have the
1182    wrong names and I'd rather define these here and keep GOFAST CYG-LOC's
1183    out of libgcc2.c.  We can't define these here if not GOFAST because then
1184    there'd be duplicate copies.  */
1185
1186 USItype
1187 float_to_usi (FLO_type arg_a)
1188 {
1189   fp_number_type a;
1190
1191   unpack_d ((FLO_union_type *) & arg_a, &a);
1192   if (iszero (&a))
1193     return 0;
1194   if (isnan (&a))
1195     return 0;
1196   /* get reasonable MAX_USI_INT... */
1197   if (isinf (&a))
1198     return a.sign ? MAX_USI_INT : 0;
1199   /* it is a negative number */
1200   if (a.sign)
1201     return 0;
1202   /* it is a number, but a small one */
1203   if (a.normal_exp < 0)
1204     return 0;
1205   if (a.normal_exp > 31)
1206     return MAX_USI_INT;
1207   else if (a.normal_exp > (FRACBITS + NGARDS))
1208     return a.fraction.ll << ((FRACBITS + NGARDS) - a.normal_exp);
1209   else
1210     return a.fraction.ll >> ((FRACBITS + NGARDS) - a.normal_exp);
1211 }
1212 #endif
1213
1214 FLO_type
1215 negate (FLO_type arg_a)
1216 {
1217   fp_number_type a;
1218
1219   unpack_d ((FLO_union_type *) & arg_a, &a);
1220   flip_sign (&a);
1221   return pack_d (&a);
1222 }
1223
1224 #ifdef FLOAT
1225
1226 SFtype
1227 __make_fp(fp_class_type class,
1228              unsigned int sign,
1229              int exp,
1230              USItype frac)
1231 {
1232   fp_number_type in;
1233
1234   in.class = class;
1235   in.sign = sign;
1236   in.normal_exp = exp;
1237   in.fraction.ll = frac;
1238   return pack_d (&in);
1239 }
1240
1241 #ifndef FLOAT_ONLY
1242
1243 /* This enables one to build an fp library that supports float but not double.
1244    Otherwise, we would get an undefined reference to __make_dp.
1245    This is needed for some 8-bit ports that can't handle well values that
1246    are 8-bytes in size, so we just don't support double for them at all.  */
1247
1248 extern DFtype __make_dp (fp_class_type, unsigned int, int, UDItype frac);
1249
1250 DFtype
1251 sf_to_df (SFtype arg_a)
1252 {
1253   fp_number_type in;
1254
1255   unpack_d ((FLO_union_type *) & arg_a, &in);
1256   return __make_dp (in.class, in.sign, in.normal_exp,
1257                     ((UDItype) in.fraction.ll) << F_D_BITOFF);
1258 }
1259
1260 #endif
1261 #endif
1262
1263 #ifndef FLOAT
1264
1265 extern SFtype __make_fp (fp_class_type, unsigned int, int, USItype);
1266
1267 DFtype
1268 __make_dp (fp_class_type class, unsigned int sign, int exp, UDItype frac)
1269 {
1270   fp_number_type in;
1271
1272   in.class = class;
1273   in.sign = sign;
1274   in.normal_exp = exp;
1275   in.fraction.ll = frac;
1276   return pack_d (&in);
1277 }
1278
1279 SFtype
1280 df_to_sf (DFtype arg_a)
1281 {
1282   fp_number_type in;
1283
1284   unpack_d ((FLO_union_type *) & arg_a, &in);
1285   return __make_fp (in.class, in.sign, in.normal_exp,
1286                     in.fraction.ll >> F_D_BITOFF);
1287 }
1288
1289 #endif