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