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