static const double TWO43 = 8796093022208.0; static const double TWO52 = 4503599627370496.0;
static const double THREEp42 = 13194139533312.0; static const double THREEp51 = 6755399441055744.0;
__ieee754_exp (double x){  static const double himark = 709.7827128933840868;  static const double lomark = -745.1332191019412221;  /* Check for usual case.  */  if (isless (x, himark) && isgreater (x, lomark))    {      static const double TWO43 = 8796093022208.0;      static const double TWO52 = 4503599627370496.0;      /* 1/ln(2).  */      static const double M_1_LN2 = 1.442695040888963387;      /* ln(2), part 1 */      static const double M_LN2_0 = .6931471805598903302;      /* ln(2), part 2 */      static const double M_LN2_1 = 5.497923018708371155e-14;      int tval, unsafe, n_i;      double x22, n, t, dely, result;      union ieee754_double ex2_u, scale_u;      fenv_t oldenv;      feholdexcept (&oldenv);      fesetround (FE_TONEAREST);      /* Calculate n.  */      if (x >= 0)	{	  n = x * M_1_LN2 + TWO52;	  n -= TWO52;	}      else	{	  n = x * M_1_LN2 - TWO52;	  n += TWO52;	}      x = x - n*M_LN2_0;      if (x >= 0)	{	  /* Calculate t/512.  */	  t = x + TWO43;	  t -= TWO43;	  x -= t;	  /* Compute tval = t.  */	  tval = (int) (t * 512.0);	  x -= __exp_deltatable[tval];	}      else	{	  /* As above, but x is negative.  */	  t = x - TWO43;	  t += TWO43;	  x -= t;	  tval = (int) (t * 512.0);	  x += __exp_deltatable[-tval];	}      /* Now, the variable x contains x + n*ln(2)_1.  */      dely = n*M_LN2_1;      /* Compute ex2 = 2^n_0 e^(t/512+delta[t]).  */      ex2_u.d = __exp_atable[tval+177];      n_i = (int)n;      /* 'unsafe' is 1 iff n_1 != 0.  */      unsafe = abs(n_i) >= -DBL_MIN_EXP - 1;      ex2_u.ieee.exponent += n_i >> unsafe;      /* Compute scale = 2^n_1.  */      scale_u.d = 1.0;      scale_u.ieee.exponent += n_i - (n_i >> unsafe);      /* Approximate e^x2 - 1, using a fourth-degree polynomial,	 with maximum error in [-2^-10-2^-28,2^-10+2^-28]	 less than 4.9e-19.  */      x22 = (((0.04166666898464281565	       * x + 0.1666666766008501610)	      * x + 0.499999999999990008)	     * x + 0.9999999999999976685) * x;      /* Allow for impact of dely.  */      x22 -= dely + dely*x22;      /* Return result.  */      fesetenv (&oldenv);      result = x22 * ex2_u.d + ex2_u.d;      if (!unsafe)	return result;      else	return result * scale_u.d;    }  /* Exceptional cases:  */  else if (isless (x, himark))    {      if (__isinf (x))	/* e^-inf == 0, with no error.  */	return 0;      else	/* Underflow */	return TWOM1000 * TWOM1000;    }  else    /* Return x, if x is a NaN or Inf; or overflow, otherwise.  */    return TWO1023*x;}