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 SFtype __attribute__ ((mode (SF)));
 160 typedef 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       sign = 0;
 464     }
 465   else
 466     {
 467       if (src->normal_exp < NORMAL_EXPMIN)
 468         {
 469           /* This number's exponent is too low to fit into the bits
 470              available in the number, so we'll store 0 in the exponent and
 471              shift the fraction to the right to make up for it.  */
 472
 473           int shift = NORMAL_EXPMIN - src->normal_exp;
 474
 475           exp = 0;
 476
 477           if (shift > FRAC_NBITS - NGARDS)
 478             {
 479               /* No point shifting, since it's more that 64 out.  */
 480               fraction = 0;
 481             }
 482           else
 483             {
 484               /* Shift by the value */
 485               fraction >>= shift;
 486             }
 487           fraction >>= NGARDS;
 488         }
 489       else if (src->normal_exp > EXPBIAS)
 490         {
 491           exp = EXPMAX;
 492           fraction = 0;
 493         }
 494       else
 495         {
 496           exp = src->normal_exp + EXPBIAS;
 497           /* IF the gard bits are the all zero, but the first, then we're
 498              half way between two numbers, choose the one which makes the
 499              lsb of the answer 0.  */
 500           if ((fraction & GARDMASK) == GARDMSB)
 501             {
 502               if (fraction & (1 << NGARDS))
 503                 fraction += GARDROUND + 1;
 504             }
 505           else
 506             {
 507               /* Add a one to the guards to round up */
 508               fraction += GARDROUND;
 509             }
 510           if (fraction >= IMPLICIT_2)
 511             {
 512               fraction >>= 1;
 513               exp += 1;
 514             }
 515           fraction >>= NGARDS;
 516         }
 517     }
 518
 519   /* We previously used bitfields to store the number, but this doesn't
 520      handle little/big endian systems conveniently, so use shifts and
 521      masks */
 522 #ifdef FLOAT_BIT_ORDER_MISMATCH
 523   dst.bits.fraction = fraction;
 524   dst.bits.exp = exp;
 525   dst.bits.sign = sign;
 526 #else
 527   dst.value_raw = fraction & ((((fractype)1) << FRACBITS) - (fractype)1);
 528   dst.value_raw |= ((fractype) (exp & ((1 << EXPBITS) - 1))) << FRACBITS;
 529   dst.value_raw |= ((fractype) (sign & 1)) << (FRACBITS | EXPBITS);
 530 #endif
 531
 532 #if defined(FLOAT_WORD_ORDER_MISMATCH) && !defined(FLOAT)
 533   {
 534     halffractype tmp = dst.words[0];
 535     dst.words[0] = dst.words[1];
 536     dst.words[1] = tmp;
 537   }
 538 #endif
 539
 540   return dst.value;
 541 }
 542 #endif
 543
 544 extern void unpack_d (FLO_union_type *, fp_number_type *);
 545
 546 #if defined(L_unpack_df) || defined(L_unpack_sf)
 547 void
 548 unpack_d (FLO_union_type * src, fp_number_type * dst)
 549 {
 550   /* We previously used bitfields to store the number, but this doesn't
 551      handle little/big endian systems conveniently, so use shifts and
 552      masks */
 553   fractype fraction;
 554   int exp;
 555   int sign;
 556
 557 #if defined(FLOAT_WORD_ORDER_MISMATCH) && !defined(FLOAT)
 558   FLO_union_type swapped;
 559
 560   swapped.words[0] = src->words[1];
 561   swapped.words[1] = src->words[0];
 562   src = &swapped;
 563 #endif
 564
 565 #ifdef FLOAT_BIT_ORDER_MISMATCH
 566   fraction = src->bits.fraction;
 567   exp = src->bits.exp;
 568   sign = src->bits.sign;
 569 #else
 570   fraction = src->value_raw & ((((fractype)1) << FRACBITS) - (fractype)1);
 571   exp = ((int)(src->value_raw >> FRACBITS)) & ((1 << EXPBITS) - 1);
 572   sign = ((int)(src->value_raw >> (FRACBITS + EXPBITS))) & 1;
 573 #endif
 574
 575   dst->sign = sign;
 576   if (exp == 0)
 577     {
 578       /* Hmm.  Looks like 0 */
 579       if (fraction == 0)
 580         {
 581           /* tastes like zero */
 582           dst->class = CLASS_ZERO;
 583         }
 584       else
 585         {
 586           /* Zero exponent with non zero fraction - it's denormalized,
 587              so there isn't a leading implicit one - we'll shift it so
 588              it gets one.  */
 589           dst->normal_exp = exp - EXPBIAS + 1;
 590           fraction <<= NGARDS;
 591
 592           dst->class = CLASS_NUMBER;
 593 #if 1
 594           while (fraction < IMPLICIT_1)
 595             {
 596               fraction <<= 1;
 597               dst->normal_exp--;
 598             }
 599 #endif
 600           dst->fraction.ll = fraction;
 601         }
 602     }
 603   else if (exp == EXPMAX)
 604     {
 605       /* Huge exponent*/
 606       if (fraction == 0)
 607         {
 608           /* Attached to a zero fraction - means infinity */
 609           dst->class = CLASS_INFINITY;
 610         }
 611       else
 612         {
 613           /* Non zero fraction, means nan */
 614           if (fraction & QUIET_NAN)
 615             {
 616               dst->class = CLASS_QNAN;
 617             }
 618           else
 619             {
 620               dst->class = CLASS_SNAN;
 621             }
 622           /* Keep the fraction part as the nan number */
 623           dst->fraction.ll = fraction;
 624         }
 625     }
 626   else
 627     {
 628       /* Nothing strange about this number */
 629       dst->normal_exp = exp - EXPBIAS;
 630       dst->class = CLASS_NUMBER;
 631       dst->fraction.ll = (fraction << NGARDS) | IMPLICIT_1;
 632     }
 633 }
 634 #endif
 635
 636 #if defined(L_addsub_sf) || defined(L_addsub_df)
 637 static fp_number_type *
 638 _fpadd_parts (fp_number_type * a,
 639               fp_number_type * b,
 640               fp_number_type * tmp)
 641 {
 642   intfrac tfraction;
 643
 644   /* Put commonly used fields in local variables.  */
 645   int a_normal_exp;
 646   int b_normal_exp;
 647   fractype a_fraction;
 648   fractype b_fraction;
 649
 650   if (isnan (a))
 651     {
 652       return a;
 653     }
 654   if (isnan (b))
 655     {
 656       return b;
 657     }
 658   if (isinf (a))
 659     {
 660       /* Adding infinities with opposite signs yields a NaN.  */
 661       if (isinf (b) && a->sign != b->sign)
 662         return nan ();
 663       return a;
 664     }
 665   if (isinf (b))
 666     {
 667       return b;
 668     }
 669   if (iszero (b))
 670     {
 671       if (iszero (a))
 672         {
 673           *tmp = *a;
 674           tmp->sign = a->sign & b->sign;
 675           return tmp;
 676         }
 677       return a;
 678     }
 679   if (iszero (a))
 680     {
 681       return b;
 682     }
 683
 684   /* Got two numbers. shift the smaller and increment the exponent till
 685      they're the same */
 686   {
 687     int diff;
 688
 689     a_normal_exp = a->normal_exp;
 690     b_normal_exp = b->normal_exp;
 691     a_fraction = a->fraction.ll;
 692     b_fraction = b->fraction.ll;
 693
 694     diff = a_normal_exp - b_normal_exp;
 695
 696     if (diff < 0)
 697       diff = -diff;
 698     if (diff < FRAC_NBITS)
 699       {
 700         /* ??? This does shifts one bit at a time.  Optimize.  */
 701         while (a_normal_exp > b_normal_exp)
 702           {
 703             b_normal_exp++;
 704             LSHIFT (b_fraction);
 705           }
 706         while (b_normal_exp > a_normal_exp)
 707           {
 708             a_normal_exp++;
 709             LSHIFT (a_fraction);
 710           }
 711       }
 712     else
 713       {
 714         /* Somethings's up.. choose the biggest */
 715         if (a_normal_exp > b_normal_exp)
 716           {
 717             b_normal_exp = a_normal_exp;
 718             b_fraction = 0;
 719           }
 720         else
 721           {
 722             a_normal_exp = b_normal_exp;
 723             a_fraction = 0;
 724           }
 725       }
 726   }
 727
 728   if (a->sign != b->sign)
 729     {
 730       if (a->sign)
 731         {
 732           tfraction = -a_fraction + b_fraction;
 733         }
 734       else
 735         {
 736           tfraction = a_fraction - b_fraction;
 737         }
 738       if (tfraction > 0)
 739         {
 740           tmp->sign = 0;
 741           tmp->normal_exp = a_normal_exp;
 742           tmp->fraction.ll = tfraction;
 743         }
 744       else
 745         {
 746           tmp->sign = 1;
 747           tmp->normal_exp = a_normal_exp;
 748           tmp->fraction.ll = -tfraction;
 749         }
 750       /* and renormalize it */
 751
 752       while (tmp->fraction.ll < IMPLICIT_1 && tmp->fraction.ll)
 753         {
 754           tmp->fraction.ll <<= 1;
 755           tmp->normal_exp--;
 756         }
 757     }
 758   else
 759     {
 760       tmp->sign = a->sign;
 761       tmp->normal_exp = a_normal_exp;
 762       tmp->fraction.ll = a_fraction + b_fraction;
 763     }
 764   tmp->class = CLASS_NUMBER;
 765   /* Now the fraction is added, we have to shift down to renormalize the
 766      number */
 767
 768   if (tmp->fraction.ll >= IMPLICIT_2)
 769     {
 770       LSHIFT (tmp->fraction.ll);
 771       tmp->normal_exp++;
 772     }
 773   return tmp;
 774
 775 }
 776
 777 FLO_type
 778 add (FLO_type arg_a, FLO_type arg_b)
 779 {
 780   fp_number_type a;
 781   fp_number_type b;
 782   fp_number_type tmp;
 783   fp_number_type *res;
 784
 785   unpack_d ((FLO_union_type *) & arg_a, &a);
 786   unpack_d ((FLO_union_type *) & arg_b, &b);
 787
 788   res = _fpadd_parts (&a, &b, &tmp);
 789
 790   return pack_d (res);
 791 }
 792
 793 FLO_type
 794 sub (FLO_type arg_a, FLO_type arg_b)
 795 {
 796   fp_number_type a;
 797   fp_number_type b;
 798   fp_number_type tmp;
 799   fp_number_type *res;
 800
 801   unpack_d ((FLO_union_type *) & arg_a, &a);
 802   unpack_d ((FLO_union_type *) & arg_b, &b);
 803
 804   b.sign ^= 1;
 805
 806   res = _fpadd_parts (&a, &b, &tmp);
 807
 808   return pack_d (res);
 809 }
 810 #endif
 811
 812 #if defined(L_mul_sf) || defined(L_mul_df)
 813 static INLINE fp_number_type *
 814 _fpmul_parts ( fp_number_type *  a,
 815                fp_number_type *  b,
 816                fp_number_type * tmp)
 817 {
 818   fractype low = 0;
 819   fractype high = 0;
 820
 821   if (isnan (a))
 822     {
 823       a->sign = a->sign != b->sign;
 824       return a;
 825     }
 826   if (isnan (b))
 827     {
 828       b->sign = a->sign != b->sign;
 829       return b;
 830     }
 831   if (isinf (a))
 832     {
 833       if (iszero (b))
 834         return nan ();
 835       a->sign = a->sign != b->sign;
 836       return a;
 837     }
 838   if (isinf (b))
 839     {
 840       if (iszero (a))
 841         {
 842           return nan ();
 843         }
 844       b->sign = a->sign != b->sign;
 845       return b;
 846     }
 847   if (iszero (a))
 848     {
 849       a->sign = a->sign != b->sign;
 850       return a;
 851     }
 852   if (iszero (b))
 853     {
 854       b->sign = a->sign != b->sign;
 855       return b;
 856     }
 857
 858   /* Calculate the mantissa by multiplying both 64bit numbers to get a
 859      128 bit number */
 860   {
 861 #if defined(NO_DI_MODE)
 862     {
 863       fractype x = a->fraction.ll;
 864       fractype ylow = b->fraction.ll;
 865       fractype yhigh = 0;
 866       int bit;
 867
 868       /* ??? This does multiplies one bit at a time.  Optimize.  */
 869       for (bit = 0; bit < FRAC_NBITS; bit++)
 870         {
 871           int carry;
 872
 873           if (x & 1)
 874             {
 875               carry = (low += ylow) < ylow;
 876               high += yhigh + carry;
 877             }
 878           yhigh <<= 1;
 879           if (ylow & FRACHIGH)
 880             {
 881               yhigh |= 1;
 882             }
 883           ylow <<= 1;
 884           x >>= 1;
 885         }
 886     }
 887 #elif defined(FLOAT)
 888     {
 889       /* Multiplying two 32 bit numbers to get a 64 bit number  on
 890         a machine with DI, so we're safe */
 891
 892       DItype answer = (DItype)(a->fraction.ll) * (DItype)(b->fraction.ll);
 893
 894       high = answer >> 32;
 895       low = answer;
 896     }
 897 #else
 898     /* Doing a 64*64 to 128 */
 899     {
 900       UDItype nl = a->fraction.ll & 0xffffffff;
 901       UDItype nh = a->fraction.ll >> 32;
 902       UDItype ml = b->fraction.ll & 0xffffffff;
 903       UDItype mh = b->fraction.ll >>32;
 904       UDItype pp_ll = ml * nl;
 905       UDItype pp_hl = mh * nl;
 906       UDItype pp_lh = ml * nh;
 907       UDItype pp_hh = mh * nh;
 908       UDItype res2 = 0;
 909       UDItype res0 = 0;
 910       UDItype ps_hh__ = pp_hl + pp_lh;
 911       if (ps_hh__ < pp_hl)
 912         res2 += 0x100000000LL;
 913       pp_hl = (ps_hh__ << 32) & 0xffffffff00000000LL;
 914       res0 = pp_ll + pp_hl;
 915       if (res0 < pp_ll)
 916         res2++;
 917       res2 += ((ps_hh__ >> 32) & 0xffffffffL) + pp_hh;
 918       high = res2;
 919       low = res0;
 920     }
 921 #endif
 922   }
 923
 924   tmp->normal_exp = a->normal_exp + b->normal_exp;
 925   tmp->sign = a->sign != b->sign;
 926 #ifdef FLOAT
 927   tmp->normal_exp += 2;         /* ??????????????? */
 928 #else
 929   tmp->normal_exp += 4;         /* ??????????????? */
 930 #endif
 931   while (high >= IMPLICIT_2)
 932     {
 933       tmp->normal_exp++;
 934       if (high & 1)
 935         {
 936           low >>= 1;
 937           low |= FRACHIGH;
 938         }
 939       high >>= 1;
 940     }
 941   while (high < IMPLICIT_1)
 942     {
 943       tmp->normal_exp--;
 944
 945       high <<= 1;
 946       if (low & FRACHIGH)
 947         high |= 1;
 948       low <<= 1;
 949     }
 950   /* rounding is tricky. if we only round if it won't make us round later. */
 951 #if 0
 952   if (low & FRACHIGH2)
 953     {
 954       if (((high & GARDMASK) != GARDMSB)
 955           && (((high + 1) & GARDMASK) == GARDMSB))
 956         {
 957           /* don't round, it gets done again later. */
 958         }
 959       else
 960         {
 961           high++;
 962         }
 963     }
 964 #endif
 965   if ((high & GARDMASK) == GARDMSB)
 966     {
 967       if (high & (1 << NGARDS))
 968         {
 969           /* half way, so round to even */
 970           high += GARDROUND + 1;
 971         }
 972       else if (low)
 973         {
 974           /* but we really weren't half way */
 975           high += GARDROUND + 1;
 976         }
 977     }
 978   tmp->fraction.ll = high;
 979   tmp->class = CLASS_NUMBER;
 980   return tmp;
 981 }
 982
 983 FLO_type
 984 multiply (FLO_type arg_a, FLO_type arg_b)
 985 {
 986   fp_number_type a;
 987   fp_number_type b;
 988   fp_number_type tmp;
 989   fp_number_type *res;
 990
 991   unpack_d ((FLO_union_type *) & arg_a, &a);
 992   unpack_d ((FLO_union_type *) & arg_b, &b);
 993
 994   res = _fpmul_parts (&a, &b, &tmp);
 995
 996   return pack_d (res);
 997 }
 998 #endif
 999
