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

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