[mirror_edk2.git] / StdLib / LibC / Math / e_pow.c

/** @file\r
  Compute the base 10 logrithm of x.\r
\r
  Copyright (c) 2010 - 2011, Intel Corporation. All rights reserved.<BR>\r
  This program and the accompanying materials are licensed and made available under\r
  the terms and conditions of the BSD License that accompanies this distribution.\r
  The full text of the license may be found at\r
  http://opensource.org/licenses/bsd-license.\r
\r
  THE PROGRAM IS DISTRIBUTED UNDER THE BSD LICENSE ON AN "AS IS" BASIS,\r
  WITHOUT WARRANTIES OR REPRESENTATIONS OF ANY KIND, EITHER EXPRESS OR IMPLIED.\r
\r
 * ====================================================\r
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.\r
 *\r
 * Developed at SunPro, a Sun Microsystems, Inc. business.\r
 * Permission to use, copy, modify, and distribute this\r
 * software is freely granted, provided that this notice\r
 * is preserved.\r
 * ====================================================\r
\r
  e_pow.c 5.1 93/09/24\r
  NetBSD: e_pow.c,v 1.13 2004/06/30 18:43:15 drochner Exp\r
**/\r
#include  <LibConfig.h>\r
#include  <sys/EfiCdefs.h>\r
\r
#if defined(_MSC_VER)           /* Handle Microsoft VC++ compiler specifics. */\r
  // C4723: potential divide by zero.\r
  #pragma warning ( disable : 4723 )\r
  // C4756: overflow in constant arithmetic\r
  #pragma warning ( disable : 4756 )\r
#endif\r
\r
/* __ieee754_pow(x,y) return x**y\r
 *\r
 *          n\r
 * Method:  Let x =  2   * (1+f)\r
 *  1. Compute and return log2(x) in two pieces:\r
 *    log2(x) = w1 + w2,\r
 *     where w1 has 53-24 = 29 bit trailing zeros.\r
 *  2. Perform y*log2(x) = n+y' by simulating multi-precision\r
 *     arithmetic, where |y'|<=0.5.\r
 *  3. Return x**y = 2**n*exp(y'*log2)\r
 *\r
 * Special cases:\r
 *  1.  (anything) ** 0  is 1\r
 *  2.  (anything) ** 1  is itself\r
 *  3.  (anything) ** NAN is NAN\r
 *  4.  NAN ** (anything except 0) is NAN\r
 *  5.  +-(|x| > 1) **  +INF is +INF\r
 *  6.  +-(|x| > 1) **  -INF is +0\r
 *  7.  +-(|x| < 1) **  +INF is +0\r
 *  8.  +-(|x| < 1) **  -INF is +INF\r
 *  9.  +-1         ** +-INF is NAN\r
 *  10. +0 ** (+anything except 0, NAN)               is +0\r
 *  11. -0 ** (+anything except 0, NAN, odd integer)  is +0\r
 *  12. +0 ** (-anything except 0, NAN)               is +INF\r
 *  13. -0 ** (-anything except 0, NAN, odd integer)  is +INF\r
 *  14. -0 ** (odd integer) = -( +0 ** (odd integer) )\r
 *  15. +INF ** (+anything except 0,NAN) is +INF\r
 *  16. +INF ** (-anything except 0,NAN) is +0\r
 *  17. -INF ** (anything)  = -0 ** (-anything)\r
 *  18. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer)\r
 *  19. (-anything except 0 and inf) ** (non-integer) is NAN\r
 *\r
 * Accuracy:\r
 *  pow(x,y) returns x**y nearly rounded. In particular\r
 *      pow(integer,integer)\r
 *  always returns the correct integer provided it is\r
 *  representable.\r
 *\r
 * Constants :\r
 * The hexadecimal values are the intended ones for the following\r
 * constants. The decimal values may be used, provided that the\r
 * compiler will convert from decimal to binary accurately enough\r
 * to produce the hexadecimal values shown.\r
 */\r
