gcc/config/fp-bit.c

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