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