+  swapped.words[0] = src->words[3];+  swapped.words[1] = src->words[2];+  swapped.words[2] = src->words[1];+  swapped.words[3] = src->words[0];+#else+  swapped.words[0] = src->words[1];+  swapped.words[1] = src->words[0];+#endif+  src = &swapped;+#endif+  +#ifdef FLOAT_BIT_ORDER_MISMATCH+  fraction = src->bits.fraction;+  exp = src->bits.exp;+  sign = src->bits.sign;+#else+# if defined TFLOAT && defined HALFFRACBITS+ {+   halffractype high, low;+   +   high = src->value_raw >> HALFSHIFT;+   low = src->value_raw & (((fractype)1 << HALFSHIFT) - 1);++   fraction = high & ((((fractype)1) << HALFFRACBITS) - 1);+   fraction <<= FRACBITS - HALFFRACBITS;+   exp = ((int)(high >> HALFFRACBITS)) & ((1 << EXPBITS) - 1);+   sign = ((int)(high >> (((HALFFRACBITS + EXPBITS))))) & 1;++   if (exp != EXPMAX && exp != 0 && low != 0)+     {+       int lowexp = ((int)(low >> HALFFRACBITS)) & ((1 << EXPBITS) - 1);+       int lowsign = ((int)(low >> (((HALFFRACBITS + EXPBITS))))) & 1;+       int shift;+       fractype xlow;++       xlow = low & ((((fractype)1) << HALFFRACBITS) - 1);+       if (lowexp)+	 xlow |= (((halffractype)1) << HALFFRACBITS);+       else+	 lowexp = 1;+       shift = (FRACBITS - HALFFRACBITS) - (exp - lowexp);+       if (shift > 0)+	 xlow <<= shift;+       else if (shift < 0)+	 xlow >>= -shift;+       if (sign == lowsign)+	 fraction += xlow;+       else if (fraction >= xlow)+	 fraction -= xlow;+       else+	 {+	   /* The high part is a power of two but the full number is lower.+	      This code will leave the implicit 1 in FRACTION, but we'd+	      have added that below anyway.  */+	   fraction = (((fractype) 1 << FRACBITS) - xlow) << 1;+	   exp--;+	 }+     }+ }+# else+  fraction = src->value_raw & ((((fractype)1) << FRACBITS) - 1);+  exp = ((int)(src->value_raw >> FRACBITS)) & ((1 << EXPBITS) - 1);+  sign = ((int)(src->value_raw >> (FRACBITS + EXPBITS))) & 1;+# endif+#endif++  dst->sign = sign;+  if (exp == 0)+    {+      /* Hmm.  Looks like 0 */+      if (fraction == 0+#ifdef NO_DENORMALS+	  || 1+#endif+	  )+	{+	  /* tastes like zero */+	  dst->class = CLASS_ZERO;+	}+      else+	{+	  /* Zero exponent with nonzero fraction - it's denormalized,+	     so there isn't a leading implicit one - we'll shift it so+	     it gets one.  */+	  dst->normal_exp = exp - EXPBIAS + 1;+	  fraction <<= NGARDS;++	  dst->class = CLASS_NUMBER;+#if 1+	  while (fraction < IMPLICIT_1)+	    {+	      fraction <<= 1;+	      dst->normal_exp--;+	    }+#endif+	  dst->fraction.ll = fraction;+	}+    }+  else if (!LARGEST_EXPONENT_IS_NORMAL (FRAC_NBITS) && exp == EXPMAX)+    {+      /* Huge exponent*/+      if (fraction == 0)+	{+	  /* Attached to a zero fraction - means infinity */+	  dst->class = CLASS_INFINITY;+	}+      else+	{+	  /* Nonzero fraction, means nan */+#ifdef QUIET_NAN_NEGATED+	  if ((fraction & QUIET_NAN) == 0)+#else+	  if (fraction & QUIET_NAN)+#endif+	    {+	      dst->class = CLASS_QNAN;+	    }+	  else+	    {+	      dst->class = CLASS_SNAN;+	    }+	  /* Keep the fraction part as the nan number */+	  dst->fraction.ll = fraction;+	}+    }+  else+    {+      /* Nothing strange about this number */+      dst->normal_exp = exp - EXPBIAS;+      dst->class = CLASS_NUMBER;+      dst->fraction.ll = (fraction << NGARDS) | IMPLICIT_1;+    }+}+#endif /* L_unpack_df || L_unpack_sf */++#if defined(L_addsub_sf) || defined(L_addsub_df) || defined(L_addsub_tf)+static fp_number_type *+_fpadd_parts (fp_number_type * a,+	      fp_number_type * b,+	      fp_number_type * tmp)+{+  intfrac tfraction;++  /* Put commonly used fields in local variables.  */+  int a_normal_exp;+  int b_normal_exp;+  fractype a_fraction;+  fractype b_fraction;++  if (isnan (a))+    {+      return a;+    }+  if (isnan (b))+    {+      return b;+    }+  if (isinf (a))+    {+      /* Adding infinities with opposite signs yields a NaN.  */+      if (isinf (b) && a->sign != b->sign)+	return nan ();+      return a;+    }+  if (isinf (b))+    {+      return b;+    }+  if (iszero (b))+    {+      if (iszero (a))+	{+	  *tmp = *a;+	  tmp->sign = a->sign & b->sign;+	  return tmp;+	}+      return a;+    }+  if (iszero (a))+    {+      return b;+    }++  /* Got two numbers. shift the smaller and increment the exponent till+     they're the same */+  {+    int diff;++    a_normal_exp = a->normal_exp;+    b_normal_exp = b->normal_exp;+    a_fraction = a->fraction.ll;+    b_fraction = b->fraction.ll;++    diff = a_normal_exp - b_normal_exp;++    if (diff < 0)+      diff = -diff;+    if (diff < FRAC_NBITS)+      {+	/* ??? This does shifts one bit at a time.  Optimize.  */+	while (a_normal_exp > b_normal_exp)+	  {+	    b_normal_exp++;+	    LSHIFT (b_fraction);+	  }+	while (b_normal_exp > a_normal_exp)+	  {+	    a_normal_exp++;+	    LSHIFT (a_fraction);+	  }+      }+    else+      {+	/* Somethings's up.. choose the biggest */+	if (a_normal_exp > b_normal_exp)+	  {+	    b_normal_exp = a_normal_exp;+	    b_fraction = 0;+	  }+	else+	  {+	    a_normal_exp = b_normal_exp;+	    a_fraction = 0;+	  }+      }+  }++  if (a->sign != b->sign)+    {+      if (a->sign)+	{+	  tfraction = -a_fraction + b_fraction;+	}+      else+	{+	  tfraction = a_fraction - b_fraction;+	}+      if (tfraction >= 0)+	{+	  tmp->sign = 0;+	  tmp->normal_exp = a_normal_exp;+	  tmp->fraction.ll = tfraction;+	}+      else+	{+	  tmp->sign = 1;+	  tmp->normal_exp = a_normal_exp;+	  tmp->fraction.ll = -tfraction;+	}+      /* and renormalize it */++      while (tmp->fraction.ll < IMPLICIT_1 && tmp->fraction.ll)+	{+	  tmp->fraction.ll <<= 1;+	  tmp->normal_exp--;+	}+    }+  else+    {+      tmp->sign = a->sign;+      tmp->normal_exp = a_normal_exp;+      tmp->fraction.ll = a_fraction + b_fraction;+    }+  tmp->class = CLASS_NUMBER;+  /* Now the fraction is added, we have to shift down to renormalize the+     number */++  if (tmp->fraction.ll >= IMPLICIT_2)+    {+      LSHIFT (tmp->fraction.ll);+      tmp->normal_exp++;+    }+  return tmp;++}++FLO_type+add (FLO_type arg_a, FLO_type arg_b)+{+  fp_number_type a;+  fp_number_type b;+  fp_number_type tmp;+  fp_number_type *res;+  FLO_union_type au, bu;++  au.value = arg_a;+  bu.value = arg_b;++  unpack_d (&au, &a);+  unpack_d (&bu, &b);++  res = _fpadd_parts (&a, &b, &tmp);++  return pack_d (res);+}++FLO_type+sub (FLO_type arg_a, FLO_type arg_b)+{+  fp_number_type a;+  fp_number_type b;+  fp_number_type tmp;+  fp_number_type *res;+  FLO_union_type au, bu;++  au.value = arg_a;+  bu.value = arg_b;++  unpack_d (&au, &a);+  unpack_d (&bu, &b);++  b.sign ^= 1;++  res = _fpadd_parts (&a, &b, &tmp);++  return pack_d (res);+}+#endif /* L_addsub_sf || L_addsub_df */++#if defined(L_mul_sf) || defined(L_mul_df) || defined(L_mul_tf)+static inline __attribute__ ((__always_inline__)) fp_number_type *+_fpmul_parts ( fp_number_type *  a,+	       fp_number_type *  b,+	       fp_number_type * tmp)+{+  fractype low = 0;+  fractype high = 0;++  if (isnan (a))+    {+      a->sign = a->sign != b->sign;+      return a;+    }+  if (isnan (b))+    {+      b->sign = a->sign != b->sign;+      return b;+    }+  if (isinf (a))+    {+      if (iszero (b))+	return nan ();+      a->sign = a->sign != b->sign;+      return a;+    }+  if (isinf (b))+    {+      if (iszero (a))+	{+	  return nan ();+	}+      b->sign = a->sign != b->sign;+      return b;+    }+  if (iszero (a))+    {+      a->sign = a->sign != b->sign;+      return a;+    }+  if (iszero (b))+    {+      b->sign = a->sign != b->sign;+      return b;+    }++  /* Calculate the mantissa by multiplying both numbers to get a+     twice-as-wide number.  */+  {+#if defined(NO_DI_MODE) || defined(TFLOAT)+    {+      fractype x = a->fraction.ll;+      fractype ylow = b->fraction.ll;+      fractype yhigh = 0;+      int bit;++      /* ??? This does multiplies one bit at a time.  Optimize.  */+      for (bit = 0; bit < FRAC_NBITS; bit++)+	{+	  int carry;++	  if (x & 1)+	    {+	      carry = (low += ylow) < ylow;+	      high += yhigh + carry;+	    }+	  yhigh <<= 1;+	  if (ylow & FRACHIGH)+	    {+	      yhigh |= 1;+	    }+	  ylow <<= 1;+	  x >>= 1;+	}+    }+#elif defined(FLOAT) +    /* Multiplying two USIs to get a UDI, we're safe.  */+    {+      UDItype answer = (UDItype)a->fraction.ll * (UDItype)b->fraction.ll;+      +      high = answer >> BITS_PER_SI;+      low = answer;+    }+#else+    /* fractype is DImode, but we need the result to be twice as wide.+       Assuming a widening multiply from DImode to TImode is not+       available, build one by hand.  */+    {+      USItype nl = a->fraction.ll;+      USItype nh = a->fraction.ll >> BITS_PER_SI;+      USItype ml = b->fraction.ll;+      USItype mh = b->fraction.ll >> BITS_PER_SI;+      UDItype pp_ll = (UDItype) ml * nl;+      UDItype pp_hl = (UDItype) mh * nl;+      UDItype pp_lh = (UDItype) ml * nh;+      UDItype pp_hh = (UDItype) mh * nh;+      UDItype res2 = 0;+      UDItype res0 = 0;+      UDItype ps_hh__ = pp_hl + pp_lh;+      if (ps_hh__ < pp_hl)+	res2 += (UDItype)1 << BITS_PER_SI;+      pp_hl = (UDItype)(USItype)ps_hh__ << BITS_PER_SI;+      res0 = pp_ll + pp_hl;+      if (res0 < pp_ll)+	res2++;+      res2 += (ps_hh__ >> BITS_PER_SI) + pp_hh;+      high = res2;+      low = res0;+    }+#endif+  }++  tmp->normal_exp = a->normal_exp + b->normal_exp+    + FRAC_NBITS - (FRACBITS + NGARDS);+  tmp->sign = a->sign != b->sign;+  while (high >= IMPLICIT_2)+    {+      tmp->normal_exp++;+      if (high & 1)+	{+	  low >>= 1;+	  low |= FRACHIGH;+	}+      high >>= 1;+    }+  while (high < IMPLICIT_1)+    {+      tmp->normal_exp--;++      high <<= 1;+      if (low & FRACHIGH)+	high |= 1;+      low <<= 1;+    }+  /* rounding is tricky. if we only round if it won't make us round later.  */+#if 0+  if (low & FRACHIGH2)+    {+      if (((high & GARDMASK) != GARDMSB)+	  && (((high + 1) & GARDMASK) == GARDMSB))+	{+	  /* don't round, it gets done again later.  */+	}+      else+	{+	  high++;+	}+    }+#endif+  if (!ROUND_TOWARDS_ZERO && (high & GARDMASK) == GARDMSB)+    {+      if (high & (1 << NGARDS))+	{+	  /* half way, so round to even */+	  high += GARDROUND + 1;+	}+      else if (low)+	{+	  /* but we really weren't half way */+	  high += GARDROUND + 1;+	}+    }+  tmp->fraction.ll = high;+  tmp->class = CLASS_NUMBER;+  return tmp;+}++FLO_type+multiply (FLO_type arg_a, FLO_type arg_b)+{+  fp_number_type a;+  fp_number_type b;+  fp_number_type tmp;+  fp_number_type *res;+  FLO_union_type au, bu;++  au.value = arg_a;+  bu.value = arg_b;++  unpack_d (&au, &a);+  unpack_d (&bu, &b);++  res = _fpmul_parts (&a, &b, &tmp);++  return pack_d (res);+}+#endif /* L_mul_sf || L_mul_df */++#if defined(L_div_sf) || defined(L_div_df) || defined(L_div_tf)+static inline __attribute__ ((__always_inline__)) fp_number_type *+_fpdiv_parts (fp_number_type * a,+	      fp_number_type * b)+{+  fractype bit;+  fractype numerator;+  fractype denominator;+  fractype quotient;++  if (isnan (a))+    {+      return a;+    }+  if (isnan (b))+    {+      return b;+    }++  a->sign = a->sign ^ b->sign;++  if (isinf (a) || iszero (a))+    {+      if (a->class == b->class)