newlib/libm/ld128/e_powl.c

   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 "../ld/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 }