1000 #if defined(L_div_sf) || defined(L_div_df)
1001 static INLINE fp_number_type *
1002 _fpdiv_parts (fp_number_type * a,
1003               fp_number_type * b,
1004               fp_number_type * tmp)
1005 {
1006   fractype bit;
1007   fractype numerator;
1008   fractype denominator;
1009   fractype quotient;
1010
1011   if (isnan (a))
1012     {
1013       return a;
1014     }
1015   if (isnan (b))
1016     {
1017       return b;
1018     }
1019
1020   a->sign = a->sign ^ b->sign;
1021
1022   if (isinf (a) || iszero (a))
1023     {
1024       if (a->class == b->class)
1025         return nan ();
1026       return a;
1027     }
1028
1029   if (isinf (b))
1030     {
1031       a->fraction.ll = 0;
1032       a->normal_exp = 0;
1033       return a;
1034     }
1035   if (iszero (b))
1036     {
1037       a->class = CLASS_INFINITY;
1038       return a;
1039     }
1040
1041   /* Calculate the mantissa by multiplying both 64bit numbers to get a
1042      128 bit number */
1043   {
1044     /* quotient =
1045        ( numerator / denominator) * 2^(numerator exponent -  denominator exponent)
1046      */
1047
1048     a->normal_exp = a->normal_exp - b->normal_exp;
1049     numerator = a->fraction.ll;
1050     denominator = b->fraction.ll;
1051
1052     if (numerator < denominator)
1053       {
1054         /* Fraction will be less than 1.0 */
1055         numerator *= 2;
1056         a->normal_exp--;
1057       }
1058     bit = IMPLICIT_1;
1059     quotient = 0;
1060     /* ??? Does divide one bit at a time.  Optimize.  */
1061     while (bit)
1062       {
1063         if (numerator >= denominator)
1064           {
1065             quotient |= bit;
1066             numerator -= denominator;
1067           }
1068         bit >>= 1;
1069         numerator *= 2;
1070       }
1071
1072     if ((quotient & GARDMASK) == GARDMSB)
1073       {
1074         if (quotient & (1 << NGARDS))
1075           {
1076             /* half way, so round to even */
1077             quotient += GARDROUND + 1;
1078           }
1079         else if (numerator)
1080           {
1081             /* but we really weren't half way, more bits exist */
1082             quotient += GARDROUND + 1;
1083           }
1084       }
1085
1086     a->fraction.ll = quotient;
1087     return (a);
1088   }
1089 }
1090
1091 FLO_type
1092 divide (FLO_type arg_a, FLO_type arg_b)
1093 {
1094   fp_number_type a;
1095   fp_number_type b;
1096   fp_number_type tmp;
1097   fp_number_type *res;
1098
1099   unpack_d ((FLO_union_type *) & arg_a, &a);
1100   unpack_d ((FLO_union_type *) & arg_b, &b);
1101
1102   res = _fpdiv_parts (&a, &b, &tmp);
1103
1104   return pack_d (res);
1105 }
1106 #endif
1107
1108 int __fpcmp_parts (fp_number_type * a, fp_number_type *b);
1109
1110 #if defined(L_fpcmp_parts_sf) || defined(L_fpcmp_parts_df)
1111 /* according to the demo, fpcmp returns a comparison with 0... thus
1112    a<b -> -1
1113    a==b -> 0
1114    a>b -> +1
1115  */
1116
1117 int
1118 __fpcmp_parts (fp_number_type * a, fp_number_type * b)
1119 {
1120 #if 0
1121   /* either nan -> unordered. Must be checked outside of this routine. */
1122   if (isnan (a) && isnan (b))
1123     {
1124       return 1;                 /* still unordered! */
1125     }
1126 #endif
1127
1128   if (isnan (a) || isnan (b))
1129     {
1130       return 1;                 /* how to indicate unordered compare? */
1131     }
1132   if (isinf (a) && isinf (b))
1133     {
1134       /* +inf > -inf, but +inf != +inf */
1135       /* b    \a| +inf(0)| -inf(1)
1136        ______\+--------+--------
1137        +inf(0)| a==b(0)| a<b(-1)
1138        -------+--------+--------
1139        -inf(1)| a>b(1) | a==b(0)
1140        -------+--------+--------
1141        So since unordered must be non zero, just line up the columns...
1142        */
1143       return b->sign - a->sign;
1144     }
1145   /* but not both... */
1146   if (isinf (a))
1147     {
1148       return a->sign ? -1 : 1;
1149     }
1150   if (isinf (b))
1151     {
1152       return b->sign ? 1 : -1;
1153     }
1154   if (iszero (a) && iszero (b))
1155     {
1156       return 0;
1157     }
1158   if (iszero (a))
1159     {
1160       return b->sign ? 1 : -1;
1161     }
1162   if (iszero (b))
1163     {
1164       return a->sign ? -1 : 1;
1165     }
1166   /* now both are "normal". */
1167   if (a->sign != b->sign)
1168     {
1169       /* opposite signs */
1170       return a->sign ? -1 : 1;
1171     }
1172   /* same sign; exponents? */
1173   if (a->normal_exp > b->normal_exp)
1174     {
1175       return a->sign ? -1 : 1;
1176     }
1177   if (a->normal_exp < b->normal_exp)
1178     {
1179       return a->sign ? 1 : -1;
1180     }
1181   /* same exponents; check size. */
1182   if (a->fraction.ll > b->fraction.ll)
1183     {
1184       return a->sign ? -1 : 1;
1185     }
1186   if (a->fraction.ll < b->fraction.ll)
1187     {
1188       return a->sign ? 1 : -1;
1189     }
1190   /* after all that, they're equal. */
1191   return 0;
1192 }
1193 #endif
1194
1195 #if defined(L_compare_sf) || defined(L_compare_df)
1196 CMPtype
1197 compare (FLO_type arg_a, FLO_type arg_b)
1198 {
1199   fp_number_type a;
1200   fp_number_type b;
1201
1202   unpack_d ((FLO_union_type *) & arg_a, &a);
1203   unpack_d ((FLO_union_type *) & arg_b, &b);
1204
1205   return __fpcmp_parts (&a, &b);
1206 }
1207 #endif
1208
1209 #ifndef US_SOFTWARE_GOFAST
1210
1211 /* These should be optimized for their specific tasks someday.  */
1212
1213 #if defined(L_eq_sf) || defined(L_eq_df)
1214 CMPtype
1215 _eq_f2 (FLO_type arg_a, FLO_type arg_b)
1216 {
1217   fp_number_type a;
1218   fp_number_type b;
1219
1220   unpack_d ((FLO_union_type *) & arg_a, &a);
1221   unpack_d ((FLO_union_type *) & arg_b, &b);
1222
1223   if (isnan (&a) || isnan (&b))
1224     return 1;                   /* false, truth == 0 */
1225
1226   return __fpcmp_parts (&a, &b) ;
1227 }
1228 #endif
1229
1230 #if defined(L_ne_sf) || defined(L_ne_df)
1231 CMPtype
1232 _ne_f2 (FLO_type arg_a, FLO_type arg_b)
1233 {
1234   fp_number_type a;
1235   fp_number_type b;
1236
1237   unpack_d ((FLO_union_type *) & arg_a, &a);
1238   unpack_d ((FLO_union_type *) & arg_b, &b);
1239
1240   if (isnan (&a) || isnan (&b))
1241     return 1;                   /* true, truth != 0 */
1242
1243   return  __fpcmp_parts (&a, &b) ;
1244 }
1245 #endif
1246
1247 #if defined(L_gt_sf) || defined(L_gt_df)
1248 CMPtype
1249 _gt_f2 (FLO_type arg_a, FLO_type arg_b)
1250 {
1251   fp_number_type a;
1252   fp_number_type b;
1253
1254   unpack_d ((FLO_union_type *) & arg_a, &a);
1255   unpack_d ((FLO_union_type *) & arg_b, &b);
1256
1257   if (isnan (&a) || isnan (&b))
1258     return -1;                  /* false, truth > 0 */
1259
1260   return __fpcmp_parts (&a, &b);
1261 }
1262 #endif
1263
1264 #if defined(L_ge_sf) || defined(L_ge_df)
1265 CMPtype
1266 _ge_f2 (FLO_type arg_a, FLO_type arg_b)
1267 {
1268   fp_number_type a;
1269   fp_number_type b;
1270
1271   unpack_d ((FLO_union_type *) & arg_a, &a);
1272   unpack_d ((FLO_union_type *) & arg_b, &b);
1273
1274   if (isnan (&a) || isnan (&b))
1275     return -1;                  /* false, truth >= 0 */
1276   return __fpcmp_parts (&a, &b) ;
1277 }
1278 #endif
1279
1280 #if defined(L_lt_sf) || defined(L_lt_df)
1281 CMPtype
1282 _lt_f2 (FLO_type arg_a, FLO_type arg_b)
1283 {
1284   fp_number_type a;
1285   fp_number_type b;
1286
1287   unpack_d ((FLO_union_type *) & arg_a, &a);
1288   unpack_d ((FLO_union_type *) & arg_b, &b);
1289
1290   if (isnan (&a) || isnan (&b))
1291     return 1;                   /* false, truth < 0 */
1292
1293   return __fpcmp_parts (&a, &b);
1294 }
1295 #endif
1296
1297 #if defined(L_le_sf) || defined(L_le_df)
1298 CMPtype
1299 _le_f2 (FLO_type arg_a, FLO_type arg_b)
1300 {
1301   fp_number_type a;
1302   fp_number_type b;
1303
1304   unpack_d ((FLO_union_type *) & arg_a, &a);
1305   unpack_d ((FLO_union_type *) & arg_b, &b);
1306
1307   if (isnan (&a) || isnan (&b))
1308     return 1;                   /* false, truth <= 0 */
1309
1310   return __fpcmp_parts (&a, &b) ;
1311 }
1312 #endif
1313
1314 #endif /* ! US_SOFTWARE_GOFAST */
1315
1316 #if defined(L_si_to_sf) || defined(L_si_to_df)
1317 FLO_type
1318 si_to_float (SItype arg_a)
1319 {
1320   fp_number_type in;
1321
1322   in.class = CLASS_NUMBER;
1323   in.sign = arg_a < 0;
1324   if (!arg_a)
1325     {
1326       in.class = CLASS_ZERO;
1327     }
1328   else
1329     {
1330       in.normal_exp = FRACBITS + NGARDS;
1331       if (in.sign)
1332         {
1333           /* Special case for minint, since there is no +ve integer
1334              representation for it */
1335           if (arg_a == 0x80000000)
1336             {
1337               return -2147483648.0;
1338             }
1339           in.fraction.ll = (-arg_a);
1340         }
1341       else
1342         in.fraction.ll = arg_a;
1343
1344       while (in.fraction.ll < (1LL << (FRACBITS + NGARDS)))
1345         {
1346           in.fraction.ll <<= 1;
1347           in.normal_exp -= 1;
1348         }
1349     }
1350   return pack_d (&in);
1351 }
1352 #endif
1353
1354 #if defined(L_sf_to_si) || defined(L_df_to_si)
1355 SItype
1356 float_to_si (FLO_type arg_a)
1357 {
1358   fp_number_type a;
1359   SItype tmp;
1360
1361   unpack_d ((FLO_union_type *) & arg_a, &a);
1362   if (iszero (&a))
1363     return 0;
1364   if (isnan (&a))
1365     return 0;
1366   /* get reasonable MAX_SI_INT... */
1367   if (isinf (&a))
1368     return a.sign ? (-MAX_SI_INT)-1 : MAX_SI_INT;
1369   /* it is a number, but a small one */
1370   if (a.normal_exp < 0)
1371     return 0;
1372   if (a.normal_exp > 30)
1373     return a.sign ? (-MAX_SI_INT)-1 : MAX_SI_INT;
1374   tmp = a.fraction.ll >> ((FRACBITS + NGARDS) - a.normal_exp);
1375   return a.sign ? (-tmp) : (tmp);
1376 }
1377 #endif
1378
1379 #if defined(L_sf_to_usi) || defined(L_df_to_usi)
1380 #ifdef US_SOFTWARE_GOFAST
1381 /* While libgcc2.c defines its own __fixunssfsi and __fixunsdfsi routines,
1382    we also define them for GOFAST because the ones in libgcc2.c have the
1383    wrong names and I'd rather define these here and keep GOFAST CYG-LOC's
1384    out of libgcc2.c.  We can't define these here if not GOFAST because then
1385    there'd be duplicate copies.  */
1386
1387 USItype
1388 float_to_usi (FLO_type arg_a)
1389 {
1390   fp_number_type a;
1391
1392   unpack_d ((FLO_union_type *) & arg_a, &a);
1393   if (iszero (&a))
1394     return 0;
1395   if (isnan (&a))
1396     return 0;
1397   /* it is a negative number */
1398   if (a.sign)
1399     return 0;
1400   /* get reasonable MAX_USI_INT... */
1401   if (isinf (&a))
1402     return MAX_USI_INT;
1403   /* it is a number, but a small one */
1404   if (a.normal_exp < 0)
1405     return 0;
1406   if (a.normal_exp > 31)
1407     return MAX_USI_INT;
1408   else if (a.normal_exp > (FRACBITS + NGARDS))
1409     return a.fraction.ll << (a.normal_exp - (FRACBITS + NGARDS));
1410   else
1411     return a.fraction.ll >> ((FRACBITS + NGARDS) - a.normal_exp);
1412 }
1413 #endif
1414 #endif
1415
1416 #if defined(L_negate_sf) || defined(L_negate_df)
1417 FLO_type
1418 negate (FLO_type arg_a)
1419 {
1420   fp_number_type a;
1421
1422   unpack_d ((FLO_union_type *) & arg_a, &a);
1423   flip_sign (&a);
1424   return pack_d (&a);
1425 }
1426 #endif
1427
1428 #ifdef FLOAT
1429
1430 #if defined(L_make_sf)
1431 SFtype
1432 __make_fp(fp_class_type class,
1433              unsigned int sign,
1434              int exp,
1435              USItype frac)
1436 {
1437   fp_number_type in;
1438
1439   in.class = class;
1440   in.sign = sign;
1441   in.normal_exp = exp;
1442   in.fraction.ll = frac;
1443   return pack_d (&in);
1444 }
1445 #endif
1446
1447 #ifndef FLOAT_ONLY
1448
1449 /* This enables one to build an fp library that supports float but not double.
1450    Otherwise, we would get an undefined reference to __make_dp.
1451    This is needed for some 8-bit ports that can't handle well values that
1452    are 8-bytes in size, so we just don't support double for them at all.  */
1453
1454 extern DFtype __make_dp (fp_class_type, unsigned int, int, UDItype frac);
1455
1456 #if defined(L_sf_to_df)
1457 DFtype
1458 sf_to_df (SFtype arg_a)
1459 {
1460   fp_number_type in;
1461
1462   unpack_d ((FLO_union_type *) & arg_a, &in);
1463   return __make_dp (in.class, in.sign, in.normal_exp,
1464                     ((UDItype) in.fraction.ll) << F_D_BITOFF);
1465 }
1466 #endif
1467
1468 #endif
1469 #endif
1470
1471 #ifndef FLOAT
1472
1473 extern SFtype __make_fp (fp_class_type, unsigned int, int, USItype);
1474
1475 #if defined(L_make_df)
1476 DFtype
1477 __make_dp (fp_class_type class, unsigned int sign, int exp, UDItype frac)
1478 {
1479   fp_number_type in;
1480
1481   in.class = class;
1482   in.sign = sign;
1483   in.normal_exp = exp;
1484   in.fraction.ll = frac;
1485   return pack_d (&in);
1486 }
1487 #endif
1488
1489 #if defined(L_df_to_sf)
1490 SFtype
1491 df_to_sf (DFtype arg_a)
1492 {
1493   fp_number_type in;
1494   USItype sffrac;
1495
1496   unpack_d ((FLO_union_type *) & arg_a, &in);
1497
1498   sffrac = in.fraction.ll >> F_D_BITOFF;
1499
1500   /* We set the lowest guard bit in SFFRAC if we discarded any non
1501      zero bits.  */
1502   if ((in.fraction.ll & (((USItype) 1 << F_D_BITOFF) - 1)) != 0)
1503     sffrac |= 1;
1504
1505   return __make_fp (in.class, in.sign, in.normal_exp, sffrac);
1506 }
1507 #endif
1508
1509 #endif
1510 #endif /* !EXTENDED_FLOAT_STUBS */