return one;   if (x == one) return one; if (x == -1.0L && iy == 0x7fff0000 && (q.parts32.w1 | q.parts32.w2 | q.parts32.w3) == 0)
__ieee754_powl (long double x, long double y)#elselong double__ieee754_powl (x, y)     long double x, y;#endif{  long double z, ax, z_h, z_l, p_h, p_l;  long double y1, t1, t2, r, s, t, u, v, w;  long double s2, s_h, s_l, t_h, t_l;  int32_t i, j, k, yisint, n;  u_int32_t ix, iy;  int32_t hx, hy;  ieee854_long_double_shape_type o, p, q;  p.value = x;  hx = p.parts32.w0;  ix = hx & 0x7fffffff;  q.value = y;  hy = q.parts32.w0;  iy = hy & 0x7fffffff;  /* y==zero: x**0 = 1 */  if ((iy | q.parts32.w1 | q.parts32.w2 | q.parts32.w3) == 0)    return one;  /* +-NaN return x+y */  if ((ix > 0x7fff0000)      || ((ix == 0x7fff0000)	  && ((p.parts32.w1 | p.parts32.w2 | p.parts32.w3) != 0))      || (iy > 0x7fff0000)      || ((iy == 0x7fff0000)	  && ((q.parts32.w1 | q.parts32.w2 | q.parts32.w3) != 0)))    return x + y;  /* determine if y is an odd int when x < 0   * yisint = 0       ... y is not an integer   * yisint = 1       ... y is an odd int   * yisint = 2       ... y is an even int   */  yisint = 0;  if (hx < 0)    {      if (iy >= 0x40700000)	/* 2^113 */	yisint = 2;		/* even integer y */      else if (iy >= 0x3fff0000)	/* 1.0 */	{	  if (__floorl (y) == y)	    {	      z = 0.5 * y;	      if (__floorl (z) == z)		yisint = 2;	      else		yisint = 1;	    }	}    }  /* special value of y */  if ((q.parts32.w1 | q.parts32.w2 | q.parts32.w3) == 0)    {      if (iy == 0x7fff0000)	/* y is +-inf */	{	  if (((ix - 0x3fff0000) | p.parts32.w1 | p.parts32.w2 | p.parts32.w3)	      == 0)	    return y - y;	/* inf**+-1 is NaN */	  else if (ix >= 0x3fff0000)	/* (|x|>1)**+-inf = inf,0 */	    return (hy >= 0) ? y : zero;	  else			/* (|x|<1)**-,+inf = inf,0 */	    return (hy < 0) ? -y : zero;	}      if (iy == 0x3fff0000)	{			/* y is  +-1 */	  if (hy < 0)	    return one / x;	  else	    return x;	}      if (hy == 0x40000000)	return x * x;		/* y is  2 */      if (hy == 0x3ffe0000)	{			/* y is  0.5 */	  if (hx >= 0)		/* x >= +0 */	    return __ieee754_sqrtl (x);	}    }  ax = fabsl (x);  /* special value of x */  if ((p.parts32.w1 | p.parts32.w2 | p.parts32.w3) == 0)    {      if (ix == 0x7fff0000 || ix == 0 || ix == 0x3fff0000)	{	  z = ax;		/*x is +-0,+-inf,+-1 */	  if (hy < 0)	    z = one / z;	/* z = (1/|x|) */	  if (hx < 0)	    {	      if (((ix - 0x3fff0000) | yisint) == 0)		{		  z = (z - z) / (z - z);	/* (-1)**non-int is NaN */		}	      else if (yisint == 1)		z = -z;		/* (x<0)**odd = -(|x|**odd) */	    }	  return z;	}    }  /* (x<0)**(non-int) is NaN */  if (((((u_int32_t) hx >> 31) - 1) | yisint) == 0)    return (x - x) / (x - x);  /* |y| is huge.     2^-16495 = 1/2 of smallest representable value.     If (1 - 1/131072)^y underflows, y > 1.4986e9 */  if (iy > 0x401d654b)    {      /* if (1 - 2^-113)^y underflows, y > 1.1873e38 */      if (iy > 0x407d654b)	{	  if (ix <= 0x3ffeffff)	    return (hy < 0) ? huge * huge : tiny * tiny;	  if (ix >= 0x3fff0000)	    return (hy > 0) ? huge * huge : tiny * tiny;	}      /* over/underflow if x is not close to one */      if (ix < 0x3ffeffff)	return (hy < 0) ? huge * huge : tiny * tiny;      if (ix > 0x3fff0000)	return (hy > 0) ? huge * huge : tiny * tiny;    }  n = 0;  /* take care subnormal number */  if (ix < 0x00010000)    {      ax *= two113;      n -= 113;      o.value = ax;      ix = o.parts32.w0;    }  n += ((ix) >> 16) - 0x3fff;  j = ix & 0x0000ffff;  /* determine interval */  ix = j | 0x3fff0000;		/* normalize ix */  if (j <= 0x3988)    k = 0;			/* |x|<sqrt(3/2) */  else if (j < 0xbb67)    k = 1;			/* |x|<sqrt(3)   */  else    {      k = 0;      n += 1;      ix -= 0x00010000;    }  o.value = ax;  o.parts32.w0 = ix;  ax = o.value;  /* compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */  u = ax - bp[k];		/* bp[0]=1.0, bp[1]=1.5 */  v = one / (ax + bp[k]);  s = u * v;  s_h = s;  o.value = s_h;  o.parts32.w3 = 0;  o.parts32.w2 &= 0xf8000000;  s_h = o.value;  /* t_h=ax+bp[k] High */  t_h = ax + bp[k];  o.value = t_h;  o.parts32.w3 = 0;  o.parts32.w2 &= 0xf8000000;  t_h = o.value;  t_l = ax - (t_h - bp[k]);  s_l = v * ((u - s_h * t_h) - s_h * t_l);  /* compute log(ax) */  s2 = s * s;  u = LN[0] + s2 * (LN[1] + s2 * (LN[2] + s2 * (LN[3] + s2 * LN[4])));  v = LD[0] + s2 * (LD[1] + s2 * (LD[2] + s2 * (LD[3] + s2 * (LD[4] + s2))));  r = s2 * s2 * u / v;  r += s_l * (s_h + s);  s2 = s_h * s_h;  t_h = 3.0 + s2 + r;  o.value = t_h;  o.parts32.w3 = 0;  o.parts32.w2 &= 0xf8000000;  t_h = o.value;  t_l = r - ((t_h - 3.0) - s2);  /* u+v = s*(1+...) */  u = s_h * t_h;  v = s_l * t_h + t_l * s;  /* 2/(3log2)*(s+...) */  p_h = u + v;  o.value = p_h;  o.parts32.w3 = 0;  o.parts32.w2 &= 0xf8000000;  p_h = o.value;  p_l = v - (p_h - u);  z_h = cp_h * p_h;		/* cp_h+cp_l = 2/(3*log2) */  z_l = cp_l * p_h + p_l * cp + dp_l[k];  /* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */  t = (long double) n;  t1 = (((z_h + z_l) + dp_h[k]) + t);  o.value = t1;  o.parts32.w3 = 0;  o.parts32.w2 &= 0xf8000000;  t1 = o.value;  t2 = z_l - (((t1 - t) - dp_h[k]) - z_h);  /* s (sign of result -ve**odd) = -1 else = 1 */  s = one;  if (((((u_int32_t) hx >> 31) - 1) | (yisint - 1)) == 0)    s = -one;			/* (-ve)**(odd int) */  /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */  y1 = y;  o.value = y1;  o.parts32.w3 = 0;  o.parts32.w2 &= 0xf8000000;  y1 = o.value;  p_l = (y - y1) * t1 + y * t2;  p_h = y1 * t1;  z = p_l + p_h;  o.value = z;  j = o.parts32.w0;  if (j >= 0x400d0000) /* z >= 16384 */    {      /* if z > 16384 */      if (((j - 0x400d0000) | o.parts32.w1 | o.parts32.w2 | o.parts32.w3) != 0)	return s * huge * huge;	/* overflow */      else	{	  if (p_l + ovt > z - p_h)	    return s * huge * huge;	/* overflow */	}    }  else if ((j & 0x7fffffff) >= 0x400d01b9)	/* z <= -16495 */    {      /* z < -16495 */      if (((j - 0xc00d01bc) | o.parts32.w1 | o.parts32.w2 | o.parts32.w3)	  != 0)	return s * tiny * tiny;	/* underflow */      else	{	  if (p_l <= z - p_h)	    return s * tiny * tiny;	/* underflow */	}    }  /* compute 2**(p_h+p_l) */  i = j & 0x7fffffff;  k = (i >> 16) - 0x3fff;  n = 0;  if (i > 0x3ffe0000)    {				/* if |z| > 0.5, set n = [z+0.5] */      n = __floorl (z + 0.5L);      t = n;      p_h -= t;    }  t = p_l + p_h;  o.value = t;  o.parts32.w3 = 0;  o.parts32.w2 &= 0xf8000000;  t = o.value;  u = t * lg2_h;  v = (p_l - (t - p_h)) * lg2 + t * lg2_l;  z = u + v;  w = v - (z - u);  /*  exp(z) */  t = z * z;  u = PN[0] + t * (PN[1] + t * (PN[2] + t * (PN[3] + t * PN[4])));  v = PD[0] + t * (PD[1] + t * (PD[2] + t * (PD[3] + t)));  t1 = z - t * u / v;  r = (z * t1) / (t1 - two) - (w + z * w);  z = one - (r - z);  o.value = z;  j = o.parts32.w0;  j += (n << 16);  if ((j >> 16) <= 0)    z = __scalbnl (z, n);	/* subnormal output */  else    {      o.parts32.w0 = j;      z = o.value;    }  return s * z;}