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 #else   /* !EXTENDED_FLOAT_STUBS, rest of file */
 154
 155
 156 typedef SFtype __attribute__ ((mode (SF)));
 157 typedef DFtype __attribute__ ((mode (DF)));
 158
 159 typedef int HItype __attribute__ ((mode (HI)));
 160 typedef int SItype __attribute__ ((mode (SI)));
 161 typedef int DItype __attribute__ ((mode (DI)));
 162
 163 /* The type of the result of a fp compare */
 164 #ifndef CMPtype
 165 #define CMPtype SItype
 166 #endif
 167
 168 typedef unsigned int UHItype __attribute__ ((mode (HI)));
 169 typedef unsigned int USItype __attribute__ ((mode (SI)));
 170 typedef unsigned int UDItype __attribute__ ((mode (DI)));
 171
 172 #define MAX_SI_INT   ((SItype) ((unsigned) (~0)>>1))
 173 #define MAX_USI_INT  ((USItype) ~0)
 174
 175
 176 #ifdef FLOAT_ONLY
 177 #define NO_DI_MODE
 178 #endif
 179
 180 #ifdef FLOAT
 181 #       define NGARDS    7L
 182 #       define GARDROUND 0x3f
 183 #       define GARDMASK  0x7f
 184 #       define GARDMSB   0x40
 185 #       define EXPBITS 8
 186 #       define EXPBIAS 127
 187 #       define FRACBITS 23
 188 #       define EXPMAX (0xff)
 189 #       define QUIET_NAN 0x100000L
 190 #       define FRAC_NBITS 32
 191 #       define FRACHIGH  0x80000000L
 192 #       define FRACHIGH2 0xc0000000L
 193 #       define pack_d __pack_f
 194 #       define unpack_d __unpack_f
 195 #       define __fpcmp_parts __fpcmp_parts_f
 196         typedef USItype fractype;
 197         typedef UHItype halffractype;
 198         typedef SFtype FLO_type;
 199         typedef SItype intfrac;
 200
 201 #else
 202 #       define PREFIXFPDP dp
 203 #       define PREFIXSFDF df
 204 #       define NGARDS 8L
 205 #       define GARDROUND 0x7f
 206 #       define GARDMASK  0xff
 207 #       define GARDMSB   0x80
 208 #       define EXPBITS 11
 209 #       define EXPBIAS 1023
 210 #       define FRACBITS 52
 211 #       define EXPMAX (0x7ff)
 212 #       define QUIET_NAN 0x8000000000000LL
 213 #       define FRAC_NBITS 64
 214 #       define FRACHIGH  0x8000000000000000LL
 215 #       define FRACHIGH2 0xc000000000000000LL
 216 #       define pack_d __pack_d
 217 #       define unpack_d __unpack_d
 218 #       define __fpcmp_parts __fpcmp_parts_d
 219         typedef UDItype fractype;
 220         typedef USItype halffractype;
 221         typedef DFtype FLO_type;
 222         typedef DItype intfrac;
 223 #endif
 224
 225 #ifdef US_SOFTWARE_GOFAST
 226 #       ifdef FLOAT
 227 #               define add              fpadd
 228 #               define sub              fpsub
 229 #               define multiply         fpmul
 230 #               define divide           fpdiv
 231 #               define compare          fpcmp
 232 #               define si_to_float      sitofp
 233 #               define float_to_si      fptosi
 234 #               define float_to_usi     fptoui
 235 #               define negate           __negsf2
 236 #               define sf_to_df         fptodp
 237 #               define dptofp           dptofp
 238 #else
 239 #               define add              dpadd
 240 #               define sub              dpsub
 241 #               define multiply         dpmul
 242 #               define divide           dpdiv
 243 #               define compare          dpcmp
 244 #               define si_to_float      litodp
 245 #               define float_to_si      dptoli
 246 #               define float_to_usi     dptoul
 247 #               define negate           __negdf2
 248 #               define df_to_sf         dptofp
 249 #endif
 250 #else
 251 #       ifdef FLOAT
 252 #               define add              __addsf3
 253 #               define sub              __subsf3
 254 #               define multiply         __mulsf3
 255 #               define divide           __divsf3
 256 #               define compare          __cmpsf2
 257 #               define _eq_f2           __eqsf2
 258 #               define _ne_f2           __nesf2
 259 #               define _gt_f2           __gtsf2
 260 #               define _ge_f2           __gesf2
 261 #               define _lt_f2           __ltsf2
 262 #               define _le_f2           __lesf2
 263 #               define si_to_float      __floatsisf
 264 #               define float_to_si      __fixsfsi
 265 #               define float_to_usi     __fixunssfsi
 266 #               define negate           __negsf2
 267 #               define sf_to_df         __extendsfdf2
 268 #else
 269 #               define add              __adddf3
 270 #               define sub              __subdf3
 271 #               define multiply         __muldf3
 272 #               define divide           __divdf3
 273 #               define compare          __cmpdf2
 274 #               define _eq_f2           __eqdf2
 275 #               define _ne_f2           __nedf2
 276 #               define _gt_f2           __gtdf2
 277 #               define _ge_f2           __gedf2
 278 #               define _lt_f2           __ltdf2
 279 #               define _le_f2           __ledf2
 280 #               define si_to_float      __floatsidf
 281 #               define float_to_si      __fixdfsi
 282 #               define float_to_usi     __fixunsdfsi
 283 #               define negate           __negdf2
 284 #               define df_to_sf         __truncdfsf2
 285 #       endif
 286 #endif
 287
 288
 289 #ifndef INLINE
 290 #define INLINE __inline__
 291 #endif
 292
 293 /* Preserve the sticky-bit when shifting fractions to the right.  */
 294 #define LSHIFT(a) { a = (a & 1) | (a >> 1); }
 295
 296 /* numeric parameters */
 297 /* F_D_BITOFF is the number of bits offset between the MSB of the mantissa
 298    of a float and of a double. Assumes there are only two float types.
 299    (double::FRAC_BITS+double::NGARDS-(float::FRAC_BITS-float::NGARDS))
 300  */
 301 #define F_D_BITOFF (52+8-(23+7))
 302
 303
 304 #define NORMAL_EXPMIN (-(EXPBIAS)+1)
 305 #define IMPLICIT_1 (1LL<<(FRACBITS+NGARDS))
 306 #define IMPLICIT_2 (1LL<<(FRACBITS+1+NGARDS))
 307
 308 /* common types */
 309
 310 typedef enum
 311 {
 312   CLASS_SNAN,
 313   CLASS_QNAN,
 314   CLASS_ZERO,
 315   CLASS_NUMBER,
 316   CLASS_INFINITY
 317 } fp_class_type;
 318
 319 typedef struct
 320 {
 321 #ifdef SMALL_MACHINE
 322   char class;
 323   unsigned char sign;
 324   short normal_exp;
 325 #else
 326   fp_class_type class;
 327   unsigned int sign;
 328   int normal_exp;
 329 #endif
 330
 331   union
 332     {
 333       fractype ll;
 334       halffractype l[2];
 335     } fraction;
 336 } fp_number_type;
 337
 338 typedef union
 339 {
 340   FLO_type value;
 341   fractype value_raw;
 342
 343 #ifndef FLOAT
 344   halffractype words[2];
 345 #endif
 346
 347 #ifdef FLOAT_BIT_ORDER_MISMATCH
 348   struct
 349     {
 350       fractype fraction:FRACBITS __attribute__ ((packed));
 351       unsigned int exp:EXPBITS __attribute__ ((packed));
 352       unsigned int sign:1 __attribute__ ((packed));
 353     }
 354   bits;
 355 #endif
 356
 357 #ifdef _DEBUG_BITFLOAT
 358   struct
 359     {
 360       unsigned int sign:1 __attribute__ ((packed));
 361       unsigned int exp:EXPBITS __attribute__ ((packed));
 362       fractype fraction:FRACBITS __attribute__ ((packed));
 363     }
 364   bits_big_endian;
 365
 366   struct
 367     {
 368       fractype fraction:FRACBITS __attribute__ ((packed));
 369       unsigned int exp:EXPBITS __attribute__ ((packed));
 370       unsigned int sign:1 __attribute__ ((packed));
 371     }
 372   bits_little_endian;
 373 #endif
 374 }
 375 FLO_union_type;
 376
 377
 378 /* end of header */
 379
 380 /* IEEE "special" number predicates */
 381
 382 #ifdef NO_NANS
 383
 384 #define nan() 0
 385 #define isnan(x) 0
 386 #define isinf(x) 0
 387 #else
 388
 389 INLINE
 390 static fp_number_type *
 391 nan ()
 392 {
 393   static fp_number_type thenan;
 394
 395   return &thenan;
 396 }
 397
 398 INLINE
 399 static int
 400 isnan ( fp_number_type *  x)
 401 {
 402   return x->class == CLASS_SNAN || x->class == CLASS_QNAN;
 403 }
 404
 405 INLINE
 406 static int
 407 isinf ( fp_number_type *  x)
 408 {
 409   return x->class == CLASS_INFINITY;
 410 }
 411
 412 #endif
 413
 414 INLINE
 415 static int
 416 iszero ( fp_number_type *  x)
 417 {
 418   return x->class == CLASS_ZERO;
 419 }
 420
 421 INLINE
 422 static void
 423 flip_sign ( fp_number_type *  x)
 424 {
 425   x->sign = !x->sign;
 426 }
 427
 428 extern FLO_type pack_d ( fp_number_type * );
 429
 430 #if defined(L_pack_df) || defined(L_pack_sf)
 431 FLO_type
 432 pack_d ( fp_number_type *  src)
 433 {
 434   FLO_union_type dst;
 435   fractype fraction = src->fraction.ll; /* wasn't unsigned before? */
 436   int sign = src->sign;
 437   int exp = 0;
 438
 439   if (isnan (src))
 440     {
 441       exp = EXPMAX;
 442       if (src->class == CLASS_QNAN || 1)
 443         {
 444           fraction |= QUIET_NAN;
 445         }
 446     }
 447   else if (isinf (src))
 448     {
 449       exp = EXPMAX;
 450       fraction = 0;
 451     }
 452   else if (iszero (src))
 453     {
 454       exp = 0;
 455       fraction = 0;
 456     }
 457   else if (fraction == 0)
 458     {
 459       exp = 0;
 460       sign = 0;
 461     }
 462   else
 463     {
 464       if (src->normal_exp < NORMAL_EXPMIN)
 465         {
 466           /* This number's exponent is too low to fit into the bits
 467              available in the number, so we'll store 0 in the exponent and
 468              shift the fraction to the right to make up for it.  */
 469
 470           int shift = NORMAL_EXPMIN - src->normal_exp;
 471
 472           exp = 0;
 473
 474           if (shift > FRAC_NBITS - NGARDS)
 475             {
 476               /* No point shifting, since it's more that 64 out.  */
 477               fraction = 0;
 478             }
 479           else
 480             {
 481               /* Shift by the value */
 482               fraction >>= shift;
 483             }
 484           fraction >>= NGARDS;
 485         }
 486       else if (src->normal_exp > EXPBIAS)
 487         {
 488           exp = EXPMAX;
 489           fraction = 0;
 490         }
 491       else
 492         {
 493           exp = src->normal_exp + EXPBIAS;
 494           /* IF the gard bits are the all zero, but the first, then we're
 495              half way between two numbers, choose the one which makes the
 496              lsb of the answer 0.  */
 497           if ((fraction & GARDMASK) == GARDMSB)
 498             {
 499               if (fraction & (1 << NGARDS))
 500                 fraction += GARDROUND + 1;
 501             }
 502           else
 503             {
 504               /* Add a one to the guards to round up */
 505               fraction += GARDROUND;
 506             }
 507           if (fraction >= IMPLICIT_2)
 508             {
 509               fraction >>= 1;
 510               exp += 1;
 511             }
 512           fraction >>= NGARDS;
 513         }
 514     }
 515
 516   /* We previously used bitfields to store the number, but this doesn't
 517      handle little/big endian systems conveniently, so use shifts and
 518      masks */
 519 #ifdef FLOAT_BIT_ORDER_MISMATCH
 520   dst.bits.fraction = fraction;
 521   dst.bits.exp = exp;
 522   dst.bits.sign = sign;
 523 #else
 524   dst.value_raw = fraction & ((((fractype)1) << FRACBITS) - (fractype)1);
 525   dst.value_raw |= ((fractype) (exp & ((1 << EXPBITS) - 1))) << FRACBITS;
 526   dst.value_raw |= ((fractype) (sign & 1)) << (FRACBITS | EXPBITS);
 527 #endif
 528
 529 #if defined(FLOAT_WORD_ORDER_MISMATCH) && !defined(FLOAT)
 530   {
 531     halffractype tmp = dst.words[0];
 532     dst.words[0] = dst.words[1];
 533     dst.words[1] = tmp;
 534   }
 535 #endif
 536
 537   return dst.value;
 538 }
 539 #endif
 540
 541 extern void unpack_d (FLO_union_type *, fp_number_type *);
 542
 543 #if defined(L_unpack_df) || defined(L_unpack_sf)
 544 void
 545 unpack_d (FLO_union_type * src, fp_number_type * dst)
 546 {
 547   /* We previously used bitfields to store the number, but this doesn't
 548      handle little/big endian systems conveniently, so use shifts and
 549      masks */
 550   fractype fraction;
 551   int exp;
 552   int sign;
 553
 554 #if defined(FLOAT_WORD_ORDER_MISMATCH) && !defined(FLOAT)
 555   FLO_union_type swapped;
 556
 557   swapped.words[0] = src->words[1];
 558   swapped.words[1] = src->words[0];
 559   src = &swapped;
 560 #endif
 561
 562 #ifdef FLOAT_BIT_ORDER_MISMATCH
 563   fraction = src->bits.fraction;
 564   exp = src->bits.exp;
 565   sign = src->bits.sign;
 566 #else
 567   fraction = src->value_raw & ((((fractype)1) << FRACBITS) - (fractype)1);
 568   exp = ((int)(src->value_raw >> FRACBITS)) & ((1 << EXPBITS) - 1);
 569   sign = ((int)(src->value_raw >> (FRACBITS + EXPBITS))) & 1;
 570 #endif
 571
 572   dst->sign = sign;
 573   if (exp == 0)
 574     {
 575       /* Hmm.  Looks like 0 */
 576       if (fraction == 0)
 577         {
 578           /* tastes like zero */
 579           dst->class = CLASS_ZERO;
 580         }
 581       else
 582         {
 583           /* Zero exponent with non zero fraction - it's denormalized,
 584              so there isn't a leading implicit one - we'll shift it so
 585              it gets one.  */
 586           dst->normal_exp = exp - EXPBIAS + 1;
 587           fraction <<= NGARDS;
 588
 589           dst->class = CLASS_NUMBER;
 590 #if 1
 591           while (fraction < IMPLICIT_1)
 592             {
 593               fraction <<= 1;
 594               dst->normal_exp--;
 595             }
 596 #endif
 597           dst->fraction.ll = fraction;
 598         }
 599     }
 600   else if (exp == EXPMAX)
 601     {
 602       /* Huge exponent*/
 603       if (fraction == 0)
 604         {
 605           /* Attached to a zero fraction - means infinity */
 606           dst->class = CLASS_INFINITY;
 607         }
 608       else
 609         {
 610           /* Non zero fraction, means nan */
 611           if (fraction & QUIET_NAN)
 612             {
 613               dst->class = CLASS_QNAN;
 614             }
 615           else
 616             {
 617               dst->class = CLASS_SNAN;
 618             }
 619           /* Keep the fraction part as the nan number */
 620           dst->fraction.ll = fraction;
 621         }
 622     }
 623   else
 624     {
 625       /* Nothing strange about this number */
 626       dst->normal_exp = exp - EXPBIAS;
 627       dst->class = CLASS_NUMBER;
 628       dst->fraction.ll = (fraction << NGARDS) | IMPLICIT_1;
 629     }
 630 }
 631 #endif
 632
 633 #if defined(L_addsub_sf) || defined(L_addsub_df)
 634 static fp_number_type *
 635 _fpadd_parts (fp_number_type * a,
 636               fp_number_type * b,
 637               fp_number_type * tmp)
 638 {
 639   intfrac tfraction;
 640
 641   /* Put commonly used fields in local variables.  */
 642   int a_normal_exp;
 643   int b_normal_exp;
 644   fractype a_fraction;
 645   fractype b_fraction;
 646
 647   if (isnan (a))
 648     {
 649       return a;
 650     }
 651   if (isnan (b))
 652     {
 653       return b;
 654     }
 655   if (isinf (a))
 656     {
 657       /* Adding infinities with opposite signs yields a NaN.  */
 658       if (isinf (b) && a->sign != b->sign)
 659         return nan ();
 660       return a;
 661     }
 662   if (isinf (b))
 663     {
 664       return b;
 665     }
 666   if (iszero (b))
 667     {
 668       if (iszero (a))
 669         {
 670           *tmp = *a;
 671           tmp->sign = a->sign & b->sign;
 672           return tmp;
 673         }
 674       return a;
 675     }
 676   if (iszero (a))
 677     {
 678       return b;
 679     }
 680
 681   /* Got two numbers. shift the smaller and increment the exponent till
 682      they're the same */
 683   {
 684     int diff;
 685
 686     a_normal_exp = a->normal_exp;
 687     b_normal_exp = b->normal_exp;
 688     a_fraction = a->fraction.ll;
 689     b_fraction = b->fraction.ll;
 690
 691     diff = a_normal_exp - b_normal_exp;
 692
 693     if (diff < 0)
 694       diff = -diff;
 695     if (diff < FRAC_NBITS)
 696       {
 697         /* ??? This does shifts one bit at a time.  Optimize.  */
 698         while (a_normal_exp > b_normal_exp)
 699           {
 700             b_normal_exp++;
 701             LSHIFT (b_fraction);
 702           }
 703         while (b_normal_exp > a_normal_exp)
 704           {
 705             a_normal_exp++;
 706             LSHIFT (a_fraction);
 707           }
 708       }
 709     else
 710       {
 711         /* Somethings's up.. choose the biggest */
 712         if (a_normal_exp > b_normal_exp)
 713           {
 714             b_normal_exp = a_normal_exp;
 715             b_fraction = 0;
 716           }
 717         else
 718           {
 719             a_normal_exp = b_normal_exp;
 720             a_fraction = 0;
 721           }
 722       }
 723   }
 724
 725   if (a->sign != b->sign)
 726     {
 727       if (a->sign)
 728         {
 729           tfraction = -a_fraction + b_fraction;
 730         }
 731       else
 732         {
 733           tfraction = a_fraction - b_fraction;
 734         }
 735       if (tfraction > 0)
 736         {
 737           tmp->sign = 0;
 738           tmp->normal_exp = a_normal_exp;
 739           tmp->fraction.ll = tfraction;
 740         }
 741       else
 742         {
 743           tmp->sign = 1;
 744           tmp->normal_exp = a_normal_exp;
 745           tmp->fraction.ll = -tfraction;
 746         }
 747       /* and renormalize it */
 748
 749       while (tmp->fraction.ll < IMPLICIT_1 && tmp->fraction.ll)
 750         {
 751           tmp->fraction.ll <<= 1;
 752           tmp->normal_exp--;
 753         }
 754     }
 755   else
 756     {
 757       tmp->sign = a->sign;
 758       tmp->normal_exp = a_normal_exp;
 759       tmp->fraction.ll = a_fraction + b_fraction;
 760     }
 761   tmp->class = CLASS_NUMBER;
 762   /* Now the fraction is added, we have to shift down to renormalize the
 763      number */
 764
 765   if (tmp->fraction.ll >= IMPLICIT_2)
 766     {
 767       LSHIFT (tmp->fraction.ll);
 768       tmp->normal_exp++;
 769     }
 770   return tmp;
 771
 772 }
 773
 774 FLO_type
 775 add (FLO_type arg_a, FLO_type arg_b)
 776 {
 777   fp_number_type a;
 778   fp_number_type b;
 779   fp_number_type tmp;
 780   fp_number_type *res;
 781
 782   unpack_d ((FLO_union_type *) & arg_a, &a);
 783   unpack_d ((FLO_union_type *) & arg_b, &b);
 784
 785   res = _fpadd_parts (&a, &b, &tmp);
 786
 787   return pack_d (res);
 788 }
 789
 790 FLO_type
 791 sub (FLO_type arg_a, FLO_type arg_b)
 792 {
 793   fp_number_type a;
 794   fp_number_type b;
 795   fp_number_type tmp;
 796   fp_number_type *res;
 797
 798   unpack_d ((FLO_union_type *) & arg_a, &a);
 799   unpack_d ((FLO_union_type *) & arg_b, &b);
 800
 801   b.sign ^= 1;
 802
 803   res = _fpadd_parts (&a, &b, &tmp);
 804
 805   return pack_d (res);
 806 }
 807 #endif
 808
 809 #if defined(L_mul_sf) || defined(L_mul_df)
 810 static INLINE fp_number_type *
 811 _fpmul_parts ( fp_number_type *  a,
 812                fp_number_type *  b,
 813                fp_number_type * tmp)
 814 {
 815   fractype low = 0;
 816   fractype high = 0;
 817
 818   if (isnan (a))
 819     {
 820       a->sign = a->sign != b->sign;
 821       return a;
 822     }
 823   if (isnan (b))
 824     {
 825       b->sign = a->sign != b->sign;
 826       return b;
 827     }
 828   if (isinf (a))
 829     {
 830       if (iszero (b))
 831         return nan ();
 832       a->sign = a->sign != b->sign;
 833       return a;
 834     }
 835   if (isinf (b))
 836     {
 837       if (iszero (a))
 838         {
 839           return nan ();
 840         }
 841       b->sign = a->sign != b->sign;
 842       return b;
 843     }
 844   if (iszero (a))
 845     {
 846       a->sign = a->sign != b->sign;
 847       return a;
 848     }
 849   if (iszero (b))
 850     {
 851       b->sign = a->sign != b->sign;
 852       return b;
 853     }
 854
 855   /* Calculate the mantissa by multiplying both 64bit numbers to get a
 856      128 bit number */
 857   {
 858 #if defined(NO_DI_MODE)
 859     {
 860       fractype x = a->fraction.ll;
 861       fractype ylow = b->fraction.ll;
 862       fractype yhigh = 0;
 863       int bit;
 864
 865       /* ??? This does multiplies one bit at a time.  Optimize.  */
 866       for (bit = 0; bit < FRAC_NBITS; bit++)
 867         {
 868           int carry;
 869
 870           if (x & 1)
 871             {
 872               carry = (low += ylow) < ylow;
 873               high += yhigh + carry;
 874             }
 875           yhigh <<= 1;
 876           if (ylow & FRACHIGH)
 877             {
 878               yhigh |= 1;
 879             }
 880           ylow <<= 1;
 881           x >>= 1;
 882         }
 883     }
 884 #elif defined(FLOAT)
 885     {
 886       /* Multiplying two 32 bit numbers to get a 64 bit number  on
 887         a machine with DI, so we're safe */
 888
 889       DItype answer = (DItype)(a->fraction.ll) * (DItype)(b->fraction.ll);
 890
 891       high = answer >> 32;
 892       low = answer;
 893     }
 894 #else
 895     /* Doing a 64*64 to 128 */
 896     {
 897       UDItype nl = a->fraction.ll & 0xffffffff;
 898       UDItype nh = a->fraction.ll >> 32;
 899       UDItype ml = b->fraction.ll & 0xffffffff;
 900       UDItype mh = b->fraction.ll >>32;
 901       UDItype pp_ll = ml * nl;
 902       UDItype pp_hl = mh * nl;
 903       UDItype pp_lh = ml * nh;
 904       UDItype pp_hh = mh * nh;
 905       UDItype res2 = 0;
 906       UDItype res0 = 0;
 907       UDItype ps_hh__ = pp_hl + pp_lh;
 908       if (ps_hh__ < pp_hl)
 909         res2 += 0x100000000LL;
 910       pp_hl = (ps_hh__ << 32) & 0xffffffff00000000LL;
 911       res0 = pp_ll + pp_hl;
 912       if (res0 < pp_ll)
 913         res2++;
 914       res2 += ((ps_hh__ >> 32) & 0xffffffffL) + pp_hh;
 915       high = res2;
 916       low = res0;
 917     }
 918 #endif
 919   }
 920
 921   tmp->normal_exp = a->normal_exp + b->normal_exp;
 922   tmp->sign = a->sign != b->sign;
 923 #ifdef FLOAT
 924   tmp->normal_exp += 2;         /* ??????????????? */
 925 #else
 926   tmp->normal_exp += 4;         /* ??????????????? */
 927 #endif
 928   while (high >= IMPLICIT_2)
 929     {
 930       tmp->normal_exp++;
 931       if (high & 1)
 932         {
 933           low >>= 1;
 934           low |= FRACHIGH;
 935         }
 936       high >>= 1;
 937     }
 938   while (high < IMPLICIT_1)
 939     {
 940       tmp->normal_exp--;
 941
 942       high <<= 1;
 943       if (low & FRACHIGH)
 944         high |= 1;
 945       low <<= 1;
 946     }
 947   /* rounding is tricky. if we only round if it won't make us round later. */
 948 #if 0
 949   if (low & FRACHIGH2)
 950     {
 951       if (((high & GARDMASK) != GARDMSB)
 952           && (((high + 1) & GARDMASK) == GARDMSB))
 953         {
 954           /* don't round, it gets done again later. */
 955         }
 956       else
 957         {
 958           high++;
 959         }
 960     }
 961 #endif
 962   if ((high & GARDMASK) == GARDMSB)
 963     {
 964       if (high & (1 << NGARDS))
 965         {
 966           /* half way, so round to even */
 967           high += GARDROUND + 1;
 968         }
 969       else if (low)
 970         {
 971           /* but we really weren't half way */
 972           high += GARDROUND + 1;
 973         }
 974     }
 975   tmp->fraction.ll = high;
 976   tmp->class = CLASS_NUMBER;
 977   return tmp;
 978 }
 979
 980 FLO_type
 981 multiply (FLO_type arg_a, FLO_type arg_b)
 982 {
 983   fp_number_type a;
 984   fp_number_type b;
 985   fp_number_type tmp;
 986   fp_number_type *res;
 987
 988   unpack_d ((FLO_union_type *) & arg_a, &a);
 989   unpack_d ((FLO_union_type *) & arg_b, &b);
 990
 991   res = _fpmul_parts (&a, &b, &tmp);
 992
 993   return pack_d (res);
 994 }
 995 #endif
 996
 997 #if defined(L_div_sf) || defined(L_div_df)
 998 static INLINE fp_number_type *
 999 _fpdiv_parts (fp_number_type * a,
1000               fp_number_type * b,
1001               fp_number_type * tmp)
1002 {
1003   fractype bit;
1004   fractype numerator;
1005   fractype denominator;
1006   fractype quotient;
1007
1008   if (isnan (a))
1009     {
1010       return a;
1011     }
1012   if (isnan (b))
1013     {
1014       return b;
1015     }
1016
1017   a->sign = a->sign ^ b->sign;
1018
1019   if (isinf (a) || iszero (a))
1020     {
1021       if (a->class == b->class)
1022         return nan ();
1023       return a;
1024     }
1025
1026   if (isinf (b))
1027     {
1028       a->fraction.ll = 0;
1029       a->normal_exp = 0;
1030       return a;
1031     }
1032   if (iszero (b))
1033     {
1034       a->class = CLASS_INFINITY;
1035       return a;
1036     }
1037
1038   /* Calculate the mantissa by multiplying both 64bit numbers to get a
1039      128 bit number */
1040   {
1041     /* quotient =
1042        ( numerator / denominator) * 2^(numerator exponent -  denominator exponent)
1043      */
1044
1045     a->normal_exp = a->normal_exp - b->normal_exp;
1046     numerator = a->fraction.ll;
1047     denominator = b->fraction.ll;
1048
1049     if (numerator < denominator)
1050       {
1051         /* Fraction will be less than 1.0 */
1052         numerator *= 2;
1053         a->normal_exp--;
1054       }
1055     bit = IMPLICIT_1;
1056     quotient = 0;
1057     /* ??? Does divide one bit at a time.  Optimize.  */
1058     while (bit)
1059       {
1060         if (numerator >= denominator)
1061           {
1062             quotient |= bit;
1063             numerator -= denominator;
1064           }
1065         bit >>= 1;
1066         numerator *= 2;
1067       }
1068
1069     if ((quotient & GARDMASK) == GARDMSB)
1070       {
1071         if (quotient & (1 << NGARDS))
1072           {
1073             /* half way, so round to even */
1074             quotient += GARDROUND + 1;
1075           }
1076         else if (numerator)
1077           {
1078             /* but we really weren't half way, more bits exist */
1079             quotient += GARDROUND + 1;
1080           }
1081       }
1082
1083     a->fraction.ll = quotient;
1084     return (a);
1085   }
1086 }
1087
1088 FLO_type
1089 divide (FLO_type arg_a, FLO_type arg_b)
1090 {
1091   fp_number_type a;
1092   fp_number_type b;
1093   fp_number_type tmp;
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, &tmp);
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 == 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 */