OSDN Git Service

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