libm: sync with upstream.
There's potential here to maybe lose some/all of builtins.cpp, but I'll
look at that separately later.
Test: treehugger
Change-Id: I2c2bc1d0753affdd214daeb09fa1ac7cd73db347
diff --git a/libm/upstream-freebsd/lib/msun/src/e_hypotl.c b/libm/upstream-freebsd/lib/msun/src/e_hypotl.c
index 9189b1f..fc43538 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_hypotl.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_hypotl.c
@@ -82,7 +82,7 @@
man_t manh, manl;
GET_LDBL_MAN(manh,manl,b);
if((manh|manl)==0) return a;
- t1=0;
+ t1=1;
SET_HIGH_WORD(t1,ESW(MAX_EXP-2)); /* t1=2^(MAX_EXP-2) */
b *= t1;
a *= t1;
diff --git a/libm/upstream-freebsd/lib/msun/src/e_powf.c b/libm/upstream-freebsd/lib/msun/src/e_powf.c
index 53f1d37..122da45 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_powf.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_powf.c
@@ -136,7 +136,7 @@
/* |y| is huge */
if(iy>0x4d000000) { /* if |y| > 2**27 */
/* over/underflow if x is not close to one */
- if(ix<0x3f7ffff8) return (hy<0)? sn*huge*huge:sn*tiny*tiny;
+ if(ix<0x3f7ffff6) return (hy<0)? sn*huge*huge:sn*tiny*tiny;
if(ix>0x3f800007) return (hy>0)? sn*huge*huge:sn*tiny*tiny;
/* now |1-x| is tiny <= 2**-20, suffice to compute
log(x) by x-x^2/2+x^3/3-x^4/4 */
diff --git a/libm/upstream-freebsd/lib/msun/src/e_sqrt.c b/libm/upstream-freebsd/lib/msun/src/e_sqrt.c
index 12fb56e..37351a4 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_sqrt.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_sqrt.c
@@ -14,6 +14,18 @@
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
+#include <float.h>
+
+#include "math.h"
+#include "math_private.h"
+
+#ifdef USE_BUILTIN_SQRT
+double
+__ieee754_sqrt(double x)
+{
+ return (__builtin_sqrt(x));
+}
+#else
/* __ieee754_sqrt(x)
* Return correctly rounded sqrt.
* ------------------------------------------
@@ -84,11 +96,6 @@
*---------------
*/
-#include <float.h>
-
-#include "math.h"
-#include "math_private.h"
-
static const double one = 1.0, tiny=1.0e-300;
double
@@ -187,6 +194,7 @@
INSERT_WORDS(z,ix0,ix1);
return z;
}
+#endif
#if (LDBL_MANT_DIG == 53)
__weak_reference(sqrt, sqrtl);
diff --git a/libm/upstream-freebsd/lib/msun/src/e_sqrtf.c b/libm/upstream-freebsd/lib/msun/src/e_sqrtf.c
index 7eba4d0..06e5d62 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_sqrtf.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_sqrtf.c
@@ -20,6 +20,13 @@
#include "math.h"
#include "math_private.h"
+#ifdef USE_BUILTIN_SQRTF
+float
+__ieee754_sqrtf(float x)
+{
+ return (__builtin_sqrtf(x));
+}
+#else
static const float one = 1.0, tiny=1.0e-30;
float
@@ -87,3 +94,4 @@
SET_FLOAT_WORD(z,ix);
return z;
}
+#endif
diff --git a/libm/upstream-freebsd/lib/msun/src/k_cospi.h b/libm/upstream-freebsd/lib/msun/src/k_cospi.h
new file mode 100644
index 0000000..32985be
--- /dev/null
+++ b/libm/upstream-freebsd/lib/msun/src/k_cospi.h
@@ -0,0 +1,44 @@
+/*-
+ * Copyright (c) 2017 Steven G. Kargl
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ * notice unmodified, this list of conditions, and the following
+ * disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ * notice, this list of conditions and the following disclaimer in the
+ * documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
+ * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
+ * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
+ * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
+ * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
+ * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+ * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+ * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
+ * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ */
+
+/*
+ * The basic kernel for x in [0,0.25]. To use the kernel for cos(x), the
+ * argument to __kernel_cospi() must be multiplied by pi.
+ */
+
+static inline double
+__kernel_cospi(double x)
+{
+ double_t hi, lo;
+
+ hi = (float)x;
+ lo = x - hi;
+ lo = lo * (pi_lo + pi_hi) + hi * pi_lo;
+ hi *= pi_hi;
+ _2sumF(hi, lo);
+ return (__kernel_cos(hi, lo));
+}
+
diff --git a/libm/upstream-freebsd/lib/msun/src/k_sincosl.h b/libm/upstream-freebsd/lib/msun/src/k_sincosl.h
index 4d4dc69..6425f14 100644
--- a/libm/upstream-freebsd/lib/msun/src/k_sincosl.h
+++ b/libm/upstream-freebsd/lib/msun/src/k_sincosl.h
@@ -76,13 +76,6 @@
#elif LDBL_MANT_DIG == 113 /* ld128 version of k_sincosl.c. */
static const long double
-C1 = 0.04166666666666666666666666666666658424671L,
-C2 = -0.001388888888888888888888888888863490893732L,
-C3 = 0.00002480158730158730158730158600795304914210L,
-C4 = -0.2755731922398589065255474947078934284324e-6L,
-C5 = 0.2087675698786809897659225313136400793948e-8L,
-C6 = -0.1147074559772972315817149986812031204775e-10L,
-C7 = 0.4779477332386808976875457937252120293400e-13L,
S1 = -0.16666666666666666666666666666666666606732416116558L,
S2 = 0.0083333333333333333333333333333331135404851288270047L,
S3 = -0.00019841269841269841269841269839935785325638310428717L,
@@ -93,15 +86,25 @@
S8 = 0.28114572543451292625024967174638477283187397621303e-14L;
static const double
-C8 = -0.1561920696721507929516718307820958119868e-15,
-C9 = 0.4110317413744594971475941557607804508039e-18,
-C10 = -0.8896592467191938803288521958313920156409e-21,
-C11 = 0.1601061435794535138244346256065192782581e-23,
S9 = -0.82206352458348947812512122163446202498005154296863e-17,
S10 = 0.19572940011906109418080609928334380560135358385256e-19,
S11 = -0.38680813379701966970673724299207480965452616911420e-22,
S12 = 0.64038150078671872796678569586315881020659912139412e-25;
+static const long double
+C1 = 4.16666666666666666666666666666666667e-02L,
+C2 = -1.38888888888888888888888888888888834e-03L,
+C3 = 2.48015873015873015873015873015446795e-05L,
+C4 = -2.75573192239858906525573190949988493e-07L,
+C5 = 2.08767569878680989792098886701451072e-09L,
+C6 = -1.14707455977297247136657111139971865e-11L,
+C7 = 4.77947733238738518870113294139830239e-14L,
+C8 = -1.56192069685858079920640872925306403e-16L,
+C9 = 4.11031762320473354032038893429515732e-19L,
+C10= -8.89679121027589608738005163931958096e-22L,
+C11= 1.61171797801314301767074036661901531e-24L,
+C12= -2.46748624357670948912574279501044295e-27L;
+
static inline void
__kernel_sincosl(long double x, long double y, int iy, long double *sn,
long double *cs)
@@ -120,12 +123,12 @@
if (iy == 0)
*sn = x + v * (S1 + z * r);
else
- *cs = x - ((z * (y / 2 - v * r) - y) - v * S1);
+ *sn = x - ((z * (y / 2 - v * r) - y) - v * S1);
hz = z / 2;
w = 1 - hz;
r = z * (C1 + z * (C2 + z * (C3 + z * (C4 + z * (C5 + z * (C6 +
- z * (C7 + z * (C8 + z * (C9 + z * (C10 + z * C11))))))))));
+ z * (C7 + z * (C8 + z * (C9 + z * (C10 + z * (C11+z*C12)))))))))));
*cs = w + (((1 - w) - hz) + (z * r - x * y));
}
diff --git a/libm/upstream-freebsd/lib/msun/src/k_sinpi.h b/libm/upstream-freebsd/lib/msun/src/k_sinpi.h
new file mode 100644
index 0000000..66152ba
--- /dev/null
+++ b/libm/upstream-freebsd/lib/msun/src/k_sinpi.h
@@ -0,0 +1,43 @@
+/*-
+ * Copyright (c) 2017 Steven G. Kargl
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ * notice unmodified, this list of conditions, and the following
+ * disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ * notice, this list of conditions and the following disclaimer in the
+ * documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
+ * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
+ * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
+ * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
+ * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
+ * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+ * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+ * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
+ * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ */
+
+/*
+ * The basic kernel for x in [0,0.25]. To use the kernel for sin(x), the
+ * argument to __kernel_sinpi() must be multiplied by pi.
+ */
+
+static inline double
+__kernel_sinpi(double x)
+{
+ double_t hi, lo;
+
+ hi = (float)x;
+ lo = x - hi;
+ lo = lo * (pi_lo + pi_hi) + hi * pi_lo;
+ hi *= pi_hi;
+ _2sumF(hi, lo);
+ return (__kernel_sin(hi, lo, 1));
+}
diff --git a/libm/upstream-freebsd/lib/msun/src/s_cexp.c b/libm/upstream-freebsd/lib/msun/src/s_cexp.c
index 2ef8ba1..a1f853e 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_cexp.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_cexp.c
@@ -30,6 +30,7 @@
__FBSDID("$FreeBSD$");
#include <complex.h>
+#include <float.h>
#include <math.h>
#include "math_private.h"
@@ -41,7 +42,7 @@
double complex
cexp(double complex z)
{
- double x, y, exp_x;
+ double c, exp_x, s, x, y;
uint32_t hx, hy, lx, ly;
x = creal(z);
@@ -55,8 +56,10 @@
return (CMPLX(exp(x), y));
EXTRACT_WORDS(hx, lx, x);
/* cexp(0 + I y) = cos(y) + I sin(y) */
- if (((hx & 0x7fffffff) | lx) == 0)
- return (CMPLX(cos(y), sin(y)));
+ if (((hx & 0x7fffffff) | lx) == 0) {
+ sincos(y, &s, &c);
+ return (CMPLX(c, s));
+ }
if (hy >= 0x7ff00000) {
if (lx != 0 || (hx & 0x7fffffff) != 0x7ff00000) {
@@ -86,6 +89,11 @@
* - x = NaN (spurious inexact exception from y)
*/
exp_x = exp(x);
- return (CMPLX(exp_x * cos(y), exp_x * sin(y)));
+ sincos(y, &s, &c);
+ return (CMPLX(exp_x * c, exp_x * s));
}
}
+
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(cexp, cexpl);
+#endif
diff --git a/libm/upstream-freebsd/lib/msun/src/s_cexpf.c b/libm/upstream-freebsd/lib/msun/src/s_cexpf.c
index b815c99..d905b74 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_cexpf.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_cexpf.c
@@ -41,7 +41,7 @@
float complex
cexpf(float complex z)
{
- float x, y, exp_x;
+ float c, exp_x, s, x, y;
uint32_t hx, hy;
x = crealf(z);
@@ -55,8 +55,10 @@
return (CMPLXF(expf(x), y));
GET_FLOAT_WORD(hx, x);
/* cexp(0 + I y) = cos(y) + I sin(y) */
- if ((hx & 0x7fffffff) == 0)
- return (CMPLXF(cosf(y), sinf(y)));
+ if ((hx & 0x7fffffff) == 0) {
+ sincosf(y, &s, &c);
+ return (CMPLXF(c, s));
+ }
if (hy >= 0x7f800000) {
if ((hx & 0x7fffffff) != 0x7f800000) {
@@ -86,6 +88,7 @@
* - x = NaN (spurious inexact exception from y)
*/
exp_x = expf(x);
- return (CMPLXF(exp_x * cosf(y), exp_x * sinf(y)));
+ sincosf(y, &s, &c);
+ return (CMPLXF(exp_x * c, exp_x * s));
}
}
diff --git a/libm/upstream-freebsd/lib/msun/src/s_cosl.c b/libm/upstream-freebsd/lib/msun/src/s_cosl.c
index 46a2e86..3d06648 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_cosl.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_cosl.c
@@ -39,12 +39,17 @@
#include <ieeefp.h>
#endif
+#include "fpmath.h"
#include "math.h"
#include "math_private.h"
#if LDBL_MANT_DIG == 64
#include "../ld80/e_rem_pio2l.h"
+static const union IEEEl2bits
+pio4u = LD80C(0xc90fdaa22168c235, -00001, 7.85398163397448309628e-01L);
+#define pio4 (pio4u.e)
#elif LDBL_MANT_DIG == 113
#include "../ld128/e_rem_pio2l.h"
+long double pio4 = 7.85398163397448309615660845819875721e-1L;
#else
#error "Unsupported long double format"
#endif
@@ -71,7 +76,7 @@
ENTERI();
/* Optimize the case where x is already within range. */
- if (z.e < M_PI_4)
+ if (z.e < pio4)
RETURNI(__kernel_cosl(z.e, 0));
e0 = __ieee754_rem_pio2l(x, y);
diff --git a/libm/upstream-freebsd/lib/msun/src/s_cospi.c b/libm/upstream-freebsd/lib/msun/src/s_cospi.c
new file mode 100644
index 0000000..860219e
--- /dev/null
+++ b/libm/upstream-freebsd/lib/msun/src/s_cospi.c
@@ -0,0 +1,152 @@
+/*-
+ * Copyright (c) 2017 Steven G. Kargl
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ * notice unmodified, this list of conditions, and the following
+ * disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ * notice, this list of conditions and the following disclaimer in the
+ * documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
+ * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
+ * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
+ * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
+ * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
+ * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+ * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+ * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
+ * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ */
+
+/**
+ * cospi(x) computes cos(pi*x) without multiplication by pi (almost). First,
+ * note that cospi(-x) = cospi(x), so the algorithm considers only |x|. The
+ * method used depends on the magnitude of x.
+ *
+ * 1. For small |x|, cospi(x) = 1 with FE_INEXACT raised where a sloppy
+ * threshold is used. The threshold is |x| < 0x1pN with N = -(P/2+M).
+ * P is the precision of the floating-point type and M = 2 to 4.
+ *
+ * 2. For |x| < 1, argument reduction is not required and sinpi(x) is
+ * computed by calling a kernel that leverages the kernels for sin(x)
+ * ans cos(x). See k_sinpi.c and k_cospi.c for details.
+ *
+ * 3. For 1 <= |x| < 0x1p(P-1), argument reduction is required where
+ * |x| = j0 + r with j0 an integer and the remainder r satisfies
+ * 0 <= r < 1. With the given domain, a simplified inline floor(x)
+ * is used. Also, note the following identity
+ *
+ * cospi(x) = cos(pi*(j0+r))
+ * = cos(pi*j0) * cos(pi*r) - sin(pi*j0) * sin(pi*r)
+ * = cos(pi*j0) * cos(pi*r)
+ * = +-cospi(r)
+ *
+ * If j0 is even, then cos(pi*j0) = 1. If j0 is odd, then cos(pi*j0) = -1.
+ * cospi(r) is then computed via an appropriate kernel.
+ *
+ * 4. For |x| >= 0x1p(P-1), |x| is integral and cospi(x) = 1.
+ *
+ * 5. Special cases:
+ *
+ * cospi(+-0) = 1.
+ * cospi(n.5) = 0 for n an integer.
+ * cospi(+-inf) = nan. Raises the "invalid" floating-point exception.
+ * cospi(nan) = nan. Raises the "invalid" floating-point exception.
+ */
+
+#include <float.h>
+#include "math.h"
+#include "math_private.h"
+
+static const double
+pi_hi = 3.1415926814079285e+00, /* 0x400921fb 0x58000000 */
+pi_lo =-2.7818135228334233e-08; /* 0xbe5dde97 0x3dcb3b3a */
+
+#include "k_cospi.h"
+#include "k_sinpi.h"
+
+volatile static const double vzero = 0;
+
+double
+cospi(double x)
+{
+ double ax, c;
+ uint32_t hx, ix, j0, lx;
+
+ EXTRACT_WORDS(hx, lx, x);
+ ix = hx & 0x7fffffff;
+ INSERT_WORDS(ax, ix, lx);
+
+ if (ix < 0x3ff00000) { /* |x| < 1 */
+ if (ix < 0x3fd00000) { /* |x| < 0.25 */
+ if (ix < 0x3e200000) { /* |x| < 0x1p-29 */
+ if ((int)ax == 0)
+ return (1);
+ }
+ return (__kernel_cospi(ax));
+ }
+
+ if (ix < 0x3fe00000) /* |x| < 0.5 */
+ c = __kernel_sinpi(0.5 - ax);
+ else if (ix < 0x3fe80000){ /* |x| < 0.75 */
+ if (ax == 0.5)
+ return (0);
+ c = -__kernel_sinpi(ax - 0.5);
+ } else
+ c = -__kernel_cospi(1 - ax);
+ return (c);
+ }
+
+ if (ix < 0x43300000) { /* 1 <= |x| < 0x1p52 */
+ /* Determine integer part of ax. */
+ j0 = ((ix >> 20) & 0x7ff) - 0x3ff;
+ if (j0 < 20) {
+ ix &= ~(0x000fffff >> j0);
+ lx = 0;
+ } else {
+ lx &= ~((uint32_t)0xffffffff >> (j0 - 20));
+ }
+ INSERT_WORDS(x, ix, lx);
+
+ ax -= x;
+ EXTRACT_WORDS(ix, lx, ax);
+
+
+ if (ix < 0x3fe00000) { /* |x| < 0.5 */
+ if (ix < 0x3fd00000) /* |x| < 0.25 */
+ c = ix == 0 ? 1 : __kernel_cospi(ax);
+ else
+ c = __kernel_sinpi(0.5 - ax);
+ } else {
+ if (ix < 0x3fe80000) { /* |x| < 0.75 */
+ if (ax == 0.5)
+ return (0);
+ c = -__kernel_sinpi(ax - 0.5);
+ } else
+ c = -__kernel_cospi(1 - ax);
+ }
+
+ if (j0 > 30)
+ x -= 0x1p30;
+ j0 = (uint32_t)x;
+ return (j0 & 1 ? -c : c);
+ }
+
+ if (ix >= 0x7f800000)
+ return (vzero / vzero);
+
+ /*
+ * |x| >= 0x1p52 is always an even integer, so return 1.
+ */
+ return (1);
+}
+
+#if LDBL_MANT_DIG == 53
+__weak_reference(cospi, cospil);
+#endif
diff --git a/libm/upstream-freebsd/lib/msun/src/s_ctanh.c b/libm/upstream-freebsd/lib/msun/src/s_ctanh.c
index 88afeb5..e5840a1 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_ctanh.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_ctanh.c
@@ -111,11 +111,13 @@
}
/*
- * ctanh(x + I NaN) = d(NaN) + I d(NaN)
- * ctanh(x +- I Inf) = dNaN + I dNaN
+ * ctanh(+-0 + i NAN) = +-0 + i NaN
+ * ctanh(+-0 +- i Inf) = +-0 + i NaN
+ * ctanh(x + i NAN) = NaN + i NaN
+ * ctanh(x +- i Inf) = NaN + i NaN
*/
if (!isfinite(y))
- return (CMPLX(y - y, y - y));
+ return (CMPLX(x ? y - y : x, y - y));
/*
* ctanh(+-huge +- I y) ~= +-1 +- I 2sin(2y)/exp(2x), using the
diff --git a/libm/upstream-freebsd/lib/msun/src/s_ctanhf.c b/libm/upstream-freebsd/lib/msun/src/s_ctanhf.c
index d2bd0b6..c46f86d 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_ctanhf.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_ctanhf.c
@@ -61,7 +61,7 @@
}
if (!isfinite(y))
- return (CMPLXF(y - y, y - y));
+ return (CMPLXF(ix ? y - y : x, y - y));
if (ix >= 0x41300000) { /* |x| >= 11 */
float exp_mx = expf(-fabsf(x));
diff --git a/libm/upstream-freebsd/lib/msun/src/s_fma.c b/libm/upstream-freebsd/lib/msun/src/s_fma.c
index 41a6424..95cffd0 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_fma.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_fma.c
@@ -35,6 +35,13 @@
#include "math_private.h"
+#ifdef USE_BUILTIN_FMA
+double
+fma(double x, double y, double z)
+{
+ return (__builtin_fma(x, y, z));
+}
+#else
/*
* A struct dd represents a floating-point number with twice the precision
* of a double. We maintain the invariant that "hi" stores the 53 high-order
@@ -284,6 +291,7 @@
else
return (add_and_denormalize(r.hi, adj, spread));
}
+#endif /* !USE_BUILTIN_FMA */
#if (LDBL_MANT_DIG == 53)
__weak_reference(fma, fmal);
diff --git a/libm/upstream-freebsd/lib/msun/src/s_fmaf.c b/libm/upstream-freebsd/lib/msun/src/s_fmaf.c
index 389cf1b..4591cc2 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_fmaf.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_fmaf.c
@@ -34,6 +34,13 @@
#include "math.h"
#include "math_private.h"
+#ifdef USE_BUILTIN_FMAF
+float
+fmaf(float x, float y, float z)
+{
+ return (__builtin_fmaf(x, y, z));
+}
+#else
/*
* Fused multiply-add: Compute x * y + z with a single rounding error.
*
@@ -69,3 +76,4 @@
SET_LOW_WORD(adjusted_result, lr + 1);
return (adjusted_result);
}
+#endif /* !USE_BUILTIN_FMAF */
diff --git a/libm/upstream-freebsd/lib/msun/src/s_fmax.c b/libm/upstream-freebsd/lib/msun/src/s_fmax.c
index 0c234bc..b53b1e6 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_fmax.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_fmax.c
@@ -34,6 +34,13 @@
#include "fpmath.h"
+#ifdef USE_BUILTIN_FMAX
+double
+fmax(double x, double y)
+{
+ return (__builtin_fmax(x, y));
+}
+#else
double
fmax(double x, double y)
{
@@ -54,6 +61,7 @@
return (x > y ? x : y);
}
+#endif
#if (LDBL_MANT_DIG == 53)
__weak_reference(fmax, fmaxl);
diff --git a/libm/upstream-freebsd/lib/msun/src/s_fmaxf.c b/libm/upstream-freebsd/lib/msun/src/s_fmaxf.c
index 8e9d1ba..8d3d14f 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_fmaxf.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_fmaxf.c
@@ -33,6 +33,13 @@
#include "fpmath.h"
+#ifdef USE_BUILTIN_FMAXF
+float
+fmaxf(float x, float y)
+{
+ return (__builtin_fmaxf(x, y));
+}
+#else
float
fmaxf(float x, float y)
{
@@ -53,3 +60,4 @@
return (x > y ? x : y);
}
+#endif
diff --git a/libm/upstream-freebsd/lib/msun/src/s_fmin.c b/libm/upstream-freebsd/lib/msun/src/s_fmin.c
index d7f24c1..53f36c1 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_fmin.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_fmin.c
@@ -34,6 +34,13 @@
#include "fpmath.h"
+#ifdef USE_BUILTIN_FMIN
+double
+fmin(double x, double y)
+{
+ return (__builtin_fmin(x, y));
+}
+#else
double
fmin(double x, double y)
{
@@ -54,6 +61,7 @@
return (x < y ? x : y);
}
+#endif
#if (LDBL_MANT_DIG == 53)
__weak_reference(fmin, fminl);
diff --git a/libm/upstream-freebsd/lib/msun/src/s_fminf.c b/libm/upstream-freebsd/lib/msun/src/s_fminf.c
index 2583167..58b6a48 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_fminf.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_fminf.c
@@ -33,6 +33,13 @@
#include "fpmath.h"
+#ifdef USE_BUILTIN_FMINF
+float
+fminf(float x, float y)
+{
+ return (__builtin_fminf(x, y));
+}
+#else
float
fminf(float x, float y)
{
@@ -53,3 +60,4 @@
return (x < y ? x : y);
}
+#endif
diff --git a/libm/upstream-freebsd/lib/msun/src/s_lround.c b/libm/upstream-freebsd/lib/msun/src/s_lround.c
index 66d9183..1dd8e69 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_lround.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_lround.c
@@ -49,9 +49,11 @@
* that everything is in range. At compile time, INRANGE(x) should reduce to
* two floating-point comparisons in the former case, or TRUE otherwise.
*/
+static const type type_min = (type)DTYPE_MIN;
+static const type type_max = (type)DTYPE_MAX;
static const type dtype_min = (type)DTYPE_MIN - 0.5;
static const type dtype_max = (type)DTYPE_MAX + 0.5;
-#define INRANGE(x) (dtype_max - (type)DTYPE_MAX != 0.5 || \
+#define INRANGE(x) (dtype_max - type_max != 0.5 || \
((x) > dtype_min && (x) < dtype_max))
dtype
diff --git a/libm/upstream-freebsd/lib/msun/src/s_scalbn.c b/libm/upstream-freebsd/lib/msun/src/s_scalbn.c
index b048b05..2d4f7a3 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_scalbn.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_scalbn.c
@@ -1,66 +1,47 @@
-/* @(#)s_scalbn.c 5.1 93/09/24 */
/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ * Copyright (c) 2005-2020 Rich Felker, et al.
*
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
+ * SPDX-License-Identifier: MIT
+ *
+ * Please see https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
+ * for all contributors to musl.
*/
-
-#include <sys/cdefs.h>
-__FBSDID("$FreeBSD$");
-
-/*
- * scalbn (double x, int n)
- * scalbn(x,n) returns x* 2**n computed by exponent
- * manipulation rather than by actually performing an
- * exponentiation or a multiplication.
- */
-
#include <float.h>
+#include <math.h>
+#include <stdint.h>
-#include "math.h"
-#include "math_private.h"
-
-static const double
-two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
-twom54 = 5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */
-huge = 1.0e+300,
-tiny = 1.0e-300;
-
-double
-scalbn (double x, int n)
+double scalbn(double x, int n)
{
- int32_t k,hx,lx;
- EXTRACT_WORDS(hx,lx,x);
- k = (hx&0x7ff00000)>>20; /* extract exponent */
- if (k==0) { /* 0 or subnormal x */
- if ((lx|(hx&0x7fffffff))==0) return x; /* +-0 */
- x *= two54;
- GET_HIGH_WORD(hx,x);
- k = ((hx&0x7ff00000)>>20) - 54;
- if (n< -50000) return tiny*x; /*underflow*/
- }
- if (k==0x7ff) return x+x; /* NaN or Inf */
- k = k+n;
- if (k > 0x7fe) return huge*copysign(huge,x); /* overflow */
- if (k > 0) /* normal result */
- {SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20)); return x;}
- if (k <= -54) {
- if (n > 50000) /* in case integer overflow in n+k */
- return huge*copysign(huge,x); /*overflow*/
- else
- return tiny*copysign(tiny,x); /*underflow*/
+ union {double f; uint64_t i;} u;
+ double_t y = x;
+
+ if (n > 1023) {
+ y *= 0x1p1023;
+ n -= 1023;
+ if (n > 1023) {
+ y *= 0x1p1023;
+ n -= 1023;
+ if (n > 1023)
+ n = 1023;
+ }
+ } else if (n < -1022) {
+ /* make sure final n < -53 to avoid double
+ rounding in the subnormal range */
+ y *= 0x1p-1022 * 0x1p53;
+ n += 1022 - 53;
+ if (n < -1022) {
+ y *= 0x1p-1022 * 0x1p53;
+ n += 1022 - 53;
+ if (n < -1022)
+ n = -1022;
+ }
}
- k += 54; /* subnormal result */
- SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20));
- return x*twom54;
+ u.i = (uint64_t)(0x3ff+n)<<52;
+ x = y * u.f;
+ return x;
}
-#if (LDBL_MANT_DIG == 53)
+#if (LDBL_MANT_DIG == 53) && !defined(scalbn)
__weak_reference(scalbn, ldexpl);
__weak_reference(scalbn, scalbnl);
#endif
diff --git a/libm/upstream-freebsd/lib/msun/src/s_scalbnf.c b/libm/upstream-freebsd/lib/msun/src/s_scalbnf.c
index 21d001c..8cf1e01 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_scalbnf.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_scalbnf.c
@@ -1,57 +1,41 @@
-/* s_scalbnf.c -- float version of s_scalbn.c.
- * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
- */
-
/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ * Copyright (c) 2005-2020 Rich Felker, et al.
*
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
+ * SPDX-License-Identifier: MIT
+ *
+ * Please see https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
+ * for all contributors to musl.
*/
+#include <math.h>
+#include <stdint.h>
-#include <sys/cdefs.h>
-__FBSDID("$FreeBSD$");
-
-#include "math.h"
-#include "math_private.h"
-
-static const float
-two25 = 3.355443200e+07, /* 0x4c000000 */
-twom25 = 2.9802322388e-08, /* 0x33000000 */
-huge = 1.0e+30,
-tiny = 1.0e-30;
-
-float
-scalbnf (float x, int n)
+float scalbnf(float x, int n)
{
- int32_t k,ix;
- GET_FLOAT_WORD(ix,x);
- k = (ix&0x7f800000)>>23; /* extract exponent */
- if (k==0) { /* 0 or subnormal x */
- if ((ix&0x7fffffff)==0) return x; /* +-0 */
- x *= two25;
- GET_FLOAT_WORD(ix,x);
- k = ((ix&0x7f800000)>>23) - 25;
- if (n< -50000) return tiny*x; /*underflow*/
- }
- if (k==0xff) return x+x; /* NaN or Inf */
- k = k+n;
- if (k > 0xfe) return huge*copysignf(huge,x); /* overflow */
- if (k > 0) /* normal result */
- {SET_FLOAT_WORD(x,(ix&0x807fffff)|(k<<23)); return x;}
- if (k <= -25) {
- if (n > 50000) /* in case integer overflow in n+k */
- return huge*copysignf(huge,x); /*overflow*/
- else
- return tiny*copysignf(tiny,x); /*underflow*/
+ union {float f; uint32_t i;} u;
+ float_t y = x;
+
+ if (n > 127) {
+ y *= 0x1p127f;
+ n -= 127;
+ if (n > 127) {
+ y *= 0x1p127f;
+ n -= 127;
+ if (n > 127)
+ n = 127;
+ }
+ } else if (n < -126) {
+ y *= 0x1p-126f * 0x1p24f;
+ n += 126 - 24;
+ if (n < -126) {
+ y *= 0x1p-126f * 0x1p24f;
+ n += 126 - 24;
+ if (n < -126)
+ n = -126;
+ }
}
- k += 25; /* subnormal result */
- SET_FLOAT_WORD(x,(ix&0x807fffff)|(k<<23));
- return x*twom25;
+ u.i = (uint32_t)(0x7f+n)<<23;
+ x = y * u.f;
+ return x;
}
__strong_reference(scalbnf, ldexpf);
diff --git a/libm/upstream-freebsd/lib/msun/src/s_scalbnl.c b/libm/upstream-freebsd/lib/msun/src/s_scalbnl.c
index 28b0cf9..6044c1b 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_scalbnl.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_scalbnl.c
@@ -1,71 +1,49 @@
-/* @(#)s_scalbn.c 5.1 93/09/24 */
/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ * Copyright (c) 2005-2020 Rich Felker, et al.
*
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
+ * SPDX-License-Identifier: MIT
+ *
+ * Please see https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
+ * for all contributors to musl.
*/
-
-#include <sys/cdefs.h>
-__FBSDID("$FreeBSD$");
-
+#include <math.h>
+#include <float.h>
+#include "math_private.h"
+#include "fpmath.h"
/*
* scalbnl (long double x, int n)
* scalbnl(x,n) returns x* 2**n computed by exponent
* manipulation rather than by actually performing an
* exponentiation or a multiplication.
*/
-
-/*
- * We assume that a long double has a 15-bit exponent. On systems
- * where long double is the same as double, scalbnl() is an alias
- * for scalbn(), so we don't use this routine.
- */
-
-#include <float.h>
-#include <math.h>
-
-#include "fpmath.h"
-
-#if LDBL_MAX_EXP != 0x4000
-#error "Unsupported long double format"
-#endif
-
-static const long double
-huge = 0x1p16000L,
-tiny = 0x1p-16000L;
-
-long double
-scalbnl (long double x, int n)
+#if (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
+long double scalbnl(long double x, int n)
{
union IEEEl2bits u;
- int k;
- u.e = x;
- k = u.bits.exp; /* extract exponent */
- if (k==0) { /* 0 or subnormal x */
- if ((u.bits.manh|u.bits.manl)==0) return x; /* +-0 */
- u.e *= 0x1p+128;
- k = u.bits.exp - 128;
- if (n< -50000) return tiny*x; /*underflow*/
- }
- if (k==0x7fff) return x+x; /* NaN or Inf */
- k = k+n;
- if (k >= 0x7fff) return huge*copysignl(huge,x); /* overflow */
- if (k > 0) /* normal result */
- {u.bits.exp = k; return u.e;}
- if (k <= -128) {
- if (n > 50000) /* in case integer overflow in n+k */
- return huge*copysign(huge,x); /*overflow*/
- else
- return tiny*copysign(tiny,x); /*underflow*/
- }
- k += 128; /* subnormal result */
- u.bits.exp = k;
- return u.e*0x1p-128;
-}
+ if (n > 16383) {
+ x *= 0x1p16383L;
+ n -= 16383;
+ if (n > 16383) {
+ x *= 0x1p16383L;
+ n -= 16383;
+ if (n > 16383)
+ n = 16383;
+ }
+ } else if (n < -16382) {
+ x *= 0x1p-16382L * 0x1p113L;
+ n += 16382 - 113;
+ if (n < -16382) {
+ x *= 0x1p-16382L * 0x1p113L;
+ n += 16382 - 113;
+ if (n < -16382)
+ n = -16382;
+ }
+ }
+ u.e = 1.0;
+ u.xbits.expsign = 0x3fff + n;
+ return x * u.e;
+}
__strong_reference(scalbnl, ldexpl);
+#endif
+
diff --git a/libm/upstream-freebsd/lib/msun/src/s_sincosl.c b/libm/upstream-freebsd/lib/msun/src/s_sincosl.c
index aef36c2..3dd3457 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_sincosl.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_sincosl.c
@@ -50,11 +50,10 @@
sincosl(long double x, long double *sn, long double *cs)
{
union IEEEl2bits z;
- int e0, sgn;
+ int e0;
long double y[2];
z.e = x;
- sgn = z.bits.sign;
z.bits.sign = 0;
ENTERV();
diff --git a/libm/upstream-freebsd/lib/msun/src/s_sinpi.c b/libm/upstream-freebsd/lib/msun/src/s_sinpi.c
new file mode 100644
index 0000000..858459a
--- /dev/null
+++ b/libm/upstream-freebsd/lib/msun/src/s_sinpi.c
@@ -0,0 +1,169 @@
+/*-
+ * Copyright (c) 2017 Steven G. Kargl
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ * notice unmodified, this list of conditions, and the following
+ * disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ * notice, this list of conditions and the following disclaimer in the
+ * documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
+ * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
+ * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
+ * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
+ * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
+ * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+ * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+ * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
+ * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ */
+
+/**
+ * sinpi(x) computes sin(pi*x) without multiplication by pi (almost). First,
+ * note that sinpi(-x) = -sinpi(x), so the algorithm considers only |x| and
+ * includes reflection symmetry by considering the sign of x on output. The
+ * method used depends on the magnitude of x.
+ *
+ * 1. For small |x|, sinpi(x) = pi * x where a sloppy threshold is used. The
+ * threshold is |x| < 0x1pN with N = -(P/2+M). P is the precision of the
+ * floating-point type and M = 2 to 4. To achieve high accuracy, pi is
+ * decomposed into high and low parts with the high part containing a
+ * number of trailing zero bits. x is also split into high and low parts.
+ *
+ * 2. For |x| < 1, argument reduction is not required and sinpi(x) is
+ * computed by calling a kernel that leverages the kernels for sin(x)
+ * ans cos(x). See k_sinpi.c and k_cospi.c for details.
+ *
+ * 3. For 1 <= |x| < 0x1p(P-1), argument reduction is required where
+ * |x| = j0 + r with j0 an integer and the remainder r satisfies
+ * 0 <= r < 1. With the given domain, a simplified inline floor(x)
+ * is used. Also, note the following identity
+ *
+ * sinpi(x) = sin(pi*(j0+r))
+ * = sin(pi*j0) * cos(pi*r) + cos(pi*j0) * sin(pi*r)
+ * = cos(pi*j0) * sin(pi*r)
+ * = +-sinpi(r)
+ *
+ * If j0 is even, then cos(pi*j0) = 1. If j0 is odd, then cos(pi*j0) = -1.
+ * sinpi(r) is then computed via an appropriate kernel.
+ *
+ * 4. For |x| >= 0x1p(P-1), |x| is integral and sinpi(x) = copysign(0,x).
+ *
+ * 5. Special cases:
+ *
+ * sinpi(+-0) = +-0
+ * sinpi(+-n) = +-0, for positive integers n.
+ * sinpi(+-inf) = nan. Raises the "invalid" floating-point exception.
+ * sinpi(nan) = nan. Raises the "invalid" floating-point exception.
+ */
+
+#include <float.h>
+#include "math.h"
+#include "math_private.h"
+
+static const double
+pi_hi = 3.1415926814079285e+00, /* 0x400921fb 0x58000000 */
+pi_lo =-2.7818135228334233e-08; /* 0xbe5dde97 0x3dcb3b3a */
+
+#include "k_cospi.h"
+#include "k_sinpi.h"
+
+volatile static const double vzero = 0;
+
+double
+sinpi(double x)
+{
+ double ax, hi, lo, s;
+ uint32_t hx, ix, j0, lx;
+
+ EXTRACT_WORDS(hx, lx, x);
+ ix = hx & 0x7fffffff;
+ INSERT_WORDS(ax, ix, lx);
+
+ if (ix < 0x3ff00000) { /* |x| < 1 */
+ if (ix < 0x3fd00000) { /* |x| < 0.25 */
+ if (ix < 0x3e200000) { /* |x| < 0x1p-29 */
+ if (x == 0)
+ return (x);
+ /*
+ * To avoid issues with subnormal values,
+ * scale the computation and rescale on
+ * return.
+ */
+ INSERT_WORDS(hi, hx, 0);
+ hi *= 0x1p53;
+ lo = x * 0x1p53 - hi;
+ s = (pi_lo + pi_hi) * lo + pi_lo * hi +
+ pi_hi * hi;
+ return (s * 0x1p-53);
+ }
+
+ s = __kernel_sinpi(ax);
+ return ((hx & 0x80000000) ? -s : s);
+ }
+
+ if (ix < 0x3fe00000) /* |x| < 0.5 */
+ s = __kernel_cospi(0.5 - ax);
+ else if (ix < 0x3fe80000) /* |x| < 0.75 */
+ s = __kernel_cospi(ax - 0.5);
+ else
+ s = __kernel_sinpi(1 - ax);
+ return ((hx & 0x80000000) ? -s : s);
+ }
+
+ if (ix < 0x43300000) { /* 1 <= |x| < 0x1p52 */
+ /* Determine integer part of ax. */
+ j0 = ((ix >> 20) & 0x7ff) - 0x3ff;
+ if (j0 < 20) {
+ ix &= ~(0x000fffff >> j0);
+ lx = 0;
+ } else {
+ lx &= ~((uint32_t)0xffffffff >> (j0 - 20));
+ }
+ INSERT_WORDS(x, ix, lx);
+
+ ax -= x;
+ EXTRACT_WORDS(ix, lx, ax);
+
+ if (ix == 0)
+ s = 0;
+ else {
+ if (ix < 0x3fe00000) { /* |x| < 0.5 */
+ if (ix < 0x3fd00000) /* |x| < 0.25 */
+ s = __kernel_sinpi(ax);
+ else
+ s = __kernel_cospi(0.5 - ax);
+ } else {
+ if (ix < 0x3fe80000) /* |x| < 0.75 */
+ s = __kernel_cospi(ax - 0.5);
+ else
+ s = __kernel_sinpi(1 - ax);
+ }
+
+ if (j0 > 30)
+ x -= 0x1p30;
+ j0 = (uint32_t)x;
+ if (j0 & 1) s = -s;
+ }
+
+ return ((hx & 0x80000000) ? -s : s);
+ }
+
+ if (ix >= 0x7f800000)
+ return (vzero / vzero);
+
+ /*
+ * |x| >= 0x1p52 is always an integer, so return +-0.
+ */
+ return (copysign(0, x));
+}
+
+#if LDBL_MANT_DIG == 53
+__weak_reference(sinpi, sinpil);
+#endif