[thirdparty/glibc.git] / sysdeps / ieee754 / ldbl-96 / s_fmal.c

/* Compute x * y + z as ternary operation.
   Copyright (C) 2010-2012 Free Software Foundation, Inc.
   This file is part of the GNU C Library.
   Contributed by Jakub Jelinek <jakub@redhat.com>, 2010.

   The GNU C Library is free software; you can redistribute it and/or
   modify it under the terms of the GNU Lesser General Public
   License as published by the Free Software Foundation; either
   version 2.1 of the License, or (at your option) any later version.

   The GNU C Library is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
   Lesser General Public License for more details.

   You should have received a copy of the GNU Lesser General Public
   License along with the GNU C Library; if not, see
   <http://www.gnu.org/licenses/>.  */

#include <float.h>
#include <math.h>
#include <fenv.h>
#include <ieee754.h>
#include <math_private.h>
#include <tininess.h>

/* This implementation uses rounding to odd to avoid problems with
   double rounding.  See a paper by Boldo and Melquiond:
   http://www.lri.fr/~melquion/doc/08-tc.pdf  */

long double
__fmal (long double x, long double y, long double z)
{
  union ieee854_long_double u, v, w;
  int adjust = 0;
  u.d = x;
  v.d = y;
  w.d = z;
  if (__builtin_expect (u.ieee.exponent + v.ieee.exponent
			>= 0x7fff + IEEE854_LONG_DOUBLE_BIAS
			   - LDBL_MANT_DIG, 0)
      || __builtin_expect (u.ieee.exponent >= 0x7fff - LDBL_MANT_DIG, 0)
      || __builtin_expect (v.ieee.exponent >= 0x7fff - LDBL_MANT_DIG, 0)
      || __builtin_expect (w.ieee.exponent >= 0x7fff - LDBL_MANT_DIG, 0)
      || __builtin_expect (u.ieee.exponent + v.ieee.exponent
			   <= IEEE854_LONG_DOUBLE_BIAS + LDBL_MANT_DIG, 0))
    {
      /* If z is Inf, but x and y are finite, the result should be
	 z rather than NaN.  */
      if (w.ieee.exponent == 0x7fff
	  && u.ieee.exponent != 0x7fff
          && v.ieee.exponent != 0x7fff)
	return (z + x) + y;
      /* If z is zero and x are y are nonzero, compute the result
	 as x * y to avoid the wrong sign of a zero result if x * y
	 underflows to 0.  */
      if (z == 0 && x != 0 && y != 0)
	return x * y;
      /* If x or y or z is Inf/NaN, or if x * y is zero, compute as
	 x * y + z.  */
      if (u.ieee.exponent == 0x7fff
	  || v.ieee.exponent == 0x7fff
	  || w.ieee.exponent == 0x7fff
	  || x == 0
	  || y == 0)
	return x * y + z;
      /* If fma will certainly overflow, compute as x * y.  */
      if (u.ieee.exponent + v.ieee.exponent
	  > 0x7fff + IEEE854_LONG_DOUBLE_BIAS)
	return x * y;
      /* If x * y is less than 1/4 of LDBL_DENORM_MIN, neither the
	 result nor whether there is underflow depends on its exact
	 value, only on its sign.  */
      if (u.ieee.exponent + v.ieee.exponent
	  < IEEE854_LONG_DOUBLE_BIAS - LDBL_MANT_DIG - 2)
	{
	  int neg = u.ieee.negative ^ v.ieee.negative;
	  long double tiny = neg ? -0x1p-16445L : 0x1p-16445L;
	  if (w.ieee.exponent >= 3)
	    return tiny + z;
	  /* Scaling up, adding TINY and scaling down produces the
	     correct result, because in round-to-nearest mode adding
	     TINY has no effect and in other modes double rounding is
	     harmless.  But it may not produce required underflow
	     exceptions.  */
	  v.d = z * 0x1p65L + tiny;
	  if (TININESS_AFTER_ROUNDING
	      ? v.ieee.exponent < 66
	      : (w.ieee.exponent == 0
		 || (w.ieee.exponent == 1
		     && w.ieee.negative != neg
		     && w.ieee.mantissa1 == 0
		     && w.ieee.mantissa0 == 0x80000000)))
	    {
	      volatile long double force_underflow = x * y;
	      (void) force_underflow;
	    }
	  return v.d * 0x1p-65L;
	}
      if (u.ieee.exponent + v.ieee.exponent
	  >= 0x7fff + IEEE854_LONG_DOUBLE_BIAS - LDBL_MANT_DIG)
	{
	  /* Compute 1p-64 times smaller result and multiply
	     at the end.  */
	  if (u.ieee.exponent > v.ieee.exponent)
	    u.ieee.exponent -= LDBL_MANT_DIG;
	  else
	    v.ieee.exponent -= LDBL_MANT_DIG;
	  /* If x + y exponent is very large and z exponent is very small,
	     it doesn't matter if we don't adjust it.  */
	  if (w.ieee.exponent > LDBL_MANT_DIG)
	    w.ieee.exponent -= LDBL_MANT_DIG;
	  adjust = 1;
	}
      else if (w.ieee.exponent >= 0x7fff - LDBL_MANT_DIG)
	{
	  /* Similarly.
	     If z exponent is very large and x and y exponents are
	     very small, it doesn't matter if we don't adjust it.  */
	  if (u.ieee.exponent > v.ieee.exponent)
	    {
	      if (u.ieee.exponent > LDBL_MANT_DIG)
		u.ieee.exponent -= LDBL_MANT_DIG;
	    }
	  else if (v.ieee.exponent > LDBL_MANT_DIG)
	    v.ieee.exponent -= LDBL_MANT_DIG;
	  w.ieee.exponent -= LDBL_MANT_DIG;
	  adjust = 1;
	}
      else if (u.ieee.exponent >= 0x7fff - LDBL_MANT_DIG)
	{
	  u.ieee.exponent -= LDBL_MANT_DIG;
	  if (v.ieee.exponent)
	    v.ieee.exponent += LDBL_MANT_DIG;
	  else
	    v.d *= 0x1p64L;
	}
      else if (v.ieee.exponent >= 0x7fff - LDBL_MANT_DIG)
	{
	  v.ieee.exponent -= LDBL_MANT_DIG;
	  if (u.ieee.exponent)
	    u.ieee.exponent += LDBL_MANT_DIG;
	  else
	    u.d *= 0x1p64L;
	}
      else /* if (u.ieee.exponent + v.ieee.exponent
		  <= IEEE854_LONG_DOUBLE_BIAS + LDBL_MANT_DIG) */
	{
	  if (u.ieee.exponent > v.ieee.exponent)
	    u.ieee.exponent += 2 * LDBL_MANT_DIG;
	  else
	    v.ieee.exponent += 2 * LDBL_MANT_DIG;
	  if (w.ieee.exponent <= 4 * LDBL_MANT_DIG + 4)
	    {
	      if (w.ieee.exponent)
		w.ieee.exponent += 2 * LDBL_MANT_DIG;
	      else
		w.d *= 0x1p128L;
	      adjust = -1;
	    }
	  /* Otherwise x * y should just affect inexact
	     and nothing else.  */
	}
      x = u.d;
      y = v.d;
      z = w.d;
    }

  /* Ensure correct sign of exact 0 + 0.  */
  if (__builtin_expect ((x == 0 || y == 0) && z == 0, 0))
    return x * y + z;

  fenv_t env;
  feholdexcept (&env);
  fesetround (FE_TONEAREST);

  /* Multiplication m1 + m2 = x * y using Dekker's algorithm.  */
#define C ((1LL << (LDBL_MANT_DIG + 1) / 2) + 1)
  long double x1 = x * C;
  long double y1 = y * C;
  long double m1 = x * y;
  x1 = (x - x1) + x1;
  y1 = (y - y1) + y1;
  long double x2 = x - x1;
  long double y2 = y - y1;
  long double m2 = (((x1 * y1 - m1) + x1 * y2) + x2 * y1) + x2 * y2;

  /* Addition a1 + a2 = z + m1 using Knuth's algorithm.  */
  long double a1 = z + m1;
  long double t1 = a1 - z;
  long double t2 = a1 - t1;
  t1 = m1 - t1;
  t2 = z - t2;
  long double a2 = t1 + t2;
  feclearexcept (FE_INEXACT);

  /* If the result is an exact zero, ensure it has the correct
     sign.  */
  if (a1 == 0 && m2 == 0)
    {
      feupdateenv (&env);
      /* Ensure that round-to-nearest value of z + m1 is not
	 reused.  */
      asm volatile ("" : "=m" (z) : "m" (z));
      return z + m1;
    }

  fesetround (FE_TOWARDZERO);
  /* Perform m2 + a2 addition with round to odd.  */
  u.d = a2 + m2;

  if (__builtin_expect (adjust == 0, 1))
    {
      if ((u.ieee.mantissa1 & 1) == 0 && u.ieee.exponent != 0x7fff)
	u.ieee.mantissa1 |= fetestexcept (FE_INEXACT) != 0;
      feupdateenv (&env);
      /* Result is a1 + u.d.  */
      return a1 + u.d;
    }
  else if (__builtin_expect (adjust > 0, 1))
    {
      if ((u.ieee.mantissa1 & 1) == 0 && u.ieee.exponent != 0x7fff)
	u.ieee.mantissa1 |= fetestexcept (FE_INEXACT) != 0;
      feupdateenv (&env);
      /* Result is a1 + u.d, scaled up.  */
      return (a1 + u.d) * 0x1p64L;
    }
  else
    {
      if ((u.ieee.mantissa1 & 1) == 0)
	u.ieee.mantissa1 |= fetestexcept (FE_INEXACT) != 0;
      v.d = a1 + u.d;
      /* Ensure the addition is not scheduled after fetestexcept call.  */
      math_force_eval (v.d);
      int j = fetestexcept (FE_INEXACT) != 0;
      feupdateenv (&env);
      /* Ensure the following computations are performed in default rounding
	 mode instead of just reusing the round to zero computation.  */
      asm volatile ("" : "=m" (u) : "m" (u));
      /* If a1 + u.d is exact, the only rounding happens during
	 scaling down.  */
      if (j == 0)
	return v.d * 0x1p-128L;
      /* If result rounded to zero is not subnormal, no double
	 rounding will occur.  */
      if (v.ieee.exponent > 128)
	return (a1 + u.d) * 0x1p-128L;
      /* If v.d * 0x1p-128L with round to zero is a subnormal above
	 or equal to LDBL_MIN / 2, then v.d * 0x1p-128L shifts mantissa
	 down just by 1 bit, which means v.ieee.mantissa1 |= j would
	 change the round bit, not sticky or guard bit.
	 v.d * 0x1p-128L never normalizes by shifting up,
	 so round bit plus sticky bit should be already enough
	 for proper rounding.  */
      if (v.ieee.exponent == 128)
	{
	  /* If the exponent would be in the normal range when
	     rounding to normal precision with unbounded exponent
	     range, the exact result is known and spurious underflows
	     must be avoided on systems detecting tininess after
	     rounding.  */
	  if (TININESS_AFTER_ROUNDING)
	    {
	      w.d = a1 + u.d;
	      if (w.ieee.exponent == 129)
		return w.d * 0x1p-128L;
	    }
	  /* v.ieee.mantissa1 & 2 is LSB bit of the result before rounding,
	     v.ieee.mantissa1 & 1 is the round bit and j is our sticky
	     bit.  */
	  w.d = 0.0L;
	  w.ieee.mantissa1 = ((v.ieee.mantissa1 & 3) << 1) | j;
	  w.ieee.negative = v.ieee.negative;
	  v.ieee.mantissa1 &= ~3U;
	  v.d *= 0x1p-128L;
	  w.d *= 0x1p-2L;
	  return v.d + w.d;
	}
      v.ieee.mantissa1 |= j;
      return v.d * 0x1p-128L;
    }
}
weak_alias (__fmal, fmal)
Commit	Line	Data
3e692e05	1	/* Compute x * y + z as ternary operation.
4842e4fe	2	Copyright (C) 2010-2012 Free Software Foundation, Inc.
3e692e05 JJ	3	This file is part of the GNU C Library.
	4	Contributed by Jakub Jelinek <jakub@redhat.com>, 2010.
	5
	6	The GNU C Library is free software; you can redistribute it and/or
	7	modify it under the terms of the GNU Lesser General Public
	8	License as published by the Free Software Foundation; either
	9	version 2.1 of the License, or (at your option) any later version.
	10
	11	The GNU C Library is distributed in the hope that it will be useful,
	12	but WITHOUT ANY WARRANTY; without even the implied warranty of
	13	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
	14	Lesser General Public License for more details.
	15
	16	You should have received a copy of the GNU Lesser General Public
59ba27a6 PE	17	License along with the GNU C Library; if not, see
59ba27a6 PE	18	<http://www.gnu.org/licenses/>. */
3e692e05 JJ	19
	20	#include <float.h>
	21	#include <math.h>
	22	#include <fenv.h>
	23	#include <ieee754.h>
