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     }
 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               fp_number_type * tmp)
1004 {
1005   fractype bit;
1006   fractype numerator;
1007   fractype denominator;
1008   fractype quotient;
1009
1010   if (isnan (a))
1011     {
1012       return a;
1013     }
1014   if (isnan (b))
1015     {
1016       return b;
1017     }
1018
1019   a->sign = a->sign ^ b->sign;
1020
1021   if (isinf (a) || iszero (a))
1022     {
1023       if (a->class == b->class)
1024         return nan ();
1025       return a;
1026     }
1027
1028   if (isinf (b))
1029     {
1030       a->fraction.ll = 0;
1031       a->normal_exp = 0;
1032       return a;
1033     }
1034   if (iszero (b))
1035     {
1036       a->class = CLASS_INFINITY;
1037       return a;
1038     }
1039
1040   /* Calculate the mantissa by multiplying both 64bit numbers to get a
1041      128 bit number */
1042   {
1043     /* quotient =
1044        ( numerator / denominator) * 2^(numerator exponent -  denominator exponent)
1045      */
1046
1047     a->normal_exp = a->normal_exp - b->normal_exp;
1048     numerator = a->fraction.ll;
1049     denominator = b->fraction.ll;
1050
1051     if (numerator < denominator)
1052       {
1053         /* Fraction will be less than 1.0 */
1054         numerator *= 2;
1055         a->normal_exp--;
1056       }
1057     bit = IMPLICIT_1;
1058     quotient = 0;
1059     /* ??? Does divide one bit at a time.  Optimize.  */
1060     while (bit)
1061       {
1062         if (numerator >= denominator)
1063           {
1064             quotient |= bit;
1065             numerator -= denominator;
1066           }
1067         bit >>= 1;
1068         numerator *= 2;
1069       }
1070
1071     if ((quotient & GARDMASK) == GARDMSB)
1072       {
1073         if (quotient & (1 << NGARDS))
1074           {
1075             /* half way, so round to even */
1076             quotient += GARDROUND + 1;
1077           }
1078         else if (numerator)
1079           {
1080             /* but we really weren't half way, more bits exist */
1081             quotient += GARDROUND + 1;
1082           }
1083       }
1084
1085     a->fraction.ll = quotient;
1086     return (a);
1087   }
1088 }
1089
1090 FLO_type
1091 divide (FLO_type arg_a, FLO_type arg_b)
1092 {
1093   fp_number_type a;
1094   fp_number_type b;
1095   fp_number_type tmp;
1096   fp_number_type *res;
1097
1098   unpack_d ((FLO_union_type *) & arg_a, &a);
1099   unpack_d ((FLO_union_type *) & arg_b, &b);
1100
1101   res = _fpdiv_parts (&a, &b, &tmp);
1102
1103   return pack_d (res);
1104 }
1105 #endif
1106
1107 int __fpcmp_parts (fp_number_type * a, fp_number_type *b);
1108
1109 #if defined(L_fpcmp_parts_sf) || defined(L_fpcmp_parts_df)
1110 /* according to the demo, fpcmp returns a comparison with 0... thus
1111    a<b -> -1
1112    a==b -> 0
1113    a>b -> +1
1114  */
1115
1116 int
1117 __fpcmp_parts (fp_number_type * a, fp_number_type * b)
1118 {
1119 #if 0
1120   /* either nan -> unordered. Must be checked outside of this routine. */
1121   if (isnan (a) && isnan (b))
1122     {
1123       return 1;                 /* still unordered! */
1124     }
1125 #endif
1126
1127   if (isnan (a) || isnan (b))
1128     {
1129       return 1;                 /* how to indicate unordered compare? */
1130     }
1131   if (isinf (a) && isinf (b))
1132     {
1133       /* +inf > -inf, but +inf != +inf */
1134       /* b    \a| +inf(0)| -inf(1)
1135        ______\+--------+--------
1136        +inf(0)| a==b(0)| a<b(-1)
1137        -------+--------+--------
1138        -inf(1)| a>b(1) | a==b(0)
1139        -------+--------+--------
1140        So since unordered must be non zero, just line up the columns...
1141        */
1142       return b->sign - a->sign;
1143     }
1144   /* but not both... */
1145   if (isinf (a))
1146     {
1147       return a->sign ? -1 : 1;
1148     }
1149   if (isinf (b))
1150     {
1151       return b->sign ? 1 : -1;
1152     }
1153   if (iszero (a) && iszero (b))
1154     {
1155       return 0;
1156     }
1157   if (iszero (a))
1158     {
1159       return b->sign ? 1 : -1;
1160     }
1161   if (iszero (b))
1162     {
1163       return a->sign ? -1 : 1;
1164     }
1165   /* now both are "normal". */
1166   if (a->sign != b->sign)
1167     {
1168       /* opposite signs */
1169       return a->sign ? -1 : 1;
1170     }
1171   /* same sign; exponents? */
1172   if (a->normal_exp > b->normal_exp)
1173     {
1174       return a->sign ? -1 : 1;
1175     }
1176   if (a->normal_exp < b->normal_exp)
1177     {
1178       return a->sign ? 1 : -1;
1179     }
1180   /* same exponents; check size. */
1181   if (a->fraction.ll > b->fraction.ll)
1182     {
1183       return a->sign ? -1 : 1;
1184     }
1185   if (a->fraction.ll < b->fraction.ll)
1186     {
1187       return a->sign ? 1 : -1;
1188     }
1189   /* after all that, they're equal. */
1190   return 0;
1191 }
1192 #endif
1193
1194 #if defined(L_compare_sf) || defined(L_compare_df)
1195 CMPtype
1196 compare (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   return __fpcmp_parts (&a, &b);
1205 }
1206 #endif
1207
1208 #ifndef US_SOFTWARE_GOFAST
1209
1210 /* These should be optimized for their specific tasks someday.  */
1211
1212 #if defined(L_eq_sf) || defined(L_eq_df)
1213 CMPtype
1214 _eq_f2 (FLO_type arg_a, FLO_type arg_b)
1215 {
1216   fp_number_type a;
1217   fp_number_type b;
1218
1219   unpack_d ((FLO_union_type *) & arg_a, &a);
1220   unpack_d ((FLO_union_type *) & arg_b, &b);
1221
1222   if (isnan (&a) || isnan (&b))
1223     return 1;                   /* false, truth == 0 */
1224
1225   return __fpcmp_parts (&a, &b) ;
1226 }
1227 #endif
1228
1229 #if defined(L_ne_sf) || defined(L_ne_df)
1230 CMPtype
1231 _ne_f2 (FLO_type arg_a, FLO_type arg_b)
1232 {
1233   fp_number_type a;
1234   fp_number_type b;
1235
1236   unpack_d ((FLO_union_type *) & arg_a, &a);
1237   unpack_d ((FLO_union_type *) & arg_b, &b);
1238
1239   if (isnan (&a) || isnan (&b))
1240     return 1;                   /* true, truth != 0 */
1241
1242   return  __fpcmp_parts (&a, &b) ;
1243 }
1244 #endif
1245
1246 #if defined(L_gt_sf) || defined(L_gt_df)
1247 CMPtype
1248 _gt_f2 (FLO_type arg_a, FLO_type arg_b)
1249 {
1250   fp_number_type a;
1251   fp_number_type b;
1252
1253   unpack_d ((FLO_union_type *) & arg_a, &a);
1254   unpack_d ((FLO_union_type *) & arg_b, &b);
1255
1256   if (isnan (&a) || isnan (&b))
1257     return -1;                  /* false, truth > 0 */
1258
1259   return __fpcmp_parts (&a, &b);
1260 }
1261 #endif
1262
1263 #if defined(L_ge_sf) || defined(L_ge_df)
1264 CMPtype
1265 _ge_f2 (FLO_type arg_a, FLO_type arg_b)
1266 {
1267   fp_number_type a;
1268   fp_number_type b;
1269
1270   unpack_d ((FLO_union_type *) & arg_a, &a);
1271   unpack_d ((FLO_union_type *) & arg_b, &b);
1272
1273   if (isnan (&a) || isnan (&b))
1274     return -1;                  /* false, truth >= 0 */
1275   return __fpcmp_parts (&a, &b) ;
1276 }
1277 #endif
1278
1279 #if defined(L_lt_sf) || defined(L_lt_df)
1280 CMPtype
1281 _lt_f2 (FLO_type arg_a, FLO_type arg_b)
1282 {
1283   fp_number_type a;
1284   fp_number_type b;
1285
1286   unpack_d ((FLO_union_type *) & arg_a, &a);
1287   unpack_d ((FLO_union_type *) & arg_b, &b);
1288
1289   if (isnan (&a) || isnan (&b))
1290     return 1;                   /* false, truth < 0 */
1291
1292   return __fpcmp_parts (&a, &b);
1293 }
1294 #endif
1295
1296 #if defined(L_le_sf) || defined(L_le_df)
1297 CMPtype
1298 _le_f2 (FLO_type arg_a, FLO_type arg_b)
1299 {
1300   fp_number_type a;
1301   fp_number_type b;
1302
1303   unpack_d ((FLO_union_type *) & arg_a, &a);
1304   unpack_d ((FLO_union_type *) & arg_b, &b);
1305
1306   if (isnan (&a) || isnan (&b))
1307     return 1;                   /* false, truth <= 0 */
1308
1309   return __fpcmp_parts (&a, &b) ;
1310 }
1311 #endif
1312
1313 #endif /* ! US_SOFTWARE_GOFAST */
1314
1315 #if defined(L_si_to_sf) || defined(L_si_to_df)
1316 FLO_type
1317 si_to_float (SItype arg_a)
1318 {
1319   fp_number_type in;
1320
1321   in.class = CLASS_NUMBER;
1322   in.sign = arg_a < 0;
1323   if (!arg_a)
1324     {
1325       in.class = CLASS_ZERO;
1326     }
1327   else
1328     {
1329       in.normal_exp = FRACBITS + NGARDS;
1330       if (in.sign)
1331         {
1332           /* Special case for minint, since there is no +ve integer
1333              representation for it */
1334           if (arg_a == 0x80000000)
1335             {
1336               return -2147483648.0;
1337             }
1338           in.fraction.ll = (-arg_a);
1339         }
1340       else
1341         in.fraction.ll = arg_a;
1342
1343       while (in.fraction.ll < (1LL << (FRACBITS + NGARDS)))
1344         {
1345           in.fraction.ll <<= 1;
1346           in.normal_exp -= 1;
1347         }
1348     }
1349   return pack_d (&in);
1350 }
1351 #endif
1352
1353 #if defined(L_sf_to_si) || defined(L_df_to_si)
1354 SItype
1355 float_to_si (FLO_type arg_a)
1356 {
1357   fp_number_type a;
1358   SItype tmp;
1359
1360   unpack_d ((FLO_union_type *) & arg_a, &a);
1361   if (iszero (&a))
1362     return 0;
1363   if (isnan (&a))
1364     return 0;
1365   /* get reasonable MAX_SI_INT... */
1366   if (isinf (&a))
1367     return a.sign ? (-MAX_SI_INT)-1 : MAX_SI_INT;
1368   /* it is a number, but a small one */
1369   if (a.normal_exp < 0)
1370     return 0;
1371   if (a.normal_exp > 30)
1372     return a.sign ? (-MAX_SI_INT)-1 : MAX_SI_INT;
1373   tmp = a.fraction.ll >> ((FRACBITS + NGARDS) - a.normal_exp);
1374   return a.sign ? (-tmp) : (tmp);
1375 }
1376 #endif
1377
1378 #if defined(L_sf_to_usi) || defined(L_df_to_usi)
1379 #ifdef US_SOFTWARE_GOFAST
1380 /* While libgcc2.c defines its own __fixunssfsi and __fixunsdfsi routines,
1381    we also define them for GOFAST because the ones in libgcc2.c have the
1382    wrong names and I'd rather define these here and keep GOFAST CYG-LOC's
1383    out of libgcc2.c.  We can't define these here if not GOFAST because then
1384    there'd be duplicate copies.  */
1385
1386 USItype
1387 float_to_usi (FLO_type arg_a)
1388 {
1389   fp_number_type a;
1390
1391   unpack_d ((FLO_union_type *) & arg_a, &a);
1392   if (iszero (&a))
1393     return 0;
1394   if (isnan (&a))
1395     return 0;
1396   /* it is a negative number */
1397   if (a.sign)
1398     return 0;
1399   /* get reasonable MAX_USI_INT... */
1400   if (isinf (&a))
1401     return MAX_USI_INT;
1402   /* it is a number, but a small one */
1403   if (a.normal_exp < 0)
1404     return 0;
1405   if (a.normal_exp > 31)
1406     return MAX_USI_INT;
1407   else if (a.normal_exp > (FRACBITS + NGARDS))
1408     return a.fraction.ll << (a.normal_exp - (FRACBITS + NGARDS));
1409   else
1410     return a.fraction.ll >> ((FRACBITS + NGARDS) - a.normal_exp);
1411 }
1412 #endif
1413 #endif
1414
1415 #if defined(L_negate_sf) || defined(L_negate_df)
1416 FLO_type
1417 negate (FLO_type arg_a)
1418 {
1419   fp_number_type a;
1420
1421   unpack_d ((FLO_union_type *) & arg_a, &a);
1422   flip_sign (&a);
1423   return pack_d (&a);
1424 }
1425 #endif
1426
1427 #ifdef FLOAT
1428
1429 #if defined(L_make_sf)
1430 SFtype
1431 __make_fp(fp_class_type class,
1432              unsigned int sign,
1433              int exp,
1434              USItype frac)
1435 {
1436   fp_number_type in;
1437
1438   in.class = class;
1439   in.sign = sign;
1440   in.normal_exp = exp;
1441   in.fraction.ll = frac;
1442   return pack_d (&in);
1443 }
1444 #endif
1445
1446 #ifndef FLOAT_ONLY
1447
1448 /* This enables one to build an fp library that supports float but not double.
1449    Otherwise, we would get an undefined reference to __make_dp.
1450    This is needed for some 8-bit ports that can't handle well values that
1451    are 8-bytes in size, so we just don't support double for them at all.  */
1452
1453 extern DFtype __make_dp (fp_class_type, unsigned int, int, UDItype frac);
1454
1455 #if defined(L_sf_to_df)
1456 DFtype
1457 sf_to_df (SFtype arg_a)
1458 {
1459   fp_number_type in;
1460
1461   unpack_d ((FLO_union_type *) & arg_a, &in);
1462   return __make_dp (in.class, in.sign, in.normal_exp,
1463                     ((UDItype) in.fraction.ll) << F_D_BITOFF);
1464 }
1465 #endif
1466
1467 #endif
1468 #endif
1469
1470 #ifndef FLOAT
1471
1472 extern SFtype __make_fp (fp_class_type, unsigned int, int, USItype);
1473
1474 #if defined(L_make_df)
1475 DFtype
1476 __make_dp (fp_class_type class, unsigned int sign, int exp, UDItype frac)
1477 {
1478   fp_number_type in;
1479
1480   in.class = class;
1481   in.sign = sign;
1482   in.normal_exp = exp;
1483   in.fraction.ll = frac;
1484   return pack_d (&in);
1485 }
1486 #endif
1487
1488 #if defined(L_df_to_sf)
1489 SFtype
1490 df_to_sf (DFtype arg_a)
1491 {
1492   fp_number_type in;
1493   USItype sffrac;
1494
1495   unpack_d ((FLO_union_type *) & arg_a, &in);
1496
1497   sffrac = in.fraction.ll >> F_D_BITOFF;
1498
1499   /* We set the lowest guard bit in SFFRAC if we discarded any non
1500      zero bits.  */
1501   if ((in.fraction.ll & (((USItype) 1 << F_D_BITOFF) - 1)) != 0)
1502     sffrac |= 1;
1503
1504   return __make_fp (in.class, in.sign, in.normal_exp, sffrac);
1505 }
1506 #endif
1507
1508 #endif
1509 #endif /* !EXTENDED_FLOAT_STUBS */