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       return a;
 489     }
 490   if (isinf (b))
 491     {
 492       return b;
 493     }
 494   if (iszero (b))
 495     {
 496       return a;
 497     }
 498   if (iszero (a))
 499     {
 500       return b;
 501     }
 502
 503   /* Got two numbers. shift the smaller and increment the exponent till
 504      they're the same */
 505   {
 506     int diff;
 507
 508     a_normal_exp = a->normal_exp;
 509     b_normal_exp = b->normal_exp;
 510     a_fraction = a->fraction.ll;
 511     b_fraction = b->fraction.ll;
 512
 513     diff = a_normal_exp - b_normal_exp;
 514
 515     if (diff < 0)
 516       diff = -diff;
 517     if (diff < FRAC_NBITS)
 518       {
 519         /* ??? This does shifts one bit at a time.  Optimize.  */
 520         while (a_normal_exp > b_normal_exp)
 521           {
 522             b_normal_exp++;
 523             LSHIFT (b_fraction);
 524           }
 525         while (b_normal_exp > a_normal_exp)
 526           {
 527             a_normal_exp++;
 528             LSHIFT (a_fraction);
 529           }
 530       }
 531     else
 532       {
 533         /* Somethings's up.. choose the biggest */
 534         if (a_normal_exp > b_normal_exp)
 535           {
 536             b_normal_exp = a_normal_exp;
 537             b_fraction = 0;
 538           }
 539         else
 540           {
 541             a_normal_exp = b_normal_exp;
 542             a_fraction = 0;
 543           }
 544       }
 545   }
 546
 547   if (a->sign != b->sign)
 548     {
 549       if (a->sign)
 550         {
 551           tfraction = -a_fraction + b_fraction;
 552         }
 553       else
 554         {
 555           tfraction = a_fraction - b_fraction;
 556         }
 557       if (tfraction > 0)
 558         {
 559           tmp->sign = 0;
 560           tmp->normal_exp = a_normal_exp;
 561           tmp->fraction.ll = tfraction;
 562         }
 563       else
 564         {
 565           tmp->sign = 1;
 566           tmp->normal_exp = a_normal_exp;
 567           tmp->fraction.ll = -tfraction;
 568         }
 569       /* and renomalize it */
 570
 571       while (tmp->fraction.ll < IMPLICIT_1 && tmp->fraction.ll)
 572         {
 573           tmp->fraction.ll <<= 1;
 574           tmp->normal_exp--;
 575         }
 576     }
 577   else
 578     {
 579       tmp->sign = a->sign;
 580       tmp->normal_exp = a_normal_exp;
 581       tmp->fraction.ll = a_fraction + b_fraction;
 582     }
 583   tmp->class = CLASS_NUMBER;
 584   /* Now the fraction is added, we have to shift down to renormalize the
 585      number */
 586
 587   if (tmp->fraction.ll >= IMPLICIT_2)
 588     {
 589       LSHIFT (tmp->fraction.ll);
 590       tmp->normal_exp++;
 591     }
 592   return tmp;
 593
 594 }
 595
 596 FLO_type
 597 add (FLO_type arg_a, FLO_type arg_b)
 598 {
 599   fp_number_type a;
 600   fp_number_type b;
 601   fp_number_type tmp;
 602   fp_number_type *res;
 603
 604   unpack_d ((FLO_union_type *) & arg_a, &a);
 605   unpack_d ((FLO_union_type *) & arg_b, &b);
 606
 607   res = _fpadd_parts (&a, &b, &tmp);
 608
 609   return pack_d (res);
 610 }
 611
 612 FLO_type
 613 sub (FLO_type arg_a, FLO_type arg_b)
 614 {
 615   fp_number_type a;
 616   fp_number_type b;
 617   fp_number_type tmp;
 618   fp_number_type *res;
 619
 620   unpack_d ((FLO_union_type *) & arg_a, &a);
 621   unpack_d ((FLO_union_type *) & arg_b, &b);
 622
 623   b.sign ^= 1;
 624
 625   res = _fpadd_parts (&a, &b, &tmp);
 626
 627   return pack_d (res);
 628 }
 629
 630 static fp_number_type *
 631 _fpmul_parts ( fp_number_type *  a,
 632                fp_number_type *  b,
 633                fp_number_type * tmp)
 634 {
 635   fractype low = 0;
 636   fractype high = 0;
 637
 638   if (isnan (a))
 639     {
 640       a->sign = a->sign != b->sign;
 641       return a;
 642     }
 643   if (isnan (b))
 644     {
 645       b->sign = a->sign != b->sign;
 646       return b;
 647     }
 648   if (isinf (a))
 649     {
 650       if (iszero (b))
 651         return nan ();
 652       a->sign = a->sign != b->sign;
 653       return a;
 654     }
 655   if (isinf (b))
 656     {
 657       if (iszero (a))
 658         {
 659           return nan ();
 660         }
 661       b->sign = a->sign != b->sign;
 662       return b;
 663     }
 664   if (iszero (a))
 665     {
 666       a->sign = a->sign != b->sign;
 667       return a;
 668     }
 669   if (iszero (b))
 670     {
 671       b->sign = a->sign != b->sign;
 672       return b;
 673     }
 674
 675   /* Calculate the mantissa by multiplying both 64bit numbers to get a
 676      128 bit number */
 677   {
 678     fractype x = a->fraction.ll;
 679     fractype ylow = b->fraction.ll;
 680     fractype yhigh = 0;
 681     int bit;
 682
 683 #if defined(NO_DI_MODE)
 684     {
 685       /* ??? This does multiplies one bit at a time.  Optimize.  */
 686       for (bit = 0; bit < FRAC_NBITS; bit++)
 687         {
 688           int carry;
 689
 690           if (x & 1)
 691             {
 692               carry = (low += ylow) < ylow;
 693               high += yhigh + carry;
 694             }
 695           yhigh <<= 1;
 696           if (ylow & FRACHIGH)
 697             {
 698               yhigh |= 1;
 699             }
 700           ylow <<= 1;
 701           x >>= 1;
 702         }
 703     }
 704 #elif defined(FLOAT)
 705     {
 706       /* Multiplying two 32 bit numbers to get a 64 bit number  on
 707         a machine with DI, so we're safe */
 708
 709       DItype answer = (DItype)(a->fraction.ll) * (DItype)(b->fraction.ll);
 710
 711       high = answer >> 32;
 712       low = answer;
 713     }
 714 #else
 715     /* Doing a 64*64 to 128 */
 716     {
 717       UDItype nl = a->fraction.ll & 0xffffffff;
 718       UDItype nh = a->fraction.ll >> 32;
 719       UDItype ml = b->fraction.ll & 0xffffffff;
 720       UDItype mh = b->fraction.ll >>32;
 721       UDItype pp_ll = ml * nl;
 722       UDItype pp_hl = mh * nl;
 723       UDItype pp_lh = ml * nh;
 724       UDItype pp_hh = mh * nh;
 725       UDItype res2 = 0;
 726       UDItype res0 = 0;
 727       UDItype ps_hh__ = pp_hl + pp_lh;
 728       if (ps_hh__ < pp_hl)
 729         res2 += 0x100000000LL;
 730       pp_hl = (ps_hh__ << 32) & 0xffffffff00000000LL;
 731       res0 = pp_ll + pp_hl;
 732       if (res0 < pp_ll)
 733         res2++;
 734       res2 += ((ps_hh__ >> 32) & 0xffffffffL) + pp_hh;
 735       high = res2;
 736       low = res0;
 737     }
 738 #endif
 739   }
 740
 741   tmp->normal_exp = a->normal_exp + b->normal_exp;
 742   tmp->sign = a->sign != b->sign;
 743 #ifdef FLOAT
 744   tmp->normal_exp += 2;         /* ??????????????? */
 745 #else
 746   tmp->normal_exp += 4;         /* ??????????????? */
 747 #endif
 748   while (high >= IMPLICIT_2)
 749     {
 750       tmp->normal_exp++;
 751       if (high & 1)
 752         {
 753           low >>= 1;
 754           low |= FRACHIGH;
 755         }
 756       high >>= 1;
 757     }
 758   while (high < IMPLICIT_1)
 759     {
 760       tmp->normal_exp--;
 761
 762       high <<= 1;
 763       if (low & FRACHIGH)
 764         high |= 1;
 765       low <<= 1;
 766     }
 767   /* rounding is tricky. if we only round if it won't make us round later. */
 768 #if 0
 769   if (low & FRACHIGH2)
 770     {
 771       if (((high & GARDMASK) != GARDMSB)
 772           && (((high + 1) & GARDMASK) == GARDMSB))
 773         {
 774           /* don't round, it gets done again later. */
 775         }
 776       else
 777         {
 778           high++;
 779         }
 780     }
 781 #endif
 782   if ((high & GARDMASK) == GARDMSB)
 783     {
 784       if (high & (1 << NGARDS))
 785         {
 786           /* half way, so round to even */
 787           high += GARDROUND + 1;
 788         }
 789       else if (low)
 790         {
 791           /* but we really weren't half way */
 792           high += GARDROUND + 1;
 793         }
 794     }
 795   tmp->fraction.ll = high;
 796   tmp->class = CLASS_NUMBER;
 797   return tmp;
 798 }
 799
 800 FLO_type
 801 multiply (FLO_type arg_a, FLO_type arg_b)
 802 {
 803   fp_number_type a;
 804   fp_number_type b;
 805   fp_number_type tmp;
 806   fp_number_type *res;
 807
 808   unpack_d ((FLO_union_type *) & arg_a, &a);
 809   unpack_d ((FLO_union_type *) & arg_b, &b);
 810
 811   res = _fpmul_parts (&a, &b, &tmp);
 812
 813   return pack_d (res);
 814 }
 815
 816 static fp_number_type *
 817 _fpdiv_parts (fp_number_type * a,
 818               fp_number_type * b,
 819               fp_number_type * tmp)
 820 {
 821   fractype low = 0;
 822   fractype high = 0;
 823   fractype r0, r1, y0, y1, bit;
 824   fractype q;
 825   fractype numerator;
 826   fractype denominator;
 827   fractype quotient;
 828   fractype remainder;
 829
 830   if (isnan (a))
 831     {
 832       return a;
 833     }
 834   if (isnan (b))
 835     {
 836       return b;
 837     }
 838   if (isinf (a) || iszero (a))
 839     {
 840       if (a->class == b->class)
 841         return nan ();
 842       return a;
 843     }
 844   a->sign = a->sign ^ b->sign;
 845
 846   if (isinf (b))
 847     {
 848       a->fraction.ll = 0;
 849       a->normal_exp = 0;
 850       return a;
 851     }
 852   if (iszero (b))
 853     {
 854       a->class = CLASS_INFINITY;
 855       return b;
 856     }
 857
 858   /* Calculate the mantissa by multiplying both 64bit numbers to get a
 859      128 bit number */
 860   {
 861     int carry;
 862     intfrac d0, d1;             /* weren't unsigned before ??? */
 863
 864     /* quotient =
 865        ( numerator / denominator) * 2^(numerator exponent -  denominator exponent)
 866      */
 867
 868     a->normal_exp = a->normal_exp - b->normal_exp;
 869     numerator = a->fraction.ll;
 870     denominator = b->fraction.ll;
 871
 872     if (numerator < denominator)
 873       {
 874         /* Fraction will be less than 1.0 */
 875         numerator *= 2;
 876         a->normal_exp--;
 877       }
 878     bit = IMPLICIT_1;
 879     quotient = 0;
 880     /* ??? Does divide one bit at a time.  Optimize.  */
 881     while (bit)
 882       {
 883         if (numerator >= denominator)
 884           {
 885             quotient |= bit;
 886             numerator -= denominator;
 887           }
 888         bit >>= 1;
 889         numerator *= 2;
 890       }
 891
 892     if ((quotient & GARDMASK) == GARDMSB)
 893       {
 894         if (quotient & (1 << NGARDS))
 895           {
 896             /* half way, so round to even */
 897             quotient += GARDROUND + 1;
 898           }
 899         else if (numerator)
 900           {
 901             /* but we really weren't half way, more bits exist */
 902             quotient += GARDROUND + 1;
 903           }
 904       }
 905
 906     a->fraction.ll = quotient;
 907     return (a);
 908   }
 909 }
 910
 911 FLO_type
 912 divide (FLO_type arg_a, FLO_type arg_b)
 913 {
 914   fp_number_type a;
 915   fp_number_type b;
 916   fp_number_type tmp;
 917   fp_number_type *res;
 918
 919   unpack_d ((FLO_union_type *) & arg_a, &a);
 920   unpack_d ((FLO_union_type *) & arg_b, &b);
 921
 922   res = _fpdiv_parts (&a, &b, &tmp);
 923
 924   return pack_d (res);
 925 }
 926
 927 /* according to the demo, fpcmp returns a comparison with 0... thus
 928    a<b -> -1
 929    a==b -> 0
 930    a>b -> +1
 931  */
 932
 933 static int
 934 _fpcmp_parts (fp_number_type * a, fp_number_type * b)
 935 {
 936 #if 0
 937   /* either nan -> unordered. Must be checked outside of this routine. */
 938   if (isnan (a) && isnan (b))
 939     {
 940       return 1;                 /* still unordered! */
 941     }
 942 #endif
 943
 944   if (isnan (a) || isnan (b))
 945     {
 946       return 1;                 /* how to indicate unordered compare? */
 947     }
 948   if (isinf (a) && isinf (b))
 949     {
 950       /* +inf > -inf, but +inf != +inf */
 951       /* b    \a| +inf(0)| -inf(1)
 952        ______\+--------+--------
 953        +inf(0)| a==b(0)| a<b(-1)
 954        -------+--------+--------
 955        -inf(1)| a>b(1) | a==b(0)
 956        -------+--------+--------
 957        So since unordered must be non zero, just line up the columns...
 958        */
 959       return b->sign - a->sign;
 960     }
 961   /* but not both... */
 962   if (isinf (a))
 963     {
 964       return a->sign ? -1 : 1;
 965     }
 966   if (isinf (b))
 967     {
 968       return b->sign ? 1 : -1;
 969     }
 970   if (iszero (a) && iszero (b))
 971     {
 972       return 0;
 973     }
 974   if (iszero (a))
 975     {
 976       return b->sign ? 1 : -1;
 977     }
 978   if (iszero (b))
 979     {
 980       return a->sign ? -1 : 1;
 981     }
 982   /* now both are "normal". */
 983   if (a->sign != b->sign)
 984     {
 985       /* opposite signs */
 986       return a->sign ? -1 : 1;
 987     }
 988   /* same sign; exponents? */
 989   if (a->normal_exp > b->normal_exp)
 990     {
 991       return a->sign ? -1 : 1;
 992     }
 993   if (a->normal_exp < b->normal_exp)
 994     {
 995       return a->sign ? 1 : -1;
 996     }
 997   /* same exponents; check size. */
 998   if (a->fraction.ll > b->fraction.ll)
 999     {
1000       return a->sign ? -1 : 1;
1001     }
1002   if (a->fraction.ll < b->fraction.ll)
1003     {
1004       return a->sign ? 1 : -1;
1005     }
1006   /* after all that, they're equal. */
1007   return 0;
1008 }
1009
1010 CMPtype
1011 compare (FLO_type arg_a, FLO_type arg_b)
1012 {
1013   fp_number_type a;
1014   fp_number_type b;
1015
1016   unpack_d ((FLO_union_type *) & arg_a, &a);
1017   unpack_d ((FLO_union_type *) & arg_b, &b);
1018
1019   return _fpcmp_parts (&a, &b);
1020 }
1021
1022 #ifndef US_SOFTWARE_GOFAST
1023
1024 /* These should be optimized for their specific tasks someday.  */
1025
1026 CMPtype
1027 _eq_f2 (FLO_type arg_a, FLO_type arg_b)
1028 {
1029   fp_number_type a;
1030   fp_number_type b;
1031
1032   unpack_d ((FLO_union_type *) & arg_a, &a);
1033   unpack_d ((FLO_union_type *) & arg_b, &b);
1034
1035   if (isnan (&a) || isnan (&b))
1036     return 1;                   /* false, truth == 0 */
1037
1038   return _fpcmp_parts (&a, &b) ;
1039 }
1040
1041 CMPtype
1042 _ne_f2 (FLO_type arg_a, FLO_type arg_b)
1043 {
1044   fp_number_type a;
1045   fp_number_type b;
1046
1047   unpack_d ((FLO_union_type *) & arg_a, &a);
1048   unpack_d ((FLO_union_type *) & arg_b, &b);
1049
1050   if (isnan (&a) || isnan (&b))
1051     return 1;                   /* true, truth != 0 */
1052
1053   return  _fpcmp_parts (&a, &b) ;
1054 }
1055
1056 CMPtype
1057 _gt_f2 (FLO_type arg_a, FLO_type arg_b)
1058 {
1059   fp_number_type a;
1060   fp_number_type b;
1061
1062   unpack_d ((FLO_union_type *) & arg_a, &a);
1063   unpack_d ((FLO_union_type *) & arg_b, &b);
1064
1065   if (isnan (&a) || isnan (&b))
1066     return -1;                  /* false, truth > 0 */
1067
1068   return _fpcmp_parts (&a, &b);
1069 }
1070
1071 CMPtype
1072 _ge_f2 (FLO_type arg_a, FLO_type arg_b)
1073 {
1074   fp_number_type a;
1075   fp_number_type b;
1076
1077   unpack_d ((FLO_union_type *) & arg_a, &a);
1078   unpack_d ((FLO_union_type *) & arg_b, &b);
1079
1080   if (isnan (&a) || isnan (&b))
1081     return -1;                  /* false, truth >= 0 */
1082   return _fpcmp_parts (&a, &b) ;
1083 }
1084
1085 CMPtype
1086 _lt_f2 (FLO_type arg_a, FLO_type arg_b)
1087 {
1088   fp_number_type a;
1089   fp_number_type b;
1090
1091   unpack_d ((FLO_union_type *) & arg_a, &a);
1092   unpack_d ((FLO_union_type *) & arg_b, &b);
1093
1094   if (isnan (&a) || isnan (&b))
1095     return 1;                   /* false, truth < 0 */
1096
1097   return _fpcmp_parts (&a, &b);
1098 }
1099
1100 CMPtype
1101 _le_f2 (FLO_type arg_a, FLO_type arg_b)
1102 {
1103   fp_number_type a;
1104   fp_number_type b;
1105
1106   unpack_d ((FLO_union_type *) & arg_a, &a);
1107   unpack_d ((FLO_union_type *) & arg_b, &b);
1108
1109   if (isnan (&a) || isnan (&b))
1110     return 1;                   /* false, truth <= 0 */
1111
1112   return _fpcmp_parts (&a, &b) ;
1113 }
1114
1115 #endif /* ! US_SOFTWARE_GOFAST */
1116
1117 FLO_type
1118 si_to_float (SItype arg_a)
1119 {
1120   fp_number_type in;
1121
1122   in.class = CLASS_NUMBER;
1123   in.sign = arg_a < 0;
1124   if (!arg_a)
1125     {
1126       in.class = CLASS_ZERO;
1127     }
1128   else
1129     {
1130       in.normal_exp = FRACBITS + NGARDS;
1131       if (in.sign)
1132         {
1133           /* Special case for minint, since there is no +ve integer
1134              representation for it */
1135           if (arg_a == 0x80000000)
1136             {
1137               return -2147483648.0;
1138             }
1139           in.fraction.ll = (-arg_a);
1140         }
1141       else
1142         in.fraction.ll = arg_a;
1143
1144       while (in.fraction.ll < (1LL << (FRACBITS + NGARDS)))
1145         {
1146           in.fraction.ll <<= 1;
1147           in.normal_exp -= 1;
1148         }
1149     }
1150   return pack_d (&in);
1151 }
1152
1153 SItype
1154 float_to_si (FLO_type arg_a)
1155 {
1156   fp_number_type a;
1157   SItype tmp;
1158
1159   unpack_d ((FLO_union_type *) & arg_a, &a);
1160   if (iszero (&a))
1161     return 0;
1162   if (isnan (&a))
1163     return 0;
1164   /* get reasonable MAX_SI_INT... */
1165   if (isinf (&a))
1166     return a.sign ? MAX_SI_INT : (-MAX_SI_INT)-1;
1167   /* it is a number, but a small one */
1168   if (a.normal_exp < 0)
1169     return 0;
1170   if (a.normal_exp > 30)
1171     return a.sign ? (-MAX_SI_INT)-1 : MAX_SI_INT;
1172   tmp = a.fraction.ll >> ((FRACBITS + NGARDS) - a.normal_exp);
1173   return a.sign ? (-tmp) : (tmp);
1174 }
1175
1176 #ifdef US_SOFTWARE_GOFAST
1177 /* While libgcc2.c defines its own __fixunssfsi and __fixunsdfsi routines,
1178    we also define them for GOFAST because the ones in libgcc2.c have the
1179    wrong names and I'd rather define these here and keep GOFAST CYG-LOC's
1180    out of libgcc2.c.  We can't define these here if not GOFAST because then
1181    there'd be duplicate copies.  */
1182
1183 USItype
1184 float_to_usi (FLO_type arg_a)
1185 {
1186   fp_number_type a;
1187
1188   unpack_d ((FLO_union_type *) & arg_a, &a);
1189   if (iszero (&a))
1190     return 0;
1191   if (isnan (&a))
1192     return 0;
1193   /* get reasonable MAX_USI_INT... */
1194   if (isinf (&a))
1195     return a.sign ? MAX_USI_INT : 0;
1196   /* it is a negative number */
1197   if (a.sign)
1198     return 0;
1199   /* it is a number, but a small one */
1200   if (a.normal_exp < 0)
1201     return 0;
1202   if (a.normal_exp > 31)
1203     return MAX_USI_INT;
1204   else if (a.normal_exp > (FRACBITS + NGARDS))
1205     return a.fraction.ll << ((FRACBITS + NGARDS) - a.normal_exp);
1206   else
1207     return a.fraction.ll >> ((FRACBITS + NGARDS) - a.normal_exp);
1208 }
1209 #endif
1210
1211 FLO_type
1212 negate (FLO_type arg_a)
1213 {
1214   fp_number_type a;
1215
1216   unpack_d ((FLO_union_type *) & arg_a, &a);
1217   flip_sign (&a);
1218   return pack_d (&a);
1219 }
1220
1221 #ifdef FLOAT
1222
1223 SFtype
1224 __make_fp(fp_class_type class,
1225              unsigned int sign,
1226              int exp,
1227              USItype frac)
1228 {
1229   fp_number_type in;
1230
1231   in.class = class;
1232   in.sign = sign;
1233   in.normal_exp = exp;
1234   in.fraction.ll = frac;
1235   return pack_d (&in);
1236 }
1237
1238 #ifndef FLOAT_ONLY
1239
1240 /* This enables one to build an fp library that supports float but not double.
1241    Otherwise, we would get an undefined reference to __make_dp.
1242    This is needed for some 8-bit ports that can't handle well values that
1243    are 8-bytes in size, so we just don't support double for them at all.  */
1244
1245 extern DFtype __make_dp (fp_class_type, unsigned int, int, UDItype frac);
1246
1247 DFtype
1248 sf_to_df (SFtype arg_a)
1249 {
1250   fp_number_type in;
1251
1252   unpack_d ((FLO_union_type *) & arg_a, &in);
1253   return __make_dp (in.class, in.sign, in.normal_exp,
1254                     ((UDItype) in.fraction.ll) << F_D_BITOFF);
1255 }
1256
1257 #endif
1258 #endif
1259
1260 #ifndef FLOAT
1261
1262 extern SFtype __make_fp (fp_class_type, unsigned int, int, USItype);
1263
1264 DFtype
1265 __make_dp (fp_class_type class, unsigned int sign, int exp, UDItype frac)
1266 {
1267   fp_number_type in;
1268
1269   in.class = class;
1270   in.sign = sign;
1271   in.normal_exp = exp;
1272   in.fraction.ll = frac;
1273   return pack_d (&in);
1274 }
1275
1276 SFtype
1277 df_to_sf (DFtype arg_a)
1278 {
1279   fp_number_type in;
1280
1281   unpack_d ((FLO_union_type *) & arg_a, &in);
1282   return __make_fp (in.class, in.sign, in.normal_exp,
1283                     in.fraction.ll >> F_D_BITOFF);
1284 }
1285
1286 #endif