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 __FBSDID("$FreeBSD$"); 64 65 #include <float.h> 66 #include <math.h> 67 68 #include "math_private.h" 69 70 static const long double bp[] = { 71 1.0L, 72 1.5L, 73 }; 74 75 /* log_2(1.5) */ 76 static const long double dp_h[] = { 77 0.0, 78 5.8496250072115607565592654282227158546448E-1L 79 }; 80 81 /* Low part of log_2(1.5) */ 82 static const long double dp_l[] = { 83 0.0, 84 1.0579781240112554492329533686862998106046E-16L 85 }; 86 87 static const long double zero = 0.0L, 88 one = 1.0L, 89 two = 2.0L, 90 two113 = 1.0384593717069655257060992658440192E34L, 91 huge = 1.0e3000L, 92 tiny = 1.0e-3000L; 93 94 /* 3/2 log x = 3 z + z^3 + z^3 (z^2 R(z^2)) 95 z = (x-1)/(x+1) 96 1 <= x <= 1.25 97 Peak relative error 2.3e-37 */ 98 static const long double LN[] = 99 { 100 -3.0779177200290054398792536829702930623200E1L, 101 6.5135778082209159921251824580292116201640E1L, 102 -4.6312921812152436921591152809994014413540E1L, 103 1.2510208195629420304615674658258363295208E1L, 104 -9.9266909031921425609179910128531667336670E-1L 105 }; 106 static const long double LD[] = 107 { 108 -5.129862866715009066465422805058933131960E1L, 109 1.452015077564081884387441590064272782044E2L, 110 -1.524043275549860505277434040464085593165E2L, 111 7.236063513651544224319663428634139768808E1L, 112 -1.494198912340228235853027849917095580053E1L 113 /* 1.0E0 */ 114 }; 115 116 /* exp(x) = 1 + x - x / (1 - 2 / (x - x^2 R(x^2))) 117 0 <= x <= 0.5 118 Peak relative error 5.7e-38 */ 119 static const long double PN[] = 120 { 121 5.081801691915377692446852383385968225675E8L, 122 9.360895299872484512023336636427675327355E6L, 123 4.213701282274196030811629773097579432957E4L, 124 5.201006511142748908655720086041570288182E1L, 125 9.088368420359444263703202925095675982530E-3L, 126 }; 127 static const long double PD[] = 128 { 129 3.049081015149226615468111430031590411682E9L, 130 1.069833887183886839966085436512368982758E8L, 131 8.259257717868875207333991924545445705394E5L, 132 1.872583833284143212651746812884298360922E3L, 133 /* 1.0E0 */ 134 }; 135 136 static const long double 137 /* ln 2 */ 138 lg2 = 6.9314718055994530941723212145817656807550E-1L, 139 lg2_h = 6.9314718055994528622676398299518041312695E-1L, 140 lg2_l = 2.3190468138462996154948554638754786504121E-17L, 141 ovt = 8.0085662595372944372e-0017L, 142 /* 2/(3*log(2)) */ 143 cp = 9.6179669392597560490661645400126142495110E-1L, 144 cp_h = 9.6179669392597555432899980587535537779331E-1L, 145 cp_l = 5.0577616648125906047157785230014751039424E-17L; 146 147 long double 148 powl(long double x, long double y) 149 { 150 long double z, ax, z_h, z_l, p_h, p_l; 151 long double yy1, t1, t2, r, s, t, u, v, w; 152 long double s2, s_h, s_l, t_h, t_l; 153 int32_t i, j, k, yisint, n; 154 u_int32_t ix, iy; 155 int32_t hx, hy; 156 ieee_quad_shape_type o, p, q; 157 158 p.value = x; 159 hx = p.parts32.mswhi; 160 ix = hx & 0x7fffffff; 161 162 q.value = y; 163 hy = q.parts32.mswhi; 164 iy = hy & 0x7fffffff; 165 166 167 /* y==zero: x**0 = 1 */ 168 if ((iy | q.parts32.mswlo | q.parts32.lswhi | q.parts32.lswlo) == 0) 169 return one; 170 171 /* 1.0**y = 1; -1.0**+-Inf = 1 */ 172 if (x == one) 173 return one; 174 if (x == -1.0L && iy == 0x7fff0000 175 && (q.parts32.mswlo | q.parts32.lswhi | q.parts32.lswlo) == 0) 176 return one; 177 178 /* +-NaN return x+y */ 179 if ((ix > 0x7fff0000) 180 || ((ix == 0x7fff0000) 181 && ((p.parts32.mswlo | p.parts32.lswhi | p.parts32.lswlo) != 0)) 182 || (iy > 0x7fff0000) 183 || ((iy == 0x7fff0000) 184 && ((q.parts32.mswlo | q.parts32.lswhi | q.parts32.lswlo) != 0))) 185 return nan_mix(x, y); 186 187 /* determine if y is an odd int when x < 0 188 * yisint = 0 ... y is not an integer 189 * yisint = 1 ... y is an odd int 190 * yisint = 2 ... y is an even int 191 */ 192 yisint = 0; 193 if (hx < 0) 194 { 195 if (iy >= 0x40700000) /* 2^113 */ 196 yisint = 2; /* even integer y */ 197 else if (iy >= 0x3fff0000) /* 1.0 */ 198 { 199 if (floorl (y) == y) 200 { 201 z = 0.5 * y; 202 if (floorl (z) == z) 203 yisint = 2; 204 else 205 yisint = 1; 206 } 207 } 208 } 209 210 /* special value of y */ 211 if ((q.parts32.mswlo | q.parts32.lswhi | q.parts32.lswlo) == 0) 212 { 213 if (iy == 0x7fff0000) /* y is +-inf */ 214 { 215 if (((ix - 0x3fff0000) | p.parts32.mswlo | p.parts32.lswhi | 216 p.parts32.lswlo) == 0) 217 return y - y; /* +-1**inf is NaN */ 218 else if (ix >= 0x3fff0000) /* (|x|>1)**+-inf = inf,0 */ 219 return (hy >= 0) ? y : zero; 220 else /* (|x|<1)**-,+inf = inf,0 */ 221 return (hy < 0) ? -y : zero; 222 } 223 if (iy == 0x3fff0000) 224 { /* y is +-1 */ 225 if (hy < 0) 226 return one / x; 227 else 228 return x; 229 } 230 if (hy == 0x40000000) 231 return x * x; /* y is 2 */ 232 if (hy == 0x3ffe0000) 233 { /* y is 0.5 */ 234 if (hx >= 0) /* x >= +0 */ 235 return sqrtl (x); 236 } 237 } 238 239 ax = fabsl (x); 240 /* special value of x */ 241 if ((p.parts32.mswlo | p.parts32.lswhi | p.parts32.lswlo) == 0) 242 { 243 if (ix == 0x7fff0000 || ix == 0 || ix == 0x3fff0000) 244 { 245 z = ax; /*x is +-0,+-inf,+-1 */ 246 if (hy < 0) 247 z = one / z; /* z = (1/|x|) */ 248 if (hx < 0) 249 { 250 if (((ix - 0x3fff0000) | yisint) == 0) 251 { 252 z = (z - z) / (z - z); /* (-1)**non-int is NaN */ 253 } 254 else if (yisint == 1) 255 z = -z; /* (x<0)**odd = -(|x|**odd) */ 256 } 257 return z; 258 } 259 } 260 261 /* (x<0)**(non-int) is NaN */ 262 if (((((u_int32_t) hx >> 31) - 1) | yisint) == 0) 263 return (x - x) / (x - x); 264 265 /* |y| is huge. 266 2^-16495 = 1/2 of smallest representable value. 267 If (1 - 1/131072)^y underflows, y > 1.4986e9 */ 268 if (iy > 0x401d654b) 269 { 270 /* if (1 - 2^-113)^y underflows, y > 1.1873e38 */ 271 if (iy > 0x407d654b) 272 { 273 if (ix <= 0x3ffeffff) 274 return (hy < 0) ? huge * huge : tiny * tiny; 275 if (ix >= 0x3fff0000) 276 return (hy > 0) ? huge * huge : tiny * tiny; 277 } 278 /* over/underflow if x is not close to one */ 279 if (ix < 0x3ffeffff) 280 return (hy < 0) ? huge * huge : tiny * tiny; 281 if (ix > 0x3fff0000) 282 return (hy > 0) ? huge * huge : tiny * tiny; 283 } 284 285 n = 0; 286 /* take care subnormal number */ 287 if (ix < 0x00010000) 288 { 289 ax *= two113; 290 n -= 113; 291 o.value = ax; 292 ix = o.parts32.mswhi; 293 } 294 n += ((ix) >> 16) - 0x3fff; 295 j = ix & 0x0000ffff; 296 /* determine interval */ 297 ix = j | 0x3fff0000; /* normalize ix */ 298 if (j <= 0x3988) 299 k = 0; /* |x|<sqrt(3/2) */ 300 else if (j < 0xbb67) 301 k = 1; /* |x|<sqrt(3) */ 302 else 303 { 304 k = 0; 305 n += 1; 306 ix -= 0x00010000; 307 } 308 309 o.value = ax; 310 o.parts32.mswhi = ix; 311 ax = o.value; 312 313 /* compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */ 314 u = ax - bp[k]; /* bp[0]=1.0, bp[1]=1.5 */ 315 v = one / (ax + bp[k]); 316 s = u * v; 317 s_h = s; 318 319 o.value = s_h; 320 o.parts32.lswlo = 0; 321 o.parts32.lswhi &= 0xf8000000; 322 s_h = o.value; 323 /* t_h=ax+bp[k] High */ 324 t_h = ax + bp[k]; 325 o.value = t_h; 326 o.parts32.lswlo = 0; 327 o.parts32.lswhi &= 0xf8000000; 328 t_h = o.value; 329 t_l = ax - (t_h - bp[k]); 330 s_l = v * ((u - s_h * t_h) - s_h * t_l); 331 /* compute log(ax) */ 332 s2 = s * s; 333 u = LN[0] + s2 * (LN[1] + s2 * (LN[2] + s2 * (LN[3] + s2 * LN[4]))); 334 v = LD[0] + s2 * (LD[1] + s2 * (LD[2] + s2 * (LD[3] + s2 * (LD[4] + s2)))); 335 r = s2 * s2 * u / v; 336 r += s_l * (s_h + s); 337 s2 = s_h * s_h; 338 t_h = 3.0 + s2 + r; 339 o.value = t_h; 340 o.parts32.lswlo = 0; 341 o.parts32.lswhi &= 0xf8000000; 342 t_h = o.value; 343 t_l = r - ((t_h - 3.0) - s2); 344 /* u+v = s*(1+...) */ 345 u = s_h * t_h; 346 v = s_l * t_h + t_l * s; 347 /* 2/(3log2)*(s+...) */ 348 p_h = u + v; 349 o.value = p_h; 350 o.parts32.lswlo = 0; 351 o.parts32.lswhi &= 0xf8000000; 352 p_h = o.value; 353 p_l = v - (p_h - u); 354 z_h = cp_h * p_h; /* cp_h+cp_l = 2/(3*log2) */ 355 z_l = cp_l * p_h + p_l * cp + dp_l[k]; 356 /* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */ 357 t = (long double) n; 358 t1 = (((z_h + z_l) + dp_h[k]) + t); 359 o.value = t1; 360 o.parts32.lswlo = 0; 361 o.parts32.lswhi &= 0xf8000000; 362 t1 = o.value; 363 t2 = z_l - (((t1 - t) - dp_h[k]) - z_h); 364 365 /* s (sign of result -ve**odd) = -1 else = 1 */ 366 s = one; 367 if (((((u_int32_t) hx >> 31) - 1) | (yisint - 1)) == 0) 368 s = -one; /* (-ve)**(odd int) */ 369 370 /* split up y into yy1+y2 and compute (yy1+y2)*(t1+t2) */ 371 yy1 = y; 372 o.value = yy1; 373 o.parts32.lswlo = 0; 374 o.parts32.lswhi &= 0xf8000000; 375 yy1 = o.value; 376 p_l = (y - yy1) * t1 + y * t2; 377 p_h = yy1 * t1; 378 z = p_l + p_h; 379 o.value = z; 380 j = o.parts32.mswhi; 381 if (j >= 0x400d0000) /* z >= 16384 */ 382 { 383 /* if z > 16384 */ 384 if (((j - 0x400d0000) | o.parts32.mswlo | o.parts32.lswhi | 385 o.parts32.lswlo) != 0) 386 return s * huge * huge; /* overflow */ 387 else 388 { 389 if (p_l + ovt > z - p_h) 390 return s * huge * huge; /* overflow */ 391 } 392 } 393 else if ((j & 0x7fffffff) >= 0x400d01b9) /* z <= -16495 */ 394 { 395 /* z < -16495 */ 396 if (((j - 0xc00d01bc) | o.parts32.mswlo | o.parts32.lswhi | 397 o.parts32.lswlo) 398 != 0) 399 return s * tiny * tiny; /* underflow */ 400 else 401 { 402 if (p_l <= z - p_h) 403 return s * tiny * tiny; /* underflow */ 404 } 405 } 406 /* compute 2**(p_h+p_l) */ 407 i = j & 0x7fffffff; 408 k = (i >> 16) - 0x3fff; 409 n = 0; 410 if (i > 0x3ffe0000) 411 { /* if |z| > 0.5, set n = [z+0.5] */ 412 n = floorl (z + 0.5L); 413 t = n; 414 p_h -= t; 415 } 416 t = p_l + p_h; 417 o.value = t; 418 o.parts32.lswlo = 0; 419 o.parts32.lswhi &= 0xf8000000; 420 t = o.value; 421 u = t * lg2_h; 422 v = (p_l - (t - p_h)) * lg2 + t * lg2_l; 423 z = u + v; 424 w = v - (z - u); 425 /* exp(z) */ 426 t = z * z; 427 u = PN[0] + t * (PN[1] + t * (PN[2] + t * (PN[3] + t * PN[4]))); 428 v = PD[0] + t * (PD[1] + t * (PD[2] + t * (PD[3] + t))); 429 t1 = z - t * u / v; 430 r = (z * t1) / (t1 - two) - (w + z * w); 431 z = one - (r - z); 432 o.value = z; 433 j = o.parts32.mswhi; 434 j += (n << 16); 435 if ((j >> 16) <= 0) 436 z = scalbnl (z, n); /* subnormal output */ 437 else 438 { 439 o.parts32.mswhi = j; 440 z = o.value; 441 } 442 return s * z; 443 } 444