1 /*
   2  * CDDL HEADER START
   3  *
   4  * The contents of this file are subject to the terms of the
   5  * Common Development and Distribution License (the "License").
   6  * You may not use this file except in compliance with the License.
   7  *
   8  * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
   9  * or http://www.opensolaris.org/os/licensing.
  10  * See the License for the specific language governing permissions
  11  * and limitations under the License.
  12  *
  13  * When distributing Covered Code, include this CDDL HEADER in each
  14  * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
  15  * If applicable, add the following below this CDDL HEADER, with the
  16  * fields enclosed by brackets "[]" replaced with your own identifying
  17  * information: Portions Copyright [yyyy] [name of copyright owner]
  18  *
  19  * CDDL HEADER END
  20  */
  21 
  22 /*
  23  * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
  24  */
  25 
  26 /*
  27  * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
  28  * Use is subject to license terms.
  29  */
  30 
  31 #pragma weak __powl = powl
  32 
  33 #include "libm.h"
  34 #include "xpg6.h"                       /* __xpg6 */
  35 #define _C99SUSv3_pow           _C99SUSv3_pow_treats_Inf_as_an_even_int
  36 
  37 #if defined(__sparc)
  38 #define i0                      0
  39 #define i1                      1
  40 #define i2                      2
  41 #define i3                      3
  42 
  43 static const long double zero = 0.0L, one = 1.0L, two = 2.0L;
  44 extern const long double _TBL_logl_hi[], _TBL_logl_lo[];
  45 static const long double two113 = 10384593717069655257060992658440192.0L,
  46         ln2hi = 6.931471805599453094172319547495844850203e-0001L,
  47         ln2lo = 1.667085920830552208890449330400379754169e-0025L,
  48         A2 = 6.666666666666666666666666666666091393804e-0001L,
  49         A3 = 4.000000000000000000000000407167070220671e-0001L,
  50         A4 = 2.857142857142857142730077490612903681164e-0001L,
  51         A5 = 2.222222222222242577702836920812882605099e-0001L,
  52         A6 = 1.818181816435493395985912667105885828356e-0001L,
  53         A7 = 1.538537835211839751112067512805496931725e-0001L,
  54         B1 = 6.666666666666666666666666666666666667787e-0001L,
  55         B2 = 3.999999999999999999999999999999848524411e-0001L,
  56         B3 = 2.857142857142857142857142865084581075070e-0001L,
  57         B4 = 2.222222222222222222222010781800643808497e-0001L,
  58         B5 = 1.818181818181818185051442171337036403674e-0001L,
  59         B6 = 1.538461538461508363540720286292008207673e-0001L,
  60         B7 = 1.333333333506731842033180638329317108428e-0001L,
  61         B8 = 1.176469984587418890634302788283946761670e-0001L,
  62         B9 = 1.053794891561452331722969901564862497132e-0001L;
  63 
  64 static long double
  65 logl_x(long double x, long double *w)
  66 {
  67         long double f, f1, v, s, z, qn, h, t;
  68         int *px = (int *)&x;
  69         int *pz = (int *)&z;
  70         int i, j, ix, n;
  71 
  72         n = 0;
  73         ix = px[i0];
  74 
  75         if (ix > 0x3ffef03f && ix < 0x3fff0820) { /* 65/63 > x > 63/65 */
  76                 f = x - one;
  77                 z = f * f;
  78 
  79                 if (((ix - 0x3fff0000) | px[i1] | px[i2] | px[i3]) == 0) {
  80                         *w = zero;
  81                         return (zero);  /* log(1)= +0 */
  82                 }
  83 
  84                 qn = one / (two + f);
  85                 s = f * qn;             /* |s|<2**-6 */
  86                 v = s * s;
  87                 h = (long double)(2.0 * (double)s);
  88                 f1 = (long double)((double)f);
  89                 t = ((two * (f - h) - h * f1) - h * (f - f1)) * qn + s * (v *
  90                     (B1 + v * (B2 + v * (B3 + v * (B4 + v * (B5 + v * (B6 + v *
  91                     (B7 + v * (B8 + v * B9)))))))));
  92                 s = (long double)((double)(h + t));
  93                 *w = t - (s - h);
  94                 return (s);
  95         }
  96 
  97         if (ix < 0x00010000) {               /* subnormal x */
  98                 x *= two113;
  99                 n = -113;
 100                 ix = px[i0];
 101         }
 102 
 103         /* LARGE_N */
 104         n += ((ix + 0x200) >> 16) - 0x3fff;
 105         ix = (ix & 0x0000ffff) | 0x3fff0000;        /* scale x to [1,2] */
 106         px[i0] = ix;
 107         i = ix + 0x200;
 108         pz[i0] = i & 0xfffffc00;
 109         pz[i1] = pz[i2] = pz[i3] = 0;
 110         qn = one / (x + z);
 111         f = x - z;
 112         s = f * qn;
 113         f1 = (long double)((double)f);
 114         h = (long double)(2.0 * (double)s);
 115         t = qn * ((two * (f - z * h) - h * f1) - h * (f - f1));
 116         j = (i >> 10) & 0x3f;
 117         v = s * s;
 118         qn = (long double)n;
 119         t += qn * ln2lo + _TBL_logl_lo[j];
 120         t += s * (v * (A2 + v * (A3 + v * (A4 + v * (A5 + v * (A6 + v *
 121             A7))))));
 122         v = qn * ln2hi + _TBL_logl_hi[j];
 123         s = h + v;
 124         t += (h - (s - v));
 125         z = (long double)((double)(s + t));
 126         *w = t - (z - s);
 127         return (z);
 128 }
 129 
 130 extern const long double _TBL_expl_hi[], _TBL_expl_lo[];
 131 static const long double
 132         invln2_32 = 4.616624130844682903551758979206054839765e+1L,
 133         ln2_32hi = 2.166084939249829091928849858592451515688e-2L,
 134         ln2_32lo = 5.209643502595475652782654157501186731779e-27L,
 135         ln2_64 = 1.083042469624914545964425189778400898568e-2L;
 136 
 137 long double
 138 powl(long double x, long double y)
 139 {
 140         long double z, ax;
 141         long double y1, y2, w1, w2;
 142         int sbx, sby, j, k, yisint, m;
 143         int hx, lx, hy, ly, ahx, ahy;
 144         int *pz = (int *)&z;
 145         int *px = (int *)&x;
 146         int *py = (int *)&y;
 147 
 148         hx = px[i0];
 149         lx = px[i1] | px[i2] | px[i3];
 150         hy = py[i0];
 151         ly = py[i1] | py[i2] | py[i3];
 152         ahx = hx & ~0x80000000;
 153         ahy = hy & ~0x80000000;
 154 
 155         if ((ahy | ly) == 0)
 156                 return (one);           /* x**+-0 = 1 */
 157         else if (hx == 0x3fff0000 && lx == 0 && (__xpg6 & _C99SUSv3_pow) != 0)
 158                 return (one);           /* C99: 1**anything = 1 */
 159         else if (ahx > 0x7fff0000 || (ahx == 0x7fff0000 && lx != 0) || ahy >
 160             0x7fff0000 || (ahy == 0x7fff0000 && ly != 0))
 161                 return (x + y);         /* +-NaN return x+y */
 162 
 163         /* includes Sun: 1**NaN = NaN */
 164         sbx = (unsigned)hx >> 31;
 165         sby = (unsigned)hy >> 31;
 166         ax = fabsl(x);
 167 
 168         /*
 169          * determine if y is an odd int when x < 0
 170          * yisint = 0 ... y is not an integer
 171          * yisint = 1 ... y is an odd int
 172          * yisint = 2 ... y is an even int
 173          */
 174         yisint = 0;
 175 
 176         if (sbx) {
 177                 if (ahy >= 0x40700000) { /* if |y|>=2**113 */
 178                         yisint = 2;                     /* even integer y */
 179                 } else if (ahy >= 0x3fff0000) {
 180                         k = (ahy >> 16) - 0x3fff; /* exponent */
 181 
 182                         if (k > 80) {
 183                                 j = ((unsigned)py[i3]) >> (112 - k);
 184 
 185                                 if ((j << (112 - k)) == py[i3])
 186                                         yisint = 2 - (j & 1);
 187                         } else if (k > 48) {
 188                                 j = ((unsigned)py[i2]) >> (80 - k);
 189 
 190                                 if ((j << (80 - k)) == py[i2])
 191                                         yisint = 2 - (j & 1);
 192                         } else if (k > 16) {
 193                                 j = ((unsigned)py[i1]) >> (48 - k);
 194 
 195                                 if ((j << (48 - k)) == py[i1])
 196                                         yisint = 2 - (j & 1);
 197                         } else if (ly == 0) {
 198                                 j = ahy >> (16 - k);
 199 
 200                                 if ((j << (16 - k)) == ahy)
 201                                         yisint = 2 - (j & 1);
 202                         }
 203                 }
 204         }
 205 
 206         /* special value of y */
 207         if (ly == 0) {
 208                 if (ahy == 0x7fff0000) {        /* y is +-inf */
 209                         if (((ahx - 0x3fff0000) | lx) == 0) {
 210                                 if ((__xpg6 & _C99SUSv3_pow) != 0)
 211                                         return (one);
 212                                 /* C99: (-1)**+-inf = 1 */
 213                                 else
 214                                         return (y - y);
 215 
 216                                 /* Sun: (+-1)**+-inf = NaN */
 217                         } else if (ahx >= 0x3fff0000) {
 218                                 /* (|x|>1)**+,-inf = inf,0 */
 219                                 return (sby == 0 ? y : zero);
 220                         } else {                /* (|x|<1)**-,+inf = inf,0 */
 221                                 return (sby != 0 ? -y : zero);
 222                         }
 223                 } else if (ahy == 0x3fff0000) { /* y is +-1 */
 224                         if (sby != 0)
 225                                 return (one / x);
 226                         else
 227                                 return (x);
 228                 } else if (hy == 0x40000000) {  /* y is 2 */
 229                         return (x * x);
 230                 } else if (hy == 0x3ffe0000) {  /* y is 0.5 */
 231                         if (!((ahx | lx) == 0 || ((ahx - 0x7fff0000) | lx) ==
 232                             0))
 233                                 return (sqrtl(x));
 234                 }
 235         }
 236 
 237         /* special value of x */
 238         if (lx == 0) {
 239                 if (ahx == 0x7fff0000 || ahx == 0 || ahx == 0x3fff0000) {
 240                         /* x is +-0,+-inf,+-1 */
 241                         z = ax;
 242 
 243                         if (sby == 1)
 244                                 z = one / z;    /* z = 1/|x| if y is negative */
 245 
 246                         if (sbx == 1) {
 247                                 if (ahx == 0x3fff0000 && yisint == 0)
 248                                         z = zero / zero;
 249                                 /* (-1)**non-int is NaN */
 250                                 else if (yisint == 1)
 251                                         z = -z; /* (x<0)**odd = -(|x|**odd) */
 252                         }
 253 
 254                         return (z);
 255                 }
 256         }
 257 
 258         /* (x<0)**(non-int) is NaN */
 259         if (sbx == 1 && yisint == 0)
 260                 return (zero / zero);   /* should be volatile */
 261 
 262         /*
 263          * Now ax is finite, y is finite
 264          * first compute log(ax) = w1+w2, with 53 bits w1
 265          */
 266         w1 = logl_x(ax, &w2);
 267 
 268         /* split up y into y1+y2 and compute (y1+y2)*(w1+w2) */
 269         if (ly == 0 || ahy >= 0x43fe0000) {
 270                 y1 = y * w1;
 271                 y2 = y * w2;
 272         } else {
 273                 y1 = (long double)((double)y);
 274                 y2 = (y - y1) * w1 + y * w2;
 275                 y1 *= w1;
 276         }
 277 
 278         z = y1 + y2;
 279         j = pz[i0];
 280 
 281         if ((unsigned)j >= 0xffff0000) {     /* NaN or -inf */
 282                 if (sbx == 1 && yisint == 1)
 283                         return (one / z);
 284                 else
 285                         return (-one / z);
 286         } else if ((j & ~0x80000000) < 0x3fc30000) {     /* |x|<2^-60 */
 287                 if (sbx == 1 && yisint == 1)
 288                         return (-one - z);
 289                 else
 290                         return (one + z);
 291         } else if (j > 0) {
 292                 if (j > 0x400d0000) {
 293                         if (sbx == 1 && yisint == 1)
 294                                 return (scalbnl(-one, 20000));
 295                         else
 296                                 return (scalbnl(one, 20000));
 297                 }
 298 
 299                 k = (int)(invln2_32 * (z + ln2_64));
 300         } else {
 301                 if ((unsigned)j > 0xc00d0000) {
 302                         if (sbx == 1 && yisint == 1)
 303                                 return (scalbnl(-one, -20000));
 304                         else
 305                                 return (scalbnl(one, -20000));
 306                 }
 307 
 308                 k = (int)(invln2_32 * (z - ln2_64));
 309         }
 310 
 311         j = k & 0x1f;
 312         m = k >> 5;
 313         {
 314                 /* rational approximation coeffs for [-(ln2)/64,(ln2)/64] */
 315                 long double t1 = 1.666666666666666666666666666660876387437e-1L,
 316                     t2 = -2.777777777777777777777707812093173478756e-3L,
 317                     t3 = 6.613756613756613482074280932874221202424e-5L,
 318                     t4 = -1.653439153392139954169609822742235851120e-6L,
 319                     t5 = 4.175314851769539751387852116610973796053e-8L;
 320                 long double t = (long double)k;
 321 
 322                 w1 = (y2 - (t * ln2_32hi - y1)) - t * ln2_32lo;
 323                 t = w1 * w1;
 324                 w2 = (w1 - t * (t1 + t * (t2 + t * (t3 + t * (t4 + t * t5))))) -
 325                     two;
 326                 z = _TBL_expl_hi[j] - ((_TBL_expl_hi[j] * (w1 + w1)) / w2 -
 327                     _TBL_expl_lo[j]);
 328         }
 329         j = m + (pz[i0] >> 16);
 330 
 331         if (j && (unsigned)j < 0x7fff)
 332                 pz[i0] += m << 16;
 333         else
 334                 z = scalbnl(z, m);
 335 
 336         if (sbx == 1 && yisint == 1)
 337                 z = -z;                 /* (-ve)**(odd int) */
 338 
 339         return (z);
 340 }
 341 #else
 342 #error Unsupported Architecture
 343 #endif /* defined(__sparc) */