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 */ 60 1.13878863475669156458e+00, /* 3FF2387A6E756238 */ 61 1.16372485877757747552e+00, /* 3FF29E9DF51FDEE1 */ 62 1.18920711500272102690e+00, /* 3FF306FE0A31B715 */ 63 1.21524735998046895524e+00, /* 3FF371A7373AA9CB */ 64 1.24185781207348400201e+00, /* 3FF3DEA64C123422 */ 65 1.26905095719173321989e+00, /* 3FF44E086061892D */ 66 1.29683955465100964055e+00, /* 3FF4BFDAD5362A27 */ 67 1.32523664315974132322e+00, /* 3FF5342B569D4F82 */ 68 1.35425554693689265129e+00, /* 3FF5AB07DD485429 */ 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 } 289