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