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