1 /*- 2 * ==================================================== 3 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 4 * 5 * Developed at SunPro, a Sun Microsystems, Inc. business. 6 * Permission to use, copy, modify, and distribute this 7 * software is freely granted, provided that this notice 8 * is preserved. 9 * ==================================================== 10 */ 11 12 /* 13 * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> 14 * 15 * Permission to use, copy, modify, and distribute this software for any 16 * purpose with or without fee is hereby granted, provided that the above 17 * copyright notice and this permission notice appear in all copies. 18 * 19 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES 20 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF 21 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR 22 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES 23 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN 24 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF 25 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. 26 */ 27 28 /* powl(x,y) return x**y 29 * 30 * n 31 * Method: Let x = 2 * (1+f) 32 * 1. Compute and return log2(x) in two pieces: 33 * log2(x) = w1 + w2, 34 * where w1 has 113-53 = 60 bit trailing zeros. 35 * 2. Perform y*log2(x) = n+y' by simulating multi-precision 36 * arithmetic, where |y'|<=0.5. 37 * 3. Return x**y = 2**n*exp(y'*log2) 38 * 39 * Special cases: 40 * 1. (anything) ** 0 is 1 41 * 2. (anything) ** 1 is itself 42 * 3. (anything) ** NAN is NAN 43 * 4. NAN ** (anything except 0) is NAN 44 * 5. +-(|x| > 1) ** +INF is +INF 45 * 6. +-(|x| > 1) ** -INF is +0 46 * 7. +-(|x| < 1) ** +INF is +0 47 * 8. +-(|x| < 1) ** -INF is +INF 48 * 9. +-1 ** +-INF is NAN 49 * 10. +0 ** (+anything except 0, NAN) is +0 50 * 11. -0 ** (+anything except 0, NAN, odd integer) is +0 51 * 12. +0 ** (-anything except 0, NAN) is +INF 52 * 13. -0 ** (-anything except 0, NAN, odd integer) is +INF 53 * 14. -0 ** (odd integer) = -( +0 ** (odd integer) ) 54 * 15. +INF ** (+anything except 0,NAN) is +INF 55 * 16. +INF ** (-anything except 0,NAN) is +0 56 * 17. -INF ** (anything) = -0 ** (-anything) 57 * 18. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer) 58 * 19. (-anything except 0 and inf) ** (non-integer) is NAN 59 * 60 */ 61 62 #include <sys/cdefs.h> 63 #include <float.h> 64 #include <math.h> 65 66 #include "math_private.h" 67 68 static const long double bp[] = { 69 1.0L, 70 1.5L, 71 }; 72 73 /* log_2(1.5) */ 74 static const long double dp_h[] = { 75 0.0, 76 5.8496250072115607565592654282227158546448E-1L 77 }; 78 79 /* Low part of log_2(1.5) */ 80 static const long double dp_l[] = { 81 0.0, 82 1.0579781240112554492329533686862998106046E-16L 83 }; 84 85 static const long double zero = 0.0L, 86 one = 1.0L, 87 two = 2.0L, 88 two113 = 1.0384593717069655257060992658440192E34L, 89 huge = 1.0e3000L, 90 tiny = 1.0e-3000L; 91 92 /* 3/2 log x = 3 z + z^3 + z^3 (z^2 R(z^2)) 93 z = (x-1)/(x+1) 94 1 <= x <= 1.25 95 Peak relative error 2.3e-37 */ 96 static const long double LN[] = 97 { 98 -3.0779177200290054398792536829702930623200E1L, 99 6.5135778082209159921251824580292116201640E1L, 100 -4.6312921812152436921591152809994014413540E1L, 101 1.2510208195629420304615674658258363295208E1L, 102 -9.9266909031921425609179910128531667336670E-1L 103 }; 104 static const long double LD[] = 105 { 106 -5.129862866715009066465422805058933131960E1L, 107 1.452015077564081884387441590064272782044E2L, 108 -1.524043275549860505277434040464085593165E2L, 109 7.236063513651544224319663428634139768808E1L, 110 -1.494198912340228235853027849917095580053E1L 111 /* 1.0E0 */ 112 }; 113 114 /* exp(x) = 1 + x - x / (1 - 2 / (x - x^2 R(x^2))) 115 0 <= x <= 0.5 116 Peak relative error 5.7e-38 */ 117 static const long double PN[] = 118 { 119 5.081801691915377692446852383385968225675E8L, 120 9.360895299872484512023336636427675327355E6L, 121 4.213701282274196030811629773097579432957E4L, 122 5.201006511142748908655720086041570288182E1L, 123 9.088368420359444263703202925095675982530E-3L, 124 }; 125 static const long double PD[] = 126 { 127 3.049081015149226615468111430031590411682E9L, 128 1.069833887183886839966085436512368982758E8L, 129 8.259257717868875207333991924545445705394E5L, 130 1.872583833284143212651746812884298360922E3L, 131 /* 1.0E0 */ 132 }; 133 134 static const long double 135 /* ln 2 */ 136 lg2 = 6.9314718055994530941723212145817656807550E-1L, 137 lg2_h = 6.9314718055994528622676398299518041312695E-1L, 138 lg2_l = 2.3190468138462996154948554638754786504121E-17L, 139 ovt = 8.0085662595372944372e-0017L, 140 /* 2/(3*log(2)) */ 141 cp = 9.6179669392597560490661645400126142495110E-1L, 142 cp_h = 9.6179669392597555432899980587535537779331E-1L, 143 cp_l = 5.0577616648125906047157785230014751039424E-17L; 144 145 long double 146 powl(long double x, long double y) 147 { 148 long double z, ax, z_h, z_l, p_h, p_l; 149 long double yy1, t1, t2, r, s, t, u, v, w; 150 long double s2, s_h, s_l, t_h, t_l; 151 int32_t i, j, k, yisint, n; 152 u_int32_t ix, iy; 153 int32_t hx, hy; 154 ieee_quad_shape_type o, p, q; 155 156 p.value = x; 157 hx = p.parts32.mswhi; 158 ix = hx & 0x7fffffff; 159 160 q.value = y; 161 hy = q.parts32.mswhi; 162 iy = hy & 0x7fffffff; 163 164 165 /* y==zero: x**0 = 1 */ 166 if ((iy | q.parts32.mswlo | q.parts32.lswhi | q.parts32.lswlo) == 0) 167 return one; 168 169 /* 1.0**y = 1; -1.0**+-Inf = 1 */ 170 if (x == one) 171 return one; 172 if (x == -1.0L && iy == 0x7fff0000 173 && (q.parts32.mswlo | q.parts32.lswhi | q.parts32.lswlo) == 0) 174 return one; 175 176 /* +-NaN return x+y */ 177 if ((ix > 0x7fff0000) 178 || ((ix == 0x7fff0000) 179 && ((p.parts32.mswlo | p.parts32.lswhi | p.parts32.lswlo) != 0)) 180 || (iy > 0x7fff0000) 181 || ((iy == 0x7fff0000) 182 && ((q.parts32.mswlo | q.parts32.lswhi | q.parts32.lswlo) != 0))) 183 return nan_mix(x, y); 184 185 /* determine if y is an odd int when x < 0 186 * yisint = 0 ... y is not an integer 187 * yisint = 1 ... y is an odd int 188 * yisint = 2 ... y is an even int 189 */ 190 yisint = 0; 191 if (hx < 0) 192 { 193 if (iy >= 0x40700000) /* 2^113 */ 194 yisint = 2; /* even integer y */ 195 else if (iy >= 0x3fff0000) /* 1.0 */ 196 { 197 if (floorl (y) == y) 198 { 199 z = 0.5 * y; 200 if (floorl (z) == z) 201 yisint = 2; 202 else 203 yisint = 1; 204 } 205 } 206 } 207 208 /* special value of y */ 209 if ((q.parts32.mswlo | q.parts32.lswhi | q.parts32.lswlo) == 0) 210 { 211 if (iy == 0x7fff0000) /* y is +-inf */ 212 { 213 if (((ix - 0x3fff0000) | p.parts32.mswlo | p.parts32.lswhi | 214 p.parts32.lswlo) == 0) 215 return y - y; /* +-1**inf is NaN */ 216 else if (ix >= 0x3fff0000) /* (|x|>1)**+-inf = inf,0 */ 217 return (hy >= 0) ? y : zero; 218 else /* (|x|<1)**-,+inf = inf,0 */ 219 return (hy < 0) ? -y : zero; 220 } 221 if (iy == 0x3fff0000) 222 { /* y is +-1 */ 223 if (hy < 0) 224 return one / x; 225 else 226 return x; 227 } 228 if (hy == 0x40000000) 229 return x * x; /* y is 2 */ 230 if (hy == 0x3ffe0000) 231 { /* y is 0.5 */ 232 if (hx >= 0) /* x >= +0 */ 233 return sqrtl (x); 234 } 235 } 236 237 ax = fabsl (x); 238 /* special value of x */ 239 if ((p.parts32.mswlo | p.parts32.lswhi | p.parts32.lswlo) == 0) 240 { 241 if (ix == 0x7fff0000 || ix == 0 || ix == 0x3fff0000) 242 { 243 z = ax; /*x is +-0,+-inf,+-1 */ 244 if (hy < 0) 245 z = one / z; /* z = (1/|x|) */ 246 if (hx < 0) 247 { 248 if (((ix - 0x3fff0000) | yisint) == 0) 249 { 250 z = (z - z) / (z - z); /* (-1)**non-int is NaN */ 251 } 252 else if (yisint == 1) 253 z = -z; /* (x<0)**odd = -(|x|**odd) */ 254 } 255 return z; 256 } 257 } 258 259 /* (x<0)**(non-int) is NaN */ 260 if (((((u_int32_t) hx >> 31) - 1) | yisint) == 0) 261 return (x - x) / (x - x); 262 263 /* |y| is huge. 264 2^-16495 = 1/2 of smallest representable value. 265 If (1 - 1/131072)^y underflows, y > 1.4986e9 */ 266 if (iy > 0x401d654b) 267 { 268 /* if (1 - 2^-113)^y underflows, y > 1.1873e38 */ 269 if (iy > 0x407d654b) 270 { 271 if (ix <= 0x3ffeffff) 272 return (hy < 0) ? huge * huge : tiny * tiny; 273 if (ix >= 0x3fff0000) 274 return (hy > 0) ? huge * huge : tiny * tiny; 275 } 276 /* over/underflow if x is not close to one */ 277 if (ix < 0x3ffeffff) 278 return (hy < 0) ? huge * huge : tiny * tiny; 279 if (ix > 0x3fff0000) 280 return (hy > 0) ? huge * huge : tiny * tiny; 281 } 282 283 n = 0; 284 /* take care subnormal number */ 285 if (ix < 0x00010000) 286 { 287 ax *= two113; 288 n -= 113; 289 o.value = ax; 290 ix = o.parts32.mswhi; 291 } 292 n += ((ix) >> 16) - 0x3fff; 293 j = ix & 0x0000ffff; 294 /* determine interval */ 295 ix = j | 0x3fff0000; /* normalize ix */ 296 if (j <= 0x3988) 297 k = 0; /* |x|<sqrt(3/2) */ 298 else if (j < 0xbb67) 299 k = 1; /* |x|<sqrt(3) */ 300 else 301 { 302 k = 0; 303 n += 1; 304 ix -= 0x00010000; 305 } 306 307 o.value = ax; 308 o.parts32.mswhi = ix; 309 ax = o.value; 310 311 /* compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */ 312 u = ax - bp[k]; /* bp[0]=1.0, bp[1]=1.5 */ 313 v = one / (ax + bp[k]); 314 s = u * v; 315 s_h = s; 316 317 o.value = s_h; 318 o.parts32.lswlo = 0; 319 o.parts32.lswhi &= 0xf8000000; 320 s_h = o.value; 321 /* t_h=ax+bp[k] High */ 322 t_h = ax + bp[k]; 323 o.value = t_h; 324 o.parts32.lswlo = 0; 325 o.parts32.lswhi &= 0xf8000000; 326 t_h = o.value; 327 t_l = ax - (t_h - bp[k]); 328 s_l = v * ((u - s_h * t_h) - s_h * t_l); 329 /* compute log(ax) */ 330 s2 = s * s; 331 u = LN[0] + s2 * (LN[1] + s2 * (LN[2] + s2 * (LN[3] + s2 * LN[4]))); 332 v = LD[0] + s2 * (LD[1] + s2 * (LD[2] + s2 * (LD[3] + s2 * (LD[4] + s2)))); 333 r = s2 * s2 * u / v; 334 r += s_l * (s_h + s); 335 s2 = s_h * s_h; 336 t_h = 3.0 + s2 + r; 337 o.value = t_h; 338 o.parts32.lswlo = 0; 339 o.parts32.lswhi &= 0xf8000000; 340 t_h = o.value; 341 t_l = r - ((t_h - 3.0) - s2); 342 /* u+v = s*(1+...) */ 343 u = s_h * t_h; 344 v = s_l * t_h + t_l * s; 345 /* 2/(3log2)*(s+...) */ 346 p_h = u + v; 347 o.value = p_h; 348 o.parts32.lswlo = 0; 349 o.parts32.lswhi &= 0xf8000000; 350 p_h = o.value; 351 p_l = v - (p_h - u); 352 z_h = cp_h * p_h; /* cp_h+cp_l = 2/(3*log2) */ 353 z_l = cp_l * p_h + p_l * cp + dp_l[k]; 354 /* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */ 355 t = (long double) n; 356 t1 = (((z_h + z_l) + dp_h[k]) + t); 357 o.value = t1; 358 o.parts32.lswlo = 0; 359 o.parts32.lswhi &= 0xf8000000; 360 t1 = o.value; 361 t2 = z_l - (((t1 - t) - dp_h[k]) - z_h); 362 363 /* s (sign of result -ve**odd) = -1 else = 1 */ 364 s = one; 365 if (((((u_int32_t) hx >> 31) - 1) | (yisint - 1)) == 0) 366 s = -one; /* (-ve)**(odd int) */ 367 368 /* split up y into yy1+y2 and compute (yy1+y2)*(t1+t2) */ 369 yy1 = y; 370 o.value = yy1; 371 o.parts32.lswlo = 0; 372 o.parts32.lswhi &= 0xf8000000; 373 yy1 = o.value; 374 p_l = (y - yy1) * t1 + y * t2; 375 p_h = yy1 * t1; 376 z = p_l + p_h; 377 o.value = z; 378 j = o.parts32.mswhi; 379 if (j >= 0x400d0000) /* z >= 16384 */ 380 { 381 /* if z > 16384 */ 382 if (((j - 0x400d0000) | o.parts32.mswlo | o.parts32.lswhi | 383 o.parts32.lswlo) != 0) 384 return s * huge * huge; /* overflow */ 385 else 386 { 387 if (p_l + ovt > z - p_h) 388 return s * huge * huge; /* overflow */ 389 } 390 } 391 else if ((j & 0x7fffffff) >= 0x400d01b9) /* z <= -16495 */ 392 { 393 /* z < -16495 */ 394 if (((j - 0xc00d01bc) | o.parts32.mswlo | o.parts32.lswhi | 395 o.parts32.lswlo) 396 != 0) 397 return s * tiny * tiny; /* underflow */ 398 else 399 { 400 if (p_l <= z - p_h) 401 return s * tiny * tiny; /* underflow */ 402 } 403 } 404 /* compute 2**(p_h+p_l) */ 405 i = j & 0x7fffffff; 406 k = (i >> 16) - 0x3fff; 407 n = 0; 408 if (i > 0x3ffe0000) 409 { /* if |z| > 0.5, set n = [z+0.5] */ 410 n = floorl (z + 0.5L); 411 t = n; 412 p_h -= t; 413 } 414 t = p_l + p_h; 415 o.value = t; 416 o.parts32.lswlo = 0; 417 o.parts32.lswhi &= 0xf8000000; 418 t = o.value; 419 u = t * lg2_h; 420 v = (p_l - (t - p_h)) * lg2 + t * lg2_l; 421 z = u + v; 422 w = v - (z - u); 423 /* exp(z) */ 424 t = z * z; 425 u = PN[0] + t * (PN[1] + t * (PN[2] + t * (PN[3] + t * PN[4]))); 426 v = PD[0] + t * (PD[1] + t * (PD[2] + t * (PD[3] + t))); 427 t1 = z - t * u / v; 428 r = (z * t1) / (t1 - two) - (w + z * w); 429 z = one - (r - z); 430 o.value = z; 431 j = o.parts32.mswhi; 432 j += (n << 16); 433 if ((j >> 16) <= 0) 434 z = scalbnl (z, n); /* subnormal output */ 435 else 436 { 437 o.parts32.mswhi = j; 438 z = o.value; 439 } 440 return s * z; 441 } 442