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