You should have received a copy of the GNU Lesser General Public
License along with this library; if not, see
- <http://www.gnu.org/licenses/>. */
+ <https://www.gnu.org/licenses/>. */
/* __ieee754_powl(x,y) return x**y
*
*/
#include <math.h>
+#include <math-barriers.h>
#include <math_private.h>
-static const long double bp[] = {
- 1.0L,
- 1.5L,
+static const _Float128 bp[] = {
+ 1,
+ L(1.5),
};
/* log_2(1.5) */
-static const long double dp_h[] = {
+static const _Float128 dp_h[] = {
0.0,
- 5.8496250072115607565592654282227158546448E-1L
+ L(5.8496250072115607565592654282227158546448E-1)
};
/* Low part of log_2(1.5) */
-static const long double dp_l[] = {
+static const _Float128 dp_l[] = {
0.0,
- 1.0579781240112554492329533686862998106046E-16L
+ L(1.0579781240112554492329533686862998106046E-16)
};
-static const long double zero = 0.0L,
- one = 1.0L,
- two = 2.0L,
- two113 = 1.0384593717069655257060992658440192E34L,
- huge = 1.0e3000L,
- tiny = 1.0e-3000L;
+static const _Float128 zero = 0,
+ one = 1,
+ two = 2,
+ two113 = L(1.0384593717069655257060992658440192E34),
+ huge = L(1.0e3000),
+ tiny = L(1.0e-3000);
/* 3/2 log x = 3 z + z^3 + z^3 (z^2 R(z^2))
z = (x-1)/(x+1)
1 <= x <= 1.25
Peak relative error 2.3e-37 */
-static const long double LN[] =
+static const _Float128 LN[] =
{
- -3.0779177200290054398792536829702930623200E1L,
- 6.5135778082209159921251824580292116201640E1L,
- -4.6312921812152436921591152809994014413540E1L,
- 1.2510208195629420304615674658258363295208E1L,
- -9.9266909031921425609179910128531667336670E-1L
+ L(-3.0779177200290054398792536829702930623200E1),
+ L(6.5135778082209159921251824580292116201640E1),
+ L(-4.6312921812152436921591152809994014413540E1),
+ L(1.2510208195629420304615674658258363295208E1),
+ L(-9.9266909031921425609179910128531667336670E-1)
};
-static const long double LD[] =
+static const _Float128 LD[] =
{
- -5.129862866715009066465422805058933131960E1L,
- 1.452015077564081884387441590064272782044E2L,
- -1.524043275549860505277434040464085593165E2L,
- 7.236063513651544224319663428634139768808E1L,
- -1.494198912340228235853027849917095580053E1L
+ L(-5.129862866715009066465422805058933131960E1),
+ L(1.452015077564081884387441590064272782044E2),
+ L(-1.524043275549860505277434040464085593165E2),
+ L(7.236063513651544224319663428634139768808E1),
+ L(-1.494198912340228235853027849917095580053E1)
/* 1.0E0 */
};
/* exp(x) = 1 + x - x / (1 - 2 / (x - x^2 R(x^2)))
0 <= x <= 0.5
Peak relative error 5.7e-38 */
-static const long double PN[] =
+static const _Float128 PN[] =
{
- 5.081801691915377692446852383385968225675E8L,
- 9.360895299872484512023336636427675327355E6L,
- 4.213701282274196030811629773097579432957E4L,
- 5.201006511142748908655720086041570288182E1L,
- 9.088368420359444263703202925095675982530E-3L,
+ L(5.081801691915377692446852383385968225675E8),
+ L(9.360895299872484512023336636427675327355E6),
+ L(4.213701282274196030811629773097579432957E4),
+ L(5.201006511142748908655720086041570288182E1),
+ L(9.088368420359444263703202925095675982530E-3),
};
-static const long double PD[] =
+static const _Float128 PD[] =
{
- 3.049081015149226615468111430031590411682E9L,
- 1.069833887183886839966085436512368982758E8L,
- 8.259257717868875207333991924545445705394E5L,
- 1.872583833284143212651746812884298360922E3L,
+ L(3.049081015149226615468111430031590411682E9),
+ L(1.069833887183886839966085436512368982758E8),
+ L(8.259257717868875207333991924545445705394E5),
+ L(1.872583833284143212651746812884298360922E3),
/* 1.0E0 */
};
-static const long double
+static const _Float128
/* ln 2 */
- lg2 = 6.9314718055994530941723212145817656807550E-1L,
- lg2_h = 6.9314718055994528622676398299518041312695E-1L,
- lg2_l = 2.3190468138462996154948554638754786504121E-17L,
- ovt = 8.0085662595372944372e-0017L,
+ lg2 = L(6.9314718055994530941723212145817656807550E-1),
+ lg2_h = L(6.9314718055994528622676398299518041312695E-1),
+ lg2_l = L(2.3190468138462996154948554638754786504121E-17),
+ ovt = L(8.0085662595372944372e-0017),
/* 2/(3*log(2)) */
- cp = 9.6179669392597560490661645400126142495110E-1L,
- cp_h = 9.6179669392597555432899980587535537779331E-1L,
- cp_l = 5.0577616648125906047157785230014751039424E-17L;
+ cp = L(9.6179669392597560490661645400126142495110E-1),
+ cp_h = L(9.6179669392597555432899980587535537779331E-1),
+ cp_l = L(5.0577616648125906047157785230014751039424E-17);
-long double
-__ieee754_powl (long double x, long double y)
+_Float128
+__ieee754_powl (_Float128 x, _Float128 y)
{
- long double z, ax, z_h, z_l, p_h, p_l;
- long double y1, t1, t2, r, s, t, u, v, w;
- long double s2, s_h, s_l, t_h, t_l;
+ _Float128 z, ax, z_h, z_l, p_h, p_l;
+ _Float128 y1, t1, t2, r, s, sgn, t, u, v, w;
+ _Float128 s2, s_h, s_l, t_h, t_l, ay;
int32_t i, j, k, yisint, n;
- u_int32_t ix, iy;
+ uint32_t ix, iy;
int32_t hx, hy;
ieee854_long_double_shape_type o, p, q;
/* y==zero: x**0 = 1 */
- if ((iy | q.parts32.w1 | q.parts32.w2 | q.parts32.w3) == 0)
+ if ((iy | q.parts32.w1 | q.parts32.w2 | q.parts32.w3) == 0
+ && !issignaling (x))
return one;
/* 1.0**y = 1; -1.0**+-Inf = 1 */
- if (x == one)
+ if (x == one && !issignaling (y))
return one;
- if (x == -1.0L && iy == 0x7fff0000
+ if (x == -1 && iy == 0x7fff0000
&& (q.parts32.w1 | q.parts32.w2 | q.parts32.w3) == 0)
return one;
yisint = 2; /* even integer y */
else if (iy >= 0x3fff0000) /* 1.0 */
{
- if (__floorl (y) == y)
+ if (floorl (y) == y)
{
z = 0.5 * y;
- if (__floorl (z) == z)
+ if (floorl (z) == z)
yisint = 2;
else
yisint = 1;
if (hy == 0x3ffe0000)
{ /* y is 0.5 */
if (hx >= 0) /* x >= +0 */
- return __ieee754_sqrtl (x);
+ return sqrtl (x);
}
}
}
/* (x<0)**(non-int) is NaN */
- if (((((u_int32_t) hx >> 31) - 1) | yisint) == 0)
+ if (((((uint32_t) hx >> 31) - 1) | yisint) == 0)
return (x - x) / (x - x);
+ /* sgn (sign of result -ve**odd) = -1 else = 1 */
+ sgn = one;
+ if (((((uint32_t) hx >> 31) - 1) | (yisint - 1)) == 0)
+ sgn = -one; /* (-ve)**(odd int) */
+
/* |y| is huge.
2^-16495 = 1/2 of smallest representable value.
If (1 - 1/131072)^y underflows, y > 1.4986e9 */
}
/* over/underflow if x is not close to one */
if (ix < 0x3ffeffff)
- return (hy < 0) ? huge * huge : tiny * tiny;
+ return (hy < 0) ? sgn * huge * huge : sgn * tiny * tiny;
if (ix > 0x3fff0000)
- return (hy > 0) ? huge * huge : tiny * tiny;
+ return (hy > 0) ? sgn * huge * huge : sgn * tiny * tiny;
}
+ ay = y > 0 ? y : -y;
+ if (ay < 0x1p-128)
+ y = y < 0 ? -0x1p-128 : 0x1p-128;
+
n = 0;
/* take care subnormal number */
if (ix < 0x00010000)
z_h = cp_h * p_h; /* cp_h+cp_l = 2/(3*log2) */
z_l = cp_l * p_h + p_l * cp + dp_l[k];
/* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */
- t = (long double) n;
+ t = (_Float128) n;
t1 = (((z_h + z_l) + dp_h[k]) + t);
o.value = t1;
o.parts32.w3 = 0;
t1 = o.value;
t2 = z_l - (((t1 - t) - dp_h[k]) - z_h);
- /* s (sign of result -ve**odd) = -1 else = 1 */
- s = one;
- if (((((u_int32_t) hx >> 31) - 1) | (yisint - 1)) == 0)
- s = -one; /* (-ve)**(odd int) */
-
/* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */
y1 = y;
o.value = y1;
{
/* if z > 16384 */
if (((j - 0x400d0000) | o.parts32.w1 | o.parts32.w2 | o.parts32.w3) != 0)
- return s * huge * huge; /* overflow */
+ return sgn * huge * huge; /* overflow */
else
{
if (p_l + ovt > z - p_h)
- return s * huge * huge; /* overflow */
+ return sgn * huge * huge; /* overflow */
}
}
else if ((j & 0x7fffffff) >= 0x400d01b9) /* z <= -16495 */
/* z < -16495 */
if (((j - 0xc00d01bc) | o.parts32.w1 | o.parts32.w2 | o.parts32.w3)
!= 0)
- return s * tiny * tiny; /* underflow */
+ return sgn * tiny * tiny; /* underflow */
else
{
if (p_l <= z - p_h)
- return s * tiny * tiny; /* underflow */
+ return sgn * tiny * tiny; /* underflow */
}
}
/* compute 2**(p_h+p_l) */
n = 0;
if (i > 0x3ffe0000)
{ /* if |z| > 0.5, set n = [z+0.5] */
- n = __floorl (z + 0.5L);
+ n = floorl (z + L(0.5));
t = n;
p_h -= t;
}
j = o.parts32.w0;
j += (n << 16);
if ((j >> 16) <= 0)
- z = __scalbnl (z, n); /* subnormal output */
+ {
+ z = __scalbnl (z, n); /* subnormal output */
+ _Float128 force_underflow = z * z;
+ math_force_eval (force_underflow);
+ }
else
{
o.parts32.w0 = j;
z = o.value;
}
- return s * z;
+ return sgn * z;
}
strong_alias (__ieee754_powl, __powl_finite)
projects / thirdparty / glibc.git / blobdiff