libquadmath/math/powq.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 /* Expansions and modifications for 128-bit long double are
  13    Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov>
  14    and are incorporated herein by permission of the author.  The author
  15    reserves the right to distribute this material elsewhere under different
  16    copying permissions.  These modifications are distributed here under
  17    the following terms:
  18
  19     This library is free software; you can redistribute it and/or
  20     modify it under the terms of the GNU Lesser General Public
  21     License as published by the Free Software Foundation; either
  22     version 2.1 of the License, or (at your option) any later version.
  23
  24     This library is distributed in the hope that it will be useful,
  25     but WITHOUT ANY WARRANTY; without even the implied warranty of
  26     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  27     Lesser General Public License for more details.
  28
  29     You should have received a copy of the GNU Lesser General Public
  30     License along with this library; if not, write to the Free Software
  31     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA */
  32
  33 /* powq(x,y) return x**y
  34  *
  35  *                    n
  36  * Method:  Let x =  2   * (1+f)
  37  *      1. Compute and return log2(x) in two pieces:
  38  *              log2(x) = w1 + w2,
  39  *         where w1 has 113-53 = 60 bit trailing zeros.
  40  *      2. Perform y*log2(x) = n+y' by simulating muti-precision
  41  *         arithmetic, where |y'|<=0.5.
  42  *      3. Return x**y = 2**n*exp(y'*log2)
  43  *
  44  * Special cases:
  45  *      1.  (anything) ** 0  is 1
  46  *      2.  (anything) ** 1  is itself
  47  *      3.  (anything) ** NAN is NAN
  48  *      4.  NAN ** (anything except 0) is NAN
  49  *      5.  +-(|x| > 1) **  +INF is +INF
  50  *      6.  +-(|x| > 1) **  -INF is +0
  51  *      7.  +-(|x| < 1) **  +INF is +0
  52  *      8.  +-(|x| < 1) **  -INF is +INF
  53  *      9.  +-1         ** +-INF is NAN
  54  *      10. +0 ** (+anything except 0, NAN)               is +0
  55  *      11. -0 ** (+anything except 0, NAN, odd integer)  is +0
  56  *      12. +0 ** (-anything except 0, NAN)               is +INF
  57  *      13. -0 ** (-anything except 0, NAN, odd integer)  is +INF
  58  *      14. -0 ** (odd integer) = -( +0 ** (odd integer) )
  59  *      15. +INF ** (+anything except 0,NAN) is +INF
  60  *      16. +INF ** (-anything except 0,NAN) is +0
  61  *      17. -INF ** (anything)  = -0 ** (-anything)
  62  *      18. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer)
  63  *      19. (-anything except 0 and inf) ** (non-integer) is NAN
  64  *
  65  */
  66
  67 #include "quadmath-imp.h"
  68
  69 static const __float128 bp[] = {
  70   1.0Q,
  71   1.5Q,
  72 };
  73
  74 /* log_2(1.5) */
  75 static const __float128 dp_h[] = {
  76   0.0,
  77   5.8496250072115607565592654282227158546448E-1Q
  78 };
  79
  80 /* Low part of log_2(1.5) */
  81 static const __float128 dp_l[] = {
  82   0.0,
  83   1.0579781240112554492329533686862998106046E-16Q
  84 };
  85
  86 static const __float128 zero = 0.0Q,
  87   one = 1.0Q,
  88   two = 2.0Q,
  89   two113 = 1.0384593717069655257060992658440192E34Q,
  90   huge = 1.0e3000Q,
  91   tiny = 1.0e-3000Q;
  92
  93 /* 3/2 log x = 3 z + z^3 + z^3 (z^2 R(z^2))
  94    z = (x-1)/(x+1)
  95    1 <= x <= 1.25
  96    Peak relative error 2.3e-37 */
  97 static const __float128 LN[] =
  98 {
  99  -3.0779177200290054398792536829702930623200E1Q,
 100   6.5135778082209159921251824580292116201640E1Q,
 101  -4.6312921812152436921591152809994014413540E1Q,
 102   1.2510208195629420304615674658258363295208E1Q,
 103  -9.9266909031921425609179910128531667336670E-1Q
 104 };
 105 static const __float128 LD[] =
 106 {
 107  -5.129862866715009066465422805058933131960E1Q,
 108   1.452015077564081884387441590064272782044E2Q,
 109  -1.524043275549860505277434040464085593165E2Q,
 110   7.236063513651544224319663428634139768808E1Q,
 111  -1.494198912340228235853027849917095580053E1Q
 112   /* 1.0E0 */
 113 };
 114
 115 /* exp(x) = 1 + x - x / (1 - 2 / (x - x^2 R(x^2)))
 116    0 <= x <= 0.5
 117    Peak relative error 5.7e-38  */
 118 static const __float128 PN[] =
 119 {
 120   5.081801691915377692446852383385968225675E8Q,
 121   9.360895299872484512023336636427675327355E6Q,
 122   4.213701282274196030811629773097579432957E4Q,
 123   5.201006511142748908655720086041570288182E1Q,
 124   9.088368420359444263703202925095675982530E-3Q,
 125 };
 126 static const __float128 PD[] =
 127 {
 128   3.049081015149226615468111430031590411682E9Q,
 129   1.069833887183886839966085436512368982758E8Q,
 130   8.259257717868875207333991924545445705394E5Q,
 131   1.872583833284143212651746812884298360922E3Q,
 132   /* 1.0E0 */
 133 };
 134
 135 static const __float128
 136   /* ln 2 */
 137   lg2 = 6.9314718055994530941723212145817656807550E-1Q,
 138   lg2_h = 6.9314718055994528622676398299518041312695E-1Q,
 139   lg2_l = 2.3190468138462996154948554638754786504121E-17Q,
 140   ovt = 8.0085662595372944372e-0017Q,
 141   /* 2/(3*log(2)) */
 142   cp = 9.6179669392597560490661645400126142495110E-1Q,
 143   cp_h = 9.6179669392597555432899980587535537779331E-1Q,
 144   cp_l = 5.0577616648125906047157785230014751039424E-17Q;
 145
 146 __float128
 147 powq (__float128 x, __float128 y)
 148 {
 149   __float128 z, ax, z_h, z_l, p_h, p_l;
 150   __float128 y1, t1, t2, r, s, t, u, v, w;
 151   __float128 s2, s_h, s_l, t_h, t_l, ay;
 152   int32_t i, j, k, yisint, n;
 153   uint32_t ix, iy;
 154   int32_t hx, hy;
 155   ieee854_float128 o, p, q;
 156
 157   p.value = x;
 158   hx = p.words32.w0;
 159   ix = hx & 0x7fffffff;
 160
 161   q.value = y;
 162   hy = q.words32.w0;
 163   iy = hy & 0x7fffffff;
 164
 165
 166   /* y==zero: x**0 = 1 */
 167   if ((iy | q.words32.w1 | q.words32.w2 | q.words32.w3) == 0)
 168     return one;
 169
 170   /* 1.0**y = 1; -1.0**+-Inf = 1 */
 171   if (x == one)
 172     return one;
 173   if (x == -1.0Q && iy == 0x7fff0000
 174       && (q.words32.w1 | q.words32.w2 | q.words32.w3) == 0)
 175     return one;
 176
 177   /* +-NaN return x+y */
 178   if ((ix > 0x7fff0000)
 179       || ((ix == 0x7fff0000)
 180           && ((p.words32.w1 | p.words32.w2 | p.words32.w3) != 0))
 181       || (iy > 0x7fff0000)
 182       || ((iy == 0x7fff0000)
 183           && ((q.words32.w1 | q.words32.w2 | q.words32.w3) != 0)))
 184     return x + y;
 185
 186   /* determine if y is an odd int when x < 0
 187    * yisint = 0       ... y is not an integer
 188    * yisint = 1       ... y is an odd int
 189    * yisint = 2       ... y is an even int
 190    */
 191   yisint = 0;
 192   if (hx < 0)
 193     {
 194       if (iy >= 0x40700000)     /* 2^113 */
 195         yisint = 2;             /* even integer y */
 196       else if (iy >= 0x3fff0000)        /* 1.0 */
 197         {
 198           if (floorq (y) == y)
 199             {
 200               z = 0.5 * y;
 201               if (floorq (z) == z)
 202                 yisint = 2;
 203               else
 204                 yisint = 1;
 205             }
 206         }
 207     }
 208
 209   /* special value of y */
 210   if ((q.words32.w1 | q.words32.w2 | q.words32.w3) == 0)
 211     {
 212       if (iy == 0x7fff0000)     /* y is +-inf */
 213         {
 214           if (((ix - 0x3fff0000) | p.words32.w1 | p.words32.w2 | p.words32.w3)
 215               == 0)
 216             return y - y;       /* +-1**inf is NaN */
 217           else if (ix >= 0x3fff0000)    /* (|x|>1)**+-inf = inf,0 */
 218             return (hy >= 0) ? y : zero;
 219           else                  /* (|x|<1)**-,+inf = inf,0 */
 220             return (hy < 0) ? -y : zero;
 221         }
 222       if (iy == 0x3fff0000)
 223         {                       /* y is  +-1 */
 224           if (hy < 0)
 225             return one / x;
 226           else
 227             return x;
 228         }
 229       if (hy == 0x40000000)
 230         return x * x;           /* y is  2 */
 231       if (hy == 0x3ffe0000)
 232         {                       /* y is  0.5 */
 233           if (hx >= 0)          /* x >= +0 */
 234             return sqrtq (x);
 235         }
 236     }
 237
 238   ax = fabsq (x);
 239   /* special value of x */
 240   if ((p.words32.w1 | p.words32.w2 | p.words32.w3) == 0)
 241     {
 242       if (ix == 0x7fff0000 || ix == 0 || ix == 0x3fff0000)
 243         {
 244           z = ax;               /*x is +-0,+-inf,+-1 */
 245           if (hy < 0)
 246             z = one / z;        /* z = (1/|x|) */
 247           if (hx < 0)
 248             {
 249               if (((ix - 0x3fff0000) | yisint) == 0)
 250                 {
 251                   z = (z - z) / (z - z);        /* (-1)**non-int is NaN */
 252                 }
 253               else if (yisint == 1)
 254                 z = -z;         /* (x<0)**odd = -(|x|**odd) */
 255             }
 256           return z;
 257         }
 258     }
 259
 260   /* (x<0)**(non-int) is NaN */
 261   if (((((uint32_t) hx >> 31) - 1) | yisint) == 0)
 262     return (x - x) / (x - x);
 263
 264   /* |y| is huge.
 265      2^-16495 = 1/2 of smallest representable value.
 266      If (1 - 1/131072)^y underflows, y > 1.4986e9 */
 267   if (iy > 0x401d654b)
 268     {
 269       /* if (1 - 2^-113)^y underflows, y > 1.1873e38 */
 270       if (iy > 0x407d654b)
 271         {
 272           if (ix <= 0x3ffeffff)
 273             return (hy < 0) ? huge * huge : tiny * tiny;
 274           if (ix >= 0x3fff0000)
 275             return (hy > 0) ? huge * huge : tiny * tiny;
 276         }
 277       /* over/underflow if x is not close to one */
 278       if (ix < 0x3ffeffff)
 279         return (hy < 0) ? huge * huge : tiny * tiny;
 280       if (ix > 0x3fff0000)
 281         return (hy > 0) ? huge * huge : tiny * tiny;
 282     }
 283
 284   ay = y > 0 ? y : -y;
 285   if (ay < 0x1p-128)
 286     y = y < 0 ? -0x1p-128 : 0x1p-128;
 287
 288   n = 0;
 289   /* take care subnormal number */
 290   if (ix < 0x00010000)
 291     {
 292       ax *= two113;
 293       n -= 113;
 294       o.value = ax;
 295       ix = o.words32.w0;
 296     }
 297   n += ((ix) >> 16) - 0x3fff;
 298   j = ix & 0x0000ffff;
 299   /* determine interval */
 300   ix = j | 0x3fff0000;          /* normalize ix */
 301   if (j <= 0x3988)
 302     k = 0;                      /* |x|<sqrt(3/2) */
 303   else if (j < 0xbb67)
 304     k = 1;                      /* |x|<sqrt(3)   */
 305   else
 306     {
 307       k = 0;
 308       n += 1;
 309       ix -= 0x00010000;
 310     }
 311
 312   o.value = ax;
 313   o.words32.w0 = ix;
 314   ax = o.value;
 315
 316   /* compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
 317   u = ax - bp[k];               /* bp[0]=1.0, bp[1]=1.5 */
 318   v = one / (ax + bp[k]);
 319   s = u * v;
 320   s_h = s;
 321
 322   o.value = s_h;
 323   o.words32.w3 = 0;
 324   o.words32.w2 &= 0xf8000000;
 325   s_h = o.value;
 326   /* t_h=ax+bp[k] High */
 327   t_h = ax + bp[k];
 328   o.value = t_h;
 329   o.words32.w3 = 0;
 330   o.words32.w2 &= 0xf8000000;
 331   t_h = o.value;
 332   t_l = ax - (t_h - bp[k]);
 333   s_l = v * ((u - s_h * t_h) - s_h * t_l);
 334   /* compute log(ax) */
 335   s2 = s * s;
 336   u = LN[0] + s2 * (LN[1] + s2 * (LN[2] + s2 * (LN[3] + s2 * LN[4])));
 337   v = LD[0] + s2 * (LD[1] + s2 * (LD[2] + s2 * (LD[3] + s2 * (LD[4] + s2))));
 338   r = s2 * s2 * u / v;
 339   r += s_l * (s_h + s);
 340   s2 = s_h * s_h;
 341   t_h = 3.0 + s2 + r;
 342   o.value = t_h;
 343   o.words32.w3 = 0;
 344   o.words32.w2 &= 0xf8000000;
 345   t_h = o.value;
 346   t_l = r - ((t_h - 3.0) - s2);
 347   /* u+v = s*(1+...) */
 348   u = s_h * t_h;
 349   v = s_l * t_h + t_l * s;
 350   /* 2/(3log2)*(s+...) */
 351   p_h = u + v;
 352   o.value = p_h;
 353   o.words32.w3 = 0;
 354   o.words32.w2 &= 0xf8000000;
 355   p_h = o.value;
 356   p_l = v - (p_h - u);
 357   z_h = cp_h * p_h;             /* cp_h+cp_l = 2/(3*log2) */
 358   z_l = cp_l * p_h + p_l * cp + dp_l[k];
 359   /* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */
 360   t = (__float128) n;
 361   t1 = (((z_h + z_l) + dp_h[k]) + t);
 362   o.value = t1;
 363   o.words32.w3 = 0;
 364   o.words32.w2 &= 0xf8000000;
 365   t1 = o.value;
 366   t2 = z_l - (((t1 - t) - dp_h[k]) - z_h);
 367
 368   /* s (sign of result -ve**odd) = -1 else = 1 */
 369   s = one;
 370   if (((((uint32_t) hx >> 31) - 1) | (yisint - 1)) == 0)
 371     s = -one;                   /* (-ve)**(odd int) */
 372
 373   /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */
 374   y1 = y;
 375   o.value = y1;
 376   o.words32.w3 = 0;
 377   o.words32.w2 &= 0xf8000000;
 378   y1 = o.value;
 379   p_l = (y - y1) * t1 + y * t2;
 380   p_h = y1 * t1;
 381   z = p_l + p_h;
 382   o.value = z;
 383   j = o.words32.w0;
 384   if (j >= 0x400d0000) /* z >= 16384 */
 385     {
 386       /* if z > 16384 */
 387       if (((j - 0x400d0000) | o.words32.w1 | o.words32.w2 | o.words32.w3) != 0)
 388         return s * huge * huge; /* overflow */
 389       else
 390         {
 391           if (p_l + ovt > z - p_h)
 392             return s * huge * huge;     /* overflow */
 393         }
 394     }
 395   else if ((j & 0x7fffffff) >= 0x400d01b9)      /* z <= -16495 */
 396     {
 397       /* z < -16495 */
 398       if (((j - 0xc00d01bc) | o.words32.w1 | o.words32.w2 | o.words32.w3)
 399           != 0)
 400         return s * tiny * tiny; /* underflow */
 401       else
 402         {
 403           if (p_l <= z - p_h)
 404             return s * tiny * tiny;     /* underflow */
 405         }
 406     }
 407   /* compute 2**(p_h+p_l) */
 408   i = j & 0x7fffffff;
 409   k = (i >> 16) - 0x3fff;
 410   n = 0;
 411   if (i > 0x3ffe0000)
 412     {                           /* if |z| > 0.5, set n = [z+0.5] */
 413       n = floorq (z + 0.5Q);
 414       t = n;
 415       p_h -= t;
 416     }
 417   t = p_l + p_h;
 418   o.value = t;
 419   o.words32.w3 = 0;
 420   o.words32.w2 &= 0xf8000000;
 421   t = o.value;
 422   u = t * lg2_h;
 423   v = (p_l - (t - p_h)) * lg2 + t * lg2_l;
 424   z = u + v;
 425   w = v - (z - u);
 426   /*  exp(z) */
 427   t = z * z;
 428   u = PN[0] + t * (PN[1] + t * (PN[2] + t * (PN[3] + t * PN[4])));
 429   v = PD[0] + t * (PD[1] + t * (PD[2] + t * (PD[3] + t)));
 430   t1 = z - t * u / v;
 431   r = (z * t1) / (t1 - two) - (w + z * w);
 432   z = one - (r - z);
 433   o.value = z;
 434   j = o.words32.w0;
 435   j += (n << 16);
 436   if ((j >> 16) <= 0)
 437     z = scalbnq (z, n); /* subnormal output */
 438   else
 439     {
 440       o.words32.w0 = j;
 441       z = o.value;
 442     }
 443   return s * z;
 444 }