149393c00Smartynas /*
249393c00Smartynas * ====================================================
349393c00Smartynas * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
449393c00Smartynas *
549393c00Smartynas * Developed at SunPro, a Sun Microsystems, Inc. business.
649393c00Smartynas * Permission to use, copy, modify, and distribute this
749393c00Smartynas * software is freely granted, provided that this notice
849393c00Smartynas * is preserved.
949393c00Smartynas * ====================================================
1049393c00Smartynas */
1149393c00Smartynas
1249393c00Smartynas /*
1349393c00Smartynas * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
1449393c00Smartynas *
1549393c00Smartynas * Permission to use, copy, modify, and distribute this software for any
1649393c00Smartynas * purpose with or without fee is hereby granted, provided that the above
1749393c00Smartynas * copyright notice and this permission notice appear in all copies.
1849393c00Smartynas *
1949393c00Smartynas * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
2049393c00Smartynas * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
2149393c00Smartynas * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
2249393c00Smartynas * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
2349393c00Smartynas * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
2449393c00Smartynas * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
2549393c00Smartynas * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
2649393c00Smartynas */
2749393c00Smartynas
2849393c00Smartynas /* powl(x,y) return x**y
2949393c00Smartynas *
3049393c00Smartynas * n
3149393c00Smartynas * Method: Let x = 2 * (1+f)
3249393c00Smartynas * 1. Compute and return log2(x) in two pieces:
3349393c00Smartynas * log2(x) = w1 + w2,
3449393c00Smartynas * where w1 has 113-53 = 60 bit trailing zeros.
3549393c00Smartynas * 2. Perform y*log2(x) = n+y' by simulating muti-precision
3649393c00Smartynas * arithmetic, where |y'|<=0.5.
3749393c00Smartynas * 3. Return x**y = 2**n*exp(y'*log2)
3849393c00Smartynas *
3949393c00Smartynas * Special cases:
4049393c00Smartynas * 1. (anything) ** 0 is 1
4149393c00Smartynas * 2. (anything) ** 1 is itself
4249393c00Smartynas * 3. (anything) ** NAN is NAN
4349393c00Smartynas * 4. NAN ** (anything except 0) is NAN
4449393c00Smartynas * 5. +-(|x| > 1) ** +INF is +INF
4549393c00Smartynas * 6. +-(|x| > 1) ** -INF is +0
4649393c00Smartynas * 7. +-(|x| < 1) ** +INF is +0
4749393c00Smartynas * 8. +-(|x| < 1) ** -INF is +INF
4849393c00Smartynas * 9. +-1 ** +-INF is NAN
4949393c00Smartynas * 10. +0 ** (+anything except 0, NAN) is +0
5049393c00Smartynas * 11. -0 ** (+anything except 0, NAN, odd integer) is +0
5149393c00Smartynas * 12. +0 ** (-anything except 0, NAN) is +INF
5249393c00Smartynas * 13. -0 ** (-anything except 0, NAN, odd integer) is +INF
5349393c00Smartynas * 14. -0 ** (odd integer) = -( +0 ** (odd integer) )
5449393c00Smartynas * 15. +INF ** (+anything except 0,NAN) is +INF
5549393c00Smartynas * 16. +INF ** (-anything except 0,NAN) is +0
5649393c00Smartynas * 17. -INF ** (anything) = -0 ** (-anything)
5749393c00Smartynas * 18. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer)
5849393c00Smartynas * 19. (-anything except 0 and inf) ** (non-integer) is NAN
5949393c00Smartynas *
6049393c00Smartynas */
6149393c00Smartynas
6249393c00Smartynas #include <math.h>
6349393c00Smartynas
6449393c00Smartynas #include "math_private.h"
6549393c00Smartynas
6649393c00Smartynas static const long double bp[] = {
6749393c00Smartynas 1.0L,
6849393c00Smartynas 1.5L,
6949393c00Smartynas };
7049393c00Smartynas
7149393c00Smartynas /* log_2(1.5) */
7249393c00Smartynas static const long double dp_h[] = {
7349393c00Smartynas 0.0,
7449393c00Smartynas 5.8496250072115607565592654282227158546448E-1L
7549393c00Smartynas };
7649393c00Smartynas
7749393c00Smartynas /* Low part of log_2(1.5) */
7849393c00Smartynas static const long double dp_l[] = {
7949393c00Smartynas 0.0,
8049393c00Smartynas 1.0579781240112554492329533686862998106046E-16L
8149393c00Smartynas };
8249393c00Smartynas
8349393c00Smartynas static const long double zero = 0.0L,
8449393c00Smartynas one = 1.0L,
8549393c00Smartynas two = 2.0L,
8649393c00Smartynas two113 = 1.0384593717069655257060992658440192E34L,
8749393c00Smartynas huge = 1.0e3000L,
8849393c00Smartynas tiny = 1.0e-3000L;
8949393c00Smartynas
9049393c00Smartynas /* 3/2 log x = 3 z + z^3 + z^3 (z^2 R(z^2))
9149393c00Smartynas z = (x-1)/(x+1)
9249393c00Smartynas 1 <= x <= 1.25
9349393c00Smartynas Peak relative error 2.3e-37 */
9449393c00Smartynas static const long double LN[] =
9549393c00Smartynas {
9649393c00Smartynas -3.0779177200290054398792536829702930623200E1L,
9749393c00Smartynas 6.5135778082209159921251824580292116201640E1L,
9849393c00Smartynas -4.6312921812152436921591152809994014413540E1L,
9949393c00Smartynas 1.2510208195629420304615674658258363295208E1L,
10049393c00Smartynas -9.9266909031921425609179910128531667336670E-1L
10149393c00Smartynas };
10249393c00Smartynas static const long double LD[] =
10349393c00Smartynas {
10449393c00Smartynas -5.129862866715009066465422805058933131960E1L,
10549393c00Smartynas 1.452015077564081884387441590064272782044E2L,
10649393c00Smartynas -1.524043275549860505277434040464085593165E2L,
10749393c00Smartynas 7.236063513651544224319663428634139768808E1L,
10849393c00Smartynas -1.494198912340228235853027849917095580053E1L
10949393c00Smartynas /* 1.0E0 */
11049393c00Smartynas };
11149393c00Smartynas
11249393c00Smartynas /* exp(x) = 1 + x - x / (1 - 2 / (x - x^2 R(x^2)))
11349393c00Smartynas 0 <= x <= 0.5
11449393c00Smartynas Peak relative error 5.7e-38 */
11549393c00Smartynas static const long double PN[] =
11649393c00Smartynas {
11749393c00Smartynas 5.081801691915377692446852383385968225675E8L,
11849393c00Smartynas 9.360895299872484512023336636427675327355E6L,
11949393c00Smartynas 4.213701282274196030811629773097579432957E4L,
12049393c00Smartynas 5.201006511142748908655720086041570288182E1L,
12149393c00Smartynas 9.088368420359444263703202925095675982530E-3L,
12249393c00Smartynas };
12349393c00Smartynas static const long double PD[] =
12449393c00Smartynas {
12549393c00Smartynas 3.049081015149226615468111430031590411682E9L,
12649393c00Smartynas 1.069833887183886839966085436512368982758E8L,
12749393c00Smartynas 8.259257717868875207333991924545445705394E5L,
12849393c00Smartynas 1.872583833284143212651746812884298360922E3L,
12949393c00Smartynas /* 1.0E0 */
13049393c00Smartynas };
13149393c00Smartynas
13249393c00Smartynas static const long double
13349393c00Smartynas /* ln 2 */
13449393c00Smartynas lg2 = 6.9314718055994530941723212145817656807550E-1L,
13549393c00Smartynas lg2_h = 6.9314718055994528622676398299518041312695E-1L,
13649393c00Smartynas lg2_l = 2.3190468138462996154948554638754786504121E-17L,
13749393c00Smartynas ovt = 8.0085662595372944372e-0017L,
13849393c00Smartynas /* 2/(3*log(2)) */
13949393c00Smartynas cp = 9.6179669392597560490661645400126142495110E-1L,
14049393c00Smartynas cp_h = 9.6179669392597555432899980587535537779331E-1L,
14149393c00Smartynas cp_l = 5.0577616648125906047157785230014751039424E-17L;
14249393c00Smartynas
14349393c00Smartynas long double
powl(long double x,long double y)14449393c00Smartynas powl(long double x, long double y)
14549393c00Smartynas {
14649393c00Smartynas long double z, ax, z_h, z_l, p_h, p_l;
14749393c00Smartynas long double yy1, t1, t2, r, s, t, u, v, w;
14849393c00Smartynas long double s2, s_h, s_l, t_h, t_l;
14949393c00Smartynas int32_t i, j, k, yisint, n;
15049393c00Smartynas u_int32_t ix, iy;
15149393c00Smartynas int32_t hx, hy;
15249393c00Smartynas ieee_quad_shape_type o, p, q;
15349393c00Smartynas
15449393c00Smartynas p.value = x;
15549393c00Smartynas hx = p.parts32.mswhi;
15649393c00Smartynas ix = hx & 0x7fffffff;
15749393c00Smartynas
15849393c00Smartynas q.value = y;
15949393c00Smartynas hy = q.parts32.mswhi;
16049393c00Smartynas iy = hy & 0x7fffffff;
16149393c00Smartynas
16249393c00Smartynas
16349393c00Smartynas /* y==zero: x**0 = 1 */
16449393c00Smartynas if ((iy | q.parts32.mswlo | q.parts32.lswhi | q.parts32.lswlo) == 0)
16549393c00Smartynas return one;
16649393c00Smartynas
16749393c00Smartynas /* 1.0**y = 1; -1.0**+-Inf = 1 */
16849393c00Smartynas if (x == one)
16949393c00Smartynas return one;
17049393c00Smartynas if (x == -1.0L && iy == 0x7fff0000
17149393c00Smartynas && (q.parts32.mswlo | q.parts32.lswhi | q.parts32.lswlo) == 0)
17249393c00Smartynas return one;
17349393c00Smartynas
17449393c00Smartynas /* +-NaN return x+y */
17549393c00Smartynas if ((ix > 0x7fff0000)
17649393c00Smartynas || ((ix == 0x7fff0000)
17749393c00Smartynas && ((p.parts32.mswlo | p.parts32.lswhi | p.parts32.lswlo) != 0))
17849393c00Smartynas || (iy > 0x7fff0000)
17949393c00Smartynas || ((iy == 0x7fff0000)
18049393c00Smartynas && ((q.parts32.mswlo | q.parts32.lswhi | q.parts32.lswlo) != 0)))
18149393c00Smartynas return x + y;
18249393c00Smartynas
18349393c00Smartynas /* determine if y is an odd int when x < 0
18449393c00Smartynas * yisint = 0 ... y is not an integer
18549393c00Smartynas * yisint = 1 ... y is an odd int
18649393c00Smartynas * yisint = 2 ... y is an even int
18749393c00Smartynas */
18849393c00Smartynas yisint = 0;
18949393c00Smartynas if (hx < 0)
19049393c00Smartynas {
19149393c00Smartynas if (iy >= 0x40700000) /* 2^113 */
19249393c00Smartynas yisint = 2; /* even integer y */
19349393c00Smartynas else if (iy >= 0x3fff0000) /* 1.0 */
19449393c00Smartynas {
19549393c00Smartynas if (floorl (y) == y)
19649393c00Smartynas {
19749393c00Smartynas z = 0.5 * y;
19849393c00Smartynas if (floorl (z) == z)
19949393c00Smartynas yisint = 2;
20049393c00Smartynas else
20149393c00Smartynas yisint = 1;
20249393c00Smartynas }
20349393c00Smartynas }
20449393c00Smartynas }
20549393c00Smartynas
20649393c00Smartynas /* special value of y */
20749393c00Smartynas if ((q.parts32.mswlo | q.parts32.lswhi | q.parts32.lswlo) == 0)
20849393c00Smartynas {
20949393c00Smartynas if (iy == 0x7fff0000) /* y is +-inf */
21049393c00Smartynas {
21149393c00Smartynas if (((ix - 0x3fff0000) | p.parts32.mswlo | p.parts32.lswhi |
21249393c00Smartynas p.parts32.lswlo) == 0)
21349393c00Smartynas return y - y; /* +-1**inf is NaN */
21449393c00Smartynas else if (ix >= 0x3fff0000) /* (|x|>1)**+-inf = inf,0 */
21549393c00Smartynas return (hy >= 0) ? y : zero;
21649393c00Smartynas else /* (|x|<1)**-,+inf = inf,0 */
21749393c00Smartynas return (hy < 0) ? -y : zero;
21849393c00Smartynas }
21949393c00Smartynas if (iy == 0x3fff0000)
22049393c00Smartynas { /* y is +-1 */
22149393c00Smartynas if (hy < 0)
22249393c00Smartynas return one / x;
22349393c00Smartynas else
22449393c00Smartynas return x;
22549393c00Smartynas }
22649393c00Smartynas if (hy == 0x40000000)
22749393c00Smartynas return x * x; /* y is 2 */
22849393c00Smartynas if (hy == 0x3ffe0000)
22949393c00Smartynas { /* y is 0.5 */
23049393c00Smartynas if (hx >= 0) /* x >= +0 */
23149393c00Smartynas return sqrtl (x);
23249393c00Smartynas }
23349393c00Smartynas }
23449393c00Smartynas
23549393c00Smartynas ax = fabsl (x);
23649393c00Smartynas /* special value of x */
23749393c00Smartynas if ((p.parts32.mswlo | p.parts32.lswhi | p.parts32.lswlo) == 0)
23849393c00Smartynas {
23949393c00Smartynas if (ix == 0x7fff0000 || ix == 0 || ix == 0x3fff0000)
24049393c00Smartynas {
24149393c00Smartynas z = ax; /*x is +-0,+-inf,+-1 */
24249393c00Smartynas if (hy < 0)
24349393c00Smartynas z = one / z; /* z = (1/|x|) */
24449393c00Smartynas if (hx < 0)
24549393c00Smartynas {
24649393c00Smartynas if (((ix - 0x3fff0000) | yisint) == 0)
24749393c00Smartynas {
24849393c00Smartynas z = (z - z) / (z - z); /* (-1)**non-int is NaN */
24949393c00Smartynas }
25049393c00Smartynas else if (yisint == 1)
25149393c00Smartynas z = -z; /* (x<0)**odd = -(|x|**odd) */
25249393c00Smartynas }
25349393c00Smartynas return z;
25449393c00Smartynas }
25549393c00Smartynas }
25649393c00Smartynas
25749393c00Smartynas /* (x<0)**(non-int) is NaN */
25849393c00Smartynas if (((((u_int32_t) hx >> 31) - 1) | yisint) == 0)
25949393c00Smartynas return (x - x) / (x - x);
26049393c00Smartynas
26149393c00Smartynas /* |y| is huge.
26249393c00Smartynas 2^-16495 = 1/2 of smallest representable value.
26349393c00Smartynas If (1 - 1/131072)^y underflows, y > 1.4986e9 */
26449393c00Smartynas if (iy > 0x401d654b)
26549393c00Smartynas {
26649393c00Smartynas /* if (1 - 2^-113)^y underflows, y > 1.1873e38 */
26749393c00Smartynas if (iy > 0x407d654b)
26849393c00Smartynas {
26949393c00Smartynas if (ix <= 0x3ffeffff)
27049393c00Smartynas return (hy < 0) ? huge * huge : tiny * tiny;
27149393c00Smartynas if (ix >= 0x3fff0000)
27249393c00Smartynas return (hy > 0) ? huge * huge : tiny * tiny;
27349393c00Smartynas }
27449393c00Smartynas /* over/underflow if x is not close to one */
27549393c00Smartynas if (ix < 0x3ffeffff)
27649393c00Smartynas return (hy < 0) ? huge * huge : tiny * tiny;
27749393c00Smartynas if (ix > 0x3fff0000)
27849393c00Smartynas return (hy > 0) ? huge * huge : tiny * tiny;
27949393c00Smartynas }
28049393c00Smartynas
28149393c00Smartynas n = 0;
28249393c00Smartynas /* take care subnormal number */
28349393c00Smartynas if (ix < 0x00010000)
28449393c00Smartynas {
28549393c00Smartynas ax *= two113;
28649393c00Smartynas n -= 113;
28749393c00Smartynas o.value = ax;
28849393c00Smartynas ix = o.parts32.mswhi;
28949393c00Smartynas }
29049393c00Smartynas n += ((ix) >> 16) - 0x3fff;
29149393c00Smartynas j = ix & 0x0000ffff;
29249393c00Smartynas /* determine interval */
29349393c00Smartynas ix = j | 0x3fff0000; /* normalize ix */
29449393c00Smartynas if (j <= 0x3988)
29549393c00Smartynas k = 0; /* |x|<sqrt(3/2) */
29649393c00Smartynas else if (j < 0xbb67)
29749393c00Smartynas k = 1; /* |x|<sqrt(3) */
29849393c00Smartynas else
29949393c00Smartynas {
30049393c00Smartynas k = 0;
30149393c00Smartynas n += 1;
30249393c00Smartynas ix -= 0x00010000;
30349393c00Smartynas }
30449393c00Smartynas
30549393c00Smartynas o.value = ax;
30649393c00Smartynas o.parts32.mswhi = ix;
30749393c00Smartynas ax = o.value;
30849393c00Smartynas
30949393c00Smartynas /* compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
31049393c00Smartynas u = ax - bp[k]; /* bp[0]=1.0, bp[1]=1.5 */
31149393c00Smartynas v = one / (ax + bp[k]);
31249393c00Smartynas s = u * v;
31349393c00Smartynas s_h = s;
31449393c00Smartynas
31549393c00Smartynas o.value = s_h;
31649393c00Smartynas o.parts32.lswlo = 0;
31749393c00Smartynas o.parts32.lswhi &= 0xf8000000;
31849393c00Smartynas s_h = o.value;
31949393c00Smartynas /* t_h=ax+bp[k] High */
32049393c00Smartynas t_h = ax + bp[k];
32149393c00Smartynas o.value = t_h;
32249393c00Smartynas o.parts32.lswlo = 0;
32349393c00Smartynas o.parts32.lswhi &= 0xf8000000;
32449393c00Smartynas t_h = o.value;
32549393c00Smartynas t_l = ax - (t_h - bp[k]);
32649393c00Smartynas s_l = v * ((u - s_h * t_h) - s_h * t_l);
32749393c00Smartynas /* compute log(ax) */
32849393c00Smartynas s2 = s * s;
32949393c00Smartynas u = LN[0] + s2 * (LN[1] + s2 * (LN[2] + s2 * (LN[3] + s2 * LN[4])));
33049393c00Smartynas v = LD[0] + s2 * (LD[1] + s2 * (LD[2] + s2 * (LD[3] + s2 * (LD[4] + s2))));
33149393c00Smartynas r = s2 * s2 * u / v;
33249393c00Smartynas r += s_l * (s_h + s);
33349393c00Smartynas s2 = s_h * s_h;
33449393c00Smartynas t_h = 3.0 + s2 + r;
33549393c00Smartynas o.value = t_h;
33649393c00Smartynas o.parts32.lswlo = 0;
33749393c00Smartynas o.parts32.lswhi &= 0xf8000000;
33849393c00Smartynas t_h = o.value;
33949393c00Smartynas t_l = r - ((t_h - 3.0) - s2);
34049393c00Smartynas /* u+v = s*(1+...) */
34149393c00Smartynas u = s_h * t_h;
34249393c00Smartynas v = s_l * t_h + t_l * s;
34349393c00Smartynas /* 2/(3log2)*(s+...) */
34449393c00Smartynas p_h = u + v;
34549393c00Smartynas o.value = p_h;
34649393c00Smartynas o.parts32.lswlo = 0;
34749393c00Smartynas o.parts32.lswhi &= 0xf8000000;
34849393c00Smartynas p_h = o.value;
34949393c00Smartynas p_l = v - (p_h - u);
35049393c00Smartynas z_h = cp_h * p_h; /* cp_h+cp_l = 2/(3*log2) */
35149393c00Smartynas z_l = cp_l * p_h + p_l * cp + dp_l[k];
35249393c00Smartynas /* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */
35349393c00Smartynas t = (long double) n;
35449393c00Smartynas t1 = (((z_h + z_l) + dp_h[k]) + t);
35549393c00Smartynas o.value = t1;
35649393c00Smartynas o.parts32.lswlo = 0;
35749393c00Smartynas o.parts32.lswhi &= 0xf8000000;
35849393c00Smartynas t1 = o.value;
35949393c00Smartynas t2 = z_l - (((t1 - t) - dp_h[k]) - z_h);
36049393c00Smartynas
36149393c00Smartynas /* s (sign of result -ve**odd) = -1 else = 1 */
36249393c00Smartynas s = one;
36349393c00Smartynas if (((((u_int32_t) hx >> 31) - 1) | (yisint - 1)) == 0)
36449393c00Smartynas s = -one; /* (-ve)**(odd int) */
36549393c00Smartynas
36649393c00Smartynas /* split up y into yy1+y2 and compute (yy1+y2)*(t1+t2) */
36749393c00Smartynas yy1 = y;
36849393c00Smartynas o.value = yy1;
36949393c00Smartynas o.parts32.lswlo = 0;
37049393c00Smartynas o.parts32.lswhi &= 0xf8000000;
37149393c00Smartynas yy1 = o.value;
37249393c00Smartynas p_l = (y - yy1) * t1 + y * t2;
37349393c00Smartynas p_h = yy1 * t1;
37449393c00Smartynas z = p_l + p_h;
37549393c00Smartynas o.value = z;
37649393c00Smartynas j = o.parts32.mswhi;
37749393c00Smartynas if (j >= 0x400d0000) /* z >= 16384 */
37849393c00Smartynas {
37949393c00Smartynas /* if z > 16384 */
38049393c00Smartynas if (((j - 0x400d0000) | o.parts32.mswlo | o.parts32.lswhi |
38149393c00Smartynas o.parts32.lswlo) != 0)
38249393c00Smartynas return s * huge * huge; /* overflow */
38349393c00Smartynas else
38449393c00Smartynas {
38549393c00Smartynas if (p_l + ovt > z - p_h)
38649393c00Smartynas return s * huge * huge; /* overflow */
38749393c00Smartynas }
38849393c00Smartynas }
38949393c00Smartynas else if ((j & 0x7fffffff) >= 0x400d01b9) /* z <= -16495 */
39049393c00Smartynas {
39149393c00Smartynas /* z < -16495 */
39249393c00Smartynas if (((j - 0xc00d01bc) | o.parts32.mswlo | o.parts32.lswhi |
39349393c00Smartynas o.parts32.lswlo)
39449393c00Smartynas != 0)
39549393c00Smartynas return s * tiny * tiny; /* underflow */
39649393c00Smartynas else
39749393c00Smartynas {
39849393c00Smartynas if (p_l <= z - p_h)
39949393c00Smartynas return s * tiny * tiny; /* underflow */
40049393c00Smartynas }
40149393c00Smartynas }
40249393c00Smartynas /* compute 2**(p_h+p_l) */
40349393c00Smartynas i = j & 0x7fffffff;
40449393c00Smartynas k = (i >> 16) - 0x3fff;
40549393c00Smartynas n = 0;
40649393c00Smartynas if (i > 0x3ffe0000)
40749393c00Smartynas { /* if |z| > 0.5, set n = [z+0.5] */
40849393c00Smartynas n = floorl (z + 0.5L);
40949393c00Smartynas t = n;
41049393c00Smartynas p_h -= t;
41149393c00Smartynas }
41249393c00Smartynas t = p_l + p_h;
41349393c00Smartynas o.value = t;
41449393c00Smartynas o.parts32.lswlo = 0;
41549393c00Smartynas o.parts32.lswhi &= 0xf8000000;
41649393c00Smartynas t = o.value;
41749393c00Smartynas u = t * lg2_h;
41849393c00Smartynas v = (p_l - (t - p_h)) * lg2 + t * lg2_l;
41949393c00Smartynas z = u + v;
42049393c00Smartynas w = v - (z - u);
42149393c00Smartynas /* exp(z) */
42249393c00Smartynas t = z * z;
42349393c00Smartynas u = PN[0] + t * (PN[1] + t * (PN[2] + t * (PN[3] + t * PN[4])));
42449393c00Smartynas v = PD[0] + t * (PD[1] + t * (PD[2] + t * (PD[3] + t)));
42549393c00Smartynas t1 = z - t * u / v;
42649393c00Smartynas r = (z * t1) / (t1 - two) - (w + z * w);
42749393c00Smartynas z = one - (r - z);
42849393c00Smartynas o.value = z;
42949393c00Smartynas j = o.parts32.mswhi;
43049393c00Smartynas j += (n << 16);
43149393c00Smartynas if ((j >> 16) <= 0)
43249393c00Smartynas z = scalbnl (z, n); /* subnormal output */
43349393c00Smartynas else
43449393c00Smartynas {
43549393c00Smartynas o.parts32.mswhi = j;
43649393c00Smartynas z = o.value;
43749393c00Smartynas }
43849393c00Smartynas return s * z;
43949393c00Smartynas }
440*2f2c0062Sguenther DEF_STD(powl);
441