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 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   if (isinf (a) || iszero (a))
 948     {
 949       if (a->class == b->class)
 950         return nan ();
 951       return a;
 952     }
 953   a->sign = a->sign ^ b->sign;
 954
 955   if (isinf (b))
 956     {
 957       a->fraction.ll = 0;
 958       a->normal_exp = 0;
 959       return a;
 960     }
 961   if (iszero (b))
 962     {
 963       a->class = CLASS_INFINITY;
 964       return b;
 965     }
 966
 967   /* Calculate the mantissa by multiplying both 64bit numbers to get a
 968      128 bit number */
 969   {
 970     int carry;
 971     intfrac d0, d1;             /* weren't unsigned before ??? */
 972
 973     /* quotient =
 974        ( numerator / denominator) * 2^(numerator exponent -  denominator exponent)
 975      */
 976
 977     a->normal_exp = a->normal_exp - b->normal_exp;
 978     numerator = a->fraction.ll;
 979     denominator = b->fraction.ll;
 980
 981     if (numerator < denominator)
 982       {
 983         /* Fraction will be less than 1.0 */
 984         numerator *= 2;
 985         a->normal_exp--;
 986       }
 987     bit = IMPLICIT_1;
 988     quotient = 0;
 989     /* ??? Does divide one bit at a time.  Optimize.  */
 990     while (bit)
 991       {
 992         if (numerator >= denominator)
 993           {
 994             quotient |= bit;
 995             numerator -= denominator;
 996           }
 997         bit >>= 1;
 998         numerator *= 2;
 999       }
1000
1001     if ((quotient & GARDMASK) == GARDMSB)
1002       {
1003         if (quotient & (1 << NGARDS))
1004           {
1005             /* half way, so round to even */
1006             quotient += GARDROUND + 1;
1007           }
1008         else if (numerator)
1009           {
1010             /* but we really weren't half way, more bits exist */
1011             quotient += GARDROUND + 1;
1012           }
1013       }
1014
1015     a->fraction.ll = quotient;
1016     return (a);
1017   }
1018 }
1019
1020 FLO_type
1021 divide (FLO_type arg_a, FLO_type arg_b)
1022 {
1023   fp_number_type a;
1024   fp_number_type b;
1025   fp_number_type tmp;
1026   fp_number_type *res;
1027
1028   unpack_d ((FLO_union_type *) & arg_a, &a);
1029   unpack_d ((FLO_union_type *) & arg_b, &b);
1030
1031   res = _fpdiv_parts (&a, &b, &tmp);
1032
1033   return pack_d (res);
1034 }
1035
1036 /* according to the demo, fpcmp returns a comparison with 0... thus
1037    a<b -> -1
1038    a==b -> 0
1039    a>b -> +1
1040  */
1041
1042 static int
1043 _fpcmp_parts (fp_number_type * a, fp_number_type * b)
1044 {
1045 #if 0
1046   /* either nan -> unordered. Must be checked outside of this routine. */
1047   if (isnan (a) && isnan (b))
1048     {
1049       return 1;                 /* still unordered! */
1050     }
1051 #endif
1052
1053   if (isnan (a) || isnan (b))
1054     {
1055       return 1;                 /* how to indicate unordered compare? */
1056     }
1057   if (isinf (a) && isinf (b))
1058     {
1059       /* +inf > -inf, but +inf != +inf */
1060       /* b    \a| +inf(0)| -inf(1)
1061        ______\+--------+--------
1062        +inf(0)| a==b(0)| a<b(-1)
1063        -------+--------+--------
1064        -inf(1)| a>b(1) | a==b(0)
1065        -------+--------+--------
1066        So since unordered must be non zero, just line up the columns...
1067        */
1068       return b->sign - a->sign;
1069     }
1070   /* but not both... */
1071   if (isinf (a))
1072     {
1073       return a->sign ? -1 : 1;
1074     }
1075   if (isinf (b))
1076     {
1077       return b->sign ? 1 : -1;
1078     }
1079   if (iszero (a) && iszero (b))
1080     {
1081       return 0;
1082     }
1083   if (iszero (a))
1084     {
1085       return b->sign ? 1 : -1;
1086     }
1087   if (iszero (b))
1088     {
1089       return a->sign ? -1 : 1;
1090     }
1091   /* now both are "normal". */
1092   if (a->sign != b->sign)
1093     {
1094       /* opposite signs */
1095       return a->sign ? -1 : 1;
1096     }
1097   /* same sign; exponents? */
1098   if (a->normal_exp > b->normal_exp)
1099     {
1100       return a->sign ? -1 : 1;
1101     }
1102   if (a->normal_exp < b->normal_exp)
1103     {
1104       return a->sign ? 1 : -1;
1105     }
1106   /* same exponents; check size. */
1107   if (a->fraction.ll > b->fraction.ll)
1108     {
1109       return a->sign ? -1 : 1;
1110     }
1111   if (a->fraction.ll < b->fraction.ll)
1112     {
1113       return a->sign ? 1 : -1;
1114     }
1115   /* after all that, they're equal. */
1116   return 0;
1117 }
1118
1119 CMPtype
1120 compare (FLO_type arg_a, FLO_type arg_b)
1121 {
1122   fp_number_type a;
1123   fp_number_type b;
1124
1125   unpack_d ((FLO_union_type *) & arg_a, &a);
1126   unpack_d ((FLO_union_type *) & arg_b, &b);
1127
1128   return _fpcmp_parts (&a, &b);
1129 }
1130
1131 #ifndef US_SOFTWARE_GOFAST
1132
1133 /* These should be optimized for their specific tasks someday.  */
1134
1135 CMPtype
1136 _eq_f2 (FLO_type arg_a, FLO_type arg_b)
1137 {
1138   fp_number_type a;
1139   fp_number_type b;
1140
1141   unpack_d ((FLO_union_type *) & arg_a, &a);
1142   unpack_d ((FLO_union_type *) & arg_b, &b);
1143
1144   if (isnan (&a) || isnan (&b))
1145     return 1;                   /* false, truth == 0 */
1146
1147   return _fpcmp_parts (&a, &b) ;
1148 }
1149
1150 CMPtype
1151 _ne_f2 (FLO_type arg_a, FLO_type arg_b)
1152 {
1153   fp_number_type a;
1154   fp_number_type b;
1155
1156   unpack_d ((FLO_union_type *) & arg_a, &a);
1157   unpack_d ((FLO_union_type *) & arg_b, &b);
1158
1159   if (isnan (&a) || isnan (&b))
1160     return 1;                   /* true, truth != 0 */
1161
1162   return  _fpcmp_parts (&a, &b) ;
1163 }
1164
1165 CMPtype
1166 _gt_f2 (FLO_type arg_a, FLO_type arg_b)
1167 {
1168   fp_number_type a;
1169   fp_number_type b;
1170
1171   unpack_d ((FLO_union_type *) & arg_a, &a);
1172   unpack_d ((FLO_union_type *) & arg_b, &b);
1173
1174   if (isnan (&a) || isnan (&b))
1175     return -1;                  /* false, truth > 0 */
1176
1177   return _fpcmp_parts (&a, &b);
1178 }
1179
1180 CMPtype
1181 _ge_f2 (FLO_type arg_a, FLO_type arg_b)
1182 {
1183   fp_number_type a;
1184   fp_number_type b;
1185
1186   unpack_d ((FLO_union_type *) & arg_a, &a);
1187   unpack_d ((FLO_union_type *) & arg_b, &b);
1188
1189   if (isnan (&a) || isnan (&b))
1190     return -1;                  /* false, truth >= 0 */
1191   return _fpcmp_parts (&a, &b) ;
1192 }
1193
1194 CMPtype
1195 _lt_f2 (FLO_type arg_a, FLO_type arg_b)
1196 {
1197   fp_number_type a;
1198   fp_number_type b;
1199
1200   unpack_d ((FLO_union_type *) & arg_a, &a);
1201   unpack_d ((FLO_union_type *) & arg_b, &b);
1202
1203   if (isnan (&a) || isnan (&b))
1204     return 1;                   /* false, truth < 0 */
1205
1206   return _fpcmp_parts (&a, &b);
1207 }
1208
1209 CMPtype
1210 _le_f2 (FLO_type arg_a, FLO_type arg_b)
1211 {
1212   fp_number_type a;
1213   fp_number_type b;
1214
1215   unpack_d ((FLO_union_type *) & arg_a, &a);
1216   unpack_d ((FLO_union_type *) & arg_b, &b);
1217
1218   if (isnan (&a) || isnan (&b))
1219     return 1;                   /* false, truth <= 0 */
1220
1221   return _fpcmp_parts (&a, &b) ;
1222 }
1223
1224 #endif /* ! US_SOFTWARE_GOFAST */
1225
1226 FLO_type
1227 si_to_float (SItype arg_a)
1228 {
1229   fp_number_type in;
1230
1231   in.class = CLASS_NUMBER;
1232   in.sign = arg_a < 0;
1233   if (!arg_a)
1234     {
1235       in.class = CLASS_ZERO;
1236     }
1237   else
1238     {
1239       in.normal_exp = FRACBITS + NGARDS;
1240       if (in.sign)
1241         {
1242           /* Special case for minint, since there is no +ve integer
1243              representation for it */
1244           if (arg_a == 0x80000000)
1245             {
1246               return -2147483648.0;
1247             }
1248           in.fraction.ll = (-arg_a);
1249         }
1250       else
1251         in.fraction.ll = arg_a;
1252
1253       while (in.fraction.ll < (1LL << (FRACBITS + NGARDS)))
1254         {
1255           in.fraction.ll <<= 1;
1256           in.normal_exp -= 1;
1257         }
1258     }
1259   return pack_d (&in);
1260 }
1261
1262 SItype
1263 float_to_si (FLO_type arg_a)
1264 {
1265   fp_number_type a;
1266   SItype tmp;
1267
1268   unpack_d ((FLO_union_type *) & arg_a, &a);
1269   if (iszero (&a))
1270     return 0;
1271   if (isnan (&a))
1272     return 0;
1273   /* get reasonable MAX_SI_INT... */
1274   if (isinf (&a))
1275     return a.sign ? MAX_SI_INT : (-MAX_SI_INT)-1;
1276   /* it is a number, but a small one */
1277   if (a.normal_exp < 0)
1278     return 0;
1279   if (a.normal_exp > 30)
1280     return a.sign ? (-MAX_SI_INT)-1 : MAX_SI_INT;
1281   tmp = a.fraction.ll >> ((FRACBITS + NGARDS) - a.normal_exp);
1282   return a.sign ? (-tmp) : (tmp);
1283 }
1284
1285 #ifdef US_SOFTWARE_GOFAST
1286 /* While libgcc2.c defines its own __fixunssfsi and __fixunsdfsi routines,
1287    we also define them for GOFAST because the ones in libgcc2.c have the
1288    wrong names and I'd rather define these here and keep GOFAST CYG-LOC's
1289    out of libgcc2.c.  We can't define these here if not GOFAST because then
1290    there'd be duplicate copies.  */
1291
1292 USItype
1293 float_to_usi (FLO_type arg_a)
1294 {
1295   fp_number_type a;
1296
1297   unpack_d ((FLO_union_type *) & arg_a, &a);
1298   if (iszero (&a))
1299     return 0;
1300   if (isnan (&a))
1301     return 0;
1302   /* get reasonable MAX_USI_INT... */
1303   if (isinf (&a))
1304     return a.sign ? MAX_USI_INT : 0;
1305   /* it is a negative number */
1306   if (a.sign)
1307     return 0;
1308   /* it is a number, but a small one */
1309   if (a.normal_exp < 0)
1310     return 0;
1311   if (a.normal_exp > 31)
1312     return MAX_USI_INT;
1313   else if (a.normal_exp > (FRACBITS + NGARDS))
1314     return a.fraction.ll << ((FRACBITS + NGARDS) - a.normal_exp);
1315   else
1316     return a.fraction.ll >> ((FRACBITS + NGARDS) - a.normal_exp);
1317 }
1318 #endif
1319
1320 FLO_type
1321 negate (FLO_type arg_a)
1322 {
1323   fp_number_type a;
1324
1325   unpack_d ((FLO_union_type *) & arg_a, &a);
1326   flip_sign (&a);
1327   return pack_d (&a);
1328 }
1329
1330 #ifdef FLOAT
1331
1332 SFtype
1333 __make_fp(fp_class_type class,
1334              unsigned int sign,
1335              int exp,
1336              USItype frac)
1337 {
1338   fp_number_type in;
1339
1340   in.class = class;
1341   in.sign = sign;
1342   in.normal_exp = exp;
1343   in.fraction.ll = frac;
1344   return pack_d (&in);
1345 }
1346
1347 #ifndef FLOAT_ONLY
1348
1349 /* This enables one to build an fp library that supports float but not double.
1350    Otherwise, we would get an undefined reference to __make_dp.
1351    This is needed for some 8-bit ports that can't handle well values that
1352    are 8-bytes in size, so we just don't support double for them at all.  */
1353
1354 extern DFtype __make_dp (fp_class_type, unsigned int, int, UDItype frac);
1355
1356 DFtype
1357 sf_to_df (SFtype arg_a)
1358 {
1359   fp_number_type in;
1360
1361   unpack_d ((FLO_union_type *) & arg_a, &in);
1362   return __make_dp (in.class, in.sign, in.normal_exp,
1363                     ((UDItype) in.fraction.ll) << F_D_BITOFF);
1364 }
1365
1366 #endif
1367 #endif
1368
1369 #ifndef FLOAT
1370
1371 extern SFtype __make_fp (fp_class_type, unsigned int, int, USItype);
1372
1373 DFtype
1374 __make_dp (fp_class_type class, unsigned int sign, int exp, UDItype frac)
1375 {
1376   fp_number_type in;
1377
1378   in.class = class;
1379   in.sign = sign;
1380   in.normal_exp = exp;
1381   in.fraction.ll = frac;
1382   return pack_d (&in);
1383 }
1384
1385 SFtype
1386 df_to_sf (DFtype arg_a)
1387 {
1388   fp_number_type in;
1389
1390   unpack_d ((FLO_union_type *) & arg_a, &in);
1391   return __make_fp (in.class, in.sign, in.normal_exp,
1392                     in.fraction.ll >> F_D_BITOFF);
1393 }
1394
1395 #endif
1396 #endif /* !EXTENDED_FLOAT_STUBS */