Print this page
11210 libm should be cstyle(1ONBLD) clean
   1 /*
   2  * CDDL HEADER START
   3  *
   4  * The contents of this file are subject to the terms of the
   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  * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
  23  */

  24 /*
  25  * Copyright 2005 Sun Microsystems, Inc.  All rights reserved.
  26  * Use is subject to license terms.
  27  */
  28 
  29 #pragma weak __powf = powf
  30 
  31 #include "libm.h"
  32 #include "xpg6.h"       /* __xpg6 */
  33 #define _C99SUSv3_pow   _C99SUSv3_pow_treats_Inf_as_an_even_int
  34 
  35 #if defined(__i386) && !defined(__amd64)
  36 extern int __swapRP(int);
  37 #endif
  38 
  39 /* INDENT OFF */
  40 static const double
  41         ln2 = 6.93147180559945286227e-01,       /* 0x3fe62e42, 0xfefa39ef */
  42         invln2 = 1.44269504088896338700e+00,    /* 0x3ff71547, 0x652b82fe */
  43         dtwo = 2.0,
  44         done = 1.0,
  45         dhalf = 0.5,
  46         d32 = 32.0,
  47         d1_32 = 0.03125,
  48         A0 = 1.999999999813723303647511146995966439250e+0000,
  49         A1 = 6.666910817935858533770138657139665608610e-0001,
  50         t0 = 2.000000000004777489262405315073203746943e+0000,
  51         t1 = 1.666663408349926379873111932994250726307e-0001;
  52 
  53 static const double S[] = {
  54         1.00000000000000000000e+00,     /* 3FF0000000000000 */
  55         1.02189714865411662714e+00,     /* 3FF059B0D3158574 */
  56         1.04427378242741375480e+00,     /* 3FF0B5586CF9890F */
  57         1.06714040067682369717e+00,     /* 3FF11301D0125B51 */
  58         1.09050773266525768967e+00,     /* 3FF172B83C7D517B */
  59         1.11438674259589243221e+00,     /* 3FF1D4873168B9AA */


  69         1.38390988196383202258e+00,     /* 3FF6247EB03A5585 */
  70         1.41421356237309514547e+00,     /* 3FF6A09E667F3BCD */
  71         1.44518080697704665027e+00,     /* 3FF71F75E8EC5F74 */
  72         1.47682614593949934623e+00,     /* 3FF7A11473EB0187 */
  73         1.50916442759342284141e+00,     /* 3FF82589994CCE13 */
  74         1.54221082540794074411e+00,     /* 3FF8ACE5422AA0DB */
  75         1.57598084510788649659e+00,     /* 3FF93737B0CDC5E5 */
  76         1.61049033194925428347e+00,     /* 3FF9C49182A3F090 */
  77         1.64575547815396494578e+00,     /* 3FFA5503B23E255D */
  78         1.68179283050742900407e+00,     /* 3FFAE89F995AD3AD */
  79         1.71861929812247793414e+00,     /* 3FFB7F76F2FB5E47 */
  80         1.75625216037329945351e+00,     /* 3FFC199BDD85529C */
  81         1.79470907500310716820e+00,     /* 3FFCB720DCEF9069 */
  82         1.83400808640934243066e+00,     /* 3FFD5818DCFBA487 */
  83         1.87416763411029996256e+00,     /* 3FFDFC97337B9B5F */
  84         1.91520656139714740007e+00,     /* 3FFEA4AFA2A490DA */
  85         1.95714412417540017941e+00,     /* 3FFF50765B6E4540 */
  86 };
  87 
  88 static const double TBL[] = {
  89         0.00000000000000000e+00,
  90         3.07716586667536873e-02,
  91         6.06246218164348399e-02,
  92         8.96121586896871380e-02,
  93         1.17783035656383456e-01,
  94         1.45182009844497889e-01,
  95         1.71850256926659228e-01,
  96         1.97825743329919868e-01,
  97         2.23143551314209765e-01,
  98         2.47836163904581269e-01,
  99         2.71933715483641758e-01,
 100         2.95464212893835898e-01,
 101         3.18453731118534589e-01,
 102         3.40926586970593193e-01,
 103         3.62905493689368475e-01,
 104         3.84411698910332056e-01,
 105         4.05465108108164385e-01,
 106         4.26084395310900088e-01,
 107         4.46287102628419530e-01,
 108         4.66089729924599239e-01,
 109         4.85507815781700824e-01,
 110         5.04556010752395312e-01,
 111         5.23248143764547868e-01,
 112         5.41597282432744409e-01,
 113         5.59615787935422659e-01,
 114         5.77315365034823613e-01,
 115         5.94707107746692776e-01,
 116         6.11801541105992941e-01,
 117         6.28608659422374094e-01,
 118         6.45137961373584701e-01,
 119         6.61398482245365016e-01,
 120         6.77398823591806143e-01,
 121 };
 122 
 123 static const float zero = 0.0F, one = 1.0F, huge = 1.0e25f, tiny = 1.0e-25f;
 124 /* INDENT ON */
 125 
 126 float
 127 powf(float x, float y) {

 128         float   fx = x, fy = y;
 129         float   fz;
 130         int     ix, iy, jx, jy, k, iw, yisint;
 131 
 132         ix = *(int *)&x;
 133         iy = *(int *)&y;
 134         jx = ix & ~0x80000000;
 135         jy = iy & ~0x80000000;
 136 
 137         if (jy == 0)
 138                 return (one);   /* x**+-0 = 1 */
 139         else if (ix == 0x3f800000 && (__xpg6 & _C99SUSv3_pow) != 0)
 140                 return (one);   /* C99: 1**anything = 1 */
 141         else if (((0x7f800000 - jx) | (0x7f800000 - jy)) < 0)
 142                 return (fx * fy);       /* at least one of x or y is NaN */
 143                                         /* includes Sun: 1**NaN = NaN */
 144         /* INDENT OFF */



 145         /*
 146          * determine if y is an odd int
 147          * yisint = 0 ... y is not an integer
 148          * yisint = 1 ... y is an odd int
 149          * yisint = 2 ... y is an even int
 150          */
 151         /* INDENT ON */
 152         yisint = 0;

 153         if (ix < 0) {
 154                 if (jy >= 0x4b800000) {
 155                         yisint = 2;     /* |y|>=2**24: y must be even */
 156                 } else if (jy >= 0x3f800000) {
 157                         k = (jy >> 23) - 0x7f;    /* exponent */
 158                         iw = jy >> (23 - k);

 159                         if ((iw << (23 - k)) == jy)
 160                                 yisint = 2 - (iw & 1);
 161                 }
 162         }
 163 
 164         /* special value of y */
 165         if ((jy & ~0x7f800000) == 0) {
 166                 if (jy == 0x7f800000) {         /* y is +-inf */
 167                         if (jx == 0x3f800000) {
 168                                 if ((__xpg6 & _C99SUSv3_pow) != 0)
 169                                         fz = one;
 170                                                 /* C99: (-1)**+-inf is 1 */
 171                                 else
 172                                         fz = fy - fy;

 173                                                 /* Sun: (+-1)**+-inf = NaN */
 174                         } else if (jx > 0x3f800000) {
 175                                                 /* (|x|>1)**+,-inf = inf,0 */
 176                                 if (iy > 0)
 177                                         fz = fy;
 178                                 else
 179                                         fz = zero;
 180                         } else {                /* (|x|<1)**-,+inf = inf,0 */
 181                                 if (iy < 0)
 182                                         fz = -fy;
 183                                 else
 184                                         fz = zero;
 185                         }

 186                         return (fz);
 187                 } else if (jy == 0x3f800000) {  /* y is +-1 */
 188                         if (iy < 0)
 189                                 fx = one / fx;  /* y is -1 */

 190                         return (fx);
 191                 } else if (iy == 0x40000000) {  /* y is 2 */
 192                         return (fx * fx);
 193                 } else if (iy == 0x3f000000) {  /* y is 0.5 */
 194                         if (jx != 0 && jx != 0x7f800000)
 195                                 return (sqrtf(x));
 196                 }
 197         }
 198 
 199         /* special value of x */
 200         if ((jx & ~0x7f800000) == 0) {
 201                 if (jx == 0x7f800000 || jx == 0 || jx == 0x3f800000) {
 202                         /* x is +-0,+-inf,-1; set fz = |x|**y */
 203                         *(int *)&fz = jx;

 204                         if (iy < 0)
 205                                 fz = one / fz;

 206                         if (ix < 0) {
 207                                 if (jx == 0x3f800000 && yisint == 0) {
 208                                         /* (-1)**non-int is NaN */
 209                                         fz = zero;
 210                                         fz /= fz;
 211                                 } else if (yisint == 1) {
 212                                         /* (x<0)**odd = -(|x|**odd) */
 213                                         fz = -fz;
 214                                 }
 215                         }

 216                         return (fz);
 217                 }
 218         }
 219 
 220         /* (x<0)**(non-int) is NaN */
 221         if (ix < 0 && yisint == 0) {
 222                 fz = zero;
 223                 return (fz / fz);
 224         }
 225 
 226         /*
 227          * compute exp(y*log(|x|))
 228          * fx = *(float *) &jx;
 229          * fz = (float) exp(((double) fy) * log((double) fx));
 230          */
 231         {
 232                 double  dx, dy, dz, ds;
 233                 int     *px = (int *)&dx, *pz = (int *)&dz, i, n, m;

 234 #if defined(__i386) && !defined(__amd64)
 235                 int     rp = __swapRP(fp_extended);
 236 #endif
 237 
 238                 fx = *(float *)&jx;
 239                 dx = (double)fx;
 240 
 241                 /* compute log(x)/ln2 */
 242                 i = px[HIWORD] + 0x4000;
 243                 n = (i >> 20) - 0x3ff;
 244                 pz[HIWORD] = i & 0xffff8000;
 245                 pz[LOWORD] = 0;
 246                 ds = (dx - dz) / (dx + dz);
 247                 i = (i >> 15) & 0x1f;
 248                 dz = ds * ds;
 249                 dy = invln2 * (TBL[i] + ds * (A0 + dz * A1));

 250                 if (n == 0)
 251                         dz = (double)fy * dy;
 252                 else
 253                         dz = (double)fy * (dy + (double)n);
 254 
 255                 /* compute exp2(dz=y*ln(x)) */
 256                 i = pz[HIWORD];

 257                 if ((i & ~0x80000000) >= 0x40640000) {   /* |z| >= 160.0 */
 258                         fz = (i > 0)? huge : tiny;

 259                         if (ix < 0 && yisint == 1)
 260                                 fz *= -fz;      /* (-ve)**(odd int) */
 261                         else
 262                                 fz *= fz;

 263 #if defined(__i386) && !defined(__amd64)
 264                         if (rp != fp_extended)
 265                                 (void) __swapRP(rp);
 266 #endif
 267                         return (fz);
 268                 }
 269 
 270                 n = (int)(d32 * dz + (i > 0 ? dhalf : -dhalf));
 271                 i = n & 0x1f;
 272                 m = n >> 5;
 273                 dy = ln2 * (dz - d1_32 * (double)n);
 274                 dx = S[i] * (done - (dtwo * dy) / (dy * (done - dy * t1) - t0));

 275                 if (m != 0)
 276                         px[HIWORD] += m << 20;

 277                 fz = (float)dx;
 278 #if defined(__i386) && !defined(__amd64)
 279                 if (rp != fp_extended)
 280                         (void) __swapRP(rp);
 281 #endif
 282         }
 283 
 284         /* end of computing exp(y*log(x)) */
 285         if (ix < 0 && yisint == 1)
 286                 fz = -fz;       /* (-ve)**(odd int) */

 287         return (fz);
 288 }
   1 /*
   2  * CDDL HEADER START
   3  *
   4  * The contents of this file are subject to the terms of the
   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 2005 Sun Microsystems, Inc.  All rights reserved.
  28  * Use is subject to license terms.
  29  */
  30 
  31 #pragma weak __powf = powf
  32 
  33 #include "libm.h"
  34 #include "xpg6.h"                       /* __xpg6 */
  35 #define _C99SUSv3_pow   _C99SUSv3_pow_treats_Inf_as_an_even_int
  36 
  37 #if defined(__i386) && !defined(__amd64)
  38 extern int __swapRP(int);
  39 #endif
  40 

  41 static const double
  42         ln2 = 6.93147180559945286227e-01,    /* 0x3fe62e42, 0xfefa39ef */
  43         invln2 = 1.44269504088896338700e+00, /* 0x3ff71547, 0x652b82fe */
  44         dtwo = 2.0,
  45         done = 1.0,
  46         dhalf = 0.5,
  47         d32 = 32.0,
  48         d1_32 = 0.03125,
  49         A0 = 1.999999999813723303647511146995966439250e+0000,
  50         A1 = 6.666910817935858533770138657139665608610e-0001,
  51         t0 = 2.000000000004777489262405315073203746943e+0000,
  52         t1 = 1.666663408349926379873111932994250726307e-0001;
  53 
  54 static const double S[] = {
  55         1.00000000000000000000e+00,     /* 3FF0000000000000 */
  56         1.02189714865411662714e+00,     /* 3FF059B0D3158574 */
  57         1.04427378242741375480e+00,     /* 3FF0B5586CF9890F */
  58         1.06714040067682369717e+00,     /* 3FF11301D0125B51 */
  59         1.09050773266525768967e+00,     /* 3FF172B83C7D517B */
  60         1.11438674259589243221e+00,     /* 3FF1D4873168B9AA */


  70         1.38390988196383202258e+00,     /* 3FF6247EB03A5585 */
  71         1.41421356237309514547e+00,     /* 3FF6A09E667F3BCD */
  72         1.44518080697704665027e+00,     /* 3FF71F75E8EC5F74 */
  73         1.47682614593949934623e+00,     /* 3FF7A11473EB0187 */
  74         1.50916442759342284141e+00,     /* 3FF82589994CCE13 */
  75         1.54221082540794074411e+00,     /* 3FF8ACE5422AA0DB */
  76         1.57598084510788649659e+00,     /* 3FF93737B0CDC5E5 */
  77         1.61049033194925428347e+00,     /* 3FF9C49182A3F090 */
  78         1.64575547815396494578e+00,     /* 3FFA5503B23E255D */
  79         1.68179283050742900407e+00,     /* 3FFAE89F995AD3AD */
  80         1.71861929812247793414e+00,     /* 3FFB7F76F2FB5E47 */
  81         1.75625216037329945351e+00,     /* 3FFC199BDD85529C */
  82         1.79470907500310716820e+00,     /* 3FFCB720DCEF9069 */
  83         1.83400808640934243066e+00,     /* 3FFD5818DCFBA487 */
  84         1.87416763411029996256e+00,     /* 3FFDFC97337B9B5F */
  85         1.91520656139714740007e+00,     /* 3FFEA4AFA2A490DA */
  86         1.95714412417540017941e+00,     /* 3FFF50765B6E4540 */
  87 };
  88 
  89 static const double TBL[] = {
  90         0.00000000000000000e+00, 3.07716586667536873e-02,
  91         6.06246218164348399e-02, 8.96121586896871380e-02,
  92         1.17783035656383456e-01, 1.45182009844497889e-01,
  93         1.71850256926659228e-01, 1.97825743329919868e-01,
  94         2.23143551314209765e-01, 2.47836163904581269e-01,
  95         2.71933715483641758e-01, 2.95464212893835898e-01,
  96         3.18453731118534589e-01, 3.40926586970593193e-01,
  97         3.62905493689368475e-01, 3.84411698910332056e-01,
  98         4.05465108108164385e-01, 4.26084395310900088e-01,
  99         4.46287102628419530e-01, 4.66089729924599239e-01,
 100         4.85507815781700824e-01, 5.04556010752395312e-01,
 101         5.23248143764547868e-01, 5.41597282432744409e-01,
 102         5.59615787935422659e-01, 5.77315365034823613e-01,
 103         5.94707107746692776e-01, 6.11801541105992941e-01,
 104         6.28608659422374094e-01, 6.45137961373584701e-01,
 105         6.61398482245365016e-01, 6.77398823591806143e-01,
















 106 };
 107 
 108 static const float zero = 0.0F, one = 1.0F, huge = 1.0e25f, tiny = 1.0e-25f;
 109 
 110 
 111 float
 112 powf(float x, float y)
 113 {
 114         float fx = x, fy = y;
 115         float fz;
 116         int ix, iy, jx, jy, k, iw, yisint;
 117 
 118         ix = *(int *)&x;
 119         iy = *(int *)&y;
 120         jx = ix & ~0x80000000;
 121         jy = iy & ~0x80000000;
 122 
 123         if (jy == 0)
 124                 return (one);           /* x**+-0 = 1 */
 125         else if (ix == 0x3f800000 && (__xpg6 & _C99SUSv3_pow) != 0)
 126                 return (one);           /* C99: 1**anything = 1 */
 127         else if (((0x7f800000 - jx) | (0x7f800000 - jy)) < 0)
 128                 return (fx * fy);       /* at least one of x or y is NaN */
 129 
 130         /*
 131          * includes Sun: 1**NaN = NaN
 132          */
 133 
 134         /*
 135          * determine if y is an odd int
 136          * yisint = 0 ... y is not an integer
 137          * yisint = 1 ... y is an odd int
 138          * yisint = 2 ... y is an even int
 139          */

 140         yisint = 0;
 141 
 142         if (ix < 0) {
 143                 if (jy >= 0x4b800000) {
 144                         yisint = 2;             /* |y|>=2**24: y must be even */
 145                 } else if (jy >= 0x3f800000) {
 146                         k = (jy >> 23) - 0x7f;    /* exponent */
 147                         iw = jy >> (23 - k);
 148 
 149                         if ((iw << (23 - k)) == jy)
 150                                 yisint = 2 - (iw & 1);
 151                 }
 152         }
 153 
 154         /* special value of y */
 155         if ((jy & ~0x7f800000) == 0) {
 156                 if (jy == 0x7f800000) { /* y is +-inf */
 157                         if (jx == 0x3f800000) {
 158                                 if ((__xpg6 & _C99SUSv3_pow) != 0)
 159                                         fz = one;
 160                                 /* C99: (-1)**+-inf is 1 */
 161                                 else
 162                                         fz = fy - fy;
 163 
 164                                 /* Sun: (+-1)**+-inf = NaN */
 165                         } else if (jx > 0x3f800000) {
 166                                 /* (|x|>1)**+,-inf = inf,0 */
 167                                 if (iy > 0)
 168                                         fz = fy;
 169                                 else
 170                                         fz = zero;
 171                         } else {        /* (|x|<1)**-,+inf = inf,0 */
 172                                 if (iy < 0)
 173                                         fz = -fy;
 174                                 else
 175                                         fz = zero;
 176                         }
 177 
 178                         return (fz);
 179                 } else if (jy == 0x3f800000) {  /* y is +-1 */
 180                         if (iy < 0)
 181                                 fx = one / fx;  /* y is -1 */
 182 
 183                         return (fx);
 184                 } else if (iy == 0x40000000) {  /* y is 2 */
 185                         return (fx * fx);
 186                 } else if (iy == 0x3f000000) {  /* y is 0.5 */
 187                         if (jx != 0 && jx != 0x7f800000)
 188                                 return (sqrtf(x));
 189                 }
 190         }
 191 
 192         /* special value of x */
 193         if ((jx & ~0x7f800000) == 0) {
 194                 if (jx == 0x7f800000 || jx == 0 || jx == 0x3f800000) {
 195                         /* x is +-0,+-inf,-1; set fz = |x|**y */
 196                         *(int *)&fz = jx;
 197 
 198                         if (iy < 0)
 199                                 fz = one / fz;
 200 
 201                         if (ix < 0) {
 202                                 if (jx == 0x3f800000 && yisint == 0) {
 203                                         /* (-1)**non-int is NaN */
 204                                         fz = zero;
 205                                         fz /= fz;
 206                                 } else if (yisint == 1) {
 207                                         /* (x<0)**odd = -(|x|**odd) */
 208                                         fz = -fz;
 209                                 }
 210                         }
 211 
 212                         return (fz);
 213                 }
 214         }
 215 
 216         /* (x<0)**(non-int) is NaN */
 217         if (ix < 0 && yisint == 0) {
 218                 fz = zero;
 219                 return (fz / fz);
 220         }
 221 
 222         /*
 223          * compute exp(y*log(|x|))
 224          * fx = *(float *) &jx;
 225          * fz = (float) exp(((double) fy) * log((double) fx));
 226          */
 227         {
 228                 double dx, dy, dz, ds;
 229                 int *px = (int *)&dx, *pz = (int *)&dz, i, n, m;
 230 
 231 #if defined(__i386) && !defined(__amd64)
 232                 int rp = __swapRP(fp_extended);
 233 #endif
 234 
 235                 fx = *(float *)&jx;
 236                 dx = (double)fx;
 237 
 238                 /* compute log(x)/ln2 */
 239                 i = px[HIWORD] + 0x4000;
 240                 n = (i >> 20) - 0x3ff;
 241                 pz[HIWORD] = i & 0xffff8000;
 242                 pz[LOWORD] = 0;
 243                 ds = (dx - dz) / (dx + dz);
 244                 i = (i >> 15) & 0x1f;
 245                 dz = ds * ds;
 246                 dy = invln2 * (TBL[i] + ds * (A0 + dz * A1));
 247 
 248                 if (n == 0)
 249                         dz = (double)fy * dy;
 250                 else
 251                         dz = (double)fy * (dy + (double)n);
 252 
 253                 /* compute exp2(dz=y*ln(x)) */
 254                 i = pz[HIWORD];
 255 
 256                 if ((i & ~0x80000000) >= 0x40640000) {   /* |z| >= 160.0 */
 257                         fz = (i > 0) ? huge : tiny;
 258 
 259                         if (ix < 0 && yisint == 1)
 260                                 fz *= -fz;      /* (-ve)**(odd int) */
 261                         else
 262                                 fz *= fz;
 263 
 264 #if defined(__i386) && !defined(__amd64)
 265                         if (rp != fp_extended)
 266                                 (void) __swapRP(rp);
 267 #endif
 268                         return (fz);
 269                 }
 270 
 271                 n = (int)(d32 * dz + (i > 0 ? dhalf : -dhalf));
 272                 i = n & 0x1f;
 273                 m = n >> 5;
 274                 dy = ln2 * (dz - d1_32 * (double)n);
 275                 dx = S[i] * (done - (dtwo * dy) / (dy * (done - dy * t1) - t0));
 276 
 277                 if (m != 0)
 278                         px[HIWORD] += m << 20;
 279 
 280                 fz = (float)dx;
 281 #if defined(__i386) && !defined(__amd64)
 282                 if (rp != fp_extended)
 283                         (void) __swapRP(rp);
 284 #endif
 285         }
 286 
 287         /* end of computing exp(y*log(x)) */
 288         if (ix < 0 && yisint == 1)
 289                 fz = -fz;               /* (-ve)**(odd int) */
 290 
 291         return (fz);
 292 }