125c28e83SPiotr Jasiukajtis /*
225c28e83SPiotr Jasiukajtis * CDDL HEADER START
325c28e83SPiotr Jasiukajtis *
425c28e83SPiotr Jasiukajtis * The contents of this file are subject to the terms of the
525c28e83SPiotr Jasiukajtis * Common Development and Distribution License (the "License").
625c28e83SPiotr Jasiukajtis * You may not use this file except in compliance with the License.
725c28e83SPiotr Jasiukajtis *
825c28e83SPiotr Jasiukajtis * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
925c28e83SPiotr Jasiukajtis * or http://www.opensolaris.org/os/licensing.
1025c28e83SPiotr Jasiukajtis * See the License for the specific language governing permissions
1125c28e83SPiotr Jasiukajtis * and limitations under the License.
1225c28e83SPiotr Jasiukajtis *
1325c28e83SPiotr Jasiukajtis * When distributing Covered Code, include this CDDL HEADER in each
1425c28e83SPiotr Jasiukajtis * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
1525c28e83SPiotr Jasiukajtis * If applicable, add the following below this CDDL HEADER, with the
1625c28e83SPiotr Jasiukajtis * fields enclosed by brackets "[]" replaced with your own identifying
1725c28e83SPiotr Jasiukajtis * information: Portions Copyright [yyyy] [name of copyright owner]
1825c28e83SPiotr Jasiukajtis *
1925c28e83SPiotr Jasiukajtis * CDDL HEADER END
2025c28e83SPiotr Jasiukajtis */
2125c28e83SPiotr Jasiukajtis
2225c28e83SPiotr Jasiukajtis /*
2325c28e83SPiotr Jasiukajtis * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
2425c28e83SPiotr Jasiukajtis */
2525c28e83SPiotr Jasiukajtis /*
2625c28e83SPiotr Jasiukajtis * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
2725c28e83SPiotr Jasiukajtis * Use is subject to license terms.
2825c28e83SPiotr Jasiukajtis */
2925c28e83SPiotr Jasiukajtis
30*ddc0e0b5SRichard Lowe #pragma weak __pow = pow
3125c28e83SPiotr Jasiukajtis
3225c28e83SPiotr Jasiukajtis /*
3325c28e83SPiotr Jasiukajtis * pow(x,y) return x**y
3425c28e83SPiotr Jasiukajtis * n
3525c28e83SPiotr Jasiukajtis * Method: Let x = 2 * (1+f)
3625c28e83SPiotr Jasiukajtis * 1. Compute and return log2(x) in two pieces:
3725c28e83SPiotr Jasiukajtis * log2(x) = w1 + w2,
3825c28e83SPiotr Jasiukajtis * where w1 has 24 bits trailing zero.
3925c28e83SPiotr Jasiukajtis * 2. Perform y*log2(x) by simulating muti-precision arithmetic
4025c28e83SPiotr Jasiukajtis * 3. Return x**y = exp2(y*log(x))
4125c28e83SPiotr Jasiukajtis *
4225c28e83SPiotr Jasiukajtis * Special cases:
4325c28e83SPiotr Jasiukajtis * 1. (anything) ** +-0 is 1
4425c28e83SPiotr Jasiukajtis * 1'. 1 ** (anything) is 1 (C99; 1 ** +-INF/NAN used to be NAN)
4525c28e83SPiotr Jasiukajtis * 2. (anything) ** 1 is itself
4625c28e83SPiotr Jasiukajtis * 3. (anything except 1) ** NAN is NAN ("except 1" is C99)
4725c28e83SPiotr Jasiukajtis * 4. NAN ** (anything except 0) is NAN
4825c28e83SPiotr Jasiukajtis * 5. +-(|x| > 1) ** +INF is +INF
4925c28e83SPiotr Jasiukajtis * 6. +-(|x| > 1) ** -INF is +0
5025c28e83SPiotr Jasiukajtis * 7. +-(|x| < 1) ** +INF is +0
5125c28e83SPiotr Jasiukajtis * 8. +-(|x| < 1) ** -INF is +INF
5225c28e83SPiotr Jasiukajtis * 9. -1 ** +-INF is 1 (C99; -1 ** +-INF used to be NAN)
5325c28e83SPiotr Jasiukajtis * 10. +0 ** (+anything except 0, NAN) is +0
5425c28e83SPiotr Jasiukajtis * 11. -0 ** (+anything except 0, NAN, odd integer) is +0
5525c28e83SPiotr Jasiukajtis * 12. +0 ** (-anything except 0, NAN) is +INF
5625c28e83SPiotr Jasiukajtis * 13. -0 ** (-anything except 0, NAN, odd integer) is +INF
5725c28e83SPiotr Jasiukajtis * 14. -0 ** (odd integer) = -( +0 ** (odd integer) )
5825c28e83SPiotr Jasiukajtis * 15. +INF ** (+anything except 0,NAN) is +INF
5925c28e83SPiotr Jasiukajtis * 16. +INF ** (-anything except 0,NAN) is +0
6025c28e83SPiotr Jasiukajtis * 17. -INF ** (anything) = -0 ** (-anything)
6125c28e83SPiotr Jasiukajtis * 18. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer)
6225c28e83SPiotr Jasiukajtis * 19. (-anything except 0 and inf) ** (non-integer) is NAN
6325c28e83SPiotr Jasiukajtis *
6425c28e83SPiotr Jasiukajtis * Accuracy:
6525c28e83SPiotr Jasiukajtis * pow(x,y) returns x**y nearly rounded. In particular
6625c28e83SPiotr Jasiukajtis * pow(integer,integer)
6725c28e83SPiotr Jasiukajtis * always returns the correct integer provided it is representable.
6825c28e83SPiotr Jasiukajtis */
6925c28e83SPiotr Jasiukajtis
7025c28e83SPiotr Jasiukajtis #include "libm.h"
7125c28e83SPiotr Jasiukajtis #include "xpg6.h" /* __xpg6 */
7225c28e83SPiotr Jasiukajtis #define _C99SUSv3_pow _C99SUSv3_pow_treats_Inf_as_an_even_int
7325c28e83SPiotr Jasiukajtis
7425c28e83SPiotr Jasiukajtis static const double zero = 0.0, one = 1.0, two = 2.0;
7525c28e83SPiotr Jasiukajtis
7625c28e83SPiotr Jasiukajtis extern const double _TBL_log2_hi[], _TBL_log2_lo[];
7725c28e83SPiotr Jasiukajtis static const double
7825c28e83SPiotr Jasiukajtis two53 = 9007199254740992.0,
7925c28e83SPiotr Jasiukajtis A1_hi = 2.8853900432586669921875,
8025c28e83SPiotr Jasiukajtis A1_lo = 3.8519259825035041963606002e-8,
8125c28e83SPiotr Jasiukajtis A1 = 2.885390081777926817222541963606002026086e+0000,
8225c28e83SPiotr Jasiukajtis A2 = 9.617966939207270828380543979852286255862e-0001,
8325c28e83SPiotr Jasiukajtis A3 = 5.770807680887875964868853124873696201995e-0001,
8425c28e83SPiotr Jasiukajtis B0_hi = 2.8853900432586669921875,
8525c28e83SPiotr Jasiukajtis B0_lo = 3.8519259822532793056374320585e-8,
8625c28e83SPiotr Jasiukajtis B0 = 2.885390081777926814720293056374320585689e+0000,
8725c28e83SPiotr Jasiukajtis B1 = 9.617966939259755138949202350396200257632e-0001,
8825c28e83SPiotr Jasiukajtis B2 = 5.770780163585687000782112776448797953382e-0001,
8925c28e83SPiotr Jasiukajtis B3 = 4.121985488948771523290174512461778354953e-0001,
9025c28e83SPiotr Jasiukajtis B4 = 3.207590534812432970433641789022666850193e-0001;
9125c28e83SPiotr Jasiukajtis
9225c28e83SPiotr Jasiukajtis static double
log2_x(double x,double * w)9325c28e83SPiotr Jasiukajtis log2_x(double x, double *w) {
9425c28e83SPiotr Jasiukajtis double f, s, z, qn, h, t;
9525c28e83SPiotr Jasiukajtis int *px = (int *) &x;
9625c28e83SPiotr Jasiukajtis int *pz = (int *) &z;
9725c28e83SPiotr Jasiukajtis int i, j, ix, n;
9825c28e83SPiotr Jasiukajtis
9925c28e83SPiotr Jasiukajtis n = 0;
10025c28e83SPiotr Jasiukajtis ix = px[HIWORD];
10125c28e83SPiotr Jasiukajtis if (ix >= 0x3fef03f1 && ix < 0x3ff08208) { /* 65/63 > x > 63/65 */
10225c28e83SPiotr Jasiukajtis double f1, v;
10325c28e83SPiotr Jasiukajtis f = x - one;
10425c28e83SPiotr Jasiukajtis if (((ix - 0x3ff00000) | px[LOWORD]) == 0) {
10525c28e83SPiotr Jasiukajtis *w = zero;
10625c28e83SPiotr Jasiukajtis return (zero); /* log2(1)= +0 */
10725c28e83SPiotr Jasiukajtis }
10825c28e83SPiotr Jasiukajtis qn = one / (two + f);
10925c28e83SPiotr Jasiukajtis s = f * qn; /* |s|<2**-6 */
11025c28e83SPiotr Jasiukajtis v = s * s;
11125c28e83SPiotr Jasiukajtis h = (double) ((float) s);
11225c28e83SPiotr Jasiukajtis f1 = (double) ((float) f);
11325c28e83SPiotr Jasiukajtis t = qn * (((f - two * h) - h * f1) - h * (f - f1));
11425c28e83SPiotr Jasiukajtis /* s = h+t */
11525c28e83SPiotr Jasiukajtis f1 = h * B0_lo + s * (v * (B1 + v * (B2 + v * (B3 + v * B4))));
11625c28e83SPiotr Jasiukajtis t = f1 + t * B0;
11725c28e83SPiotr Jasiukajtis h *= B0_hi;
11825c28e83SPiotr Jasiukajtis s = (double) ((float) (h + t));
11925c28e83SPiotr Jasiukajtis *w = t - (s - h);
12025c28e83SPiotr Jasiukajtis return (s);
12125c28e83SPiotr Jasiukajtis }
12225c28e83SPiotr Jasiukajtis if (ix < 0x00100000) { /* subnormal x */
12325c28e83SPiotr Jasiukajtis x *= two53;
12425c28e83SPiotr Jasiukajtis n = -53;
12525c28e83SPiotr Jasiukajtis ix = px[HIWORD];
12625c28e83SPiotr Jasiukajtis }
12725c28e83SPiotr Jasiukajtis /* LARGE N */
12825c28e83SPiotr Jasiukajtis n += ((ix + 0x1000) >> 20) - 0x3ff;
12925c28e83SPiotr Jasiukajtis ix = (ix & 0x000fffff) | 0x3ff00000; /* scale x to [1,2] */
13025c28e83SPiotr Jasiukajtis px[HIWORD] = ix;
13125c28e83SPiotr Jasiukajtis i = ix + 0x1000;
13225c28e83SPiotr Jasiukajtis pz[HIWORD] = i & 0xffffe000;
13325c28e83SPiotr Jasiukajtis pz[LOWORD] = 0;
13425c28e83SPiotr Jasiukajtis qn = one / (x + z);
13525c28e83SPiotr Jasiukajtis f = x - z;
13625c28e83SPiotr Jasiukajtis s = f * qn;
13725c28e83SPiotr Jasiukajtis h = (double) ((float) s);
13825c28e83SPiotr Jasiukajtis t = qn * ((f - (h + h) * z) - h * f);
13925c28e83SPiotr Jasiukajtis j = (i >> 13) & 0x7f;
14025c28e83SPiotr Jasiukajtis f = s * s;
14125c28e83SPiotr Jasiukajtis t = t * A1 + h * A1_lo;
14225c28e83SPiotr Jasiukajtis t += (s * f) * (A2 + f * A3);
14325c28e83SPiotr Jasiukajtis qn = h * A1_hi;
14425c28e83SPiotr Jasiukajtis s = n + _TBL_log2_hi[j];
14525c28e83SPiotr Jasiukajtis h = qn + s;
14625c28e83SPiotr Jasiukajtis t += _TBL_log2_lo[j] - ((h - s) - qn);
14725c28e83SPiotr Jasiukajtis f = (double) ((float) (h + t));
14825c28e83SPiotr Jasiukajtis *w = t - (f - h);
14925c28e83SPiotr Jasiukajtis return (f);
15025c28e83SPiotr Jasiukajtis }
15125c28e83SPiotr Jasiukajtis
15225c28e83SPiotr Jasiukajtis extern const double _TBL_exp2_hi[], _TBL_exp2_lo[];
15325c28e83SPiotr Jasiukajtis static const double /* poly app of 2^x-1 on [-1e-10,2^-7+1e-10] */
15425c28e83SPiotr Jasiukajtis E1 = 6.931471805599453100674958533810346197328e-0001,
15525c28e83SPiotr Jasiukajtis E2 = 2.402265069587779347846769151717493815979e-0001,
15625c28e83SPiotr Jasiukajtis E3 = 5.550410866475410512631124892773937864699e-0002,
15725c28e83SPiotr Jasiukajtis E4 = 9.618143209991026824853712740162451423355e-0003,
15825c28e83SPiotr Jasiukajtis E5 = 1.333357676549940345096774122231849082991e-0003;
15925c28e83SPiotr Jasiukajtis
16025c28e83SPiotr Jasiukajtis double
pow(double x,double y)16125c28e83SPiotr Jasiukajtis pow(double x, double y) {
16225c28e83SPiotr Jasiukajtis double z, ax;
16325c28e83SPiotr Jasiukajtis double y1, y2, w1, w2;
16425c28e83SPiotr Jasiukajtis int sbx, sby, j, k, yisint;
16525c28e83SPiotr Jasiukajtis int hx, hy, ahx, ahy;
16625c28e83SPiotr Jasiukajtis unsigned lx, ly;
16725c28e83SPiotr Jasiukajtis int *pz = (int *) &z;
16825c28e83SPiotr Jasiukajtis
16925c28e83SPiotr Jasiukajtis hx = ((int *) &x)[HIWORD];
17025c28e83SPiotr Jasiukajtis lx = ((unsigned *) &x)[LOWORD];
17125c28e83SPiotr Jasiukajtis hy = ((int *) &y)[HIWORD];
17225c28e83SPiotr Jasiukajtis ly = ((unsigned *) &y)[LOWORD];
17325c28e83SPiotr Jasiukajtis ahx = hx & ~0x80000000;
17425c28e83SPiotr Jasiukajtis ahy = hy & ~0x80000000;
17525c28e83SPiotr Jasiukajtis if ((ahy | ly) == 0) { /* y==zero */
17625c28e83SPiotr Jasiukajtis if ((ahx | lx) == 0)
17725c28e83SPiotr Jasiukajtis z = _SVID_libm_err(x, y, 20); /* +-0**+-0 */
17825c28e83SPiotr Jasiukajtis else if ((ahx | (((lx | -lx) >> 31) & 1)) > 0x7ff00000)
17925c28e83SPiotr Jasiukajtis z = _SVID_libm_err(x, y, 42); /* NaN**+-0 */
18025c28e83SPiotr Jasiukajtis else
18125c28e83SPiotr Jasiukajtis z = one; /* x**+-0 = 1 */
18225c28e83SPiotr Jasiukajtis return (z);
18325c28e83SPiotr Jasiukajtis } else if (hx == 0x3ff00000 && lx == 0 &&
18425c28e83SPiotr Jasiukajtis (__xpg6 & _C99SUSv3_pow) != 0)
18525c28e83SPiotr Jasiukajtis return (one); /* C99: 1**anything = 1 */
18625c28e83SPiotr Jasiukajtis else if (ahx > 0x7ff00000 || (ahx == 0x7ff00000 && lx != 0) ||
18725c28e83SPiotr Jasiukajtis ahy > 0x7ff00000 || (ahy == 0x7ff00000 && ly != 0))
18825c28e83SPiotr Jasiukajtis return (x * y); /* +-NaN return x*y; + -> * for Cheetah */
18925c28e83SPiotr Jasiukajtis /* includes Sun: 1**NaN = NaN */
19025c28e83SPiotr Jasiukajtis sbx = (unsigned) hx >> 31;
19125c28e83SPiotr Jasiukajtis sby = (unsigned) hy >> 31;
19225c28e83SPiotr Jasiukajtis ax = fabs(x);
19325c28e83SPiotr Jasiukajtis
19425c28e83SPiotr Jasiukajtis /*
19525c28e83SPiotr Jasiukajtis * determine if y is an odd int when x < 0
19625c28e83SPiotr Jasiukajtis * yisint = 0 ... y is not an integer
19725c28e83SPiotr Jasiukajtis * yisint = 1 ... y is an odd int
19825c28e83SPiotr Jasiukajtis * yisint = 2 ... y is an even int
19925c28e83SPiotr Jasiukajtis */
20025c28e83SPiotr Jasiukajtis yisint = 0;
20125c28e83SPiotr Jasiukajtis if (sbx) {
20225c28e83SPiotr Jasiukajtis if (ahy >= 0x43400000)
20325c28e83SPiotr Jasiukajtis yisint = 2; /* even integer y */
20425c28e83SPiotr Jasiukajtis else if (ahy >= 0x3ff00000) {
20525c28e83SPiotr Jasiukajtis k = (ahy >> 20) - 0x3ff; /* exponent */
20625c28e83SPiotr Jasiukajtis if (k > 20) {
20725c28e83SPiotr Jasiukajtis j = ly >> (52 - k);
20825c28e83SPiotr Jasiukajtis if ((j << (52 - k)) == ly)
20925c28e83SPiotr Jasiukajtis yisint = 2 - (j & 1);
21025c28e83SPiotr Jasiukajtis } else if (ly == 0) {
21125c28e83SPiotr Jasiukajtis j = ahy >> (20 - k);
21225c28e83SPiotr Jasiukajtis if ((j << (20 - k)) == ahy)
21325c28e83SPiotr Jasiukajtis yisint = 2 - (j & 1);
21425c28e83SPiotr Jasiukajtis }
21525c28e83SPiotr Jasiukajtis }
21625c28e83SPiotr Jasiukajtis }
21725c28e83SPiotr Jasiukajtis /* special value of y */
21825c28e83SPiotr Jasiukajtis if (ly == 0) {
21925c28e83SPiotr Jasiukajtis if (ahy == 0x7ff00000) { /* y is +-inf */
22025c28e83SPiotr Jasiukajtis if (((ahx - 0x3ff00000) | lx) == 0) {
22125c28e83SPiotr Jasiukajtis if ((__xpg6 & _C99SUSv3_pow) != 0)
22225c28e83SPiotr Jasiukajtis return (one);
22325c28e83SPiotr Jasiukajtis /* C99: (-1)**+-inf = 1 */
22425c28e83SPiotr Jasiukajtis else
22525c28e83SPiotr Jasiukajtis return (y - y);
22625c28e83SPiotr Jasiukajtis /* Sun: (+-1)**+-inf = NaN */
22725c28e83SPiotr Jasiukajtis } else if (ahx >= 0x3ff00000)
22825c28e83SPiotr Jasiukajtis /* (|x|>1)**+,-inf = inf,0 */
22925c28e83SPiotr Jasiukajtis return (sby == 0 ? y : zero);
23025c28e83SPiotr Jasiukajtis else /* (|x|<1)**-,+inf = inf,0 */
23125c28e83SPiotr Jasiukajtis return (sby != 0 ? -y : zero);
23225c28e83SPiotr Jasiukajtis }
23325c28e83SPiotr Jasiukajtis if (ahy == 0x3ff00000) { /* y is +-1 */
23425c28e83SPiotr Jasiukajtis if (sby != 0) { /* y is -1 */
23525c28e83SPiotr Jasiukajtis if (x == zero) /* divided by zero */
23625c28e83SPiotr Jasiukajtis return (_SVID_libm_err(x, y, 23));
23725c28e83SPiotr Jasiukajtis else if (ahx < 0x40000 || ((ahx - 0x40000) |
23825c28e83SPiotr Jasiukajtis lx) == 0) /* overflow */
23925c28e83SPiotr Jasiukajtis return (_SVID_libm_err(x, y, 21));
24025c28e83SPiotr Jasiukajtis else
24125c28e83SPiotr Jasiukajtis return (one / x);
24225c28e83SPiotr Jasiukajtis } else
24325c28e83SPiotr Jasiukajtis return (x);
24425c28e83SPiotr Jasiukajtis }
24525c28e83SPiotr Jasiukajtis if (hy == 0x40000000) { /* y is 2 */
24625c28e83SPiotr Jasiukajtis if (ahx >= 0x5ff00000 && ahx < 0x7ff00000)
24725c28e83SPiotr Jasiukajtis return (_SVID_libm_err(x, y, 21));
24825c28e83SPiotr Jasiukajtis /* x*x overflow */
24925c28e83SPiotr Jasiukajtis else if ((ahx < 0x1e56a09e && (ahx | lx) != 0) ||
25025c28e83SPiotr Jasiukajtis (ahx == 0x1e56a09e && lx < 0x667f3bcd))
25125c28e83SPiotr Jasiukajtis return (_SVID_libm_err(x, y, 22));
25225c28e83SPiotr Jasiukajtis /* x*x underflow */
25325c28e83SPiotr Jasiukajtis else
25425c28e83SPiotr Jasiukajtis return (x * x);
25525c28e83SPiotr Jasiukajtis }
25625c28e83SPiotr Jasiukajtis if (hy == 0x3fe00000) {
25725c28e83SPiotr Jasiukajtis if (!((ahx | lx) == 0 || ((ahx - 0x7ff00000) | lx) ==
25825c28e83SPiotr Jasiukajtis 0 || sbx == 1))
25925c28e83SPiotr Jasiukajtis return (sqrt(x)); /* y is 0.5 and x > 0 */
26025c28e83SPiotr Jasiukajtis }
26125c28e83SPiotr Jasiukajtis }
26225c28e83SPiotr Jasiukajtis /* special value of x */
26325c28e83SPiotr Jasiukajtis if (lx == 0) {
26425c28e83SPiotr Jasiukajtis if (ahx == 0x7ff00000 || ahx == 0 || ahx == 0x3ff00000) {
26525c28e83SPiotr Jasiukajtis /* x is +-0,+-inf,-1 */
26625c28e83SPiotr Jasiukajtis z = ax;
26725c28e83SPiotr Jasiukajtis if (sby == 1) {
26825c28e83SPiotr Jasiukajtis z = one / z; /* z = |x|**y */
26925c28e83SPiotr Jasiukajtis if (ahx == 0)
27025c28e83SPiotr Jasiukajtis return (_SVID_libm_err(x, y, 23));
27125c28e83SPiotr Jasiukajtis }
27225c28e83SPiotr Jasiukajtis if (sbx == 1) {
27325c28e83SPiotr Jasiukajtis if (ahx == 0x3ff00000 && yisint == 0)
27425c28e83SPiotr Jasiukajtis z = _SVID_libm_err(x, y, 24);
27525c28e83SPiotr Jasiukajtis /* neg**non-integral is NaN + invalid */
27625c28e83SPiotr Jasiukajtis else if (yisint == 1)
27725c28e83SPiotr Jasiukajtis z = -z; /* (x<0)**odd = -(|x|**odd) */
27825c28e83SPiotr Jasiukajtis }
27925c28e83SPiotr Jasiukajtis return (z);
28025c28e83SPiotr Jasiukajtis }
28125c28e83SPiotr Jasiukajtis }
28225c28e83SPiotr Jasiukajtis /* (x<0)**(non-int) is NaN */
28325c28e83SPiotr Jasiukajtis if (sbx == 1 && yisint == 0)
28425c28e83SPiotr Jasiukajtis return (_SVID_libm_err(x, y, 24));
28525c28e83SPiotr Jasiukajtis /* Now ax is finite, y is finite */
28625c28e83SPiotr Jasiukajtis /* first compute log2(ax) = w1+w2, with 24 bits w1 */
28725c28e83SPiotr Jasiukajtis w1 = log2_x(ax, &w2);
28825c28e83SPiotr Jasiukajtis
28925c28e83SPiotr Jasiukajtis /* split up y into y1+y2 and compute (y1+y2)*(w1+w2) */
29025c28e83SPiotr Jasiukajtis if (((ly & 0x07ffffff) == 0) || ahy >= 0x47e00000 ||
29125c28e83SPiotr Jasiukajtis ahy <= 0x38100000) {
29225c28e83SPiotr Jasiukajtis /* no need to split if y is short or too large or too small */
29325c28e83SPiotr Jasiukajtis y1 = y * w1;
29425c28e83SPiotr Jasiukajtis y2 = y * w2;
29525c28e83SPiotr Jasiukajtis } else {
29625c28e83SPiotr Jasiukajtis y1 = (double) ((float) y);
29725c28e83SPiotr Jasiukajtis y2 = (y - y1) * w1 + y * w2;
29825c28e83SPiotr Jasiukajtis y1 *= w1;
29925c28e83SPiotr Jasiukajtis }
30025c28e83SPiotr Jasiukajtis z = y1 + y2;
30125c28e83SPiotr Jasiukajtis j = pz[HIWORD];
30225c28e83SPiotr Jasiukajtis if (j >= 0x40900000) { /* z >= 1024 */
30325c28e83SPiotr Jasiukajtis if (!(j == 0x40900000 && pz[LOWORD] == 0)) /* z > 1024 */
30425c28e83SPiotr Jasiukajtis return (_SVID_libm_err(x, y, 21)); /* overflow */
30525c28e83SPiotr Jasiukajtis else {
30625c28e83SPiotr Jasiukajtis w2 = y1 - z;
30725c28e83SPiotr Jasiukajtis w2 += y2;
30825c28e83SPiotr Jasiukajtis /* rounded to inf */
30925c28e83SPiotr Jasiukajtis if (w2 >= -8.008566259537296567160e-17)
31025c28e83SPiotr Jasiukajtis return (_SVID_libm_err(x, y, 21));
31125c28e83SPiotr Jasiukajtis /* overflow */
31225c28e83SPiotr Jasiukajtis }
31325c28e83SPiotr Jasiukajtis } else if ((j & ~0x80000000) >= 0x4090cc00) { /* z <= -1075 */
31425c28e83SPiotr Jasiukajtis if (!(j == 0xc090cc00 && pz[LOWORD] == 0)) /* z < -1075 */
31525c28e83SPiotr Jasiukajtis return (_SVID_libm_err(x, y, 22)); /* underflow */
31625c28e83SPiotr Jasiukajtis else {
31725c28e83SPiotr Jasiukajtis w2 = y1 - z;
31825c28e83SPiotr Jasiukajtis w2 += y2;
31925c28e83SPiotr Jasiukajtis if (w2 <= zero) /* underflow */
32025c28e83SPiotr Jasiukajtis return (_SVID_libm_err(x, y, 22));
32125c28e83SPiotr Jasiukajtis }
32225c28e83SPiotr Jasiukajtis }
32325c28e83SPiotr Jasiukajtis /*
32425c28e83SPiotr Jasiukajtis * compute 2**(k+f[j]+g)
32525c28e83SPiotr Jasiukajtis */
32625c28e83SPiotr Jasiukajtis k = (int) (z * 64.0 + (((hy ^ (ahx - 0x3ff00000)) > 0) ? 0.5 : -0.5));
32725c28e83SPiotr Jasiukajtis j = k & 63;
32825c28e83SPiotr Jasiukajtis w1 = y2 - ((double) k * 0.015625 - y1);
32925c28e83SPiotr Jasiukajtis w2 = _TBL_exp2_hi[j];
33025c28e83SPiotr Jasiukajtis z = _TBL_exp2_lo[j] + (w2 * w1) * (E1 + w1 * (E2 + w1 * (E3 + w1 *
33125c28e83SPiotr Jasiukajtis (E4 + w1 * E5))));
33225c28e83SPiotr Jasiukajtis z += w2;
33325c28e83SPiotr Jasiukajtis k >>= 6;
33425c28e83SPiotr Jasiukajtis if (k < -1021)
33525c28e83SPiotr Jasiukajtis z = scalbn(z, k);
33625c28e83SPiotr Jasiukajtis else /* subnormal output */
33725c28e83SPiotr Jasiukajtis pz[HIWORD] += k << 20;
33825c28e83SPiotr Jasiukajtis if (sbx == 1 && yisint == 1)
33925c28e83SPiotr Jasiukajtis z = -z; /* (-ve)**(odd int) */
34025c28e83SPiotr Jasiukajtis return (z);
34125c28e83SPiotr Jasiukajtis }
342