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