1/* This is a software floating point library which can be used
2   for targets without hardware floating point.
3   Copyright (C) 1994-2015 Free Software Foundation, Inc.
4
5This file is part of GCC.
6
7GCC is free software; you can redistribute it and/or modify it under
8the terms of the GNU General Public License as published by the Free
9Software Foundation; either version 3, or (at your option) any later
10version.
11
12GCC is distributed in the hope that it will be useful, but WITHOUT ANY
13WARRANTY; without even the implied warranty of MERCHANTABILITY or
14FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
15for more details.
16
17Under Section 7 of GPL version 3, you are granted additional
18permissions described in the GCC Runtime Library Exception, version
193.1, as published by the Free Software Foundation.
20
21You should have received a copy of the GNU General Public License and
22a copy of the GCC Runtime Library Exception along with this program;
23see 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
73extern void abort (void);
74void __extendsfxf2 (void) { abort(); }
75void __extenddfxf2 (void) { abort(); }
76void __truncxfdf2 (void) { abort(); }
77void __truncxfsf2 (void) { abort(); }
78void __fixxfsi (void) { abort(); }
79void __floatsixf (void) { abort(); }
80void __addxf3 (void) { abort(); }
81void __subxf3 (void) { abort(); }
82void __mulxf3 (void) { abort(); }
83void __divxf3 (void) { abort(); }
84void __negxf2 (void) { abort(); }
85void __eqxf2 (void) { abort(); }
86void __nexf2 (void) { abort(); }
87void __gtxf2 (void) { abort(); }
88void __gexf2 (void) { abort(); }
89void __lexf2 (void) { abort(); }
90void __ltxf2 (void) { abort(); }
91
92void __extendsftf2 (void) { abort(); }
93void __extenddftf2 (void) { abort(); }
94void __trunctfdf2 (void) { abort(); }
95void __trunctfsf2 (void) { abort(); }
96void __fixtfsi (void) { abort(); }
97void __floatsitf (void) { abort(); }
98void __addtf3 (void) { abort(); }
99void __subtf3 (void) { abort(); }
100void __multf3 (void) { abort(); }
101void __divtf3 (void) { abort(); }
102void __negtf2 (void) { abort(); }
103void __eqtf2 (void) { abort(); }
104void __netf2 (void) { abort(); }
105void __gttf2 (void) { abort(); }
106void __getf2 (void) { abort(); }
107void __letf2 (void) { abort(); }
108void __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
121const fp_number_type __thenan_sf = { CLASS_SNAN, 0, 0, {(fractype) 0} };
122#elif defined L_thenan_df
123const fp_number_type __thenan_df = { CLASS_SNAN, 0, 0, {(fractype) 0} };
124#elif defined L_thenan_tf
125const fp_number_type __thenan_tf = { CLASS_SNAN, 0, 0, {(fractype) 0} };
126#elif defined TFLOAT
127extern const fp_number_type __thenan_tf;
128#elif defined FLOAT
129extern const fp_number_type __thenan_sf;
130#else
131extern const fp_number_type __thenan_df;
132#endif
133
134INLINE
135static const fp_number_type *
136makenan (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
147INLINE
148static int
149isnan (const fp_number_type *x)
150{
151  return __builtin_expect (x->class == CLASS_SNAN || x->class == CLASS_QNAN,
152			   0);
153}
154
155INLINE
156static int
157isinf (const fp_number_type *  x)
158{
159  return __builtin_expect (x->class == CLASS_INFINITY, 0);
160}
161
162#endif /* NO_NANS */
163
164INLINE
165static int
166iszero (const fp_number_type *  x)
167{
168  return x->class == CLASS_ZERO;
169}
170
171INLINE
172static void
173flip_sign ( fp_number_type *  x)
174{
175  x->sign = !x->sign;
176}
177
178/* Count leading zeroes in N.  */
179INLINE
180static int
181clzusi (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
194extern FLO_type pack_d (const fp_number_type * );
195
196#if defined(L_pack_df) || defined(L_pack_sf) || defined(L_pack_tf)
197FLO_type
198pack_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)
428void
429unpack_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)
580static 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
720FLO_type
721add (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
740FLO_type
741sub (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)
764static 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
926FLO_type
927multiply (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)
948static 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
1045FLO_type
1046divide (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
1073int
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)
1152CMPtype
1153compare (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)
1172CMPtype
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)
1193CMPtype
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)
1214CMPtype
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)
1235CMPtype
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)
1255CMPtype
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)
1276CMPtype
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)
1297CMPtype
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)
1315FLO_type
1316si_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)
1357FLO_type
1358usi_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)
1393SItype
1394float_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)
1421USItype
1422float_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)
1453FLO_type
1454negate (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)
1470SFtype
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)
1494DFtype
1495sf_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)
1509TFtype
1510sf_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
1528extern SFtype __make_fp (fp_class_type, unsigned int, int, USItype);
1529
1530#if defined(L_make_df)
1531DFtype
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)
1545SFtype
1546df_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)
1568TFtype
1569df_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)
1584TFtype
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)
1601DFtype
1602tf_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)
1623SFtype
1624tf_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