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 * Copyright 2006 Sun Microsystems, Inc. All rights reserved. 27 * Use is subject to license terms. 28 */ 29 30 #pragma weak __powl = powl 31 32 #include "libm.h" 33 #include "xpg6.h" /* __xpg6 */ 34 #define _C99SUSv3_pow _C99SUSv3_pow_treats_Inf_as_an_even_int 35 36 #if defined(__sparc) 37 #define i0 0 38 #define i1 1 39 #define i2 2 40 #define i3 3 41 42 static const long double zero = 0.0L, one = 1.0L, two = 2.0L; 43 44 extern const long double _TBL_logl_hi[], _TBL_logl_lo[]; 45 46 static const long double 47 two113 = 10384593717069655257060992658440192.0L, 48 ln2hi = 6.931471805599453094172319547495844850203e-0001L, 49 ln2lo = 1.667085920830552208890449330400379754169e-0025L, 50 A2 = 6.666666666666666666666666666666091393804e-0001L, 51 A3 = 4.000000000000000000000000407167070220671e-0001L, 52 A4 = 2.857142857142857142730077490612903681164e-0001L, 53 A5 = 2.222222222222242577702836920812882605099e-0001L, 54 A6 = 1.818181816435493395985912667105885828356e-0001L, 55 A7 = 1.538537835211839751112067512805496931725e-0001L, 56 B1 = 6.666666666666666666666666666666666667787e-0001L, 57 B2 = 3.999999999999999999999999999999848524411e-0001L, 58 B3 = 2.857142857142857142857142865084581075070e-0001L, 59 B4 = 2.222222222222222222222010781800643808497e-0001L, 60 B5 = 1.818181818181818185051442171337036403674e-0001L, 61 B6 = 1.538461538461508363540720286292008207673e-0001L, 62 B7 = 1.333333333506731842033180638329317108428e-0001L, 63 B8 = 1.176469984587418890634302788283946761670e-0001L, 64 B9 = 1.053794891561452331722969901564862497132e-0001L; 65 66 static long double 67 logl_x(long double x, long double *w) { 68 long double f, f1, v, s, z, qn, h, t; 69 int *px = (int *) &x; 70 int *pz = (int *) &z; 71 int i, j, ix, n; 72 73 n = 0; 74 ix = px[i0]; 75 if (ix > 0x3ffef03f && ix < 0x3fff0820) { /* 65/63 > x > 63/65 */ 76 f = x - one; 77 z = f * f; 78 if (((ix - 0x3fff0000) | px[i1] | px[i2] | px[i3]) == 0) { 79 *w = zero; 80 return (zero); /* log(1)= +0 */ 81 } 82 qn = one / (two + f); 83 s = f * qn; /* |s|<2**-6 */ 84 v = s * s; 85 h = (long double) (2.0 * (double) s); 86 f1 = (long double) ((double) f); 87 t = ((two * (f - h) - h * f1) - h * (f - f1)) * qn + 88 s * (v * (B1 + v * (B2 + v * (B3 + v * (B4 + 89 v * (B5 + v * (B6 + v * (B7 + v * (B8 + v * B9))))))))); 90 s = (long double) ((double) (h + t)); 91 *w = t - (s - h); 92 return (s); 93 } 94 if (ix < 0x00010000) { /* subnormal x */ 95 x *= two113; 96 n = -113; 97 ix = px[i0]; 98 } 99 /* LARGE_N */ 100 n += ((ix + 0x200) >> 16) - 0x3fff; 101 ix = (ix & 0x0000ffff) | 0x3fff0000; /* scale x to [1,2] */ 102 px[i0] = ix; 103 i = ix + 0x200; 104 pz[i0] = i & 0xfffffc00; 105 pz[i1] = pz[i2] = pz[i3] = 0; 106 qn = one / (x + z); 107 f = x - z; 108 s = f * qn; 109 f1 = (long double) ((double) f); 110 h = (long double) (2.0 * (double) s); 111 t = qn * ((two * (f - z * h) - h * f1) - h * (f - f1)); 112 j = (i >> 10) & 0x3f; 113 v = s * s; 114 qn = (long double) n; 115 t += qn * ln2lo + _TBL_logl_lo[j]; 116 t += s * (v * (A2 + v * (A3 + v * (A4 + v * (A5 + v * (A6 + 117 v * A7)))))); 118 v = qn * ln2hi + _TBL_logl_hi[j]; 119 s = h + v; 120 t += (h - (s - v)); 121 z = (long double) ((double) (s + t)); 122 *w = t - (z - s); 123 return (z); 124 } 125 126 extern const long double _TBL_expl_hi[], _TBL_expl_lo[]; 127 static const long double 128 invln2_32 = 4.616624130844682903551758979206054839765e+1L, 129 ln2_32hi = 2.166084939249829091928849858592451515688e-2L, 130 ln2_32lo = 5.209643502595475652782654157501186731779e-27L, 131 ln2_64 = 1.083042469624914545964425189778400898568e-2L; 132 133 long double 134 powl(long double x, long double y) { 135 long double z, ax; 136 long double y1, y2, w1, w2; 137 int sbx, sby, j, k, yisint, m; 138 int hx, lx, hy, ly, ahx, ahy; 139 int *pz = (int *) &z; 140 int *px = (int *) &x; 141 int *py = (int *) &y; 142 143 hx = px[i0]; 144 lx = px[i1] | px[i2] | px[i3]; 145 hy = py[i0]; 146 ly = py[i1] | py[i2] | py[i3]; 147 ahx = hx & ~0x80000000; 148 ahy = hy & ~0x80000000; 149 150 if ((ahy | ly) == 0) 151 return (one); /* x**+-0 = 1 */ 152 else if (hx == 0x3fff0000 && lx == 0 && 153 (__xpg6 & _C99SUSv3_pow) != 0) 154 return (one); /* C99: 1**anything = 1 */ 155 else if (ahx > 0x7fff0000 || (ahx == 0x7fff0000 && lx != 0) || 156 ahy > 0x7fff0000 || (ahy == 0x7fff0000 && ly != 0)) 157 return (x + y); /* +-NaN return x+y */ 158 /* includes Sun: 1**NaN = NaN */ 159 sbx = (unsigned) hx >> 31; 160 sby = (unsigned) hy >> 31; 161 ax = fabsl(x); 162 /* 163 * determine if y is an odd int when x < 0 164 * yisint = 0 ... y is not an integer 165 * yisint = 1 ... y is an odd int 166 * yisint = 2 ... y is an even int 167 */ 168 yisint = 0; 169 if (sbx) { 170 if (ahy >= 0x40700000) /* if |y|>=2**113 */ 171 yisint = 2; /* even integer y */ 172 else if (ahy >= 0x3fff0000) { 173 k = (ahy >> 16) - 0x3fff; /* exponent */ 174 if (k > 80) { 175 j = ((unsigned) py[i3]) >> (112 - k); 176 if ((j << (112 - k)) == py[i3]) 177 yisint = 2 - (j & 1); 178 } else if (k > 48) { 179 j = ((unsigned) py[i2]) >> (80 - k); 180 if ((j << (80 - k)) == py[i2]) 181 yisint = 2 - (j & 1); 182 } else if (k > 16) { 183 j = ((unsigned) py[i1]) >> (48 - k); 184 if ((j << (48 - k)) == py[i1]) 185 yisint = 2 - (j & 1); 186 } else if (ly == 0) { 187 j = ahy >> (16 - k); 188 if ((j << (16 - k)) == ahy) 189 yisint = 2 - (j & 1); 190 } 191 } 192 } 193 194 /* special value of y */ 195 if (ly == 0) { 196 if (ahy == 0x7fff0000) { /* y is +-inf */ 197 if (((ahx - 0x3fff0000) | lx) == 0) { 198 if ((__xpg6 & _C99SUSv3_pow) != 0) 199 return (one); 200 /* C99: (-1)**+-inf = 1 */ 201 else 202 return (y - y); 203 /* Sun: (+-1)**+-inf = NaN */ 204 } else if (ahx >= 0x3fff0000) 205 /* (|x|>1)**+,-inf = inf,0 */ 206 return (sby == 0 ? y : zero); 207 else /* (|x|<1)**-,+inf = inf,0 */ 208 return (sby != 0 ? -y : zero); 209 } else if (ahy == 0x3fff0000) { /* y is +-1 */ 210 if (sby != 0) 211 return (one / x); 212 else 213 return (x); 214 } else if (hy == 0x40000000) /* y is 2 */ 215 return (x * x); 216 else if (hy == 0x3ffe0000) { /* y is 0.5 */ 217 if (!((ahx | lx) == 0 || ((ahx - 0x7fff0000) | lx) == 218 0)) 219 return (sqrtl(x)); 220 } 221 } 222 223 /* special value of x */ 224 if (lx == 0) { 225 if (ahx == 0x7fff0000 || ahx == 0 || ahx == 0x3fff0000) { 226 /* x is +-0,+-inf,+-1 */ 227 z = ax; 228 if (sby == 1) 229 z = one / z; /* z = 1/|x| if y is negative */ 230 if (sbx == 1) { 231 if (ahx == 0x3fff0000 && yisint == 0) 232 z = zero / zero; 233 /* (-1)**non-int is NaN */ 234 else if (yisint == 1) 235 z = -z; /* (x<0)**odd = -(|x|**odd) */ 236 } 237 return (z); 238 } 239 } 240 241 /* (x<0)**(non-int) is NaN */ 242 if (sbx == 1 && yisint == 0) 243 return (zero / zero); /* should be volatile */ 244 245 /* Now ax is finite, y is finite */ 246 /* first compute log(ax) = w1+w2, with 53 bits w1 */ 247 w1 = logl_x(ax, &w2); 248 249 /* split up y into y1+y2 and compute (y1+y2)*(w1+w2) */ 250 if (ly == 0 || ahy >= 0x43fe0000) { 251 y1 = y * w1; 252 y2 = y * w2; 253 } else { 254 y1 = (long double) ((double) y); 255 y2 = (y - y1) * w1 + y * w2; 256 y1 *= w1; 257 } 258 z = y1 + y2; 259 j = pz[i0]; 260 if ((unsigned) j >= 0xffff0000) { /* NaN or -inf */ 261 if (sbx == 1 && yisint == 1) 262 return (one / z); 263 else 264 return (-one / z); 265 } else if ((j & ~0x80000000) < 0x3fc30000) { /* |x|<2^-60 */ 266 if (sbx == 1 && yisint == 1) 267 return (-one - z); 268 else 269 return (one + z); 270 } else if (j > 0) { 271 if (j > 0x400d0000) { 272 if (sbx == 1 && yisint == 1) 273 return (scalbnl(-one, 20000)); 274 else 275 return (scalbnl(one, 20000)); 276 } 277 k = (int) (invln2_32 * (z + ln2_64)); 278 } else { 279 if ((unsigned) j > 0xc00d0000) { 280 if (sbx == 1 && yisint == 1) 281 return (scalbnl(-one, -20000)); 282 else 283 return (scalbnl(one, -20000)); 284 } 285 k = (int) (invln2_32 * (z - ln2_64)); 286 } 287 j = k & 0x1f; 288 m = k >> 5; 289 { 290 /* rational approximation coeffs for [-(ln2)/64,(ln2)/64] */ 291 long double 292 t1 = 1.666666666666666666666666666660876387437e-1L, 293 t2 = -2.777777777777777777777707812093173478756e-3L, 294 t3 = 6.613756613756613482074280932874221202424e-5L, 295 t4 = -1.653439153392139954169609822742235851120e-6L, 296 t5 = 4.175314851769539751387852116610973796053e-8L; 297 long double t = (long double) k; 298 299 w1 = (y2 - (t * ln2_32hi - y1)) - t * ln2_32lo; 300 t = w1 * w1; 301 w2 = (w1 - t * (t1 + t * (t2 + t * (t3 + t * (t4 + t * t5))))) - 302 two; 303 z = _TBL_expl_hi[j] - ((_TBL_expl_hi[j] * (w1 + w1)) / w2 - 304 _TBL_expl_lo[j]); 305 } 306 j = m + (pz[i0] >> 16); 307 if (j && (unsigned) j < 0x7fff) 308 pz[i0] += m << 16; 309 else 310 z = scalbnl(z, m); 311 312 if (sbx == 1 && yisint == 1) 313 z = -z; /* (-ve)**(odd int) */ 314 return (z); 315 } 316 #else 317 #error Unsupported Architecture 318 #endif /* defined(__sparc) */ 319