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-1998, 2000 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_unord_sf
  82 #define L_unord_df
  83 #define L_si_to_sf
  84 #define L_si_to_df
  85 #define L_sf_to_si
  86 #define L_df_to_si
  87 #define L_f_to_usi
  88 #define L_df_to_usi
  89 #define L_negate_sf
  90 #define L_negate_df
  91 #define L_make_sf
  92 #define L_make_df
  93 #define L_sf_to_df
  94 #define L_df_to_sf
  95 #ifdef FLOAT
  96 #define L_thenan_sf
  97 #else
  98 #define L_thenan_df
  99 #endif
 100 #endif
 101
 102 /* The following macros can be defined to change the behaviour of this file:
 103    FLOAT: Implement a `float', aka SFmode, fp library.  If this is not
 104      defined, then this file implements a `double', aka DFmode, fp library.
 105    FLOAT_ONLY: Used with FLOAT, to implement a `float' only library, i.e.
 106      don't include float->double conversion which requires the double library.
 107      This is useful only for machines which can't support doubles, e.g. some
 108      8-bit processors.
 109    CMPtype: Specify the type that floating point compares should return.
 110      This defaults to SItype, aka int.
 111    US_SOFTWARE_GOFAST: This makes all entry points use the same names as the
 112      US Software goFast library.  If this is not defined, the entry points use
 113      the same names as libgcc1.c.
 114    _DEBUG_BITFLOAT: This makes debugging the code a little easier, by adding
 115      two integers to the FLO_union_type.
 116    NO_NANS: Disable nan and infinity handling
 117    SMALL_MACHINE: Useful when operations on QIs and HIs are faster
 118      than on an SI */
 119
 120 /* We don't currently support extended floats (long doubles) on machines
 121    without hardware to deal with them.
 122
 123    These stubs are just to keep the linker from complaining about unresolved
 124    references which can be pulled in from libio & libstdc++, even if the
 125    user isn't using long doubles.  However, they may generate an unresolved
 126    external to abort if abort is not used by the function, and the stubs
 127    are referenced from within libc, since libgcc goes before and after the
 128    system library.  */
 129
 130 #ifdef EXTENDED_FLOAT_STUBS
 131 __truncxfsf2 (){ abort(); }
 132 __extendsfxf2 (){ abort(); }
 133 __addxf3 (){ abort(); }
 134 __divxf3 (){ abort(); }
 135 __eqxf2 (){ abort(); }
 136 __extenddfxf2 (){ abort(); }
 137 __gtxf2 (){ abort(); }
 138 __lexf2 (){ abort(); }
 139 __ltxf2 (){ abort(); }
 140 __mulxf3 (){ abort(); }
 141 __negxf2 (){ abort(); }
 142 __nexf2 (){ abort(); }
 143 __subxf3 (){ abort(); }
 144 __truncxfdf2 (){ abort(); }
 145
 146 __trunctfsf2 (){ abort(); }
 147 __extendsftf2 (){ abort(); }
 148 __addtf3 (){ abort(); }
 149 __divtf3 (){ abort(); }
 150 __eqtf2 (){ abort(); }
 151 __extenddftf2 (){ abort(); }
 152 __gttf2 (){ abort(); }
 153 __letf2 (){ abort(); }
 154 __lttf2 (){ abort(); }
 155 __multf3 (){ abort(); }
 156 __negtf2 (){ abort(); }
 157 __netf2 (){ abort(); }
 158 __subtf3 (){ abort(); }
 159 __trunctfdf2 (){ abort(); }
 160 __gexf2 (){ abort(); }
 161 __fixxfsi (){ abort(); }
 162 __floatsixf (){ abort(); }
 163 #else   /* !EXTENDED_FLOAT_STUBS, rest of file */
 164
 165
 166 typedef float SFtype __attribute__ ((mode (SF)));
 167 typedef float DFtype __attribute__ ((mode (DF)));
 168
 169 typedef int HItype __attribute__ ((mode (HI)));
 170 typedef int SItype __attribute__ ((mode (SI)));
 171 typedef int DItype __attribute__ ((mode (DI)));
 172
 173 /* The type of the result of a fp compare */
 174 #ifndef CMPtype
 175 #define CMPtype SItype
 176 #endif
 177
 178 typedef unsigned int UHItype __attribute__ ((mode (HI)));
 179 typedef unsigned int USItype __attribute__ ((mode (SI)));
 180 typedef unsigned int UDItype __attribute__ ((mode (DI)));
 181
 182 #define MAX_SI_INT   ((SItype) ((unsigned) (~0)>>1))
 183 #define MAX_USI_INT  ((USItype) ~0)
 184
 185
 186 #ifdef FLOAT_ONLY
 187 #define NO_DI_MODE
 188 #endif
 189
 190 #ifdef FLOAT
 191 #       define NGARDS    7L
 192 #       define GARDROUND 0x3f
 193 #       define GARDMASK  0x7f
 194 #       define GARDMSB   0x40
 195 #       define EXPBITS 8
 196 #       define EXPBIAS 127
 197 #       define FRACBITS 23
 198 #       define EXPMAX (0xff)
 199 #       define QUIET_NAN 0x100000L
 200 #       define FRAC_NBITS 32
 201 #       define FRACHIGH  0x80000000L
 202 #       define FRACHIGH2 0xc0000000L
 203 #       define pack_d __pack_f
 204 #       define unpack_d __unpack_f
 205 #       define __fpcmp_parts __fpcmp_parts_f
 206         typedef USItype fractype;
 207         typedef UHItype halffractype;
 208         typedef SFtype FLO_type;
 209         typedef SItype intfrac;
 210
 211 #else
 212 #       define PREFIXFPDP dp
 213 #       define PREFIXSFDF df
 214 #       define NGARDS 8L
 215 #       define GARDROUND 0x7f
 216 #       define GARDMASK  0xff
 217 #       define GARDMSB   0x80
 218 #       define EXPBITS 11
 219 #       define EXPBIAS 1023
 220 #       define FRACBITS 52
 221 #       define EXPMAX (0x7ff)
 222 #       define QUIET_NAN 0x8000000000000LL
 223 #       define FRAC_NBITS 64
 224 #       define FRACHIGH  0x8000000000000000LL
 225 #       define FRACHIGH2 0xc000000000000000LL
 226 #       define pack_d __pack_d
 227 #       define unpack_d __unpack_d
 228 #       define __fpcmp_parts __fpcmp_parts_d
 229         typedef UDItype fractype;
 230         typedef USItype halffractype;
 231         typedef DFtype FLO_type;
 232         typedef DItype intfrac;
 233 #endif
 234
 235 #ifdef US_SOFTWARE_GOFAST
 236 #       ifdef FLOAT
 237 #               define add              fpadd
 238 #               define sub              fpsub
 239 #               define multiply         fpmul
 240 #               define divide           fpdiv
 241 #               define compare          fpcmp
 242 #               define si_to_float      sitofp
 243 #               define float_to_si      fptosi
 244 #               define float_to_usi     fptoui
 245 #               define negate           __negsf2
 246 #               define sf_to_df         fptodp
 247 #               define dptofp           dptofp
 248 #else
 249 #               define add              dpadd
 250 #               define sub              dpsub
 251 #               define multiply         dpmul
 252 #               define divide           dpdiv
 253 #               define compare          dpcmp
 254 #               define si_to_float      litodp
 255 #               define float_to_si      dptoli
 256 #               define float_to_usi     dptoul
 257 #               define negate           __negdf2
 258 #               define df_to_sf         dptofp
 259 #endif
 260 #else
 261 #       ifdef FLOAT
 262 #               define add              __addsf3
 263 #               define sub              __subsf3
 264 #               define multiply         __mulsf3
 265 #               define divide           __divsf3
 266 #               define compare          __cmpsf2
 267 #               define _eq_f2           __eqsf2
 268 #               define _ne_f2           __nesf2
 269 #               define _gt_f2           __gtsf2
 270 #               define _ge_f2           __gesf2
 271 #               define _lt_f2           __ltsf2
 272 #               define _le_f2           __lesf2
 273 #               define _unord_f2        __unordsf2
 274 #               define si_to_float      __floatsisf
 275 #               define float_to_si      __fixsfsi
 276 #               define float_to_usi     __fixunssfsi
 277 #               define negate           __negsf2
 278 #               define sf_to_df         __extendsfdf2
 279 #else
 280 #               define add              __adddf3
 281 #               define sub              __subdf3
 282 #               define multiply         __muldf3
 283 #               define divide           __divdf3
 284 #               define compare          __cmpdf2
 285 #               define _eq_f2           __eqdf2
 286 #               define _ne_f2           __nedf2
 287 #               define _gt_f2           __gtdf2
 288 #               define _ge_f2           __gedf2
 289 #               define _lt_f2           __ltdf2
 290 #               define _le_f2           __ledf2
 291 #               define _unord_f2        __unorddf2
 292 #               define si_to_float      __floatsidf
 293 #               define float_to_si      __fixdfsi
 294 #               define float_to_usi     __fixunsdfsi
 295 #               define negate           __negdf2
 296 #               define df_to_sf         __truncdfsf2
 297 #       endif
 298 #endif
 299
 300
 301 #ifndef INLINE
 302 #define INLINE __inline__
 303 #endif
 304
 305 /* Preserve the sticky-bit when shifting fractions to the right.  */
 306 #define LSHIFT(a) { a = (a & 1) | (a >> 1); }
 307
 308 /* numeric parameters */
 309 /* F_D_BITOFF is the number of bits offset between the MSB of the mantissa
 310    of a float and of a double. Assumes there are only two float types.
 311    (double::FRAC_BITS+double::NGARDS-(float::FRAC_BITS-float::NGARDS))
 312  */
 313 #define F_D_BITOFF (52+8-(23+7))
 314
 315
 316 #define NORMAL_EXPMIN (-(EXPBIAS)+1)
 317 #define IMPLICIT_1 (1LL<<(FRACBITS+NGARDS))
 318 #define IMPLICIT_2 (1LL<<(FRACBITS+1+NGARDS))
 319
 320 /* common types */
 321
 322 typedef enum
 323 {
 324   CLASS_SNAN,
 325   CLASS_QNAN,
 326   CLASS_ZERO,
 327   CLASS_NUMBER,
 328   CLASS_INFINITY
 329 } fp_class_type;
 330
 331 typedef struct
 332 {
 333 #ifdef SMALL_MACHINE
 334   char class;
 335   unsigned char sign;
 336   short normal_exp;
 337 #else
 338   fp_class_type class;
 339   unsigned int sign;
 340   int normal_exp;
 341 #endif
 342
 343   union
 344     {
 345       fractype ll;
 346       halffractype l[2];
 347     } fraction;
 348 } fp_number_type;
 349
 350 typedef union
 351 {
 352   FLO_type value;
 353   fractype value_raw;
 354
 355 #ifndef FLOAT
 356   halffractype words[2];
 357 #endif
 358
 359 #ifdef FLOAT_BIT_ORDER_MISMATCH
 360   struct
 361     {
 362       fractype fraction:FRACBITS __attribute__ ((packed));
 363       unsigned int exp:EXPBITS __attribute__ ((packed));
 364       unsigned int sign:1 __attribute__ ((packed));
 365     }
 366   bits;
 367 #endif
 368
 369 #ifdef _DEBUG_BITFLOAT
 370   struct
 371     {
 372       unsigned int sign:1 __attribute__ ((packed));
 373       unsigned int exp:EXPBITS __attribute__ ((packed));
 374       fractype fraction:FRACBITS __attribute__ ((packed));
 375     }
 376   bits_big_endian;
 377
 378   struct
 379     {
 380       fractype fraction:FRACBITS __attribute__ ((packed));
 381       unsigned int exp:EXPBITS __attribute__ ((packed));
 382       unsigned int sign:1 __attribute__ ((packed));
 383     }
 384   bits_little_endian;
 385 #endif
 386 }
 387 FLO_union_type;
 388
 389
 390 /* end of header */
 391
 392 /* IEEE "special" number predicates */
 393
 394 #ifdef NO_NANS
 395
 396 #define nan() 0
 397 #define isnan(x) 0
 398 #define isinf(x) 0
 399 #else
 400
 401 #if   defined L_thenan_sf
 402 const fp_number_type __thenan_sf = { CLASS_SNAN, 0, 0, (fractype) 0 };
 403 #elif defined L_thenan_df
 404 const fp_number_type __thenan_df = { CLASS_SNAN, 0, 0, (fractype) 0 };
 405 #elif defined FLOAT
 406 extern const fp_number_type __thenan_sf;
 407 #else
 408 extern const fp_number_type __thenan_df;
 409 #endif
 410
 411 INLINE
 412 static fp_number_type *
 413 nan ()
 414 {
 415   /* Discard the const qualifier... */
 416 #ifdef FLOAT
 417   return (fp_number_type *) (& __thenan_sf);
 418 #else
 419   return (fp_number_type *) (& __thenan_df);
 420 #endif
 421 }
 422
 423 INLINE
 424 static int
 425 isnan ( fp_number_type *  x)
 426 {
 427   return x->class == CLASS_SNAN || x->class == CLASS_QNAN;
 428 }
 429
 430 INLINE
 431 static int
 432 isinf ( fp_number_type *  x)
 433 {
 434   return x->class == CLASS_INFINITY;
 435 }
 436
 437 #endif
 438
 439 INLINE
 440 static int
 441 iszero ( fp_number_type *  x)
 442 {
 443   return x->class == CLASS_ZERO;
 444 }
 445
 446 INLINE
 447 static void
 448 flip_sign ( fp_number_type *  x)
 449 {
 450   x->sign = !x->sign;
 451 }
 452
 453 extern FLO_type pack_d ( fp_number_type * );
 454
 455 #if defined(L_pack_df) || defined(L_pack_sf)
 456 FLO_type
 457 pack_d ( fp_number_type *  src)
 458 {
 459   FLO_union_type dst;
 460   fractype fraction = src->fraction.ll; /* wasn't unsigned before? */
 461   int sign = src->sign;
 462   int exp = 0;
 463
 464   if (isnan (src))
 465     {
 466       exp = EXPMAX;
 467       if (src->class == CLASS_QNAN || 1)
 468         {
 469           fraction |= QUIET_NAN;
 470         }
 471     }
 472   else if (isinf (src))
 473     {
 474       exp = EXPMAX;
 475       fraction = 0;
 476     }
 477   else if (iszero (src))
 478     {
 479       exp = 0;
 480       fraction = 0;
 481     }
 482   else if (fraction == 0)
 483     {
 484       exp = 0;
 485     }
 486   else
 487     {
 488       if (src->normal_exp < NORMAL_EXPMIN)
 489         {
 490           /* This number's exponent is too low to fit into the bits
 491              available in the number, so we'll store 0 in the exponent and
 492              shift the fraction to the right to make up for it.  */
 493
 494           int shift = NORMAL_EXPMIN - src->normal_exp;
 495
 496           exp = 0;
 497
 498           if (shift > FRAC_NBITS - NGARDS)
 499             {
 500               /* No point shifting, since it's more that 64 out.  */
 501               fraction = 0;
 502             }
 503           else
 504             {
 505               int lowbit = (fraction & ((1 << shift) - 1)) ? 1 : 0;
 506               fraction = (fraction >> shift) | lowbit;
 507             }
 508           if ((fraction & GARDMASK) == GARDMSB)
 509             {
 510               if ((fraction & (1 << NGARDS)))
 511                 fraction += GARDROUND + 1;
 512             }
 513           else
 514             {
 515               /* Add to the guards to round up.  */
 516               fraction += GARDROUND;
 517             }
 518           /* Perhaps the rounding means we now need to change the
 519              exponent.  */
 520           if (fraction >= IMPLICIT_2)
 521             {
 522               fraction >>= 1;
 523               exp += 1;
 524             }
 525           fraction >>= NGARDS;
 526         }
 527       else if (src->normal_exp > EXPBIAS)
 528         {
 529           exp = EXPMAX;
 530           fraction = 0;
 531         }
 532       else
 533         {
 534           exp = src->normal_exp + EXPBIAS;
 535           /* IF the gard bits are the all zero, but the first, then we're
 536              half way between two numbers, choose the one which makes the
 537              lsb of the answer 0.  */
 538           if ((fraction & GARDMASK) == GARDMSB)
 539             {
 540               if (fraction & (1 << NGARDS))
 541                 fraction += GARDROUND + 1;
 542             }
 543           else
 544             {
 545               /* Add a one to the guards to round up */
 546               fraction += GARDROUND;
 547             }
 548           if (fraction >= IMPLICIT_2)
 549             {
 550               fraction >>= 1;
 551               exp += 1;
 552             }
 553           fraction >>= NGARDS;
 554         }
 555     }
 556
 557   /* We previously used bitfields to store the number, but this doesn't
 558      handle little/big endian systems conveniently, so use shifts and
 559      masks */
 560 #ifdef FLOAT_BIT_ORDER_MISMATCH
 561   dst.bits.fraction = fraction;
 562   dst.bits.exp = exp;
 563   dst.bits.sign = sign;
 564 #else
 565   dst.value_raw = fraction & ((((fractype)1) << FRACBITS) - (fractype)1);
 566   dst.value_raw |= ((fractype) (exp & ((1 << EXPBITS) - 1))) << FRACBITS;
 567   dst.value_raw |= ((fractype) (sign & 1)) << (FRACBITS | EXPBITS);
 568 #endif
 569
 570 #if defined(FLOAT_WORD_ORDER_MISMATCH) && !defined(FLOAT)
 571   {
 572     halffractype tmp = dst.words[0];
 573     dst.words[0] = dst.words[1];
 574     dst.words[1] = tmp;
 575   }
 576 #endif
 577
 578   return dst.value;
 579 }
 580 #endif
 581
 582 extern void unpack_d (FLO_union_type *, fp_number_type *);
 583
 584 #if defined(L_unpack_df) || defined(L_unpack_sf)
 585 void
 586 unpack_d (FLO_union_type * src, fp_number_type * dst)
 587 {
 588   /* We previously used bitfields to store the number, but this doesn't
 589      handle little/big endian systems conveniently, so use shifts and
 590      masks */
 591   fractype fraction;
 592   int exp;
 593   int sign;
 594
 595 #if defined(FLOAT_WORD_ORDER_MISMATCH) && !defined(FLOAT)
 596   FLO_union_type swapped;
 597
 598   swapped.words[0] = src->words[1];
 599   swapped.words[1] = src->words[0];
 600   src = &swapped;
 601 #endif
 602
 603 #ifdef FLOAT_BIT_ORDER_MISMATCH
 604   fraction = src->bits.fraction;
 605   exp = src->bits.exp;
 606   sign = src->bits.sign;
 607 #else
 608   fraction = src->value_raw & ((((fractype)1) << FRACBITS) - (fractype)1);
 609   exp = ((int)(src->value_raw >> FRACBITS)) & ((1 << EXPBITS) - 1);
 610   sign = ((int)(src->value_raw >> (FRACBITS + EXPBITS))) & 1;
 611 #endif
 612
 613   dst->sign = sign;
 614   if (exp == 0)
 615     {
 616       /* Hmm.  Looks like 0 */
 617       if (fraction == 0)
 618         {
 619           /* tastes like zero */
 620           dst->class = CLASS_ZERO;
 621         }
 622       else
 623         {
 624           /* Zero exponent with non zero fraction - it's denormalized,
 625              so there isn't a leading implicit one - we'll shift it so
 626              it gets one.  */
 627           dst->normal_exp = exp - EXPBIAS + 1;
 628           fraction <<= NGARDS;
 629
 630           dst->class = CLASS_NUMBER;
 631 #if 1
 632           while (fraction < IMPLICIT_1)
 633             {
 634               fraction <<= 1;
 635               dst->normal_exp--;
 636             }
 637 #endif
 638           dst->fraction.ll = fraction;
 639         }
 640     }
 641   else if (exp == EXPMAX)
 642     {
 643       /* Huge exponent*/
 644       if (fraction == 0)
 645         {
 646           /* Attached to a zero fraction - means infinity */
 647           dst->class = CLASS_INFINITY;
 648         }
 649       else
 650         {
 651           /* Non zero fraction, means nan */
 652           if (fraction & QUIET_NAN)
 653             {
 654               dst->class = CLASS_QNAN;
 655             }
 656           else
 657             {
 658               dst->class = CLASS_SNAN;
 659             }
 660           /* Keep the fraction part as the nan number */
 661           dst->fraction.ll = fraction;
 662         }
 663     }
 664   else
 665     {
 666       /* Nothing strange about this number */
 667       dst->normal_exp = exp - EXPBIAS;
 668       dst->class = CLASS_NUMBER;
 669       dst->fraction.ll = (fraction << NGARDS) | IMPLICIT_1;
 670     }
 671 }
 672 #endif
 673
 674 #if defined(L_addsub_sf) || defined(L_addsub_df)
 675 static fp_number_type *
 676 _fpadd_parts (fp_number_type * a,
 677               fp_number_type * b,
 678               fp_number_type * tmp)
 679 {
 680   intfrac tfraction;
 681
 682   /* Put commonly used fields in local variables.  */
 683   int a_normal_exp;
 684   int b_normal_exp;
 685   fractype a_fraction;
 686   fractype b_fraction;
 687
 688   if (isnan (a))
 689     {
 690       return a;
 691     }
 692   if (isnan (b))
 693     {
 694       return b;
 695     }
 696   if (isinf (a))
 697     {
 698       /* Adding infinities with opposite signs yields a NaN.  */
 699       if (isinf (b) && a->sign != b->sign)
 700         return nan ();
 701       return a;
 702     }
 703   if (isinf (b))
 704     {
 705       return b;
 706     }
 707   if (iszero (b))
 708     {
 709       if (iszero (a))
 710         {
 711           *tmp = *a;
 712           tmp->sign = a->sign & b->sign;
 713           return tmp;
 714         }
 715       return a;
 716     }
 717   if (iszero (a))
 718     {
 719       return b;
 720     }
 721
 722   /* Got two numbers. shift the smaller and increment the exponent till
 723      they're the same */
 724   {
 725     int diff;
 726
 727     a_normal_exp = a->normal_exp;
 728     b_normal_exp = b->normal_exp;
 729     a_fraction = a->fraction.ll;
 730     b_fraction = b->fraction.ll;
 731
 732     diff = a_normal_exp - b_normal_exp;
 733
 734     if (diff < 0)
 735       diff = -diff;
 736     if (diff < FRAC_NBITS)
 737       {
 738         /* ??? This does shifts one bit at a time.  Optimize.  */
 739         while (a_normal_exp > b_normal_exp)
 740           {
 741             b_normal_exp++;
 742             LSHIFT (b_fraction);
 743           }
 744         while (b_normal_exp > a_normal_exp)
 745           {
 746             a_normal_exp++;
 747             LSHIFT (a_fraction);
 748           }
 749       }
 750     else
 751       {
 752         /* Somethings's up.. choose the biggest */
 753         if (a_normal_exp > b_normal_exp)
 754           {
 755             b_normal_exp = a_normal_exp;
 756             b_fraction = 0;
 757           }
 758         else
 759           {
 760             a_normal_exp = b_normal_exp;
 761             a_fraction = 0;
 762           }
 763       }
 764   }
 765
 766   if (a->sign != b->sign)
 767     {
 768       if (a->sign)
 769         {
 770           tfraction = -a_fraction + b_fraction;
 771         }
 772       else
 773         {
 774           tfraction = a_fraction - b_fraction;
 775         }
 776       if (tfraction >= 0)
 777         {
 778           tmp->sign = 0;
 779           tmp->normal_exp = a_normal_exp;
 780           tmp->fraction.ll = tfraction;
 781         }
 782       else
 783         {
 784           tmp->sign = 1;
 785           tmp->normal_exp = a_normal_exp;
 786           tmp->fraction.ll = -tfraction;
 787         }
 788       /* and renormalize it */
 789
 790       while (tmp->fraction.ll < IMPLICIT_1 && tmp->fraction.ll)
 791         {
 792           tmp->fraction.ll <<= 1;
 793           tmp->normal_exp--;
 794         }
 795     }
 796   else
 797     {
 798       tmp->sign = a->sign;
 799       tmp->normal_exp = a_normal_exp;
 800       tmp->fraction.ll = a_fraction + b_fraction;
 801     }
 802   tmp->class = CLASS_NUMBER;
 803   /* Now the fraction is added, we have to shift down to renormalize the
 804      number */
 805
 806   if (tmp->fraction.ll >= IMPLICIT_2)
 807     {
 808       LSHIFT (tmp->fraction.ll);
 809       tmp->normal_exp++;
 810     }
 811   return tmp;
 812
 813 }
 814
 815 FLO_type
 816 add (FLO_type arg_a, FLO_type arg_b)
 817 {
 818   fp_number_type a;
 819   fp_number_type b;
 820   fp_number_type tmp;
 821   fp_number_type *res;
 822   FLO_union_type au, bu;
 823
 824   au.value = arg_a;
 825   bu.value = arg_b;
 826
 827   unpack_d (&au, &a);
 828   unpack_d (&bu, &b);
 829
 830   res = _fpadd_parts (&a, &b, &tmp);
 831
 832   return pack_d (res);
 833 }
 834
 835 FLO_type
 836 sub (FLO_type arg_a, FLO_type arg_b)
 837 {
 838   fp_number_type a;
 839   fp_number_type b;
 840   fp_number_type tmp;
 841   fp_number_type *res;
 842   FLO_union_type au, bu;
 843
 844   au.value = arg_a;
 845   bu.value = arg_b;
 846
 847   unpack_d (&au, &a);
 848   unpack_d (&bu, &b);
 849
 850   b.sign ^= 1;
 851
 852   res = _fpadd_parts (&a, &b, &tmp);
 853
 854   return pack_d (res);
 855 }
 856 #endif
 857
 858 #if defined(L_mul_sf) || defined(L_mul_df)
 859 static INLINE fp_number_type *
 860 _fpmul_parts ( fp_number_type *  a,
 861                fp_number_type *  b,
 862                fp_number_type * tmp)
 863 {
 864   fractype low = 0;
 865   fractype high = 0;
 866
 867   if (isnan (a))
 868     {
 869       a->sign = a->sign != b->sign;
 870       return a;
 871     }
 872   if (isnan (b))
 873     {
 874       b->sign = a->sign != b->sign;
 875       return b;
 876     }
 877   if (isinf (a))
 878     {
 879       if (iszero (b))
 880         return nan ();
 881       a->sign = a->sign != b->sign;
 882       return a;
 883     }
 884   if (isinf (b))
 885     {
 886       if (iszero (a))
 887         {
 888           return nan ();
 889         }
 890       b->sign = a->sign != b->sign;
 891       return b;
 892     }
 893   if (iszero (a))
 894     {
 895       a->sign = a->sign != b->sign;
 896       return a;
 897     }
 898   if (iszero (b))
 899     {
 900       b->sign = a->sign != b->sign;
 901       return b;
 902     }
 903
 904   /* Calculate the mantissa by multiplying both 64bit numbers to get a
 905      128 bit number */
 906   {
 907 #if defined(NO_DI_MODE)
 908     {
 909       fractype x = a->fraction.ll;
 910       fractype ylow = b->fraction.ll;
 911       fractype yhigh = 0;
 912       int bit;
 913
 914       /* ??? This does multiplies one bit at a time.  Optimize.  */
 915       for (bit = 0; bit < FRAC_NBITS; bit++)
 916         {
 917           int carry;
 918
 919           if (x & 1)
 920             {
 921               carry = (low += ylow) < ylow;
 922               high += yhigh + carry;
 923             }
 924           yhigh <<= 1;
 925           if (ylow & FRACHIGH)
 926             {
 927               yhigh |= 1;
 928             }
 929           ylow <<= 1;
 930           x >>= 1;
 931         }
 932     }
 933 #elif defined(FLOAT)
 934     {
 935       /* Multiplying two 32 bit numbers to get a 64 bit number  on
 936         a machine with DI, so we're safe */
 937
 938       DItype answer = (DItype)(a->fraction.ll) * (DItype)(b->fraction.ll);
 939
 940       high = answer >> 32;
 941       low = answer;
 942     }
 943 #else
 944     /* Doing a 64*64 to 128 */
 945     {
 946       UDItype nl = a->fraction.ll & 0xffffffff;
 947       UDItype nh = a->fraction.ll >> 32;
 948       UDItype ml = b->fraction.ll & 0xffffffff;
 949       UDItype mh = b->fraction.ll >>32;
 950       UDItype pp_ll = ml * nl;
 951       UDItype pp_hl = mh * nl;
 952       UDItype pp_lh = ml * nh;
 953       UDItype pp_hh = mh * nh;
 954       UDItype res2 = 0;
 955       UDItype res0 = 0;
 956       UDItype ps_hh__ = pp_hl + pp_lh;
 957       if (ps_hh__ < pp_hl)
 958         res2 += 0x100000000LL;
 959       pp_hl = (ps_hh__ << 32) & 0xffffffff00000000LL;
 960       res0 = pp_ll + pp_hl;
 961       if (res0 < pp_ll)
 962         res2++;
 963       res2 += ((ps_hh__ >> 32) & 0xffffffffL) + pp_hh;
 964       high = res2;
 965       low = res0;
 966     }
 967 #endif
 968   }
 969
 970   tmp->normal_exp = a->normal_exp + b->normal_exp;
 971   tmp->sign = a->sign != b->sign;
 972 #ifdef FLOAT
 973   tmp->normal_exp += 2;         /* ??????????????? */
 974 #else
 975   tmp->normal_exp += 4;         /* ??????????????? */
 976 #endif
 977   while (high >= IMPLICIT_2)
 978     {
 979       tmp->normal_exp++;
 980       if (high & 1)
 981         {
 982           low >>= 1;
 983           low |= FRACHIGH;
 984         }
 985       high >>= 1;
 986     }
 987   while (high < IMPLICIT_1)
 988     {
 989       tmp->normal_exp--;
 990
 991       high <<= 1;
 992       if (low & FRACHIGH)
 993         high |= 1;
 994       low <<= 1;
 995     }
 996   /* rounding is tricky. if we only round if it won't make us round later. */
 997 #if 0
 998   if (low & FRACHIGH2)
 999     {
1000       if (((high & GARDMASK) != GARDMSB)
1001           && (((high + 1) & GARDMASK) == GARDMSB))
1002         {
1003           /* don't round, it gets done again later. */
1004         }
1005       else
1006         {
1007           high++;
1008         }
1009     }
1010 #endif
1011   if ((high & GARDMASK) == GARDMSB)
1012     {
1013       if (high & (1 << NGARDS))
1014         {
1015           /* half way, so round to even */
1016           high += GARDROUND + 1;
1017         }
1018       else if (low)
1019         {
1020           /* but we really weren't half way */
1021           high += GARDROUND + 1;
1022         }
1023     }
1024   tmp->fraction.ll = high;
1025   tmp->class = CLASS_NUMBER;
1026   return tmp;
1027 }
1028
1029 FLO_type
1030 multiply (FLO_type arg_a, FLO_type arg_b)
1031 {
1032   fp_number_type a;
1033   fp_number_type b;
1034   fp_number_type tmp;
1035   fp_number_type *res;
1036   FLO_union_type au, bu;
1037
1038   au.value = arg_a;
1039   bu.value = arg_b;
1040
1041   unpack_d (&au, &a);
1042   unpack_d (&bu, &b);
1043
1044   res = _fpmul_parts (&a, &b, &tmp);
1045
1046   return pack_d (res);
1047 }
1048 #endif
1049
1050 #if defined(L_div_sf) || defined(L_div_df)
1051 static INLINE fp_number_type *
1052 _fpdiv_parts (fp_number_type * a,
1053               fp_number_type * b)
1054 {
1055   fractype bit;
1056   fractype numerator;
1057   fractype denominator;
1058   fractype quotient;
1059
1060   if (isnan (a))
1061     {
1062       return a;
1063     }
1064   if (isnan (b))
1065     {
1066       return b;
1067     }
1068
1069   a->sign = a->sign ^ b->sign;
1070
1071   if (isinf (a) || iszero (a))
1072     {
1073       if (a->class == b->class)
1074         return nan ();
1075       return a;
1076     }
1077
1078   if (isinf (b))
1079     {
1080       a->fraction.ll = 0;
1081       a->normal_exp = 0;
1082       return a;
1083     }
1084   if (iszero (b))
1085     {
1086       a->class = CLASS_INFINITY;
1087       return a;
1088     }
1089
1090   /* Calculate the mantissa by multiplying both 64bit numbers to get a
1091      128 bit number */
1092   {
1093     /* quotient =
1094        ( numerator / denominator) * 2^(numerator exponent -  denominator exponent)
1095      */
1096
1097     a->normal_exp = a->normal_exp - b->normal_exp;
1098     numerator = a->fraction.ll;
1099     denominator = b->fraction.ll;
1100
1101     if (numerator < denominator)
1102       {
1103         /* Fraction will be less than 1.0 */
1104         numerator *= 2;
1105         a->normal_exp--;
1106       }
1107     bit = IMPLICIT_1;
1108     quotient = 0;
1109     /* ??? Does divide one bit at a time.  Optimize.  */
1110     while (bit)
1111       {
1112         if (numerator >= denominator)
1113           {
1114             quotient |= bit;
1115             numerator -= denominator;
1116           }
1117         bit >>= 1;
1118         numerator *= 2;
1119       }
1120
1121     if ((quotient & GARDMASK) == GARDMSB)
1122       {
1123         if (quotient & (1 << NGARDS))
1124           {
1125             /* half way, so round to even */
1126             quotient += GARDROUND + 1;
1127           }
1128         else if (numerator)
1129           {
1130             /* but we really weren't half way, more bits exist */
1131             quotient += GARDROUND + 1;
1132           }
1133       }
1134
1135     a->fraction.ll = quotient;
1136     return (a);
1137   }
1138 }
1139
1140 FLO_type
1141 divide (FLO_type arg_a, FLO_type arg_b)
1142 {
1143   fp_number_type a;
1144   fp_number_type b;
1145   fp_number_type *res;
1146   FLO_union_type au, bu;
1147
1148   au.value = arg_a;
1149   bu.value = arg_b;
1150
1151   unpack_d (&au, &a);
1152   unpack_d (&bu, &b);
1153
1154   res = _fpdiv_parts (&a, &b);
1155
1156   return pack_d (res);
1157 }
1158 #endif
1159
1160 int __fpcmp_parts (fp_number_type * a, fp_number_type *b);
1161
1162 #if defined(L_fpcmp_parts_sf) || defined(L_fpcmp_parts_df)
1163 /* according to the demo, fpcmp returns a comparison with 0... thus
1164    a<b -> -1
1165    a==b -> 0
1166    a>b -> +1
1167  */
1168
1169 int
1170 __fpcmp_parts (fp_number_type * a, fp_number_type * b)
1171 {
1172 #if 0
1173   /* either nan -> unordered. Must be checked outside of this routine. */
1174   if (isnan (a) && isnan (b))
1175     {
1176       return 1;                 /* still unordered! */
1177     }
1178 #endif
1179
1180   if (isnan (a) || isnan (b))
1181     {
1182       return 1;                 /* how to indicate unordered compare? */
1183     }
1184   if (isinf (a) && isinf (b))
1185     {
1186       /* +inf > -inf, but +inf != +inf */
1187       /* b    \a| +inf(0)| -inf(1)
1188        ______\+--------+--------
1189        +inf(0)| a==b(0)| a<b(-1)
1190        -------+--------+--------
1191        -inf(1)| a>b(1) | a==b(0)
1192        -------+--------+--------
1193        So since unordered must be non zero, just line up the columns...
1194        */
1195       return b->sign - a->sign;
1196     }
1197   /* but not both... */
1198   if (isinf (a))
1199     {
1200       return a->sign ? -1 : 1;
1201     }
1202   if (isinf (b))
1203     {
1204       return b->sign ? 1 : -1;
1205     }
1206   if (iszero (a) && iszero (b))
1207     {
1208       return 0;
1209     }
1210   if (iszero (a))
1211     {
1212       return b->sign ? 1 : -1;
1213     }
1214   if (iszero (b))
1215     {
1216       return a->sign ? -1 : 1;
1217     }
1218   /* now both are "normal". */
1219   if (a->sign != b->sign)
1220     {
1221       /* opposite signs */
1222       return a->sign ? -1 : 1;
1223     }
1224   /* same sign; exponents? */
1225   if (a->normal_exp > b->normal_exp)
1226     {
1227       return a->sign ? -1 : 1;
1228     }
1229   if (a->normal_exp < b->normal_exp)
1230     {
1231       return a->sign ? 1 : -1;
1232     }
1233   /* same exponents; check size. */
1234   if (a->fraction.ll > b->fraction.ll)
1235     {
1236       return a->sign ? -1 : 1;
1237     }
1238   if (a->fraction.ll < b->fraction.ll)
1239     {
1240       return a->sign ? 1 : -1;
1241     }
1242   /* after all that, they're equal. */
1243   return 0;
1244 }
1245 #endif
1246
1247 #if defined(L_compare_sf) || defined(L_compare_df)
1248 CMPtype
1249 compare (FLO_type arg_a, FLO_type arg_b)
1250 {
1251   fp_number_type a;
1252   fp_number_type b;
1253   FLO_union_type au, bu;
1254
1255   au.value = arg_a;
1256   bu.value = arg_b;
1257
1258   unpack_d (&au, &a);
1259   unpack_d (&bu, &b);
1260
1261   return __fpcmp_parts (&a, &b);
1262 }
1263 #endif
1264
1265 #ifndef US_SOFTWARE_GOFAST
1266
1267 /* These should be optimized for their specific tasks someday.  */
1268
1269 #if defined(L_eq_sf) || defined(L_eq_df)
1270 CMPtype
1271 _eq_f2 (FLO_type arg_a, FLO_type arg_b)
1272 {
1273   fp_number_type a;
1274   fp_number_type b;
1275   FLO_union_type au, bu;
1276
1277   au.value = arg_a;
1278   bu.value = arg_b;
1279
1280   unpack_d (&au, &a);
1281   unpack_d (&bu, &b);
1282
1283   if (isnan (&a) || isnan (&b))
1284     return 1;                   /* false, truth == 0 */
1285
1286   return __fpcmp_parts (&a, &b) ;
1287 }
1288 #endif
1289
1290 #if defined(L_ne_sf) || defined(L_ne_df)
1291 CMPtype
1292 _ne_f2 (FLO_type arg_a, FLO_type arg_b)
1293 {
1294   fp_number_type a;
1295   fp_number_type b;
1296   FLO_union_type au, bu;
1297
1298   au.value = arg_a;
1299   bu.value = arg_b;
1300
1301   unpack_d (&au, &a);
1302   unpack_d (&bu, &b);
1303
1304   if (isnan (&a) || isnan (&b))
1305     return 1;                   /* true, truth != 0 */
1306
1307   return  __fpcmp_parts (&a, &b) ;
1308 }
1309 #endif
1310
1311 #if defined(L_gt_sf) || defined(L_gt_df)
1312 CMPtype
1313 _gt_f2 (FLO_type arg_a, FLO_type arg_b)
1314 {
1315   fp_number_type a;
1316   fp_number_type b;
1317   FLO_union_type au, bu;
1318
1319   au.value = arg_a;
1320   bu.value = arg_b;
1321
1322   unpack_d (&au, &a);
1323   unpack_d (&bu, &b);
1324
1325   if (isnan (&a) || isnan (&b))
1326     return -1;                  /* false, truth > 0 */
1327
1328   return __fpcmp_parts (&a, &b);
1329 }
1330 #endif
1331
1332 #if defined(L_ge_sf) || defined(L_ge_df)
1333 CMPtype
1334 _ge_f2 (FLO_type arg_a, FLO_type arg_b)
1335 {
1336   fp_number_type a;
1337   fp_number_type b;
1338   FLO_union_type au, bu;
1339
1340   au.value = arg_a;
1341   bu.value = arg_b;
1342
1343   unpack_d (&au, &a);
1344   unpack_d (&bu, &b);
1345
1346   if (isnan (&a) || isnan (&b))
1347     return -1;                  /* false, truth >= 0 */
1348   return __fpcmp_parts (&a, &b) ;
1349 }
1350 #endif
1351
1352 #if defined(L_lt_sf) || defined(L_lt_df)
1353 CMPtype
1354 _lt_f2 (FLO_type arg_a, FLO_type arg_b)
1355 {
1356   fp_number_type a;
1357   fp_number_type b;
1358   FLO_union_type au, bu;
1359
1360   au.value = arg_a;
1361   bu.value = arg_b;
1362
1363   unpack_d (&au, &a);
1364   unpack_d (&bu, &b);
1365
1366   if (isnan (&a) || isnan (&b))
1367     return 1;                   /* false, truth < 0 */
1368
1369   return __fpcmp_parts (&a, &b);
1370 }
1371 #endif
1372
1373 #if defined(L_le_sf) || defined(L_le_df)
1374 CMPtype
1375 _le_f2 (FLO_type arg_a, FLO_type arg_b)
1376 {
1377   fp_number_type a;
1378   fp_number_type b;
1379   FLO_union_type au, bu;
1380
1381   au.value = arg_a;
1382   bu.value = arg_b;
1383
1384   unpack_d (&au, &a);
1385   unpack_d (&bu, &b);
1386
1387   if (isnan (&a) || isnan (&b))
1388     return 1;                   /* false, truth <= 0 */
1389
1390   return __fpcmp_parts (&a, &b) ;
1391 }
1392 #endif
1393
1394 #if defined(L_unord_sf) || defined(L_unord_df)
1395 CMPtype
1396 _unord_f2 (FLO_type arg_a, FLO_type arg_b)
1397 {
1398   fp_number_type a;
1399   fp_number_type b;
1400   FLO_union_type au, bu;
1401
1402   au.value = arg_a;
1403   bu.value = arg_b;
1404
1405   unpack_d (&au, &a);
1406   unpack_d (&bu, &b);
1407
1408   return (isnan (&a) || isnan (&b));
1409 }
1410 #endif
1411
1412 #endif /* ! US_SOFTWARE_GOFAST */
1413
1414 #if defined(L_si_to_sf) || defined(L_si_to_df)
1415 FLO_type
1416 si_to_float (SItype arg_a)
1417 {
1418   fp_number_type in;
1419
1420   in.class = CLASS_NUMBER;
1421   in.sign = arg_a < 0;
1422   if (!arg_a)
1423     {
1424       in.class = CLASS_ZERO;
1425     }
1426   else
1427     {
1428       in.normal_exp = FRACBITS + NGARDS;
1429       if (in.sign)
1430         {
1431           /* Special case for minint, since there is no +ve integer
1432              representation for it */
1433           if (arg_a == (SItype) 0x80000000)
1434             {
1435               return -2147483648.0;
1436             }
1437           in.fraction.ll = (-arg_a);
1438         }
1439       else
1440         in.fraction.ll = arg_a;
1441
1442       while (in.fraction.ll < (1LL << (FRACBITS + NGARDS)))
1443         {
1444           in.fraction.ll <<= 1;
1445           in.normal_exp -= 1;
1446         }
1447     }
1448   return pack_d (&in);
1449 }
1450 #endif
1451
1452 #if defined(L_sf_to_si) || defined(L_df_to_si)
1453 SItype
1454 float_to_si (FLO_type arg_a)
1455 {
1456   fp_number_type a;
1457   SItype tmp;
1458   FLO_union_type au;
1459
1460   au.value = arg_a;
1461   unpack_d (&au, &a);
1462
1463   if (iszero (&a))
1464     return 0;
1465   if (isnan (&a))
1466     return 0;
1467   /* get reasonable MAX_SI_INT... */
1468   if (isinf (&a))
1469     return a.sign ? (-MAX_SI_INT)-1 : MAX_SI_INT;
1470   /* it is a number, but a small one */
1471   if (a.normal_exp < 0)
1472     return 0;
1473   if (a.normal_exp > 30)
1474     return a.sign ? (-MAX_SI_INT)-1 : MAX_SI_INT;
1475   tmp = a.fraction.ll >> ((FRACBITS + NGARDS) - a.normal_exp);
1476   return a.sign ? (-tmp) : (tmp);
1477 }
1478 #endif
1479
1480 #if defined(L_sf_to_usi) || defined(L_df_to_usi)
1481 #ifdef US_SOFTWARE_GOFAST
1482 /* While libgcc2.c defines its own __fixunssfsi and __fixunsdfsi routines,
1483    we also define them for GOFAST because the ones in libgcc2.c have the
1484    wrong names and I'd rather define these here and keep GOFAST CYG-LOC's
1485    out of libgcc2.c.  We can't define these here if not GOFAST because then
1486    there'd be duplicate copies.  */
1487
1488 USItype
1489 float_to_usi (FLO_type arg_a)
1490 {
1491   fp_number_type a;
1492   FLO_union_type au;
1493
1494   au.value = arg_a;
1495   unpack_d (&au, &a);
1496
1497   if (iszero (&a))
1498     return 0;
1499   if (isnan (&a))
1500     return 0;
1501   /* it is a negative number */
1502   if (a.sign)
1503     return 0;
1504   /* get reasonable MAX_USI_INT... */
1505   if (isinf (&a))
1506     return MAX_USI_INT;
1507   /* it is a number, but a small one */
1508   if (a.normal_exp < 0)
1509     return 0;
1510   if (a.normal_exp > 31)
1511     return MAX_USI_INT;
1512   else if (a.normal_exp > (FRACBITS + NGARDS))
1513     return a.fraction.ll << (a.normal_exp - (FRACBITS + NGARDS));
1514   else
1515     return a.fraction.ll >> ((FRACBITS + NGARDS) - a.normal_exp);
1516 }
1517 #endif
1518 #endif
1519
1520 #if defined(L_negate_sf) || defined(L_negate_df)
1521 FLO_type
1522 negate (FLO_type arg_a)
1523 {
1524   fp_number_type a;
1525   FLO_union_type au;
1526
1527   au.value = arg_a;
1528   unpack_d (&au, &a);
1529
1530   flip_sign (&a);
1531   return pack_d (&a);
1532 }
1533 #endif
1534
1535 #ifdef FLOAT
1536
1537 #if defined(L_make_sf)
1538 SFtype
1539 __make_fp(fp_class_type class,
1540              unsigned int sign,
1541              int exp,
1542              USItype frac)
1543 {
1544   fp_number_type in;
1545
1546   in.class = class;
1547   in.sign = sign;
1548   in.normal_exp = exp;
1549   in.fraction.ll = frac;
1550   return pack_d (&in);
1551 }
1552 #endif
1553
1554 #ifndef FLOAT_ONLY
1555
1556 /* This enables one to build an fp library that supports float but not double.
1557    Otherwise, we would get an undefined reference to __make_dp.
1558    This is needed for some 8-bit ports that can't handle well values that
1559    are 8-bytes in size, so we just don't support double for them at all.  */
1560
1561 extern DFtype __make_dp (fp_class_type, unsigned int, int, UDItype frac);
1562
1563 #if defined(L_sf_to_df)
1564 DFtype
1565 sf_to_df (SFtype arg_a)
1566 {
1567   fp_number_type in;
1568   FLO_union_type au;
1569
1570   au.value = arg_a;
1571   unpack_d (&au, &in);
1572
1573   return __make_dp (in.class, in.sign, in.normal_exp,
1574                     ((UDItype) in.fraction.ll) << F_D_BITOFF);
1575 }
1576 #endif
1577
1578 #endif
1579 #endif
1580
1581 #ifndef FLOAT
1582
1583 extern SFtype __make_fp (fp_class_type, unsigned int, int, USItype);
1584
1585 #if defined(L_make_df)
1586 DFtype
1587 __make_dp (fp_class_type class, unsigned int sign, int exp, UDItype frac)
1588 {
1589   fp_number_type in;
1590
1591   in.class = class;
1592   in.sign = sign;
1593   in.normal_exp = exp;
1594   in.fraction.ll = frac;
1595   return pack_d (&in);
1596 }
1597 #endif
1598
1599 #if defined(L_df_to_sf)
1600 SFtype
1601 df_to_sf (DFtype arg_a)
1602 {
1603   fp_number_type in;
1604   USItype sffrac;
1605   FLO_union_type au;
1606
1607   au.value = arg_a;
1608   unpack_d (&au, &in);
1609
1610   sffrac = in.fraction.ll >> F_D_BITOFF;
1611
1612   /* We set the lowest guard bit in SFFRAC if we discarded any non
1613      zero bits.  */
1614   if ((in.fraction.ll & (((USItype) 1 << F_D_BITOFF) - 1)) != 0)
1615     sffrac |= 1;
1616
1617   return __make_fp (in.class, in.sign, in.normal_exp, sffrac);
1618 }
1619 #endif
1620
1621 #endif
1622 #endif /* !EXTENDED_FLOAT_STUBS */