4842e4fe	24	#include <math_private.h>
ef82f4da	25	#include <tininess.h>
3e692e05 JJ	26
	27	/* This implementation uses rounding to odd to avoid problems with
	28	double rounding. See a paper by Boldo and Melquiond:
	29	http://www.lri.fr/~melquion/doc/08-tc.pdf */
	30
	31	long double
	32	__fmal (long double x, long double y, long double z)
	33	{
	34	union ieee854_long_double u, v, w;
	35	int adjust = 0;
	36	u.d = x;
	37	v.d = y;
	38	w.d = z;
	39	if (__builtin_expect (u.ieee.exponent + v.ieee.exponent
	40	>= 0x7fff + IEEE854_LONG_DOUBLE_BIAS
	41	- LDBL_MANT_DIG, 0)
	42	\|\| __builtin_expect (u.ieee.exponent >= 0x7fff - LDBL_MANT_DIG, 0)
	43	\|\| __builtin_expect (v.ieee.exponent >= 0x7fff - LDBL_MANT_DIG, 0)
	44	\|\| __builtin_expect (w.ieee.exponent >= 0x7fff - LDBL_MANT_DIG, 0)
	45	\|\| __builtin_expect (u.ieee.exponent + v.ieee.exponent
	46	<= IEEE854_LONG_DOUBLE_BIAS + LDBL_MANT_DIG, 0))
	47	{
	48	/* If z is Inf, but x and y are finite, the result should be
	49	z rather than NaN. */
	50	if (w.ieee.exponent == 0x7fff
	51	&& u.ieee.exponent != 0x7fff
	52	&& v.ieee.exponent != 0x7fff)
	53	return (z + x) + y;
bec749fd JM	54	/* If z is zero and x are y are nonzero, compute the result
	55	as x * y to avoid the wrong sign of a zero result if x * y
	56	underflows to 0. */
	57	if (z == 0 && x != 0 && y != 0)
	58	return x * y;
a0c2940d JM	59	/* If x or y or z is Inf/NaN, or if x * y is zero, compute as
a0c2940d JM	60	x * y + z. */
3e692e05 JJ	61	if (u.ieee.exponent == 0x7fff
	62	\|\| v.ieee.exponent == 0x7fff
	63	\|\| w.ieee.exponent == 0x7fff
473611b2 JM	64	\|\| x == 0
473611b2 JM	65	\|\| y == 0)
3e692e05	66	return x * y + z;
a0c2940d JM	67	/* If fma will certainly overflow, compute as x * y. */
	68	if (u.ieee.exponent + v.ieee.exponent
	69	> 0x7fff + IEEE854_LONG_DOUBLE_BIAS)
	70	return x * y;
473611b2 JM	71	/* If x * y is less than 1/4 of LDBL_DENORM_MIN, neither the
	72	result nor whether there is underflow depends on its exact
	73	value, only on its sign. */
	74	if (u.ieee.exponent + v.ieee.exponent
	75	< IEEE854_LONG_DOUBLE_BIAS - LDBL_MANT_DIG - 2)
	76	{
	77	int neg = u.ieee.negative ^ v.ieee.negative;
	78	long double tiny = neg ? -0x1p-16445L : 0x1p-16445L;
	79	if (w.ieee.exponent >= 3)
	80	return tiny + z;
	81	/* Scaling up, adding TINY and scaling down produces the
	82	correct result, because in round-to-nearest mode adding
	83	TINY has no effect and in other modes double rounding is
	84	harmless. But it may not produce required underflow
	85	exceptions. */
	86	v.d = z * 0x1p65L + tiny;
	87	if (TININESS_AFTER_ROUNDING
	88	? v.ieee.exponent < 66
	89	: (w.ieee.exponent == 0
	90	\|\| (w.ieee.exponent == 1
	91	&& w.ieee.negative != neg
	92	&& w.ieee.mantissa1 == 0
	93	&& w.ieee.mantissa0 == 0x80000000)))
	94	{
	95	volatile long double force_underflow = x * y;
	96	(void) force_underflow;
	97	}
	98	return v.d * 0x1p-65L;
	99	}
3e692e05 JJ	100	if (u.ieee.exponent + v.ieee.exponent
	101	>= 0x7fff + IEEE854_LONG_DOUBLE_BIAS - LDBL_MANT_DIG)
	102	{
	103	/* Compute 1p-64 times smaller result and multiply
	104	at the end. */
	105	if (u.ieee.exponent > v.ieee.exponent)
	106	u.ieee.exponent -= LDBL_MANT_DIG;
	107	else
	108	v.ieee.exponent -= LDBL_MANT_DIG;
	109	/* If x + y exponent is very large and z exponent is very small,
	110	it doesn't matter if we don't adjust it. */
	111	if (w.ieee.exponent > LDBL_MANT_DIG)
	112	w.ieee.exponent -= LDBL_MANT_DIG;
	113	adjust = 1;
	114	}
	115	else if (w.ieee.exponent >= 0x7fff - LDBL_MANT_DIG)
	116	{
	117	/* Similarly.
	118	If z exponent is very large and x and y exponents are
	119	very small, it doesn't matter if we don't adjust it. */
	120	if (u.ieee.exponent > v.ieee.exponent)
	121	{
	122	if (u.ieee.exponent > LDBL_MANT_DIG)
	123	u.ieee.exponent -= LDBL_MANT_DIG;
	124	}
	125	else if (v.ieee.exponent > LDBL_MANT_DIG)
	126	v.ieee.exponent -= LDBL_MANT_DIG;
	127	w.ieee.exponent -= LDBL_MANT_DIG;
	128	adjust = 1;
	129	}
	130	else if (u.ieee.exponent >= 0x7fff - LDBL_MANT_DIG)
	131	{
	132	u.ieee.exponent -= LDBL_MANT_DIG;
	133	if (v.ieee.exponent)
	134	v.ieee.exponent += LDBL_MANT_DIG;
	135	else
	136	v.d *= 0x1p64L;
	137	}
	138	else if (v.ieee.exponent >= 0x7fff - LDBL_MANT_DIG)
	139	{
	140	v.ieee.exponent -= LDBL_MANT_DIG;
	141	if (u.ieee.exponent)
	142	u.ieee.exponent += LDBL_MANT_DIG;
	143	else
	144	u.d *= 0x1p64L;
	145	}
	146	else /* if (u.ieee.exponent + v.ieee.exponent
	147	<= IEEE854_LONG_DOUBLE_BIAS + LDBL_MANT_DIG) */
	148	{
	149	if (u.ieee.exponent > v.ieee.exponent)
	150	u.ieee.exponent += 2 * LDBL_MANT_DIG;
	151	else
	152	v.ieee.exponent += 2 * LDBL_MANT_DIG;
	153	if (w.ieee.exponent <= 4 * LDBL_MANT_DIG + 4)
	154	{
	155	if (w.ieee.exponent)
	156	w.ieee.exponent += 2 * LDBL_MANT_DIG;
	157	else
	158	w.d *= 0x1p128L;
	159	adjust = -1;
	160	}
	161	/* Otherwise x * y should just affect inexact
	162	and nothing else. */
	163	}
164	x = u.d;
165	y = v.d;
166	z = w.d;
167	}
8ec5b013 JM	168
	169	/* Ensure correct sign of exact 0 + 0. */
	170	if (__builtin_expect ((x == 0 \|\| y == 0) && z == 0, 0))
	171	return x * y + z;
	172
5b5b04d6 JM	173	fenv_t env;
	174	feholdexcept (&env);
	175	fesetround (FE_TONEAREST);
	176
3e692e05 JJ	177	/* Multiplication m1 + m2 = x * y using Dekker's algorithm. */
	178	#define C ((1LL << (LDBL_MANT_DIG + 1) / 2) + 1)
	179	long double x1 = x * C;
	180	long double y1 = y * C;
	181	long double m1 = x * y;
	182	x1 = (x - x1) + x1;
	183	y1 = (y - y1) + y1;
	184	long double x2 = x - x1;
	185	long double y2 = y - y1;
	186	long double m2 = (((x1 * y1 - m1) + x1 * y2) + x2 * y1) + x2 * y2;
	187
	188	/* Addition a1 + a2 = z + m1 using Knuth's algorithm. */
	189	long double a1 = z + m1;
	190	long double t1 = a1 - z;
	191	long double t2 = a1 - t1;
	192	t1 = m1 - t1;
	193	t2 = z - t2;
	194	long double a2 = t1 + t2;
5b5b04d6 JM	195	feclearexcept (FE_INEXACT);
	196
	197	/* If the result is an exact zero, ensure it has the correct
	198	sign. */
	199	if (a1 == 0 && m2 == 0)
	200	{
	201	feupdateenv (&env);
	202	/* Ensure that round-to-nearest value of z + m1 is not
	203	reused. */
	204	asm volatile ("" : "=m" (z) : "m" (z));
	205	return z + m1;
	206	}
3e692e05	207
3e692e05 JJ	208	fesetround (FE_TOWARDZERO);
	209	/* Perform m2 + a2 addition with round to odd. */
	210	u.d = a2 + m2;
	211
	212	if (__builtin_expect (adjust == 0, 1))
	213	{
	214	if ((u.ieee.mantissa1 & 1) == 0 && u.ieee.exponent != 0x7fff)
	215	u.ieee.mantissa1 \|= fetestexcept (FE_INEXACT) != 0;
	216	feupdateenv (&env);
	217	/* Result is a1 + u.d. */
	218	return a1 + u.d;
	219	}
	220	else if (__builtin_expect (adjust > 0, 1))
	221	{
	222	if ((u.ieee.mantissa1 & 1) == 0 && u.ieee.exponent != 0x7fff)
	223	u.ieee.mantissa1 \|= fetestexcept (FE_INEXACT) != 0;
	224	feupdateenv (&env);
	225	/* Result is a1 + u.d, scaled up. */
	226	return (a1 + u.d) * 0x1p64L;
	227	}
	228	else
	229	{
	230	if ((u.ieee.mantissa1 & 1) == 0)
	231	u.ieee.mantissa1 \|= fetestexcept (FE_INEXACT) != 0;
	232	v.d = a1 + u.d;
4842e4fe JM	233	/* Ensure the addition is not scheduled after fetestexcept call. */
4842e4fe JM	234	math_force_eval (v.d);
3e692e05 JJ	235	int j = fetestexcept (FE_INEXACT) != 0;
	236	feupdateenv (&env);
	237	/* Ensure the following computations are performed in default rounding
	238	mode instead of just reusing the round to zero computation. */
	239	asm volatile ("" : "=m" (u) : "m" (u));
	240	/* If a1 + u.d is exact, the only rounding happens during
	241	scaling down. */
	242	if (j == 0)
	243	return v.d * 0x1p-128L;
	244	/* If result rounded to zero is not subnormal, no double
	245	rounding will occur. */
	246	if (v.ieee.exponent > 128)
	247	return (a1 + u.d) * 0x1p-128L;
	248	/* If v.d * 0x1p-128L with round to zero is a subnormal above
	249	or equal to LDBL_MIN / 2, then v.d * 0x1p-128L shifts mantissa
	250	down just by 1 bit, which means v.ieee.mantissa1 \|= j would
	251	change the round bit, not sticky or guard bit.
	252	v.d * 0x1p-128L never normalizes by shifting up,
	253	so round bit plus sticky bit should be already enough
	254	for proper rounding. */
	255	if (v.ieee.exponent == 128)
	256	{
ef82f4da JM	257	/* If the exponent would be in the normal range when
	258	rounding to normal precision with unbounded exponent
	259	range, the exact result is known and spurious underflows
	260	must be avoided on systems detecting tininess after
	261	rounding. */
	262	if (TININESS_AFTER_ROUNDING)
	263	{
	264	w.d = a1 + u.d;
	265	if (w.ieee.exponent == 129)
	266	return w.d * 0x1p-128L;
	267	}
3e692e05 JJ	268	/* v.ieee.mantissa1 & 2 is LSB bit of the result before rounding,
3e692e05 JJ	269	v.ieee.mantissa1 & 1 is the round bit and j is our sticky
8627a232 JM	270	bit. */
	271	w.d = 0.0L;
	272	w.ieee.mantissa1 = ((v.ieee.mantissa1 & 3) << 1) \| j;
	273	w.ieee.negative = v.ieee.negative;
	274	v.ieee.mantissa1 &= ~3U;
	275	v.d *= 0x1p-128L;
	276	w.d *= 0x1p-2L;
	277	return v.d + w.d;
3e692e05 JJ	278	}
	279	v.ieee.mantissa1 \|= j;
	280	return v.d * 0x1p-128L;
	281	}
	282	}
	283	weak_alias (__fmal, fmal)