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 /* __ieee754_powl(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;
 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   n = 0;
 285   /* take care subnormal number */
 286   if (ix < 0x00010000)
 287     {
 288       ax *= two113;
 289       n -= 113;
 290       o.value = ax;
 291       ix = o.words32.w0;
 292     }
 293   n += ((ix) >> 16) - 0x3fff;
 294   j = ix & 0x0000ffff;
 295   /* determine interval */
 296   ix = j | 0x3fff0000;          /* normalize ix */
 297   if (j <= 0x3988)
 298     k = 0;                      /* |x|<sqrt(3/2) */
 299   else if (j < 0xbb67)
 300     k = 1;                      /* |x|<sqrt(3)   */
 301   else
 302     {
 303       k = 0;
 304       n += 1;
 305       ix -= 0x00010000;
 306     }
 307
 308   o.value = ax;
 309   o.words32.w0 = ix;
 310   ax = o.value;
 311
 312   /* compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
 313   u = ax - bp[k];               /* bp[0]=1.0, bp[1]=1.5 */
 314   v = one / (ax + bp[k]);
 315   s = u * v;
 316   s_h = s;
 317
 318   o.value = s_h;
 319   o.words32.w3 = 0;
 320   o.words32.w2 &= 0xf8000000;
 321   s_h = o.value;
 322   /* t_h=ax+bp[k] High */
 323   t_h = ax + bp[k];
 324   o.value = t_h;
 325   o.words32.w3 = 0;
 326   o.words32.w2 &= 0xf8000000;
 327   t_h = o.value;
 328   t_l = ax - (t_h - bp[k]);
 329   s_l = v * ((u - s_h * t_h) - s_h * t_l);
 330   /* compute log(ax) */
 331   s2 = s * s;
 332   u = LN[0] + s2 * (LN[1] + s2 * (LN[2] + s2 * (LN[3] + s2 * LN[4])));
 333   v = LD[0] + s2 * (LD[1] + s2 * (LD[2] + s2 * (LD[3] + s2 * (LD[4] + s2))));
 334   r = s2 * s2 * u / v;
 335   r += s_l * (s_h + s);
 336   s2 = s_h * s_h;
 337   t_h = 3.0 + s2 + r;
 338   o.value = t_h;
 339   o.words32.w3 = 0;
 340   o.words32.w2 &= 0xf8000000;
 341   t_h = o.value;
 342   t_l = r - ((t_h - 3.0) - s2);
 343   /* u+v = s*(1+...) */
 344   u = s_h * t_h;
 345   v = s_l * t_h + t_l * s;
 346   /* 2/(3log2)*(s+...) */
 347   p_h = u + v;
 348   o.value = p_h;
 349   o.words32.w3 = 0;
 350   o.words32.w2 &= 0xf8000000;
 351   p_h = o.value;
 352   p_l = v - (p_h - u);
 353   z_h = cp_h * p_h;             /* cp_h+cp_l = 2/(3*log2) */
 354   z_l = cp_l * p_h + p_l * cp + dp_l[k];
 355   /* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */
 356   t = (__float128) n;
 357   t1 = (((z_h + z_l) + dp_h[k]) + t);
 358   o.value = t1;
 359   o.words32.w3 = 0;
 360   o.words32.w2 &= 0xf8000000;
 361   t1 = o.value;
 362   t2 = z_l - (((t1 - t) - dp_h[k]) - z_h);
 363
 364   /* s (sign of result -ve**odd) = -1 else = 1 */
 365   s = one;
 366   if (((((uint32_t) hx >> 31) - 1) | (yisint - 1)) == 0)
 367     s = -one;                   /* (-ve)**(odd int) */
 368
 369   /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */
 370   y1 = y;
 371   o.value = y1;
 372   o.words32.w3 = 0;
 373   o.words32.w2 &= 0xf8000000;
 374   y1 = o.value;
 375   p_l = (y - y1) * t1 + y * t2;
 376   p_h = y1 * t1;
 377   z = p_l + p_h;
 378   o.value = z;
 379   j = o.words32.w0;
 380   if (j >= 0x400d0000) /* z >= 16384 */
 381     {
 382       /* if z > 16384 */
 383       if (((j - 0x400d0000) | o.words32.w1 | o.words32.w2 | o.words32.w3) != 0)
 384         return s * huge * huge; /* overflow */
 385       else
 386         {
 387           if (p_l + ovt > z - p_h)
 388             return s * huge * huge;     /* overflow */
 389         }
 390     }
 391   else if ((j & 0x7fffffff) >= 0x400d01b9)      /* z <= -16495 */
 392     {
 393       /* z < -16495 */
 394       if (((j - 0xc00d01bc) | o.words32.w1 | o.words32.w2 | o.words32.w3)
 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 = floorq (z + 0.5Q);
 410       t = n;
 411       p_h -= t;
 412     }
 413   t = p_l + p_h;
 414   o.value = t;
 415   o.words32.w3 = 0;
 416   o.words32.w2 &= 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.words32.w0;
 431   j += (n << 16);
 432   if ((j >> 16) <= 0)
 433     z = scalbnq (z, n); /* subnormal output */
 434   else
 435     {
 436       o.words32.w0 = j;
 437       z = o.value;
 438     }
 439   return s * z;
 440 }