[thirdparty/glibc.git] / sysdeps / ieee754 / ldbl-128ibm / e_powl.c

/*
 * ====================================================
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 *
 * Developed at SunPro, a Sun Microsystems, Inc. business.
 * Permission to use, copy, modify, and distribute this
 * software is freely granted, provided that this notice
 * is preserved.
 * ====================================================
 */

/* Expansions and modifications for 128-bit long double are
   Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov>
   and are incorporated herein by permission of the author.  The author
   reserves the right to distribute this material elsewhere under different
   copying permissions.  These modifications are distributed here under
   the following terms:

    This 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.

    This 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 this library; if not, see
    <https://www.gnu.org/licenses/>.  */

/* __ieee754_powl(x,y) return x**y
 *
 *		      n
 * Method:  Let x =  2   * (1+f)
 *	1. Compute and return log2(x) in two pieces:
 *		log2(x) = w1 + w2,
 *	   where w1 has 113-53 = 60 bit trailing zeros.
 *	2. Perform y*log2(x) = n+y' by simulating muti-precision
 *	   arithmetic, where |y'|<=0.5.
 *	3. Return x**y = 2**n*exp(y'*log2)
 *
 * Special cases:
 *	1.  (anything) ** 0  is 1
 *	2.  (anything) ** 1  is itself
 *	3.  (anything) ** NAN is NAN
 *	4.  NAN ** (anything except 0) is NAN
 *	5.  +-(|x| > 1) **  +INF is +INF
 *	6.  +-(|x| > 1) **  -INF is +0
 *	7.  +-(|x| < 1) **  +INF is +0
 *	8.  +-(|x| < 1) **  -INF is +INF
 *	9.  +-1         ** +-INF is NAN
 *	10. +0 ** (+anything except 0, NAN)               is +0
 *	11. -0 ** (+anything except 0, NAN, odd integer)  is +0
 *	12. +0 ** (-anything except 0, NAN)               is +INF
 *	13. -0 ** (-anything except 0, NAN, odd integer)  is +INF
 *	14. -0 ** (odd integer) = -( +0 ** (odd integer) )
 *	15. +INF ** (+anything except 0,NAN) is +INF
 *	16. +INF ** (-anything except 0,NAN) is +0
 *	17. -INF ** (anything)  = -0 ** (-anything)
 *	18. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer)
 *	19. (-anything except 0 and inf) ** (non-integer) is NAN
 *
 */

#include <math.h>
#include <math_private.h>
#include <math-underflow.h>
#include <libm-alias-finite.h>

static const long double bp[] = {
  1.0L,
  1.5L,
};

/* log_2(1.5) */
static const long double dp_h[] = {
  0.0,
  5.8496250072115607565592654282227158546448E-1L
};

/* Low part of log_2(1.5) */
static const long double dp_l[] = {
  0.0,
  1.0579781240112554492329533686862998106046E-16L
};

static const long double zero = 0.0L,
  one = 1.0L,
  two = 2.0L,
  two113 = 1.0384593717069655257060992658440192E34L,
  huge = 1.0e300L,
  tiny = 1.0e-300L;

/* 3/2 log x = 3 z + z^3 + z^3 (z^2 R(z^2))
   z = (x-1)/(x+1)
   1 <= x <= 1.25
   Peak relative error 2.3e-37 */
static const long double LN[] =
{
 -3.0779177200290054398792536829702930623200E1L,
  6.5135778082209159921251824580292116201640E1L,
 -4.6312921812152436921591152809994014413540E1L,
  1.2510208195629420304615674658258363295208E1L,
 -9.9266909031921425609179910128531667336670E-1L
};
static const long double LD[] =
{
 -5.129862866715009066465422805058933131960E1L,
  1.452015077564081884387441590064272782044E2L,
 -1.524043275549860505277434040464085593165E2L,
  7.236063513651544224319663428634139768808E1L,
 -1.494198912340228235853027849917095580053E1L
  /* 1.0E0 */
};

