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