\r
#include "math.h"\r
#include "math_private.h"\r
#include  <errno.h>\r
\r
static const double\r
bp[] = {1.0, 1.5,},\r
dp_h[] = { 0.0, 5.84962487220764160156e-01,}, /* 0x3FE2B803, 0x40000000 */\r
dp_l[] = { 0.0, 1.35003920212974897128e-08,}, /* 0x3E4CFDEB, 0x43CFD006 */\r
zero    =  0.0,\r
one =  1.0,\r
two =  2.0,\r
two53 =  9007199254740992.0,  /* 0x43400000, 0x00000000 */\r
huge  =  1.0e300,\r
tiny    =  1.0e-300,\r
  /* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */\r
L1  =  5.99999999999994648725e-01, /* 0x3FE33333, 0x33333303 */\r
L2  =  4.28571428578550184252e-01, /* 0x3FDB6DB6, 0xDB6FABFF */\r
L3  =  3.33333329818377432918e-01, /* 0x3FD55555, 0x518F264D */\r
L4  =  2.72728123808534006489e-01, /* 0x3FD17460, 0xA91D4101 */\r
L5  =  2.30660745775561754067e-01, /* 0x3FCD864A, 0x93C9DB65 */\r
L6  =  2.06975017800338417784e-01, /* 0x3FCA7E28, 0x4A454EEF */\r
P1   =  1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */\r
P2   = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */\r
P3   =  6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */\r
P4   = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */\r
P5   =  4.13813679705723846039e-08, /* 0x3E663769, 0x72BEA4D0 */\r
lg2  =  6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */\r
lg2_h  =  6.93147182464599609375e-01, /* 0x3FE62E43, 0x00000000 */\r
lg2_l  = -1.90465429995776804525e-09, /* 0xBE205C61, 0x0CA86C39 */\r
ovt =  8.0085662595372944372e-0017, /* -(1024-log2(ovfl+.5ulp)) */\r
cp    =  9.61796693925975554329e-01, /* 0x3FEEC709, 0xDC3A03FD =2/(3ln2) */\r
cp_h  =  9.61796700954437255859e-01, /* 0x3FEEC709, 0xE0000000 =(float)cp */\r
cp_l  = -7.02846165095275826516e-09, /* 0xBE3E2FE0, 0x145B01F5 =tail of cp_h*/\r
ivln2    =  1.44269504088896338700e+00, /* 0x3FF71547, 0x652B82FE =1/ln2 */\r
ivln2_h  =  1.44269502162933349609e+00, /* 0x3FF71547, 0x60000000 =24b 1/ln2*/\r
ivln2_l  =  1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/\r
\r
double\r
__ieee754_pow(double x, double y)\r
{\r
  double z,ax,z_h,z_l,p_h,p_l;\r
  double y1,t1,t2,r,s,t,u,v,w;\r
  int32_t i,j,k,yisint,n;\r
  int32_t hx,hy,ix,iy;\r
  u_int32_t lx,ly;\r
\r
  EXTRACT_WORDS(hx,lx,x);\r
  EXTRACT_WORDS(hy,ly,y);\r
  ix = hx&0x7fffffff;  iy = hy&0x7fffffff;\r
\r
    /* y==zero: x**0 = 1 */\r
  if((iy|ly)==0) return one;\r
\r
    /* +-NaN return x+y */\r
  if(ix > 0x7ff00000 || ((ix==0x7ff00000)&&(lx!=0)) ||\r
     iy > 0x7ff00000 || ((iy==0x7ff00000)&&(ly!=0)))\r
    return x+y;\r
\r
    /* determine if y is an odd int when x < 0\r
     * yisint = 0 ... y is not an integer\r
     * yisint = 1 ... y is an odd int\r
     * yisint = 2 ... y is an even int\r
     */\r
  yisint  = 0;\r
  if(hx<0) {\r
      if(iy>=0x43400000) yisint = 2; /* even integer y */\r
      else if(iy>=0x3ff00000) {\r
    k = (iy>>20)-0x3ff;    /* exponent */\r
    if(k>20) {\r
        j = ly>>(52-k);\r
        if((u_int32_t)(j<<(52-k))==ly) yisint = 2-(j&1);\r
    } else if(ly==0) {\r
        j = iy>>(20-k);\r
        if((j<<(20-k))==iy) yisint = 2-(j&1);\r
    }\r
      }\r
  }\r
\r
    /* special value of y */\r
  if(ly==0) {\r
      if (iy==0x7ff00000) { /* y is +-inf */\r
          if(((ix-0x3ff00000)|lx)==0)\r
        return  y - y;  /* inf**+-1 is NaN */\r
          else if (ix >= 0x3ff00000)/* (|x|>1)**+-inf = inf,0 */\r
        return (hy>=0)? y: zero;\r
          else      /* (|x|<1)**-,+inf = inf,0 */\r
        return (hy<0)?-y: zero;\r
      }\r
      if(iy==0x3ff00000) {  /* y is  +-1 */\r
    if(hy<0) return one/x; else return x;\r
      }\r
      if(hy==0x40000000) return x*x; /* y is  2 */\r
      if(hy==0x3fe00000) {  /* y is  0.5 */\r
    if(hx>=0) /* x >= +0 */\r
    return __ieee754_sqrt(x);\r
      }\r
  }\r
\r
  ax   = fabs(x);\r
    /* special value of x */\r
  if(lx==0) {\r
      if(ix==0x7ff00000||ix==0||ix==0x3ff00000){\r
    z = ax;     /*x is +-0,+-inf,+-1*/\r
    if(hy<0) z = one/z; /* z = (1/|x|) */\r
    if(hx<0) {\r
        if(((ix-0x3ff00000)|yisint)==0) {\r
      z = (z-z)/(z-z); /* (-1)**non-int is NaN */\r
        } else if(yisint==1)\r
      z = -z;   /* (x<0)**odd = -(|x|**odd) */\r
    }\r
    return z;\r
      }\r
  }\r
\r
  n = (hx>>31)+1;\r
\r
    /* (x<0)**(non-int) is NaN */\r
    if((n|yisint)==0) {\r
      errno = EDOM;\r
      return (x-x)/(x-x);\r
    }\r
\r
  s = one; /* s (sign of result -ve**odd) = -1 else = 1 */\r
  if((n|(yisint-1))==0) s = -one;/* (-ve)**(odd int) */\r
\r
    /* |y| is huge */\r
  if(iy>0x41e00000) { /* if |y| > 2**31 */\r
      if(iy>0x43f00000){  /* if |y| > 2**64, must o/uflow */\r
    if(ix<=0x3fefffff) return (hy<0)? huge*huge:tiny*tiny;\r
    if(ix>=0x3ff00000) return (hy>0)? huge*huge:tiny*tiny;\r
      }\r
  /* over/underflow if x is not close to one */\r
      if(ix<0x3fefffff) return (hy<0)? s*huge*huge:s*tiny*tiny;\r
      if(ix>0x3ff00000) return (hy>0)? s*huge*huge:s*tiny*tiny;\r
  /* now |1-x| is tiny <= 2**-20, suffice to compute\r
     log(x) by x-x^2/2+x^3/3-x^4/4 */\r
      t = ax-one;   /* t has 20 trailing zeros */\r
      w = (t*t)*(0.5-t*(0.3333333333333333333333-t*0.25));\r
      u = ivln2_h*t;  /* ivln2_h has 21 sig. bits */\r
      v = t*ivln2_l-w*ivln2;\r
      t1 = u+v;\r
      SET_LOW_WORD(t1,0);\r
      t2 = v-(t1-u);\r
  } else {\r
      double ss,s2,s_h,s_l,t_h,t_l;\r
      n = 0;\r
  /* take care subnormal number */\r
      if(ix<0x00100000)\r
    {ax *= two53; n -= 53; GET_HIGH_WORD(ix,ax); }\r
      n  += ((ix)>>20)-0x3ff;\r
      j  = ix&0x000fffff;\r
  /* determine interval */\r
      ix = j|0x3ff00000;    /* normalize ix */\r
      if(j<=0x3988E) k=0;   /* |x|<sqrt(3/2) */\r
      else if(j<0xBB67A) k=1; /* |x|<sqrt(3)   */\r
      else {k=0;n+=1;ix -= 0x00100000;}\r
      SET_HIGH_WORD(ax,ix);\r
\r
  /* compute ss = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */\r
      u = ax-bp[k];   /* bp[0]=1.0, bp[1]=1.5 */\r
      v = one/(ax+bp[k]);\r
      ss = u*v;\r
      s_h = ss;\r
      SET_LOW_WORD(s_h,0);\r
  /* t_h=ax+bp[k] High */\r
      t_h = zero;\r
      SET_HIGH_WORD(t_h,((ix>>1)|0x20000000)+0x00080000+(k<<18));\r
      t_l = ax - (t_h-bp[k]);\r
      s_l = v*((u-s_h*t_h)-s_h*t_l);\r
  /* compute log(ax) */\r
      s2 = ss*ss;\r
      r = s2*s2*(L1+s2*(L2+s2*(L3+s2*(L4+s2*(L5+s2*L6)))));\r
      r += s_l*(s_h+ss);\r
      s2  = s_h*s_h;\r
      t_h = 3.0+s2+r;\r
      SET_LOW_WORD(t_h,0);\r
      t_l = r-((t_h-3.0)-s2);\r
  /* u+v = ss*(1+...) */\r
      u = s_h*t_h;\r
      v = s_l*t_h+t_l*ss;\r
  /* 2/(3log2)*(ss+...) */\r
      p_h = u+v;\r
      SET_LOW_WORD(p_h,0);\r
      p_l = v-(p_h-u);\r
      z_h = cp_h*p_h;   /* cp_h+cp_l = 2/(3*log2) */\r
      z_l = cp_l*p_h+p_l*cp+dp_l[k];\r
  /* log2(ax) = (ss+..)*2/(3*log2) = n + dp_h + z_h + z_l */\r
      t = (double)n;\r
      t1 = (((z_h+z_l)+dp_h[k])+t);\r
      SET_LOW_WORD(t1,0);\r
      t2 = z_l-(((t1-t)-dp_h[k])-z_h);\r
  }\r
\r
    /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */\r
  y1  = y;\r
  SET_LOW_WORD(y1,0);\r
  p_l = (y-y1)*t1+y*t2;\r
  p_h = y1*t1;\r
  z = p_l+p_h;\r
  EXTRACT_WORDS(j,i,z);\r
  if (j>=0x40900000) {        /* z >= 1024 */\r
      if(((j-0x40900000)|i)!=0)     /* if z > 1024 */\r
    return s*huge*huge;     /* overflow */\r
      else {\r
    if(p_l+ovt>z-p_h) return s*huge*huge; /* overflow */\r
      }\r
  } else if((j&0x7fffffff)>=0x4090cc00 ) {  /* z <= -1075 */\r
      if(((j-0xc090cc00)|i)!=0)     /* z < -1075 */\r
    return s*tiny*tiny;   /* underflow */\r
      else {\r
    if(p_l<=z-p_h) return s*tiny*tiny;  /* underflow */\r
      }\r
  }\r
    /*\r
     * compute 2**(p_h+p_l)\r
     */\r
  i = j&0x7fffffff;\r
  k = (i>>20)-0x3ff;\r
  n = 0;\r
  if(i>0x3fe00000) {    /* if |z| > 0.5, set n = [z+0.5] */\r
      n = j+(0x00100000>>(k+1));\r
      k = ((n&0x7fffffff)>>20)-0x3ff; /* new k for n */\r
      t = zero;\r
      SET_HIGH_WORD(t,n&~(0x000fffff>>k));\r
      n = ((n&0x000fffff)|0x00100000)>>(20-k);\r
      if(j<0) n = -n;\r
      p_h -= t;\r
  }\r
  t = p_l+p_h;\r
  SET_LOW_WORD(t,0);\r
  u = t*lg2_h;\r
  v = (p_l-(t-p_h))*lg2+t*lg2_l;\r
  z = u+v;\r
  w = v-(z-u);\r
  t  = z*z;\r
  t1  = z - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));\r
  r  = (z*t1)/(t1-two)-(w+z*w);\r
  z  = one-(r-z);\r
  GET_HIGH_WORD(j,z);\r
  j += (n<<20);\r
  if((j>>20)<=0) z = scalbn(z,n); /* subnormal output */\r
  else SET_HIGH_WORD(z,j);\r
  return s*z;\r
}\r
Commit	Line	Data
2aa62f2b	1	/** @file\r
	2	Compute the base 10 logrithm of x.\r
	3	\r
	4	Copyright (c) 2010 - 2011, Intel Corporation. All rights reserved.<BR>\r
	5	This program and the accompanying materials are licensed and made available under\r
	6	the terms and conditions of the BSD License that accompanies this distribution.\r
	7	The full text of the license may be found at\r
	8	http://opensource.org/licenses/bsd-license.\r
	9	\r
	10	THE PROGRAM IS DISTRIBUTED UNDER THE BSD LICENSE ON AN "AS IS" BASIS,\r
	11	WITHOUT WARRANTIES OR REPRESENTATIONS OF ANY KIND, EITHER EXPRESS OR IMPLIED.\r
	12	\r
	13	* ====================================================\r
	14	* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.\r
	15	*\r
	16	* Developed at SunPro, a Sun Microsystems, Inc. business.\r
	17	* Permission to use, copy, modify, and distribute this\r
	18	* software is freely granted, provided that this notice\r
	19	* is preserved.\r
	20	* ====================================================\r
	21	\r
	22	e_pow.c 5.1 93/09/24\r
	23	NetBSD: e_pow.c,v 1.13 2004/06/30 18:43:15 drochner Exp\r
	24	**/\r
	25	#include <LibConfig.h>\r
	26	#include <sys/EfiCdefs.h>\r
	27	\r
	28	#if defined(_MSC_VER) /* Handle Microsoft VC++ compiler specifics. */\r
	29	// C4723: potential divide by zero.\r
	30	#pragma warning ( disable : 4723 )\r
	31	// C4756: overflow in constant arithmetic\r
	32	#pragma warning ( disable : 4756 )\r
	33	#endif\r
	34	\r
	35	/* __ieee754_pow(x,y) return x**y\r
	36	*\r
	37	* n\r
	38	* Method: Let x = 2 * (1+f)\r
	39	* 1. Compute and return log2(x) in two pieces:\r
	40	* log2(x) = w1 + w2,\r
	41	* where w1 has 53-24 = 29 bit trailing zeros.\r
	42	* 2. Perform y*log2(x) = n+y' by simulating multi-precision\r
	43	* arithmetic, where \|y'\|<=0.5.\r
	44	* 3. Return xy = 2nexp(y'log2)\r
	45	*\r
	46	* Special cases:\r
	47	* 1. (anything) ** 0 is 1\r
	48	* 2. (anything) ** 1 is itself\r
	49	* 3. (anything) ** NAN is NAN\r
	50	* 4. NAN ** (anything except 0) is NAN\r
	51	* 5. +-(\|x\| > 1) ** +INF is +INF\r
	52	* 6. +-(\|x\| > 1) ** -INF is +0\r
	53	* 7. +-(\|x\| < 1) ** +INF is +0\r
	54	* 8. +-(\|x\| < 1) ** -INF is +INF\r
	55	* 9. +-1 ** +-INF is NAN\r
	56	* 10. +0 ** (+anything except 0, NAN) is +0\r
	57	* 11. -0 ** (+anything except 0, NAN, odd integer) is +0\r
	58	* 12. +0 ** (-anything except 0, NAN) is +INF\r
	59	* 13. -0 ** (-anything except 0, NAN, odd integer) is +INF\r
	60	* 14. -0 (odd integer) = -( +0 (odd integer) )\r
	61	* 15. +INF ** (+anything except 0,NAN) is +INF\r
	62	* 16. +INF ** (-anything except 0,NAN) is +0\r
	63	* 17. -INF (anything) = -0 (-anything)\r
	64	* 18. (-anything) (integer) is (-1)(integer)(+anything*integer)\r
65	* 19. (-anything except 0 and inf) ** (non-integer) is NAN\r
66	*\r
67	* Accuracy:\r
68	* pow(x,y) returns x**y nearly rounded. In particular\r
69	* pow(integer,integer)\r
70	* always returns the correct integer provided it is\r
71	* representable.\r
72	*\r
73	* Constants :\r
74	* The hexadecimal values are the intended ones for the following\r
75	* constants. The decimal values may be used, provided that the\r
76	* compiler will convert from decimal to binary accurately enough\r
77	* to produce the hexadecimal values shown.\r
78	*/\r
79	\r
80	#include "math.h"\r
81	#include "math_private.h"\r
82	#include <errno.h>\r
83	\r
84	static const double\r
85	bp[] = {1.0, 1.5,},\r
86	dp_h[] = { 0.0, 5.84962487220764160156e-01,}, /* 0x3FE2B803, 0x40000000 */\r
87	dp_l[] = { 0.0, 1.35003920212974897128e-08,}, /* 0x3E4CFDEB, 0x43CFD006 */\r
88	zero = 0.0,\r
89	one = 1.0,\r
90	two = 2.0,\r
91	two53 = 9007199254740992.0, /* 0x43400000, 0x00000000 */\r
92	huge = 1.0e300,\r
93	tiny = 1.0e-300,\r
94	/* poly coefs for (3/2)(log(x)-2s-2/3s*3 /\r
95	L1 = 5.99999999999994648725e-01, /* 0x3FE33333, 0x33333303 */\r
96	L2 = 4.28571428578550184252e-01, /* 0x3FDB6DB6, 0xDB6FABFF */\r
97	L3 = 3.33333329818377432918e-01, /* 0x3FD55555, 0x518F264D */\r
98	L4 = 2.72728123808534006489e-01, /* 0x3FD17460, 0xA91D4101 */\r
99	L5 = 2.30660745775561754067e-01, /* 0x3FCD864A, 0x93C9DB65 */\r
100	L6 = 2.06975017800338417784e-01, /* 0x3FCA7E28, 0x4A454EEF */\r
101	P1 = 1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */\r
102	P2 = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */\r
103	P3 = 6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */\r
104	P4 = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */\r
105	P5 = 4.13813679705723846039e-08, /* 0x3E663769, 0x72BEA4D0 */\r
106	lg2 = 6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */\r
107	lg2_h = 6.93147182464599609375e-01, /* 0x3FE62E43, 0x00000000 */\r
108	lg2_l = -1.90465429995776804525e-09, /* 0xBE205C61, 0x0CA86C39 */\r
109	ovt = 8.0085662595372944372e-0017, /* -(1024-log2(ovfl+.5ulp)) */\r
110	cp = 9.61796693925975554329e-01, /* 0x3FEEC709, 0xDC3A03FD =2/(3ln2) */\r
111	cp_h = 9.61796700954437255859e-01, /* 0x3FEEC709, 0xE0000000 =(float)cp */\r
112	cp_l = -7.02846165095275826516e-09, /* 0xBE3E2FE0, 0x145B01F5 =tail of cp_h*/\r
113	ivln2 = 1.44269504088896338700e+00, /* 0x3FF71547, 0x652B82FE =1/ln2 */\r
114	ivln2_h = 1.44269502162933349609e+00, /* 0x3FF71547, 0x60000000 =24b 1/ln2*/\r
115	ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/\r
116	\r
117	double\r
118	__ieee754_pow(double x, double y)\r
119	{\r
120	double z,ax,z_h,z_l,p_h,p_l;\r
121	double y1,t1,t2,r,s,t,u,v,w;\r
122	int32_t i,j,k,yisint,n;\r
123	int32_t hx,hy,ix,iy;\r
124	u_int32_t lx,ly;\r
125	\r
126	EXTRACT_WORDS(hx,lx,x);\r
127	EXTRACT_WORDS(hy,ly,y);\r
128	ix = hx&0x7fffffff; iy = hy&0x7fffffff;\r
129	\r
130	/* y==zero: x*0 = 1 /\r
131	if((iy\|ly)==0) return one;\r
132	\r
133	/* +-NaN return x+y */\r
134	if(ix > 0x7ff00000 \|\| ((ix==0x7ff00000)&&(lx!=0)) \|\|\r
135	iy > 0x7ff00000 \|\| ((iy==0x7ff00000)&&(ly!=0)))\r
136	return x+y;\r
137	\r
138	/* determine if y is an odd int when x < 0\r
139	* yisint = 0 ... y is not an integer\r
140	* yisint = 1 ... y is an odd int\r
141	* yisint = 2 ... y is an even int\r
142	*/\r
143	yisint = 0;\r
144	if(hx<0) {\r
145	if(iy>=0x43400000) yisint = 2; /* even integer y */\r
146	else if(iy>=0x3ff00000) {\r
147	k = (iy>>20)-0x3ff; /* exponent */\r
148	if(k>20) {\r
149	j = ly>>(52-k);\r
150	if((u_int32_t)(j<<(52-k))==ly) yisint = 2-(j&1);\r
151	} else if(ly==0) {\r
152	j = iy>>(20-k);\r
153	if((j<<(20-k))==iy) yisint = 2-(j&1);\r
154	}\r
155	}\r
156	}\r
157	\r
158	/* special value of y */\r
159	if(ly==0) {\r
160	if (iy==0x7ff00000) { /* y is +-inf */\r
161	if(((ix-0x3ff00000)\|lx)==0)\r
162	return y - y; /* inf*+-1 is NaN /\r
163	else if (ix >= 0x3ff00000)/* (\|x\|>1)*+-inf = inf,0 /\r
164	return (hy>=0)? y: zero;\r
165	else /* (\|x\|<1)*-,+inf = inf,0 /\r
166	return (hy<0)?-y: zero;\r
167	}\r
168	if(iy==0x3ff00000) { /* y is +-1 */\r
169	if(hy<0) return one/x; else return x;\r
170	}\r
171	if(hy==0x40000000) return xx; / y is 2 */\r
172	if(hy==0x3fe00000) { /* y is 0.5 */\r
173	if(hx>=0) /* x >= +0 */\r
174	return __ieee754_sqrt(x);\r
175	}\r
176	}\r
177	\r
178	ax = fabs(x);\r
179	/* special value of x */\r
180	if(lx==0) {\r
181	if(ix==0x7ff00000\|\|ix==0\|\|ix==0x3ff00000){\r
182	z = ax; /x is +-0,+-inf,+-1/\r
183	if(hy<0) z = one/z; /* z = (1/\|x\|) */\r
184	if(hx<0) {\r
185	if(((ix-0x3ff00000)\|yisint)==0) {\r
186	z = (z-z)/(z-z); /* (-1)*non-int is NaN /\r
187	} else if(yisint==1)\r
188	z = -z; /* (x<0)odd = -(\|x\|odd) */\r
189	}\r
190	return z;\r
191	}\r
192	}\r
193	\r
194	n = (hx>>31)+1;\r
195	\r
196	/* (x<0)*(non-int) is NaN /\r
197	if((n\|yisint)==0) {\r
198	errno = EDOM;\r
199	return (x-x)/(x-x);\r
200	}\r
201	\r
202	s = one; /* s (sign of result -ve*odd) = -1 else = 1 /\r
203	if((n\|(yisint-1))==0) s = -one;/* (-ve)*(odd int) /\r
204	\r
205	/* \|y\| is huge */\r
206	if(iy>0x41e00000) { /* if \|y\| > 2*31 /\r
207	if(iy>0x43f00000){ /* if \|y\| > 2*64, must o/uflow /\r
208	if(ix<=0x3fefffff) return (hy<0)? hugehuge:tinytiny;\r
209	if(ix>=0x3ff00000) return (hy>0)? hugehuge:tinytiny;\r
210	}\r
211	/* over/underflow if x is not close to one */\r
212	if(ix<0x3fefffff) return (hy<0)? shugehuge:stinytiny;\r
213	if(ix>0x3ff00000) return (hy>0)? shugehuge:stinytiny;\r
214	/* now \|1-x\| is tiny <= 2**-20, suffice to compute\r
215	log(x) by x-x^2/2+x^3/3-x^4/4 */\r
216	t = ax-one; /* t has 20 trailing zeros */\r
217	w = (tt)(0.5-t(0.3333333333333333333333-t0.25));\r
218	u = ivln2_ht; / ivln2_h has 21 sig. bits */\r
219	v = tivln2_l-wivln2;\r
220	t1 = u+v;\r
221	SET_LOW_WORD(t1,0);\r
222	t2 = v-(t1-u);\r
223	} else {\r
224	double ss,s2,s_h,s_l,t_h,t_l;\r
225	n = 0;\r
226	/* take care subnormal number */\r
227	if(ix<0x00100000)\r
228	{ax *= two53; n -= 53; GET_HIGH_WORD(ix,ax); }\r
229	n += ((ix)>>20)-0x3ff;\r
230	j = ix&0x000fffff;\r
231	/* determine interval */\r
232	ix = j\|0x3ff00000; /* normalize ix */\r
233	if(j<=0x3988E) k=0; /* \|x\|<sqrt(3/2) */\r
234	else if(j<0xBB67A) k=1; /* \|x\|<sqrt(3) */\r
235	else {k=0;n+=1;ix -= 0x00100000;}\r
236	SET_HIGH_WORD(ax,ix);\r
237	\r
238	/* compute ss = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */\r
239	u = ax-bp[k]; /* bp[0]=1.0, bp[1]=1.5 */\r
240	v = one/(ax+bp[k]);\r
241	ss = u*v;\r
242	s_h = ss;\r
243	SET_LOW_WORD(s_h,0);\r
244	/* t_h=ax+bp[k] High */\r
245	t_h = zero;\r
246	SET_HIGH_WORD(t_h,((ix>>1)\|0x20000000)+0x00080000+(k<<18));\r
247	t_l = ax - (t_h-bp[k]);\r
248	s_l = v((u-s_ht_h)-s_h*t_l);\r
249	/* compute log(ax) */\r
250	s2 = ss*ss;\r
251	r = s2s2(L1+s2(L2+s2(L3+s2(L4+s2(L5+s2*L6)))));\r
252	r += s_l*(s_h+ss);\r
253	s2 = s_h*s_h;\r
254	t_h = 3.0+s2+r;\r
255	SET_LOW_WORD(t_h,0);\r
256	t_l = r-((t_h-3.0)-s2);\r
257	/* u+v = ss(1+...) /\r
258	u = s_h*t_h;\r
259	v = s_lt_h+t_lss;\r
260	/* 2/(3log2)(ss+...) /\r
261	p_h = u+v;\r
262	SET_LOW_WORD(p_h,0);\r
263	p_l = v-(p_h-u);\r
264	z_h = cp_hp_h; / cp_h+cp_l = 2/(3log2) /\r
265	z_l = cp_lp_h+p_lcp+dp_l[k];\r
266	/* log2(ax) = (ss+..)2/(3log2) = n + dp_h + z_h + z_l */\r
267	t = (double)n;\r
268	t1 = (((z_h+z_l)+dp_h[k])+t);\r
269	SET_LOW_WORD(t1,0);\r
270	t2 = z_l-(((t1-t)-dp_h[k])-z_h);\r
271	}\r
272	\r
273	/* split up y into y1+y2 and compute (y1+y2)(t1+t2) /\r
274	y1 = y;\r
275	SET_LOW_WORD(y1,0);\r
276	p_l = (y-y1)t1+yt2;\r
277	p_h = y1*t1;\r
278	z = p_l+p_h;\r
279	EXTRACT_WORDS(j,i,z);\r
280	if (j>=0x40900000) { /* z >= 1024 */\r
281	if(((j-0x40900000)\|i)!=0) /* if z > 1024 */\r
282	return shugehuge; /* overflow */\r
283	else {\r
284	if(p_l+ovt>z-p_h) return shugehuge; /* overflow */\r
285	}\r
286	} else if((j&0x7fffffff)>=0x4090cc00 ) { /* z <= -1075 */\r
287	if(((j-0xc090cc00)\|i)!=0) /* z < -1075 */\r
288	return stinytiny; /* underflow */\r
289	else {\r
290	if(p_l<=z-p_h) return stinytiny; /* underflow */\r
291	}\r
292	}\r
293	/*\r
294	* compute 2**(p_h+p_l)\r
295	*/\r
296	i = j&0x7fffffff;\r
297	k = (i>>20)-0x3ff;\r
298	n = 0;\r
299	if(i>0x3fe00000) { /* if \|z\| > 0.5, set n = [z+0.5] */\r
300	n = j+(0x00100000>>(k+1));\r
301	k = ((n&0x7fffffff)>>20)-0x3ff; /* new k for n */\r
302	t = zero;\r
303	SET_HIGH_WORD(t,n&~(0x000fffff>>k));\r
304	n = ((n&0x000fffff)\|0x00100000)>>(20-k);\r
305	if(j<0) n = -n;\r
306	p_h -= t;\r
307	}\r
308	t = p_l+p_h;\r
309	SET_LOW_WORD(t,0);\r
310	u = t*lg2_h;\r
311	v = (p_l-(t-p_h))lg2+tlg2_l;\r
312	z = u+v;\r
313	w = v-(z-u);\r
314	t = z*z;\r
315	t1 = z - t(P1+t(P2+t(P3+t(P4+t*P5))));\r
316	r = (zt1)/(t1-two)-(w+zw);\r
317	z = one-(r-z);\r
318	GET_HIGH_WORD(j,z);\r
319	j += (n<<20);\r
320	if((j>>20)<=0) z = scalbn(z,n); /* subnormal output */\r
321	else SET_HIGH_WORD(z,j);\r
322	return s*z;\r
323	}\r