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