Print this page
11210 libm should be cstyle(1ONBLD) clean


   5  * Common Development and Distribution License (the "License").
   6  * You may not use this file except in compliance with the License.
   7  *
   8  * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
   9  * or http://www.opensolaris.org/os/licensing.
  10  * See the License for the specific language governing permissions
  11  * and limitations under the License.
  12  *
  13  * When distributing Covered Code, include this CDDL HEADER in each
  14  * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
  15  * If applicable, add the following below this CDDL HEADER, with the
  16  * fields enclosed by brackets "[]" replaced with your own identifying
  17  * information: Portions Copyright [yyyy] [name of copyright owner]
  18  *
  19  * CDDL HEADER END
  20  */
  21 
  22 /*
  23  * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
  24  */

  25 /*
  26  * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
  27  * Use is subject to license terms.
  28  */
  29 
  30 #pragma weak __pow = pow
  31 
  32 /*
  33  * pow(x,y) return x**y
  34  *                    n
  35  * Method:  Let x =  2   * (1+f)
  36  *      1. Compute and return log2(x) in two pieces:
  37  *              log2(x) = w1 + w2,
  38  *         where w1 has 24 bits trailing zero.
  39  *      2. Perform y*log2(x) by simulating muti-precision arithmetic
  40  *      3. Return x**y = exp2(y*log(x))
  41  *
  42  * Special cases:
  43  *      1.  (anything) ** +-0 is 1
  44  *      1'. 1 ** (anything)   is 1      (C99; 1 ** +-INF/NAN used to be NAN)


  54  *      11. -0 ** (+anything except 0, NAN, odd integer)  is +0
  55  *      12. +0 ** (-anything except 0, NAN)               is +INF
  56  *      13. -0 ** (-anything except 0, NAN, odd integer)  is +INF
  57  *      14. -0 ** (odd integer) = -( +0 ** (odd integer) )
  58  *      15. +INF ** (+anything except 0,NAN) is +INF
  59  *      16. +INF ** (-anything except 0,NAN) is +0
  60  *      17. -INF ** (anything)  = -0 ** (-anything)
  61  *      18. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer)
  62  *      19. (-anything except 0 and inf) ** (non-integer) is NAN
  63  *
  64  * Accuracy:
  65  *      pow(x,y) returns x**y nearly rounded. In particular
  66  *                      pow(integer,integer)
  67  *      always returns the correct integer provided it is representable.
  68  */
  69 
  70 #include "libm.h"
  71 #include "xpg6.h"       /* __xpg6 */
  72 #define _C99SUSv3_pow   _C99SUSv3_pow_treats_Inf_as_an_even_int
  73 
  74 static const double zero = 0.0, one = 1.0, two = 2.0;


  75 
  76 extern const double _TBL_log2_hi[], _TBL_log2_lo[];
  77 static const double
  78         two53 = 9007199254740992.0,
  79         A1_hi = 2.8853900432586669921875,
  80         A1_lo = 3.8519259825035041963606002e-8,
  81         A1 = 2.885390081777926817222541963606002026086e+0000,
  82         A2 = 9.617966939207270828380543979852286255862e-0001,
  83         A3 = 5.770807680887875964868853124873696201995e-0001,
  84         B0_hi = 2.8853900432586669921875,
  85         B0_lo = 3.8519259822532793056374320585e-8,
  86         B0 = 2.885390081777926814720293056374320585689e+0000,
  87         B1 = 9.617966939259755138949202350396200257632e-0001,
  88         B2 = 5.770780163585687000782112776448797953382e-0001,
  89         B3 = 4.121985488948771523290174512461778354953e-0001,
  90         B4 = 3.207590534812432970433641789022666850193e-0001;
  91 
  92 static double
  93 log2_x(double x, double *w) {

  94         double f, s, z, qn, h, t;
  95         int *px = (int *) &x;
  96         int *pz = (int *) &z;
  97         int i, j, ix, n;
  98 
  99         n = 0;
 100         ix = px[HIWORD];

 101         if (ix >= 0x3fef03f1 && ix < 0x3ff08208) {        /* 65/63 > x > 63/65 */
 102                 double f1, v;

 103                 f = x - one;

 104                 if (((ix - 0x3ff00000) | px[LOWORD]) == 0) {
 105                         *w = zero;
 106                         return (zero);          /* log2(1)= +0 */
 107                 }

 108                 qn = one / (two + f);
 109                 s = f * qn;                             /* |s|<2**-6 */
 110                 v = s * s;
 111                 h = (double) ((float) s);
 112                 f1 = (double) ((float) f);
 113                 t = qn * (((f - two * h) - h * f1) - h * (f - f1));
 114                                                                 /* s = h+t */
 115                 f1 = h * B0_lo + s * (v * (B1 + v * (B2 + v * (B3 + v * B4))));
 116                 t = f1 + t * B0;
 117                 h *= B0_hi;
 118                 s = (double) ((float) (h + t));
 119                 *w = t - (s - h);
 120                 return (s);
 121         }

 122         if (ix < 0x00100000) {                               /* subnormal x */
 123                 x *= two53;
 124                 n = -53;
 125                 ix = px[HIWORD];
 126         }

 127         /* LARGE N */
 128         n += ((ix + 0x1000) >> 20) - 0x3ff;
 129         ix = (ix & 0x000fffff) | 0x3ff00000;                /* scale x to [1,2] */
 130         px[HIWORD] = ix;
 131         i = ix + 0x1000;
 132         pz[HIWORD] = i & 0xffffe000;
 133         pz[LOWORD] = 0;
 134         qn = one / (x + z);
 135         f = x - z;
 136         s = f * qn;
 137         h = (double) ((float) s);
 138         t = qn * ((f - (h + h) * z) - h * f);
 139         j = (i >> 13) & 0x7f;
 140         f = s * s;
 141         t = t * A1 + h * A1_lo;
 142         t += (s * f) * (A2 + f * A3);
 143         qn = h * A1_hi;
 144         s = n + _TBL_log2_hi[j];
 145         h = qn + s;
 146         t += _TBL_log2_lo[j] - ((h - s) - qn);
 147         f = (double) ((float) (h + t));
 148         *w = t - (f - h);
 149         return (f);
 150 }
 151 
 152 extern const double _TBL_exp2_hi[], _TBL_exp2_lo[];
 153 static const double             /* poly app of 2^x-1 on [-1e-10,2^-7+1e-10] */
 154         E1 = 6.931471805599453100674958533810346197328e-0001,
 155         E2 = 2.402265069587779347846769151717493815979e-0001,
 156         E3 = 5.550410866475410512631124892773937864699e-0002,
 157         E4 = 9.618143209991026824853712740162451423355e-0003,
 158         E5 = 1.333357676549940345096774122231849082991e-0003;
 159 
 160 double
 161 pow(double x, double y) {

 162         double z, ax;
 163         double y1, y2, w1, w2;
 164         int sbx, sby, j, k, yisint;
 165         int hx, hy, ahx, ahy;
 166         unsigned lx, ly;
 167         int *pz = (int *) &z;
 168 
 169         hx = ((int *) &x)[HIWORD];
 170         lx = ((unsigned *) &x)[LOWORD];
 171         hy = ((int *) &y)[HIWORD];
 172         ly = ((unsigned *) &y)[LOWORD];
 173         ahx = hx & ~0x80000000;
 174         ahy = hy & ~0x80000000;

 175         if ((ahy | ly) == 0) {  /* y==zero  */
 176                 if ((ahx | lx) == 0)
 177                         z = _SVID_libm_err(x, y, 20);   /* +-0**+-0 */
 178                 else if ((ahx | (((lx | -lx) >> 31) & 1)) > 0x7ff00000)
 179                         z = _SVID_libm_err(x, y, 42);   /* NaN**+-0 */
 180                 else
 181                         z = one;                        /* x**+-0 = 1 */

 182                 return (z);
 183         } else if (hx == 0x3ff00000 && lx == 0 &&
 184                 (__xpg6 & _C99SUSv3_pow) != 0)
 185                 return (one);                   /* C99: 1**anything = 1 */
 186         else if (ahx > 0x7ff00000 || (ahx == 0x7ff00000 && lx != 0) ||
 187                 ahy > 0x7ff00000 || (ahy == 0x7ff00000 && ly != 0))
 188                 return (x * y); /* +-NaN return x*y; + -> * for Cheetah */


 189                                 /* includes Sun: 1**NaN = NaN */
 190         sbx = (unsigned) hx >> 31;
 191         sby = (unsigned) hy >> 31;
 192         ax = fabs(x);
 193 
 194         /*
 195          * determine if y is an odd int when x < 0
 196          * yisint = 0 ... y is not an integer
 197          * yisint = 1 ... y is an odd int
 198          * yisint = 2 ... y is an even int
 199          */
 200         yisint = 0;

 201         if (sbx) {
 202                 if (ahy >= 0x43400000)
 203                         yisint = 2;             /* even integer y */
 204                 else if (ahy >= 0x3ff00000) {
 205                         k = (ahy >> 20) - 0x3ff;  /* exponent */

 206                         if (k > 20) {
 207                                 j = ly >> (52 - k);

 208                                 if ((j << (52 - k)) == ly)
 209                                         yisint = 2 - (j & 1);
 210                         } else if (ly == 0) {
 211                                 j = ahy >> (20 - k);

 212                                 if ((j << (20 - k)) == ahy)
 213                                         yisint = 2 - (j & 1);
 214                         }
 215                 }
 216         }

 217         /* special value of y */
 218         if (ly == 0) {
 219                 if (ahy == 0x7ff00000) {        /* y is +-inf */
 220                         if (((ahx - 0x3ff00000) | lx) == 0) {
 221                                 if ((__xpg6 & _C99SUSv3_pow) != 0)
 222                                         return (one);
 223                                                 /* C99: (-1)**+-inf = 1 */
 224                                 else
 225                                         return (y - y);

 226                                                 /* Sun: (+-1)**+-inf = NaN */
 227                         } else if (ahx >= 0x3ff00000)
 228                                                 /* (|x|>1)**+,-inf = inf,0 */
 229                                 return (sby == 0 ? y : zero);
 230                         else                    /* (|x|<1)**-,+inf = inf,0 */
 231                                 return (sby != 0 ? -y : zero);
 232                 }


 233                 if (ahy == 0x3ff00000) {        /* y is  +-1 */
 234                         if (sby != 0) { /* y is -1 */
 235                                 if (x == zero)  /* divided by zero */
 236                                         return (_SVID_libm_err(x, y, 23));
 237                                 else if (ahx < 0x40000 || ((ahx - 0x40000) |
 238                                         lx) == 0)       /* overflow */
 239                                         return (_SVID_libm_err(x, y, 21));
 240                                 else
 241                                         return (one / x);
 242                         } else
 243                                 return (x);
 244                 }


 245                 if (hy == 0x40000000) {         /* y is  2 */
 246                         if (ahx >= 0x5ff00000 && ahx < 0x7ff00000)
 247                                 return (_SVID_libm_err(x, y, 21));
 248                                                         /* x*x overflow */
 249                         else if ((ahx < 0x1e56a09e && (ahx | lx) != 0) ||
 250                                 (ahx == 0x1e56a09e && lx < 0x667f3bcd))
 251                                 return (_SVID_libm_err(x, y, 22));
 252                                                         /* x*x underflow */
 253                         else
 254                                 return (x * x);
 255                 }

 256                 if (hy == 0x3fe00000) {
 257                         if (!((ahx | lx) == 0 || ((ahx - 0x7ff00000) | lx) ==
 258                                 0 || sbx == 1))
 259                                 return (sqrt(x));       /* y is 0.5 and x > 0 */
 260                 }
 261         }

 262         /* special value of x */
 263         if (lx == 0) {
 264                 if (ahx == 0x7ff00000 || ahx == 0 || ahx == 0x3ff00000) {
 265                         /* x is +-0,+-inf,-1 */
 266                         z = ax;

 267                         if (sby == 1) {
 268                                 z = one / z;    /* z = |x|**y */

 269                                 if (ahx == 0)
 270                                         return (_SVID_libm_err(x, y, 23));
 271                         }

 272                         if (sbx == 1) {
 273                                 if (ahx == 0x3ff00000 && yisint == 0)
 274                                         z = _SVID_libm_err(x, y, 24);
 275                                         /* neg**non-integral is NaN + invalid */
 276                                 else if (yisint == 1)
 277                                         z = -z; /* (x<0)**odd = -(|x|**odd) */
 278                         }

 279                         return (z);
 280                 }
 281         }

 282         /* (x<0)**(non-int) is NaN */
 283         if (sbx == 1 && yisint == 0)
 284                 return (_SVID_libm_err(x, y, 24));
 285         /* Now ax is finite, y is finite */
 286         /* first compute log2(ax) = w1+w2, with 24 bits w1 */



 287         w1 = log2_x(ax, &w2);
 288 
 289         /* split up y into y1+y2 and compute (y1+y2)*(w1+w2) */
 290         if (((ly & 0x07ffffff) == 0) || ahy >= 0x47e00000 ||
 291                 ahy <= 0x38100000) {
 292                 /* no need to split if y is short or too large or too small */
 293                 y1 = y * w1;
 294                 y2 = y * w2;
 295         } else {
 296                 y1 = (double) ((float) y);
 297                 y2 = (y - y1) * w1 + y * w2;
 298                 y1 *= w1;
 299         }

 300         z = y1 + y2;
 301         j = pz[HIWORD];

 302         if (j >= 0x40900000) {                               /* z >= 1024 */
 303                 if (!(j == 0x40900000 && pz[LOWORD] == 0))      /* z > 1024 */
 304                         return (_SVID_libm_err(x, y, 21));      /* overflow */
 305                 else {
 306                         w2 = y1 - z;
 307                         w2 += y2;

 308                                                         /* rounded to inf */
 309                         if (w2 >= -8.008566259537296567160e-17)
 310                                 return (_SVID_libm_err(x, y, 21));

 311                                                                 /* overflow */
 312                 }
 313         } else if ((j & ~0x80000000) >= 0x4090cc00) {    /* z <= -1075 */
 314                 if (!(j == 0xc090cc00 && pz[LOWORD] == 0))      /* z < -1075 */
 315                         return (_SVID_libm_err(x, y, 22));      /* underflow */
 316                 else {
 317                         w2 = y1 - z;
 318                         w2 += y2;

 319                         if (w2 <= zero)                      /* underflow */
 320                                 return (_SVID_libm_err(x, y, 22));
 321                 }
 322         }

 323         /*
 324          * compute 2**(k+f[j]+g)
 325          */
 326         k = (int) (z * 64.0 + (((hy ^ (ahx - 0x3ff00000)) > 0) ? 0.5 : -0.5));
 327         j = k & 63;
 328         w1 = y2 - ((double) k * 0.015625 - y1);
 329         w2 = _TBL_exp2_hi[j];
 330         z = _TBL_exp2_lo[j] + (w2 * w1) * (E1 + w1 * (E2 + w1 * (E3 + w1 *
 331                 (E4 + w1 * E5))));
 332         z += w2;
 333         k >>= 6;

 334         if (k < -1021)
 335                 z = scalbn(z, k);
 336         else                    /* subnormal output */
 337                 pz[HIWORD] += k << 20;

 338         if (sbx == 1 && yisint == 1)
 339                 z = -z;         /* (-ve)**(odd int) */

 340         return (z);
 341 }


   5  * Common Development and Distribution License (the "License").
   6  * You may not use this file except in compliance with the License.
   7  *
   8  * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
   9  * or http://www.opensolaris.org/os/licensing.
  10  * See the License for the specific language governing permissions
  11  * and limitations under the License.
  12  *
  13  * When distributing Covered Code, include this CDDL HEADER in each
  14  * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
  15  * If applicable, add the following below this CDDL HEADER, with the
  16  * fields enclosed by brackets "[]" replaced with your own identifying
  17  * information: Portions Copyright [yyyy] [name of copyright owner]
  18  *
  19  * CDDL HEADER END
  20  */
  21 
  22 /*
  23  * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
  24  */
  25 
  26 /*
  27  * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
  28  * Use is subject to license terms.
  29  */
  30 
  31 #pragma weak __pow = pow
  32 
  33 /*
  34  * pow(x,y) return x**y
  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 24 bits trailing zero.
  40  *      2. Perform y*log2(x) by simulating muti-precision arithmetic
  41  *      3. Return x**y = exp2(y*log(x))
  42  *
  43  * Special cases:
  44  *      1.  (anything) ** +-0 is 1
  45  *      1'. 1 ** (anything)   is 1      (C99; 1 ** +-INF/NAN used to be NAN)


  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  * Accuracy:
  66  *      pow(x,y) returns x**y nearly rounded. In particular
  67  *                      pow(integer,integer)
  68  *      always returns the correct integer provided it is representable.
  69  */
  70 
  71 #include "libm.h"
  72 #include "xpg6.h"                       /* __xpg6 */
  73 #define _C99SUSv3_pow   _C99SUSv3_pow_treats_Inf_as_an_even_int
  74 
  75 static const double zero = 0.0,
  76         one = 1.0,
  77         two = 2.0;
  78 
  79 extern const double _TBL_log2_hi[], _TBL_log2_lo[];
  80 
  81 static const double two53 = 9007199254740992.0,
  82         A1_hi = 2.8853900432586669921875,
  83         A1_lo = 3.8519259825035041963606002e-8,
  84         A1 = 2.885390081777926817222541963606002026086e+0000,
  85         A2 = 9.617966939207270828380543979852286255862e-0001,
  86         A3 = 5.770807680887875964868853124873696201995e-0001,
  87         B0_hi = 2.8853900432586669921875,
  88         B0_lo = 3.8519259822532793056374320585e-8,
  89         B0 = 2.885390081777926814720293056374320585689e+0000,
  90         B1 = 9.617966939259755138949202350396200257632e-0001,
  91         B2 = 5.770780163585687000782112776448797953382e-0001,
  92         B3 = 4.121985488948771523290174512461778354953e-0001,
  93         B4 = 3.207590534812432970433641789022666850193e-0001;
  94 
  95 static double
  96 log2_x(double x, double *w)
  97 {
  98         double f, s, z, qn, h, t;
  99         int *px = (int *)&x;
 100         int *pz = (int *)&z;
 101         int i, j, ix, n;
 102 
 103         n = 0;
 104         ix = px[HIWORD];
 105 
 106         if (ix >= 0x3fef03f1 && ix < 0x3ff08208) {        /* 65/63 > x > 63/65 */
 107                 double f1, v;
 108 
 109                 f = x - one;
 110 
 111                 if (((ix - 0x3ff00000) | px[LOWORD]) == 0) {
 112                         *w = zero;
 113                         return (zero);  /* log2(1)= +0 */
 114                 }
 115 
 116                 qn = one / (two + f);
 117                 s = f * qn;             /* |s|<2**-6 */
 118                 v = s * s;
 119                 h = (double)((float)s);
 120                 f1 = (double)((float)f);
 121                 t = qn * (((f - two * h) - h * f1) - h * (f - f1));
 122                 /* s = h+t */
 123                 f1 = h * B0_lo + s * (v * (B1 + v * (B2 + v * (B3 + v * B4))));
 124                 t = f1 + t * B0;
 125                 h *= B0_hi;
 126                 s = (double)((float)(h + t));
 127                 *w = t - (s - h);
 128                 return (s);
 129         }
 130 
 131         if (ix < 0x00100000) {               /* subnormal x */
 132                 x *= two53;
 133                 n = -53;
 134                 ix = px[HIWORD];
 135         }
 136 
 137         /* LARGE N */
 138         n += ((ix + 0x1000) >> 20) - 0x3ff;
 139         ix = (ix & 0x000fffff) | 0x3ff00000;        /* scale x to [1,2] */
 140         px[HIWORD] = ix;
 141         i = ix + 0x1000;
 142         pz[HIWORD] = i & 0xffffe000;
 143         pz[LOWORD] = 0;
 144         qn = one / (x + z);
 145         f = x - z;
 146         s = f * qn;
 147         h = (double)((float)s);
 148         t = qn * ((f - (h + h) * z) - h * f);
 149         j = (i >> 13) & 0x7f;
 150         f = s * s;
 151         t = t * A1 + h * A1_lo;
 152         t += (s * f) * (A2 + f * A3);
 153         qn = h * A1_hi;
 154         s = n + _TBL_log2_hi[j];
 155         h = qn + s;
 156         t += _TBL_log2_lo[j] - ((h - s) - qn);
 157         f = (double)((float)(h + t));
 158         *w = t - (f - h);
 159         return (f);
 160 }
 161 
 162 extern const double _TBL_exp2_hi[], _TBL_exp2_lo[];
 163 static const double             /* poly app of 2^x-1 on [-1e-10,2^-7+1e-10] */
 164         E1 = 6.931471805599453100674958533810346197328e-0001,
 165         E2 = 2.402265069587779347846769151717493815979e-0001,
 166         E3 = 5.550410866475410512631124892773937864699e-0002,
 167         E4 = 9.618143209991026824853712740162451423355e-0003,
 168         E5 = 1.333357676549940345096774122231849082991e-0003;
 169 
 170 double
 171 pow(double x, double y)
 172 {
 173         double z, ax;
 174         double y1, y2, w1, w2;
 175         int sbx, sby, j, k, yisint;
 176         int hx, hy, ahx, ahy;
 177         unsigned lx, ly;
 178         int *pz = (int *)&z;
 179 
 180         hx = ((int *)&x)[HIWORD];
 181         lx = ((unsigned *)&x)[LOWORD];
 182         hy = ((int *)&y)[HIWORD];
 183         ly = ((unsigned *)&y)[LOWORD];
 184         ahx = hx & ~0x80000000;
 185         ahy = hy & ~0x80000000;
 186 
 187         if ((ahy | ly) == 0) {                          /* y==zero  */
 188                 if ((ahx | lx) == 0)
 189                         z = _SVID_libm_err(x, y, 20);   /* +-0**+-0 */
 190                 else if ((ahx | (((lx | -lx) >> 31) & 1)) > 0x7ff00000)
 191                         z = _SVID_libm_err(x, y, 42);   /* NaN**+-0 */
 192                 else
 193                         z = one;                        /* x**+-0 = 1 */
 194 
 195                 return (z);
 196         } else if (hx == 0x3ff00000 && lx == 0 && (__xpg6 & _C99SUSv3_pow) !=
 197             0) {
 198                 return (one);           /* C99: 1**anything = 1 */
 199         } else if (ahx > 0x7ff00000 || (ahx == 0x7ff00000 && lx != 0) || ahy >
 200             0x7ff00000 || (ahy == 0x7ff00000 && ly != 0)) {
 201                 return (x * y); /* +-NaN return x*y; + -> * for Cheetah */
 202         }
 203 
 204         /* includes Sun: 1**NaN = NaN */
 205         sbx = (unsigned)hx >> 31;
 206         sby = (unsigned)hy >> 31;
 207         ax = fabs(x);
 208 
 209         /*
 210          * determine if y is an odd int when x < 0
 211          * yisint = 0 ... y is not an integer
 212          * yisint = 1 ... y is an odd int
 213          * yisint = 2 ... y is an even int
 214          */
 215         yisint = 0;
 216 
 217         if (sbx) {
 218                 if (ahy >= 0x43400000) {
 219                         yisint = 2;                     /* even integer y */
 220                 } else if (ahy >= 0x3ff00000) {
 221                         k = (ahy >> 20) - 0x3ff;  /* exponent */
 222 
 223                         if (k > 20) {
 224                                 j = ly >> (52 - k);
 225 
 226                                 if ((j << (52 - k)) == ly)
 227                                         yisint = 2 - (j & 1);
 228                         } else if (ly == 0) {
 229                                 j = ahy >> (20 - k);
 230 
 231                                 if ((j << (20 - k)) == ahy)
 232                                         yisint = 2 - (j & 1);
 233                         }
 234                 }
 235         }
 236 
 237         /* special value of y */
 238         if (ly == 0) {
 239                 if (ahy == 0x7ff00000) {        /* y is +-inf */
 240                         if (((ahx - 0x3ff00000) | lx) == 0) {
 241                                 if ((__xpg6 & _C99SUSv3_pow) != 0)
 242                                         return (one);
 243                                 /* C99: (-1)**+-inf = 1 */
 244                                 else
 245                                         return (y - y);
 246 
 247                                 /* Sun: (+-1)**+-inf = NaN */
 248                         } else if (ahx >= 0x3ff00000) {
 249                                 /* (|x|>1)**+,-inf = inf,0 */
 250                                 return (sby == 0 ? y : zero);
 251                         } else {        /* (|x|<1)**-,+inf = inf,0 */
 252                                 return (sby != 0 ? -y : zero);
 253                         }
 254                 }
 255 
 256                 if (ahy == 0x3ff00000) {        /* y is  +-1 */
 257                         if (sby != 0) {         /* y is -1 */
 258                                 if (x == zero)  /* divided by zero */
 259                                         return (_SVID_libm_err(x, y, 23));
 260                                 else if (ahx < 0x40000 || ((ahx - 0x40000) |
 261                                     lx) == 0)   /* overflow */
 262                                         return (_SVID_libm_err(x, y, 21));
 263                                 else
 264                                         return (one / x);
 265                         } else {
 266                                 return (x);
 267                         }
 268                 }
 269 
 270                 if (hy == 0x40000000) { /* y is  2 */
 271                         if (ahx >= 0x5ff00000 && ahx < 0x7ff00000)
 272                                 return (_SVID_libm_err(x, y, 21));
 273                         /* x*x overflow */
 274                         else if ((ahx < 0x1e56a09e && (ahx | lx) != 0) ||
 275                             (ahx == 0x1e56a09e && lx < 0x667f3bcd))
 276                                 return (_SVID_libm_err(x, y, 22));
 277                         /* x*x underflow */
 278                         else
 279                                 return (x * x);
 280                 }
 281 
 282                 if (hy == 0x3fe00000) {
 283                         if (!((ahx | lx) == 0 || ((ahx - 0x7ff00000) | lx) ==
 284                             0 || sbx == 1))
 285                                 return (sqrt(x));       /* y is 0.5 and x > 0 */
 286                 }
 287         }
 288 
 289         /* special value of x */
 290         if (lx == 0) {
 291                 if (ahx == 0x7ff00000 || ahx == 0 || ahx == 0x3ff00000) {
 292                         /* x is +-0,+-inf,-1 */
 293                         z = ax;
 294 
 295                         if (sby == 1) {
 296                                 z = one / z;    /* z = |x|**y */
 297 
 298                                 if (ahx == 0)
 299                                         return (_SVID_libm_err(x, y, 23));
 300                         }
 301 
 302                         if (sbx == 1) {
 303                                 if (ahx == 0x3ff00000 && yisint == 0)
 304                                         z = _SVID_libm_err(x, y, 24);
 305                                 /* neg**non-integral is NaN + invalid */
 306                                 else if (yisint == 1)
 307                                         z = -z; /* (x<0)**odd = -(|x|**odd) */
 308                         }
 309 
 310                         return (z);
 311                 }
 312         }
 313 
 314         /* (x<0)**(non-int) is NaN */
 315         if (sbx == 1 && yisint == 0)
 316                 return (_SVID_libm_err(x, y, 24));
 317 
 318         /*
 319          * Now ax is finite, y is finite
 320          * first compute log2(ax) = w1+w2, with 24 bits w1
 321          */
 322         w1 = log2_x(ax, &w2);
 323 
 324         /* split up y into y1+y2 and compute (y1+y2)*(w1+w2) */
 325         if (((ly & 0x07ffffff) == 0) || ahy >= 0x47e00000 || ahy <=
 326             0x38100000) {
 327                 /* no need to split if y is short or too large or too small */
 328                 y1 = y * w1;
 329                 y2 = y * w2;
 330         } else {
 331                 y1 = (double)((float)y);
 332                 y2 = (y - y1) * w1 + y * w2;
 333                 y1 *= w1;
 334         }
 335 
 336         z = y1 + y2;
 337         j = pz[HIWORD];
 338 
 339         if (j >= 0x40900000) {                                       /* z >= 1024 */
 340                 if (!(j == 0x40900000 && pz[LOWORD] == 0)) {    /* z > 1024 */
 341                         return (_SVID_libm_err(x, y, 21));      /* overflow */
 342                 } else {
 343                         w2 = y1 - z;
 344                         w2 += y2;
 345 
 346                         /* rounded to inf */
 347                         if (w2 >= -8.008566259537296567160e-17)
 348                                 return (_SVID_libm_err(x, y, 21));
 349 
 350                         /* overflow */
 351                 }
 352         } else if ((j & ~0x80000000) >= 0x4090cc00) {            /* z <= -1075 */
 353                 if (!(j == 0xc090cc00 && pz[LOWORD] == 0)) {    /* z < -1075 */
 354                         return (_SVID_libm_err(x, y, 22));      /* underflow */
 355                 } else {
 356                         w2 = y1 - z;
 357                         w2 += y2;
 358 
 359                         if (w2 <= zero)      /* underflow */
 360                                 return (_SVID_libm_err(x, y, 22));
 361                 }
 362         }
 363 
 364         /*
 365          * compute 2**(k+f[j]+g)
 366          */
 367         k = (int)(z * 64.0 + (((hy ^ (ahx - 0x3ff00000)) > 0) ? 0.5 : -0.5));
 368         j = k & 63;
 369         w1 = y2 - ((double)k * 0.015625 - y1);
 370         w2 = _TBL_exp2_hi[j];
 371         z = _TBL_exp2_lo[j] + (w2 * w1) * (E1 + w1 * (E2 + w1 * (E3 + w1 * (E4 +
 372             w1 * E5))));
 373         z += w2;
 374         k >>= 6;
 375 
 376         if (k < -1021)
 377                 z = scalbn(z, k);
 378         else                            /* subnormal output */
 379                 pz[HIWORD] += k << 20;
 380 
 381         if (sbx == 1 && yisint == 1)
 382                 z = -z;                 /* (-ve)**(odd int) */
 383 
 384         return (z);
 385 }