/* exp(x) = 1 + x - x / (1 - 2 / (x - x^2 R(x^2)))
   0 <= x <= 0.5
   Peak relative error 5.7e-38  */
static const long double PN[] =
{
  5.081801691915377692446852383385968225675E8L,
  9.360895299872484512023336636427675327355E6L,
  4.213701282274196030811629773097579432957E4L,
  5.201006511142748908655720086041570288182E1L,
  9.088368420359444263703202925095675982530E-3L,
};
static const long double PD[] =
{
  3.049081015149226615468111430031590411682E9L,
  1.069833887183886839966085436512368982758E8L,
  8.259257717868875207333991924545445705394E5L,
  1.872583833284143212651746812884298360922E3L,
  /* 1.0E0 */
};

static const long double
  /* ln 2 */
  lg2 = 6.9314718055994530941723212145817656807550E-1L,
  lg2_h = 6.9314718055994528622676398299518041312695E-1L,
  lg2_l = 2.3190468138462996154948554638754786504121E-17L,
  ovt = 8.0085662595372944372e-0017L,
  /* 2/(3*log(2)) */
  cp = 9.6179669392597560490661645400126142495110E-1L,
  cp_h = 9.6179669392597555432899980587535537779331E-1L,
  cp_l = 5.0577616648125906047157785230014751039424E-17L;

