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