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