long double
__ieee754_powl (long double x, long double y)
{
  long double z, ax, z_h, z_l, p_h, p_l;
  long double y1, t1, t2, r, s, sgn, t, u, v, w;
  long double s2, s_h, s_l, t_h, t_l, ay;
  int32_t i, j, k, yisint, n;
  uint32_t ix, iy;
  int32_t hx, hy, hax;
  double ohi, xhi, xlo, yhi, ylo;
  uint32_t lx, ly, lj;

  ldbl_unpack (x, &xhi, &xlo);
  EXTRACT_WORDS (hx, lx, xhi);
  ix = hx & 0x7fffffff;

  ldbl_unpack (y, &yhi, &ylo);
  EXTRACT_WORDS (hy, ly, yhi);
  iy = hy & 0x7fffffff;

  /* y==zero: x**0 = 1 */
  if ((iy | ly) == 0 && !issignaling (x))
    return one;

  /* 1.0**y = 1; -1.0**+-Inf = 1 */
  if (x == one && !issignaling (y))
    return one;
  if (x == -1.0L && ((iy - 0x7ff00000) | ly) == 0)
    return one;

  /* +-NaN return x+y */
  if ((ix >= 0x7ff00000 && ((ix - 0x7ff00000) | lx) != 0)
      || (iy >= 0x7ff00000 && ((iy - 0x7ff00000) | ly) != 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)
    {
      uint32_t low_ye;

      GET_HIGH_WORD (low_ye, ylo);
      if ((low_ye & 0x7fffffff) >= 0x43400000)	/* Low part >= 2^53 */
	yisint = 2;		/* even integer y */
      else if (iy >= 0x3ff00000)	/* 1.0 */
	{
	  if (floorl (y) == y)
	    {
	      z = 0.5 * y;
	      if (floorl (z) == z)
		yisint = 2;
	      else
		yisint = 1;
	    }
	}
    }

  ax = fabsl (x);

  /* special value of y */
  if (ly == 0)
    {
      if (iy == 0x7ff00000)	/* y is +-inf */
	{
	  if (ax > one)
	    /* (|x|>1)**+-inf = inf,0 */
	    return (hy >= 0) ? y : zero;
	  else
	    /* (|x|<1)**-,+inf = inf,0 */
	    return (hy < 0) ? -y : zero;
	}
      if (ylo == 0.0)
	{
	  if (iy == 0x3ff00000)
	    {			/* y is  +-1 */
	      if (hy < 0)
		return one / x;
	      else
		return x;
	    }
	  if (hy == 0x40000000)
	    return x * x;		/* y is  2 */
	  if (hy == 0x3fe00000)
	    {			/* y is  0.5 */
	      if (hx >= 0)		/* x >= +0 */
		return sqrtl (x);
	    }
	}
    }

  /* special value of x */
  if (lx == 0)
    {
      if (ix == 0x7ff00000 || ix == 0 || (ix == 0x3ff00000 && xlo == 0.0))
	{
	  z = ax;		/*x is +-0,+-inf,+-1 */
	  if (hy < 0)
	    z = one / z;	/* z = (1/|x|) */
	  if (hx < 0)
	    {
	      if (((ix - 0x3ff00000) | 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 (((((uint32_t) hx >> 31) - 1) | yisint) == 0)
    return (x - x) / (x - x);

  /* sgn (sign of result -ve**odd) = -1 else = 1 */
  sgn = one;
  if (((((uint32_t) hx >> 31) - 1) | (yisint - 1)) == 0)
    sgn = -one;			/* (-ve)**(odd int) */

  /* |y| is huge.
     2^-16495 = 1/2 of smallest representable value.
     If (1 - 1/131072)^y underflows, y > 1.4986e9 */
  if (iy > 0x41d654b0)
    {
      /* if (1 - 2^-113)^y underflows, y > 1.1873e38 */
      if (iy > 0x47d654b0)
	{
	  if (ix <= 0x3fefffff)
	    return (hy < 0) ? sgn * huge * huge : sgn * tiny * tiny;
	  if (ix >= 0x3ff00000)
	    return (hy > 0) ? sgn * huge * huge : sgn * tiny * tiny;
	}
      /* over/underflow if x is not close to one */
      if (ix < 0x3fefffff)
	return (hy < 0) ? sgn * huge * huge : sgn * tiny * tiny;
      if (ix > 0x3ff00000)
	return (hy > 0) ? sgn * huge * huge : sgn * tiny * tiny;
    }

  ay = y > 0 ? y : -y;
  if (ay < 0x1p-117)
    y = y < 0 ? -0x1p-117 : 0x1p-117;

  n = 0;
  /* take care subnormal number */
  if (ix < 0x00100000)
    {
      ax *= two113;
      n -= 113;
      ohi = ldbl_high (ax);
      GET_HIGH_WORD (ix, ohi);
    }
  n += ((ix) >> 20) - 0x3ff;
  j = ix & 0x000fffff;
  /* determine interval */
  ix = j | 0x3ff00000;		/* normalize ix */
  if (j <= 0x39880)
    k = 0;			/* |x|<sqrt(3/2) */
  else if (j < 0xbb670)
    k = 1;			/* |x|<sqrt(3)   */
  else
    {
      k = 0;
      n += 1;
      ix -= 0x00100000;
    }

  ohi = ldbl_high (ax);
  GET_HIGH_WORD (hax, ohi);
  ax = __scalbnl (ax, ((int) ((ix - hax) * 2)) >> 21);

  /* 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 = ldbl_high (s);

  /* t_h=ax+bp[k] High */
  t_h = ax + bp[k];
  t_h = ldbl_high (t_h);
  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;
  t_h = ldbl_high (t_h);
  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;
  p_h = ldbl_high (p_h);
  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);
  t1 = ldbl_high (t1);
  t2 = z_l - (((t1 - t) - dp_h[k]) - z_h);

  /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */
  y1 = ldbl_high (y);
  p_l = (y - y1) * t1 + y * t2;
  p_h = y1 * t1;
  z = p_l + p_h;
  ohi = ldbl_high (z);
  EXTRACT_WORDS (j, lj, ohi);
  if (j >= 0x40d00000) /* z >= 16384 */
    {
      /* if z > 16384 */
      if (((j - 0x40d00000) | lj) != 0)
	return sgn * huge * huge;	/* overflow */
      else
	{
	  if (p_l + ovt > z - p_h)
	    return sgn * huge * huge;	/* overflow */
	}
    }
  else if ((j & 0x7fffffff) >= 0x40d01b90)	/* z <= -16495 */
    {
      /* z < -16495 */
      if (((j - 0xc0d01bc0) | lj) != 0)
	return sgn * tiny * tiny;	/* underflow */
      else
	{
	  if (p_l <= z - p_h)
	    return sgn * tiny * tiny;	/* underflow */
	}
    }
  /* compute 2**(p_h+p_l) */
  i = j & 0x7fffffff;
  k = (i >> 20) - 0x3ff;
  n = 0;
  if (i > 0x3fe00000)
    {				/* if |z| > 0.5, set n = [z+0.5] */
      n = floorl (z + 0.5L);
      t = n;
      p_h -= t;
    }
  t = p_l + p_h;
  t = ldbl_high (t);
  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);
  z = __scalbnl (sgn * z, n);
  math_check_force_underflow (z);
  return z;
}
libm_alias_finite (__ieee754_powl, __powl)
Commit	Line	Data
f964490f RM	1	/*
	2	* ====================================================
	3	* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
	4	*
	5	* Developed at SunPro, a Sun Microsystems, Inc. business.
	6	* Permission to use, copy, modify, and distribute this
	7	* software is freely granted, provided that this notice
	8	* is preserved.
	9	* ====================================================
	10	*/
	11
	12	/* Expansions and modifications for 128-bit long double are
	13	Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov>
	14	and are incorporated herein by permission of the author. The author
	15	reserves the right to distribute this material elsewhere under different
	16	copying permissions. These modifications are distributed here under
	17	the following terms:
	18
	19	This library is free software; you can redistribute it and/or
	20	modify it under the terms of the GNU Lesser General Public
	21	License as published by the Free Software Foundation; either
	22	version 2.1 of the License, or (at your option) any later version.
	23
	24	This library is distributed in the hope that it will be useful,
	25	but WITHOUT ANY WARRANTY; without even the implied warranty of
	26	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
	27	Lesser General Public License for more details.
	28
	29	You should have received a copy of the GNU Lesser General Public
59ba27a6	30	License along with this library; if not, see
5a82c748	31	<https://www.gnu.org/licenses/>. */
f964490f RM	32
	33	/* __ieee754_powl(x,y) return x**y
	34	*
	35	* n
	36	* Method: Let x = 2 * (1+f)
	37	* 1. Compute and return log2(x) in two pieces:
	38	* log2(x) = w1 + w2,
	39	* where w1 has 113-53 = 60 bit trailing zeros.
	40	* 2. Perform y*log2(x) = n+y' by simulating muti-precision
	41	* arithmetic, where \|y'\|<=0.5.
	42	* 3. Return xy = 2nexp(y'log2)
	43	*
	44	* Special cases:
	45	* 1. (anything) ** 0 is 1
	46	* 2. (anything) ** 1 is itself
	47	* 3. (anything) ** NAN is NAN
	48	* 4. NAN ** (anything except 0) is NAN
	49	* 5. +-(\|x\| > 1) ** +INF is +INF
	50	* 6. +-(\|x\| > 1) ** -INF is +0
	51	* 7. +-(\|x\| < 1) ** +INF is +0
	52	* 8. +-(\|x\| < 1) ** -INF is +INF
	53	* 9. +-1 ** +-INF is NAN
	54	* 10. +0 ** (+anything except 0, NAN) is +0
	55	* 11. -0 ** (+anything except 0, NAN, odd integer) is +0
	56	* 12. +0 ** (-anything except 0, NAN) is +INF
	57	* 13. -0 ** (-anything except 0, NAN, odd integer) is +INF
	58	* 14. -0 (odd integer) = -( +0 (odd integer) )
	59	* 15. +INF ** (+anything except 0,NAN) is +INF
	60	* 16. +INF ** (-anything except 0,NAN) is +0
	61	* 17. -INF (anything) = -0 (-anything)
	62	* 18. (-anything) (integer) is (-1)(integer)(+anything*integer)
	63	* 19. (-anything except 0 and inf) ** (non-integer) is NAN
	64	*
	65	*/
	66
1ed0291c RH	67	#include <math.h>
1ed0291c RH	68	#include <math_private.h>
8f5b00d3	69	#include <math-underflow.h>
220622dd	70	#include <libm-alias-finite.h>
f964490f RM	71
	72	static const long double bp[] = {
	73	1.0L,
	74	1.5L,
	75	};
	76
	77	/* log_2(1.5) */
	78	static const long double dp_h[] = {
	79	0.0,
	80	5.8496250072115607565592654282227158546448E-1L
	81	};
	82
	83	/* Low part of log_2(1.5) */
	84	static const long double dp_l[] = {
	85	0.0,
	86	1.0579781240112554492329533686862998106046E-16L
	87	};
	88
	89	static const long double zero = 0.0L,
	90	one = 1.0L,
	91	two = 2.0L,
	92	two113 = 1.0384593717069655257060992658440192E34L,
dcb33988 AS	93	huge = 1.0e300L,
dcb33988 AS	94	tiny = 1.0e-300L;
f964490f RM	95
	96	/* 3/2 log x = 3 z + z^3 + z^3 (z^2 R(z^2))
	97	z = (x-1)/(x+1)
	98	1 <= x <= 1.25
	99	Peak relative error 2.3e-37 */
	100	static const long double LN[] =
	101	{
	102	-3.0779177200290054398792536829702930623200E1L,
	103	6.5135778082209159921251824580292116201640E1L,
	104	-4.6312921812152436921591152809994014413540E1L,
	105	1.2510208195629420304615674658258363295208E1L,
	106	-9.9266909031921425609179910128531667336670E-1L
	107	};
	108	static const long double LD[] =
	109	{
	110	-5.129862866715009066465422805058933131960E1L,
	111	1.452015077564081884387441590064272782044E2L,
	112	-1.524043275549860505277434040464085593165E2L,
	113	7.236063513651544224319663428634139768808E1L,
	114	-1.494198912340228235853027849917095580053E1L
	115	/* 1.0E0 */
	116	};
	117
	118	/* exp(x) = 1 + x - x / (1 - 2 / (x - x^2 R(x^2)))
	119	0 <= x <= 0.5
	120	Peak relative error 5.7e-38 */
	121	static const long double PN[] =
	122	{
	123	5.081801691915377692446852383385968225675E8L,
	124	9.360895299872484512023336636427675327355E6L,
	125	4.213701282274196030811629773097579432957E4L,
	126	5.201006511142748908655720086041570288182E1L,
	127	9.088368420359444263703202925095675982530E-3L,
	128	};
	129	static const long double PD[] =
	130	{
	131	3.049081015149226615468111430031590411682E9L,
	132	1.069833887183886839966085436512368982758E8L,
	133	8.259257717868875207333991924545445705394E5L,
	134	1.872583833284143212651746812884298360922E3L,
	135	/* 1.0E0 */
	136	};
	137
	138	static const long double
	139	/* ln 2 */
	140	lg2 = 6.9314718055994530941723212145817656807550E-1L,
	141	lg2_h = 6.9314718055994528622676398299518041312695E-1L,
	142	lg2_l = 2.3190468138462996154948554638754786504121E-17L,
	143	ovt = 8.0085662595372944372e-0017L,
	144	/* 2/(3log(2)) /
	145	cp = 9.6179669392597560490661645400126142495110E-1L,
	146	cp_h = 9.6179669392597555432899980587535537779331E-1L,
	147	cp_l = 5.0577616648125906047157785230014751039424E-17L;
	148
f964490f RM	149	long double
f964490f RM	150	__ieee754_powl (long double x, long double y)
f964490f RM	151	{
f964490f RM	152	long double z, ax, z_h, z_l, p_h, p_l;
c091488e	153	long double y1, t1, t2, r, s, sgn, t, u, v, w;
ef1e0867	154	long double s2, s_h, s_l, t_h, t_l, ay;
f964490f	155	int32_t i, j, k, yisint, n;
765714ca AM	156	uint32_t ix, iy;
	157	int32_t hx, hy, hax;
	158	double ohi, xhi, xlo, yhi, ylo;
	159	uint32_t lx, ly, lj;
f964490f	160
765714ca AM	161	ldbl_unpack (x, &xhi, &xlo);
765714ca AM	162	EXTRACT_WORDS (hx, lx, xhi);
f964490f RM	163	ix = hx & 0x7fffffff;
f964490f RM	164
765714ca AM	165	ldbl_unpack (y, &yhi, &ylo);
765714ca AM	166	EXTRACT_WORDS (hy, ly, yhi);
f964490f RM	167	iy = hy & 0x7fffffff;
f964490f RM	168
f964490f	169	/* y==zero: x*0 = 1 /
90ab295a	170	if ((iy \| ly) == 0 && !issignaling (x))
f964490f RM	171	return one;
	172
	173	/* 1.0y = 1; -1.0+-Inf = 1 */
90ab295a	174	if (x == one && !issignaling (y))
f964490f	175	return one;
765714ca	176	if (x == -1.0L && ((iy - 0x7ff00000) \| ly) == 0)
f964490f RM	177	return one;
	178
	179	/* +-NaN return x+y */
765714ca AM	180	if ((ix >= 0x7ff00000 && ((ix - 0x7ff00000) \| lx) != 0)
765714ca AM	181	\|\| (iy >= 0x7ff00000 && ((iy - 0x7ff00000) \| ly) != 0))
f964490f RM	182	return x + y;
	183
	184	/* determine if y is an odd int when x < 0
	185	* yisint = 0 ... y is not an integer
	186	* yisint = 1 ... y is an odd int
	187	* yisint = 2 ... y is an even int
	188	*/
	189	yisint = 0;
	190	if (hx < 0)
	191	{
765714ca AM	192	uint32_t low_ye;
	193
	194	GET_HIGH_WORD (low_ye, ylo);
	195	if ((low_ye & 0x7fffffff) >= 0x43400000) /* Low part >= 2^53 */
f964490f RM	196	yisint = 2; /* even integer y */
	197	else if (iy >= 0x3ff00000) /* 1.0 */
	198	{
e44acb20	199	if (floorl (y) == y)
f964490f RM	200	{
f964490f RM	201	z = 0.5 * y;
e44acb20	202	if (floorl (z) == z)
f964490f RM	203	yisint = 2;
	204	else
	205	yisint = 1;
	206	}
	207	}
	208	}
	209
765714ca AM	210	ax = fabsl (x);
765714ca AM	211
f964490f	212	/* special value of y */
765714ca	213	if (ly == 0)
f964490f	214	{
765714ca	215	if (iy == 0x7ff00000) /* y is +-inf */
f964490f	216	{
765714ca	217	if (ax > one)
f964490f RM	218	/* (\|x\|>1)*+-inf = inf,0 /
	219	return (hy >= 0) ? y : zero;
	220	else
	221	/* (\|x\|<1)*-,+inf = inf,0 /
	222	return (hy < 0) ? -y : zero;
	223	}
765714ca AM	224	if (ylo == 0.0)
	225	{
	226	if (iy == 0x3ff00000)
	227	{ /* y is +-1 */
	228	if (hy < 0)
	229	return one / x;
	230	else
	231	return x;
	232	}
	233	if (hy == 0x40000000)
	234	return x * x; /* y is 2 */
	235	if (hy == 0x3fe00000)
	236	{ /* y is 0.5 */
	237	if (hx >= 0) /* x >= +0 */
f67a8147	238	return sqrtl (x);
765714ca	239	}
f964490f RM	240	}
	241	}
	242
f964490f	243	/* special value of x */
765714ca	244	if (lx == 0)
f964490f	245	{
765714ca	246	if (ix == 0x7ff00000 \|\| ix == 0 \|\| (ix == 0x3ff00000 && xlo == 0.0))
f964490f RM	247	{
	248	z = ax; /x is +-0,+-inf,+-1 /
	249	if (hy < 0)
	250	z = one / z; /* z = (1/\|x\|) */
	251	if (hx < 0)
	252	{
	253	if (((ix - 0x3ff00000) \| yisint) == 0)
	254	{
	255	z = (z - z) / (z - z); /* (-1)*non-int is NaN /
	256	}
	257	else if (yisint == 1)
	258	z = -z; /* (x<0)odd = -(\|x\|odd) */
	259	}
	260	return z;
	261	}
	262	}
	263
	264	/* (x<0)*(non-int) is NaN /
24ab7723	265	if (((((uint32_t) hx >> 31) - 1) \| yisint) == 0)
f964490f RM	266	return (x - x) / (x - x);
f964490f RM	267
c091488e JM	268	/* sgn (sign of result -ve*odd) = -1 else = 1 /
c091488e JM	269	sgn = one;
24ab7723	270	if (((((uint32_t) hx >> 31) - 1) \| (yisint - 1)) == 0)
c091488e JM	271	sgn = -one; /* (-ve)*(odd int) /
c091488e JM	272
f964490f RM	273	/* \|y\| is huge.
	274	2^-16495 = 1/2 of smallest representable value.
	275	If (1 - 1/131072)^y underflows, y > 1.4986e9 */
	276	if (iy > 0x41d654b0)
	277	{
	278	/* if (1 - 2^-113)^y underflows, y > 1.1873e38 */
	279	if (iy > 0x47d654b0)
	280	{
	281	if (ix <= 0x3fefffff)
c091488e	282	return (hy < 0) ? sgn * huge * huge : sgn * tiny * tiny;
f964490f	283	if (ix >= 0x3ff00000)
c091488e	284	return (hy > 0) ? sgn * huge * huge : sgn * tiny * tiny;
f964490f RM	285	}
	286	/* over/underflow if x is not close to one */
	287	if (ix < 0x3fefffff)
c091488e	288	return (hy < 0) ? sgn * huge * huge : sgn * tiny * tiny;
f964490f	289	if (ix > 0x3ff00000)
c091488e	290	return (hy > 0) ? sgn * huge * huge : sgn * tiny * tiny;
f964490f RM	291	}
f964490f RM	292
ef1e0867 JM	293	ay = y > 0 ? y : -y;
	294	if (ay < 0x1p-117)
	295	y = y < 0 ? -0x1p-117 : 0x1p-117;
	296
f964490f RM	297	n = 0;
	298	/* take care subnormal number */
	299	if (ix < 0x00100000)
	300	{
	301	ax *= two113;
	302	n -= 113;
765714ca AM	303	ohi = ldbl_high (ax);
765714ca AM	304	GET_HIGH_WORD (ix, ohi);
f964490f RM	305	}
	306	n += ((ix) >> 20) - 0x3ff;
	307	j = ix & 0x000fffff;
	308	/* determine interval */
	309	ix = j \| 0x3ff00000; /* normalize ix */
	310	if (j <= 0x39880)
	311	k = 0; /* \|x\|<sqrt(3/2) */
	312	else if (j < 0xbb670)
	313	k = 1; /* \|x\|<sqrt(3) */
	314	else
	315	{
	316	k = 0;
	317	n += 1;
	318	ix -= 0x00100000;
	319	}
	320
765714ca AM	321	ohi = ldbl_high (ax);
	322	GET_HIGH_WORD (hax, ohi);
	323	ax = __scalbnl (ax, ((int) ((ix - hax) * 2)) >> 21);
f964490f RM	324
	325	/* compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
	326	u = ax - bp[k]; /* bp[0]=1.0, bp[1]=1.5 */
	327	v = one / (ax + bp[k]);
	328	s = u * v;
765714ca	329	s_h = ldbl_high (s);
f964490f	330
f964490f RM	331	/* t_h=ax+bp[k] High */
f964490f RM	332	t_h = ax + bp[k];
765714ca	333	t_h = ldbl_high (t_h);
f964490f RM	334	t_l = ax - (t_h - bp[k]);
	335	s_l = v * ((u - s_h * t_h) - s_h * t_l);
	336	/* compute log(ax) */
	337	s2 = s * s;
	338	u = LN[0] + s2 * (LN[1] + s2 * (LN[2] + s2 * (LN[3] + s2 * LN[4])));
	339	v = LD[0] + s2 * (LD[1] + s2 * (LD[2] + s2 * (LD[3] + s2 * (LD[4] + s2))));
	340	r = s2 * s2 * u / v;
	341	r += s_l * (s_h + s);
	342	s2 = s_h * s_h;
	343	t_h = 3.0 + s2 + r;
765714ca	344	t_h = ldbl_high (t_h);
f964490f RM	345	t_l = r - ((t_h - 3.0) - s2);
	346	/* u+v = s(1+...) /
	347	u = s_h * t_h;
	348	v = s_l * t_h + t_l * s;
	349	/* 2/(3log2)(s+...) /
	350	p_h = u + v;
765714ca	351	p_h = ldbl_high (p_h);
f964490f RM	352	p_l = v - (p_h - u);
	353	z_h = cp_h * p_h; /* cp_h+cp_l = 2/(3log2) /
	354	z_l = cp_l * p_h + p_l * cp + dp_l[k];
	355	/* log2(ax) = (s+..)2/(3log2) = n + dp_h + z_h + z_l */
	356	t = (long double) n;
	357	t1 = (((z_h + z_l) + dp_h[k]) + t);
765714ca	358	t1 = ldbl_high (t1);
f964490f RM	359	t2 = z_l - (((t1 - t) - dp_h[k]) - z_h);
f964490f RM	360
f964490f	361	/* split up y into y1+y2 and compute (y1+y2)(t1+t2) /
765714ca	362	y1 = ldbl_high (y);
f964490f RM	363	p_l = (y - y1) * t1 + y * t2;
	364	p_h = y1 * t1;
	365	z = p_l + p_h;
765714ca AM	366	ohi = ldbl_high (z);
765714ca AM	367	EXTRACT_WORDS (j, lj, ohi);
f964490f RM	368	if (j >= 0x40d00000) /* z >= 16384 */
	369	{
	370	/* if z > 16384 */
765714ca	371	if (((j - 0x40d00000) \| lj) != 0)
c091488e	372	return sgn * huge * huge; /* overflow */
f964490f RM	373	else
	374	{
	375	if (p_l + ovt > z - p_h)
c091488e	376	return sgn * huge * huge; /* overflow */
f964490f RM	377	}
	378	}
	379	else if ((j & 0x7fffffff) >= 0x40d01b90) /* z <= -16495 */
	380	{
	381	/* z < -16495 */
765714ca	382	if (((j - 0xc0d01bc0) \| lj) != 0)
c091488e	383	return sgn * tiny * tiny; /* underflow */
f964490f RM	384	else
	385	{
	386	if (p_l <= z - p_h)
c091488e	387	return sgn * tiny * tiny; /* underflow */
f964490f RM	388	}
	389	}
	390	/* compute 2*(p_h+p_l) /
	391	i = j & 0x7fffffff;
	392	k = (i >> 20) - 0x3ff;
	393	n = 0;
	394	if (i > 0x3fe00000)
	395	{ /* if \|z\| > 0.5, set n = [z+0.5] */
e44acb20	396	n = floorl (z + 0.5L);
f964490f RM	397	t = n;
	398	p_h -= t;
	399	}
	400	t = p_l + p_h;
765714ca	401	t = ldbl_high (t);
f964490f RM	402	u = t * lg2_h;
	403	v = (p_l - (t - p_h)) * lg2 + t * lg2_l;
	404	z = u + v;
	405	w = v - (z - u);
	406	/* exp(z) */
	407	t = z * z;
	408	u = PN[0] + t * (PN[1] + t * (PN[2] + t * (PN[3] + t * PN[4])));
	409	v = PD[0] + t * (PD[1] + t * (PD[2] + t * (PD[3] + t)));
	410	t1 = z - t * u / v;
	411	r = (z * t1) / (t1 - two) - (w + z * w);
	412	z = one - (r - z);
c091488e JM	413	z = __scalbnl (sgn * z, n);
	414	math_check_force_underflow (z);
	415	return z;
f964490f	416	}
220622dd	417	libm_alias_finite (__ieee754_powl, __powl)