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