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