summaryrefslogtreecommitdiffstats
path: root/lib/msun/src
diff options
context:
space:
mode:
Diffstat (limited to 'lib/msun/src')
-rw-r--r--lib/msun/src/catrig.c9
-rw-r--r--lib/msun/src/catrigl.c417
-rw-r--r--lib/msun/src/e_asin.c2
-rw-r--r--lib/msun/src/e_coshl.c2
-rw-r--r--lib/msun/src/e_lgammaf_r.c2
-rw-r--r--lib/msun/src/e_sinhl.c5
-rw-r--r--lib/msun/src/k_sincos.h52
-rw-r--r--lib/msun/src/k_sincosf.h43
-rw-r--r--lib/msun/src/k_sincosl.h134
-rw-r--r--lib/msun/src/math.h3
-rw-r--r--lib/msun/src/math_private.h14
-rw-r--r--lib/msun/src/s_fabs.c5
-rw-r--r--lib/msun/src/s_fmax.c5
-rw-r--r--lib/msun/src/s_fmin.c5
-rw-r--r--lib/msun/src/s_logbl.c5
-rw-r--r--lib/msun/src/s_scalbn.c12
-rw-r--r--lib/msun/src/s_scalbnf.c11
-rw-r--r--lib/msun/src/s_scalbnl.c12
-rw-r--r--lib/msun/src/s_sincos.c80
-rw-r--r--lib/msun/src/s_sincosf.c126
-rw-r--r--lib/msun/src/s_sincosl.c105
-rw-r--r--lib/msun/src/s_tanhl.c2
22 files changed, 1023 insertions, 28 deletions
diff --git a/lib/msun/src/catrig.c b/lib/msun/src/catrig.c
index f392862..ebfe356 100644
--- a/lib/msun/src/catrig.c
+++ b/lib/msun/src/catrig.c
@@ -637,3 +637,12 @@ catan(double complex z)
return (CMPLX(cimag(w), creal(w)));
}
+
+#if LDBL_MANT_DIG == 53
+__weak_reference(cacosh, cacoshl);
+__weak_reference(cacos, cacosl);
+__weak_reference(casinh, casinhl);
+__weak_reference(casin, casinl);
+__weak_reference(catanh, catanhl);
+__weak_reference(catan, catanl);
+#endif
diff --git a/lib/msun/src/catrigl.c b/lib/msun/src/catrigl.c
new file mode 100644
index 0000000..b798570
--- /dev/null
+++ b/lib/msun/src/catrigl.c
@@ -0,0 +1,417 @@
+/*-
+ * Copyright (c) 2012 Stephen Montgomery-Smith <stephen@FreeBSD.ORG>
+ * Copyright (c) 2017 Mahdi Mokhtari <mmokhi@FreeBSD.org>
+ * 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, 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 AND CONTRIBUTORS ``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 OR CONTRIBUTORS 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 algorithm is very close to that in "Implementing the complex arcsine
+ * and arccosine functions using exception handling" by T. E. Hull, Thomas F.
+ * Fairgrieve, and Ping Tak Peter Tang, published in ACM Transactions on
+ * Mathematical Software, Volume 23 Issue 3, 1997, Pages 299-335,
+ * http://dl.acm.org/citation.cfm?id=275324.
+ *
+ * See catrig.c for complete comments.
+ *
+ * XXX comments were removed automatically, and even short ones on the right
+ * of statements were removed (all of them), contrary to normal style. Only
+ * a few comments on the right of declarations remain.
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+#include <complex.h>
+#include <float.h>
+
+#include "invtrig.h"
+#include "math.h"
+#include "math_private.h"
+
+#undef isinf
+#define isinf(x) (fabsl(x) == INFINITY)
+#undef isnan
+#define isnan(x) ((x) != (x))
+#define raise_inexact() do { volatile float junk __unused = 1 + tiny; } while(0)
+#undef signbit
+#define signbit(x) (__builtin_signbitl(x))
+
+#if LDBL_MAX_EXP != 0x4000
+#error "Unsupported long double format"
+#endif
+
+static const long double
+A_crossover = 10,
+B_crossover = 0.6417,
+FOUR_SQRT_MIN = 0x1p-8189L,
+HALF_MAX = 0x1p16383L,
+QUARTER_SQRT_MAX = 0x1p8189L,
+RECIP_EPSILON = 1 / LDBL_EPSILON,
+SQRT_MIN = 0x1p-8191L;
+
+#if LDBL_MANT_DIG == 64
+static const union IEEEl2bits
+um_e = LD80C(0xadf85458a2bb4a9b, 1, 2.71828182845904523536e+0L),
+um_ln2 = LD80C(0xb17217f7d1cf79ac, -1, 6.93147180559945309417e-1L);
+#define m_e um_e.e
+#define m_ln2 um_ln2.e
+static const long double
+/* The next 2 literals for non-i386. Misrounding them on i386 is harmless. */
+SQRT_3_EPSILON = 5.70316273435758915310e-10, /* 0x9cc470a0490973e8.0p-94 */
+SQRT_6_EPSILON = 8.06549008734932771664e-10; /* 0xddb3d742c265539e.0p-94 */
+#elif LDBL_MANT_DIG == 113
+static const long double
+m_e = 2.71828182845904523536028747135266250e0L, /* 0x15bf0a8b1457695355fb8ac404e7a.0p-111 */
+m_ln2 = 6.93147180559945309417232121458176568e-1L, /* 0x162e42fefa39ef35793c7673007e6.0p-113 */
+SQRT_3_EPSILON = 2.40370335797945490975336727199878124e-17, /* 0x1bb67ae8584caa73b25742d7078b8.0p-168 */
+SQRT_6_EPSILON = 3.39934988877629587239082586223300391e-17; /* 0x13988e1409212e7d0321914321a55.0p-167 */
+#else
+#error "Unsupported long double format"
+#endif
+
+static const volatile float
+tiny = 0x1p-100;
+
+static long double complex clog_for_large_values(long double complex z);
+
+static inline long double
+f(long double a, long double b, long double hypot_a_b)
+{
+ if (b < 0)
+ return ((hypot_a_b - b) / 2);
+ if (b == 0)
+ return (a / 2);
+ return (a * a / (hypot_a_b + b) / 2);
+}
+
+static inline void
+do_hard_work(long double x, long double y, long double *rx, int *B_is_usable,
+ long double *B, long double *sqrt_A2my2, long double *new_y)
+{
+ long double R, S, A;
+ long double Am1, Amy;
+
+ R = hypotl(x, y + 1);
+ S = hypotl(x, y - 1);
+
+ A = (R + S) / 2;
+ if (A < 1)
+ A = 1;
+
+ if (A < A_crossover) {
+ if (y == 1 && x < LDBL_EPSILON * LDBL_EPSILON / 128) {
+ *rx = sqrtl(x);
+ } else if (x >= LDBL_EPSILON * fabsl(y - 1)) {
+ Am1 = f(x, 1 + y, R) + f(x, 1 - y, S);
+ *rx = log1pl(Am1 + sqrtl(Am1 * (A + 1)));
+ } else if (y < 1) {
+ *rx = x / sqrtl((1 - y) * (1 + y));
+ } else {
+ *rx = log1pl((y - 1) + sqrtl((y - 1) * (y + 1)));
+ }
+ } else {
+ *rx = logl(A + sqrtl(A * A - 1));
+ }
+
+ *new_y = y;
+
+ if (y < FOUR_SQRT_MIN) {
+ *B_is_usable = 0;
+ *sqrt_A2my2 = A * (2 / LDBL_EPSILON);
+ *new_y = y * (2 / LDBL_EPSILON);
+ return;
+ }
+
+ *B = y / A;
+ *B_is_usable = 1;
+
+ if (*B > B_crossover) {
+ *B_is_usable = 0;
+ if (y == 1 && x < LDBL_EPSILON / 128) {
+ *sqrt_A2my2 = sqrtl(x) * sqrtl((A + y) / 2);
+ } else if (x >= LDBL_EPSILON * fabsl(y - 1)) {
+ Amy = f(x, y + 1, R) + f(x, y - 1, S);
+ *sqrt_A2my2 = sqrtl(Amy * (A + y));
+ } else if (y > 1) {
+ *sqrt_A2my2 = x * (4 / LDBL_EPSILON / LDBL_EPSILON) * y /
+ sqrtl((y + 1) * (y - 1));
+ *new_y = y * (4 / LDBL_EPSILON / LDBL_EPSILON);
+ } else {
+ *sqrt_A2my2 = sqrtl((1 - y) * (1 + y));
+ }
+ }
+}
+
+long double complex
+casinhl(long double complex z)
+{
+ long double x, y, ax, ay, rx, ry, B, sqrt_A2my2, new_y;
+ int B_is_usable;
+ long double complex w;
+
+ x = creall(z);
+ y = cimagl(z);
+ ax = fabsl(x);
+ ay = fabsl(y);
+
+ if (isnan(x) || isnan(y)) {
+ if (isinf(x))
+ return (CMPLXL(x, y + y));
+ if (isinf(y))
+ return (CMPLXL(y, x + x));
+ if (y == 0)
+ return (CMPLXL(x + x, y));
+ return (CMPLXL(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
+ }
+
+ if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
+ if (signbit(x) == 0)
+ w = clog_for_large_values(z) + m_ln2;
+ else
+ w = clog_for_large_values(-z) + m_ln2;
+ return (CMPLXL(copysignl(creall(w), x),
+ copysignl(cimagl(w), y)));
+ }
+
+ if (x == 0 && y == 0)
+ return (z);
+
+ raise_inexact();
+
+ if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4)
+ return (z);
+
+ do_hard_work(ax, ay, &rx, &B_is_usable, &B, &sqrt_A2my2, &new_y);
+ if (B_is_usable)
+ ry = asinl(B);
+ else
+ ry = atan2l(new_y, sqrt_A2my2);
+ return (CMPLXL(copysignl(rx, x), copysignl(ry, y)));
+}
+
+long double complex
+casinl(long double complex z)
+{
+ long double complex w;
+
+ w = casinhl(CMPLXL(cimagl(z), creall(z)));
+ return (CMPLXL(cimagl(w), creall(w)));
+}
+
+long double complex
+cacosl(long double complex z)
+{
+ long double x, y, ax, ay, rx, ry, B, sqrt_A2mx2, new_x;
+ int sx, sy;
+ int B_is_usable;
+ long double complex w;
+
+ x = creall(z);
+ y = cimagl(z);
+ sx = signbit(x);
+ sy = signbit(y);
+ ax = fabsl(x);
+ ay = fabsl(y);
+
+ if (isnan(x) || isnan(y)) {
+ if (isinf(x))
+ return (CMPLXL(y + y, -INFINITY));
+ if (isinf(y))
+ return (CMPLXL(x + x, -y));
+ if (x == 0)
+ return (CMPLXL(pio2_hi + pio2_lo, y + y));
+ return (CMPLXL(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
+ }
+
+ if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
+ w = clog_for_large_values(z);
+ rx = fabsl(cimagl(w));
+ ry = creall(w) + m_ln2;
+ if (sy == 0)
+ ry = -ry;
+ return (CMPLXL(rx, ry));
+ }
+
+ if (x == 1 && y == 0)
+ return (CMPLXL(0, -y));
+
+ raise_inexact();
+
+ if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4)
+ return (CMPLXL(pio2_hi - (x - pio2_lo), -y));
+
+ do_hard_work(ay, ax, &ry, &B_is_usable, &B, &sqrt_A2mx2, &new_x);
+ if (B_is_usable) {
+ if (sx == 0)
+ rx = acosl(B);
+ else
+ rx = acosl(-B);
+ } else {
+ if (sx == 0)
+ rx = atan2l(sqrt_A2mx2, new_x);
+ else
+ rx = atan2l(sqrt_A2mx2, -new_x);
+ }
+ if (sy == 0)
+ ry = -ry;
+ return (CMPLXL(rx, ry));
+}
+
+long double complex
+cacoshl(long double complex z)
+{
+ long double complex w;
+ long double rx, ry;
+
+ w = cacosl(z);
+ rx = creall(w);
+ ry = cimagl(w);
+ if (isnan(rx) && isnan(ry))
+ return (CMPLXL(ry, rx));
+ if (isnan(rx))
+ return (CMPLXL(fabsl(ry), rx));
+ if (isnan(ry))
+ return (CMPLXL(ry, ry));
+ return (CMPLXL(fabsl(ry), copysignl(rx, cimagl(z))));
+}
+
+static long double complex
+clog_for_large_values(long double complex z)
+{
+ long double x, y;
+ long double ax, ay, t;
+
+ x = creall(z);
+ y = cimagl(z);
+ ax = fabsl(x);
+ ay = fabsl(y);
+ if (ax < ay) {
+ t = ax;
+ ax = ay;
+ ay = t;
+ }
+
+ if (ax > HALF_MAX)
+ return (CMPLXL(logl(hypotl(x / m_e, y / m_e)) + 1,
+ atan2l(y, x)));
+
+ if (ax > QUARTER_SQRT_MAX || ay < SQRT_MIN)
+ return (CMPLXL(logl(hypotl(x, y)), atan2l(y, x)));
+
+ return (CMPLXL(logl(ax * ax + ay * ay) / 2, atan2l(y, x)));
+}
+
+static inline long double
+sum_squares(long double x, long double y)
+{
+
+ if (y < SQRT_MIN)
+ return (x * x);
+
+ return (x * x + y * y);
+}
+
+static inline long double
+real_part_reciprocal(long double x, long double y)
+{
+ long double scale;
+ uint16_t hx, hy;
+ int16_t ix, iy;
+
+ GET_LDBL_EXPSIGN(hx, x);
+ ix = hx & 0x7fff;
+ GET_LDBL_EXPSIGN(hy, y);
+ iy = hy & 0x7fff;
+#define BIAS (LDBL_MAX_EXP - 1)
+#define CUTOFF (LDBL_MANT_DIG / 2 + 1)
+ if (ix - iy >= CUTOFF || isinf(x))
+ return (1 / x);
+ if (iy - ix >= CUTOFF)
+ return (x / y / y);
+ if (ix <= BIAS + LDBL_MAX_EXP / 2 - CUTOFF)
+ return (x / (x * x + y * y));
+ scale = 1;
+ SET_LDBL_EXPSIGN(scale, 0x7fff - ix);
+ x *= scale;
+ y *= scale;
+ return (x / (x * x + y * y) * scale);
+}
+
+long double complex
+catanhl(long double complex z)
+{
+ long double x, y, ax, ay, rx, ry;
+
+ x = creall(z);
+ y = cimagl(z);
+ ax = fabsl(x);
+ ay = fabsl(y);
+
+ if (y == 0 && ax <= 1)
+ return (CMPLXL(atanhl(x), y));
+
+ if (x == 0)
+ return (CMPLXL(x, atanl(y)));
+
+ if (isnan(x) || isnan(y)) {
+ if (isinf(x))
+ return (CMPLXL(copysignl(0, x), y + y));
+ if (isinf(y))
+ return (CMPLXL(copysignl(0, x),
+ copysignl(pio2_hi + pio2_lo, y)));
+ return (CMPLXL(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
+ }
+
+ if (ax > RECIP_EPSILON || ay > RECIP_EPSILON)
+ return (CMPLXL(real_part_reciprocal(x, y),
+ copysignl(pio2_hi + pio2_lo, y)));
+
+ if (ax < SQRT_3_EPSILON / 2 && ay < SQRT_3_EPSILON / 2) {
+ raise_inexact();
+ return (z);
+ }
+
+ if (ax == 1 && ay < LDBL_EPSILON)
+ rx = (m_ln2 - logl(ay)) / 2;
+ else
+ rx = log1pl(4 * ax / sum_squares(ax - 1, ay)) / 4;
+
+ if (ax == 1)
+ ry = atan2l(2, -ay) / 2;
+ else if (ay < LDBL_EPSILON)
+ ry = atan2l(2 * ay, (1 - ax) * (1 + ax)) / 2;
+ else
+ ry = atan2l(2 * ay, (1 - ax) * (1 + ax) - ay * ay) / 2;
+
+ return (CMPLXL(copysignl(rx, x), copysignl(ry, y)));
+}
+
+long double complex
+catanl(long double complex z)
+{
+ long double complex w;
+
+ w = catanhl(CMPLXL(cimagl(z), creall(z)));
+ return (CMPLXL(cimagl(w), creall(w)));
+}
diff --git a/lib/msun/src/e_asin.c b/lib/msun/src/e_asin.c
index 27de207..931b270 100644
--- a/lib/msun/src/e_asin.c
+++ b/lib/msun/src/e_asin.c
@@ -6,7 +6,7 @@
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
+ * software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
diff --git a/lib/msun/src/e_coshl.c b/lib/msun/src/e_coshl.c
index 0a21277..4e3b283 100644
--- a/lib/msun/src/e_coshl.c
+++ b/lib/msun/src/e_coshl.c
@@ -86,7 +86,9 @@ long double
coshl(long double x)
{
long double hi,lo,x2,x4;
+#if LDBL_MANT_DIG == 113
double dx2;
+#endif
uint16_t ix;
GET_LDBL_EXPSIGN(ix,x);
diff --git a/lib/msun/src/e_lgammaf_r.c b/lib/msun/src/e_lgammaf_r.c
index 9084e18..48346c33 100644
--- a/lib/msun/src/e_lgammaf_r.c
+++ b/lib/msun/src/e_lgammaf_r.c
@@ -122,7 +122,7 @@ sin_pif(float x)
float
__ieee754_lgammaf_r(float x, int *signgamp)
{
- float nadj,p,p1,p2,p3,q,r,t,w,y,z;
+ float nadj,p,p1,p2,q,r,t,w,y,z;
int32_t hx;
int i,ix;
diff --git a/lib/msun/src/e_sinhl.c b/lib/msun/src/e_sinhl.c
index ce7e333..38d3df1 100644
--- a/lib/msun/src/e_sinhl.c
+++ b/lib/msun/src/e_sinhl.c
@@ -85,7 +85,10 @@ long double
sinhl(long double x)
{
long double hi,lo,x2,x4;
- double dx2,s;
+#if LDBL_MANT_DIG == 113
+ double dx2;
+#endif
+ double s;
int16_t ix,jx;
GET_LDBL_EXPSIGN(jx,x);
diff --git a/lib/msun/src/k_sincos.h b/lib/msun/src/k_sincos.h
new file mode 100644
index 0000000..6f03be2
--- /dev/null
+++ b/lib/msun/src/k_sincos.h
@@ -0,0 +1,52 @@
+/*-
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ *
+ * k_sin.c and k_cos.c merged by Steven G. Kargl.
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+static const double
+S1 = -1.66666666666666324348e-01, /* 0xBFC55555, 0x55555549 */
+S2 = 8.33333333332248946124e-03, /* 0x3F811111, 0x1110F8A6 */
+S3 = -1.98412698298579493134e-04, /* 0xBF2A01A0, 0x19C161D5 */
+S4 = 2.75573137070700676789e-06, /* 0x3EC71DE3, 0x57B1FE7D */
+S5 = -2.50507602534068634195e-08, /* 0xBE5AE5E6, 0x8A2B9CEB */
+S6 = 1.58969099521155010221e-10; /* 0x3DE5D93A, 0x5ACFD57C */
+
+static const double
+C1 = 4.16666666666666019037e-02, /* 0x3FA55555, 0x5555554C */
+C2 = -1.38888888888741095749e-03, /* 0xBF56C16C, 0x16C15177 */
+C3 = 2.48015872894767294178e-05, /* 0x3EFA01A0, 0x19CB1590 */
+C4 = -2.75573143513906633035e-07, /* 0xBE927E4F, 0x809C52AD */
+C5 = 2.08757232129817482790e-09, /* 0x3E21EE9E, 0xBDB4B1C4 */
+C6 = -1.13596475577881948265e-11; /* 0xBDA8FAE9, 0xBE8838D4 */
+
+static inline void
+__kernel_sincos(double x, double y, int iy, double *sn, double *cs)
+{
+ double hz, r, v, w, z;
+
+ z = x * x;
+ w = z * z;
+ r = S2 + z * (S3 + z * S4) + z * w * (S5 + z * S6);
+ v = z * x;
+
+ if (iy == 0)
+ *sn = x + v * (S1 + z * r);
+ else
+ *sn = x - ((z * (y / 2 - v * r) - y) - v * S1);
+
+ r = z * (C1 + z * (C2 + z * C3)) + w * w * (C4 + z * (C5 + z * C6));
+ hz = z / 2;
+ w = 1 - hz;
+ *cs = w + (((1 - w) - hz) + (z * r - x * y));
+}
diff --git a/lib/msun/src/k_sincosf.h b/lib/msun/src/k_sincosf.h
new file mode 100644
index 0000000..073986d
--- /dev/null
+++ b/lib/msun/src/k_sincosf.h
@@ -0,0 +1,43 @@
+/*-
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * 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.
+ * ====================================================
+ *
+ * k_sinf.c and k_cosf.c merged by Steven G. Kargl.
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/* |sin(x)/x - s(x)| < 2**-37.5 (~[-4.89e-12, 4.824e-12]). */
+static const double
+S1 = -0x15555554cbac77.0p-55, /* -0.166666666416265235595 */
+S2 = 0x111110896efbb2.0p-59, /* 0.0083333293858894631756 */
+S3 = -0x1a00f9e2cae774.0p-65, /* -0.000198393348360966317347 */
+S4 = 0x16cd878c3b46a7.0p-71; /* 0.0000027183114939898219064 */
+
+/* |cos(x) - c(x)| < 2**-34.1 (~[-5.37e-11, 5.295e-11]). */
+static const double
+C0 = -0x1ffffffd0c5e81.0p-54, /* -0.499999997251031003120 */
+C1 = 0x155553e1053a42.0p-57, /* 0.0416666233237390631894 */
+C2 = -0x16c087e80f1e27.0p-62, /* -0.00138867637746099294692 */
+C3 = 0x199342e0ee5069.0p-68; /* 0.0000243904487962774090654 */
+
+static inline void
+__kernel_sincosdf(double x, float *sn, float *cs)
+{
+ double r, s, w, z;
+
+ z = x * x;
+ w = z * z;
+ r = S3 + z * S4;
+ s = z * x;
+ *sn = (x + s * (S1 + z * S2)) + s * w * r;
+ r = C2 + z * C3;
+ *cs = ((1 + z * C0) + w * C1) + (w * z) * r;
+}
diff --git a/lib/msun/src/k_sincosl.h b/lib/msun/src/k_sincosl.h
new file mode 100644
index 0000000..3a609bc
--- /dev/null
+++ b/lib/msun/src/k_sincosl.h
@@ -0,0 +1,134 @@
+/*-
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ * Copyright (c) 2008 Steven G. Kargl, David Schultz, Bruce D. Evans.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ *
+ * k_sinl.c and k_cosl.c merged by Steven G. Kargl
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+#if LDBL_MANT_DIG == 64 /* ld80 version of k_sincosl.c. */
+
+#if defined(__amd64__) || defined(__i386__)
+/* Long double constants are slow on these arches, and broken on i386. */
+static const volatile double
+C1hi = 0.041666666666666664, /* 0x15555555555555.0p-57 */
+C1lo = 2.2598839032744733e-18, /* 0x14d80000000000.0p-111 */
+S1hi = -0.16666666666666666, /* -0x15555555555555.0p-55 */
+S1lo = -9.2563760475949941e-18; /* -0x15580000000000.0p-109 */
+#define S1 ((long double)S1hi + S1lo)
+#define C1 ((long double)C1hi + C1lo)
+#else
+static const long double
+C1 = 0.0416666666666666666136L; /* 0xaaaaaaaaaaaaaa9b.0p-68 */
+S1 = -0.166666666666666666671L, /* -0xaaaaaaaaaaaaaaab.0p-66 */
+#endif
+
+static const double
+C2 = -0.0013888888888888874, /* -0x16c16c16c16c10.0p-62 */
+C3 = 0.000024801587301571716, /* 0x1a01a01a018e22.0p-68 */
+C4 = -0.00000027557319215507120, /* -0x127e4fb7602f22.0p-74 */
+C5 = 0.0000000020876754400407278, /* 0x11eed8caaeccf1.0p-81 */
+C6 = -1.1470297442401303e-11, /* -0x19393412bd1529.0p-89 */
+C7 = 4.7383039476436467e-14, /* 0x1aac9d9af5c43e.0p-97 */
+S2 = 0.0083333333333333332, /* 0x11111111111111.0p-59 */
+S3 = -0.00019841269841269427, /* -0x1a01a01a019f81.0p-65 */
+S4 = 0.0000027557319223597490, /* 0x171de3a55560f7.0p-71 */
+S5 = -0.000000025052108218074604, /* -0x1ae64564f16cad.0p-78 */
+S6 = 1.6059006598854211e-10, /* 0x161242b90243b5.0p-85 */
+S7 = -7.6429779983024564e-13, /* -0x1ae42ebd1b2e00.0p-93 */
+S8 = 2.6174587166648325e-15; /* 0x179372ea0b3f64.0p-101 */
+
+static inline void
+__kernel_sincosl(long double x, long double y, int iy, long double *sn,
+ long double *cs)
+{
+ long double hz, r, v, w, z;
+
+ z = x * x;
+ v = z * x;
+ /*
+ * XXX Replace Horner scheme with an algorithm suitable for CPUs
+ * with more complex pipelines.
+ */
+ r = S2 + z * (S3 + z * (S4 + z * (S5 + z * (S6 + z * (S7 + z * S8)))));
+
+ if (iy == 0)
+ *sn = x + v * (S1 + z * r);
+ else
+ *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))))));
+ *cs = w + (((1 - w) - hz) + (z * r - x * y));
+}
+
+#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,
+S4 = 0.27557319223985890652557316053039946268333231205686e-5L,
+S5 = -0.25052108385441718775048214826384312253862930064745e-7L,
+S6 = 0.16059043836821614596571832194524392581082444805729e-9L,
+S7 = -0.76471637318198151807063387954939213287488216303768e-12L,
+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 inline void
+__kernel_sincosl(long double x, long double y, int iy, long double *sn,
+ long double *cs)
+{
+ long double hz, r, v, w, z;
+
+ z = x * x;
+ v = z * x;
+ /*
+ * XXX Replace Horner scheme with an algorithm suitable for CPUs
+ * with more complex pipelines.
+ */
+ r = S2 + z * (S3 + z * (S4 + z * (S5 + z * (S6 + z * (S7 + z * (S8 +
+ z * (S9 + z * (S10 + z * (S11 + z * S12)))))))));
+
+ if (iy == 0)
+ *sn = x + v * (S1 + z * r);
+ else
+ *cs = 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))))))))));
+
+ *cs = w + (((1 - w) - hz) + (z * r - x * y));
+}
+#else
+#error "Unsupported long double format"
+#endif
diff --git a/lib/msun/src/math.h b/lib/msun/src/math.h
index 2214c07..a8f4554 100644
--- a/lib/msun/src/math.h
+++ b/lib/msun/src/math.h
@@ -500,6 +500,9 @@ long double truncl(long double);
#if __BSD_VISIBLE
long double lgammal_r(long double, int *);
+void sincos(double, double *, double *);
+void sincosf(float, float *, float *);
+void sincosl(long double, long double *, long double *);
#endif
__END_DECLS
diff --git a/lib/msun/src/math_private.h b/lib/msun/src/math_private.h
index afaf201..7487902 100644
--- a/lib/msun/src/math_private.h
+++ b/lib/msun/src/math_private.h
@@ -306,9 +306,21 @@ do { \
fpsetprec(__oprec); \
RETURNF(__retval); \
} while (0)
+#define ENTERV() \
+ fp_prec_t __oprec; \
+ \
+ if ((__oprec = fpgetprec()) != FP_PE) \
+ fpsetprec(FP_PE)
+#define RETURNV() do { \
+ if (__oprec != FP_PE) \
+ fpsetprec(__oprec); \
+ return; \
+} while (0)
#else
-#define ENTERI(x)
+#define ENTERI()
#define RETURNI(x) RETURNF(x)
+#define ENTERV()
+#define RETURNV() return
#endif
/* Default return statement if hack*_t() is not used. */
diff --git a/lib/msun/src/s_fabs.c b/lib/msun/src/s_fabs.c
index 15529e5..48aab25 100644
--- a/lib/msun/src/s_fabs.c
+++ b/lib/msun/src/s_fabs.c
@@ -10,9 +10,8 @@
* ====================================================
*/
-#ifndef lint
-static char rcsid[] = "$FreeBSD$";
-#endif
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
/*
* fabs(x) returns the absolute value of x.
diff --git a/lib/msun/src/s_fmax.c b/lib/msun/src/s_fmax.c
index b51b865..77e27be 100644
--- a/lib/msun/src/s_fmax.c
+++ b/lib/msun/src/s_fmax.c
@@ -27,6 +27,7 @@
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
+#include <float.h>
#include <math.h>
#include "fpmath.h"
@@ -51,3 +52,7 @@ fmax(double x, double y)
return (x > y ? x : y);
}
+
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(fmax, fmaxl);
+#endif
diff --git a/lib/msun/src/s_fmin.c b/lib/msun/src/s_fmin.c
index 3500c84..5545ff2 100644
--- a/lib/msun/src/s_fmin.c
+++ b/lib/msun/src/s_fmin.c
@@ -27,6 +27,7 @@
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
+#include <float.h>
#include <math.h>
#include "fpmath.h"
@@ -51,3 +52,7 @@ fmin(double x, double y)
return (x < y ? x : y);
}
+
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(fmin, fminl);
+#endif
diff --git a/lib/msun/src/s_logbl.c b/lib/msun/src/s_logbl.c
index 7e88e36..ee1a91f 100644
--- a/lib/msun/src/s_logbl.c
+++ b/lib/msun/src/s_logbl.c
@@ -10,9 +10,8 @@
* ====================================================
*/
-#ifndef lint
-static char rcsid[] = "$FreeBSD$";
-#endif
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
#include <float.h>
#include <limits.h>
diff --git a/lib/msun/src/s_scalbn.c b/lib/msun/src/s_scalbn.c
index e7efaab..b048b05 100644
--- a/lib/msun/src/s_scalbn.c
+++ b/lib/msun/src/s_scalbn.c
@@ -10,9 +10,8 @@
* ====================================================
*/
-#ifndef lint
-static char rcsid[] = "$FreeBSD$";
-#endif
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
/*
* scalbn (double x, int n)
@@ -21,7 +20,6 @@ static char rcsid[] = "$FreeBSD$";
* exponentiation or a multiplication.
*/
-#include <sys/cdefs.h>
#include <float.h>
#include "math.h"
@@ -51,10 +49,12 @@ scalbn (double x, int 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 (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*/
+ else
+ return tiny*copysign(tiny,x); /*underflow*/
+ }
k += 54; /* subnormal result */
SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20));
return x*twom54;
diff --git a/lib/msun/src/s_scalbnf.c b/lib/msun/src/s_scalbnf.c
index 7666c74..21d001c 100644
--- a/lib/msun/src/s_scalbnf.c
+++ b/lib/msun/src/s_scalbnf.c
@@ -13,11 +13,8 @@
* ====================================================
*/
-#ifndef lint
-static char rcsid[] = "$FreeBSD$";
-#endif
-
#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
#include "math.h"
#include "math_private.h"
@@ -46,10 +43,12 @@ scalbnf (float x, int 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 (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*/
+ else
+ return tiny*copysignf(tiny,x); /*underflow*/
+ }
k += 25; /* subnormal result */
SET_FLOAT_WORD(x,(ix&0x807fffff)|(k<<23));
return x*twom25;
diff --git a/lib/msun/src/s_scalbnl.c b/lib/msun/src/s_scalbnl.c
index fc89f8d..28b0cf9 100644
--- a/lib/msun/src/s_scalbnl.c
+++ b/lib/msun/src/s_scalbnl.c
@@ -10,9 +10,8 @@
* ====================================================
*/
-#ifndef lint
-static char rcsid[] = "$FreeBSD$";
-#endif
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
/*
* scalbnl (long double x, int n)
@@ -27,7 +26,6 @@ static char rcsid[] = "$FreeBSD$";
* for scalbn(), so we don't use this routine.
*/
-#include <sys/cdefs.h>
#include <float.h>
#include <math.h>
@@ -59,10 +57,12 @@ scalbnl (long double x, int 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 (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*/
+ else
+ return tiny*copysign(tiny,x); /*underflow*/
+ }
k += 128; /* subnormal result */
u.bits.exp = k;
return u.e*0x1p-128;
diff --git a/lib/msun/src/s_sincos.c b/lib/msun/src/s_sincos.c
new file mode 100644
index 0000000..85e8d74
--- /dev/null
+++ b/lib/msun/src/s_sincos.c
@@ -0,0 +1,80 @@
+/*-
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * 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.
+ * ====================================================
+ *
+ * s_sin.c and s_cos.c merged by Steven G. Kargl. Descriptions of the
+ * algorithms are contained in the original files.
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+#include <float.h>
+
+#include "math.h"
+#define INLINE_REM_PIO2
+#include "math_private.h"
+#include "e_rem_pio2.c"
+#include "k_sincos.h"
+
+void
+sincos(double x, double *sn, double *cs)
+{
+ double y[2];
+ int32_t n, ix;
+
+ /* High word of x. */
+ GET_HIGH_WORD(ix, x);
+
+ /* |x| ~< pi/4 */
+ ix &= 0x7fffffff;
+ if (ix <= 0x3fe921fb) {
+ if (ix < 0x3e400000) { /* |x| < 2**-27 */
+ if ((int)x == 0) { /* Generate inexact. */
+ *sn = x;
+ *cs = 1;
+ return;
+ }
+ }
+ __kernel_sincos(x, 0, 0, sn, cs);
+ return;
+ }
+
+ /* If x = Inf or NaN, then sin(x) = NaN and cos(x) = NaN. */
+ if (ix >= 0x7ff00000) {
+ *sn = x - x;
+ *cs = x - x;
+ return;
+ }
+
+ /* Argument reduction. */
+ n = __ieee754_rem_pio2(x, y);
+
+ switch(n & 3) {
+ case 0:
+ __kernel_sincos(y[0], y[1], 1, sn, cs);
+ break;
+ case 1:
+ __kernel_sincos(y[0], y[1], 1, cs, sn);
+ *cs = -*cs;
+ break;
+ case 2:
+ __kernel_sincos(y[0], y[1], 1, sn, cs);
+ *sn = -*sn;
+ *cs = -*cs;
+ break;
+ default:
+ __kernel_sincos(y[0], y[1], 1, cs, sn);
+ *sn = -*sn;
+ }
+}
+
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(sincos, sincosl);
+#endif
diff --git a/lib/msun/src/s_sincosf.c b/lib/msun/src/s_sincosf.c
new file mode 100644
index 0000000..755ff05
--- /dev/null
+++ b/lib/msun/src/s_sincosf.c
@@ -0,0 +1,126 @@
+/*-
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * 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.
+ * ====================================================
+ */
+
+/* s_sincosf.c -- float version of s_sincos.c.
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ * Optimized by Bruce D. Evans.
+ * Merged s_sinf.c and s_cosf.c by Steven G. Kargl.
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+#include <float.h>
+
+#include "math.h"
+#define INLINE_REM_PIO2F
+#include "math_private.h"
+#include "e_rem_pio2f.c"
+#include "k_sincosf.h"
+
+/* Small multiples of pi/2 rounded to double precision. */
+static const double
+p1pio2 = 1*M_PI_2, /* 0x3FF921FB, 0x54442D18 */
+p2pio2 = 2*M_PI_2, /* 0x400921FB, 0x54442D18 */
+p3pio2 = 3*M_PI_2, /* 0x4012D97C, 0x7F3321D2 */
+p4pio2 = 4*M_PI_2; /* 0x401921FB, 0x54442D18 */
+
+void
+sincosf(float x, float *sn, float *cs)
+{
+ float c, s;
+ double y;
+ int32_t n, hx, ix;
+
+ GET_FLOAT_WORD(hx, x);
+ ix = hx & 0x7fffffff;
+
+ if (ix <= 0x3f490fda) { /* |x| ~<= pi/4 */
+ if (ix < 0x39800000) { /* |x| < 2**-12 */
+ if ((int)x == 0) {
+ *sn = x; /* x with inexact if x != 0 */
+ *cs = 1;
+ return;
+ }
+ }
+ __kernel_sincosdf(x, sn, cs);
+ return;
+ }
+
+ if (ix <= 0x407b53d1) { /* |x| ~<= 5*pi/4 */
+ if (ix <= 0x4016cbe3) { /* |x| ~<= 3pi/4 */
+ if (hx > 0) {
+ __kernel_sincosdf(x - p1pio2, cs, sn);
+ *cs = -*cs;
+ } else {
+ __kernel_sincosdf(x + p1pio2, cs, sn);
+ *sn = -*sn;
+ }
+ } else {
+ if (hx > 0)
+ __kernel_sincosdf(x - p2pio2, sn, cs);
+ else
+ __kernel_sincosdf(x + p2pio2, sn, cs);
+ *sn = -*sn;
+ *cs = -*cs;
+ }
+ return;
+ }
+
+ if (ix <= 0x40e231d5) { /* |x| ~<= 9*pi/4 */
+ if (ix <= 0x40afeddf) { /* |x| ~<= 7*pi/4 */
+ if (hx > 0) {
+ __kernel_sincosdf(x - p3pio2, cs, sn);
+ *sn = -*sn;
+ } else {
+ __kernel_sincosdf(x + p3pio2, cs, sn);
+ *cs = -*cs;
+ }
+ } else {
+ if (hx > 0)
+ __kernel_sincosdf(x - p4pio2, sn, cs);
+ else
+ __kernel_sincosdf(x + p4pio2, sn, cs);
+ }
+ return;
+ }
+
+ /* If x = Inf or NaN, then sin(x) = NaN and cos(x) = NaN. */
+ if (ix >= 0x7f800000) {
+ *sn = x - x;
+ *cs = x - x;
+ return;
+ }
+
+ /* Argument reduction. */
+ n = __ieee754_rem_pio2f(x, &y);
+ __kernel_sincosdf(y, &s, &c);
+
+ switch(n & 3) {
+ case 0:
+ *sn = s;
+ *cs = c;
+ break;
+ case 1:
+ *sn = c;
+ *cs = -s;
+ break;
+ case 2:
+ *sn = -s;
+ *cs = -c;
+ break;
+ default:
+ *sn = -c;
+ *cs = s;
+ }
+}
+
+
diff --git a/lib/msun/src/s_sincosl.c b/lib/msun/src/s_sincosl.c
new file mode 100644
index 0000000..aef36c2
--- /dev/null
+++ b/lib/msun/src/s_sincosl.c
@@ -0,0 +1,105 @@
+/*-
+ * Copyright (c) 2007, 2010-2013 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.
+ *
+ * s_sinl.c and s_cosl.c merged by Steven G. Kargl.
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+#include <float.h>
+#ifdef __i386__
+#include <ieeefp.h>
+#endif
+
+#include "math.h"
+#include "math_private.h"
+#include "k_sincosl.h"
+
+#if LDBL_MANT_DIG == 64
+#include "../ld80/e_rem_pio2l.h"
+#elif LDBL_MANT_DIG == 113
+#include "../ld128/e_rem_pio2l.h"
+#else
+#error "Unsupported long double format"
+#endif
+
+void
+sincosl(long double x, long double *sn, long double *cs)
+{
+ union IEEEl2bits z;
+ int e0, sgn;
+ long double y[2];
+
+ z.e = x;
+ sgn = z.bits.sign;
+ z.bits.sign = 0;
+
+ ENTERV();
+
+ /* Optimize the case where x is already within range. */
+ if (z.e < M_PI_4) {
+ /*
+ * If x = +-0 or x is a subnormal number, then sin(x) = x and
+ * cos(x) = 1.
+ */
+ if (z.bits.exp == 0) {
+ *sn = x;
+ *cs = 1;
+ } else
+ __kernel_sincosl(x, 0, 0, sn, cs);
+ RETURNV();
+ }
+
+ /* If x = NaN or Inf, then sin(x) and cos(x) are NaN. */
+ if (z.bits.exp == 32767) {
+ *sn = x - x;
+ *cs = x - x;
+ RETURNV();
+ }
+
+ /* Range reduction. */
+ e0 = __ieee754_rem_pio2l(x, y);
+
+ switch (e0 & 3) {
+ case 0:
+ __kernel_sincosl(y[0], y[1], 1, sn, cs);
+ break;
+ case 1:
+ __kernel_sincosl(y[0], y[1], 1, cs, sn);
+ *cs = -*cs;
+ break;
+ case 2:
+ __kernel_sincosl(y[0], y[1], 1, sn, cs);
+ *sn = -*sn;
+ *cs = -*cs;
+ break;
+ default:
+ __kernel_sincosl(y[0], y[1], 1, cs, sn);
+ *sn = -*sn;
+ }
+
+ RETURNV();
+}
diff --git a/lib/msun/src/s_tanhl.c b/lib/msun/src/s_tanhl.c
index 886158b..b753186 100644
--- a/lib/msun/src/s_tanhl.c
+++ b/lib/msun/src/s_tanhl.c
@@ -113,7 +113,9 @@ long double
tanhl(long double x)
{
long double hi,lo,s,x2,x4,z;
+#if LDBL_MANT_DIG == 113
double dx2;
+#endif
int16_t jx,ix;
GET_LDBL_EXPSIGN(jx,x);
OpenPOWER on IntegriCloud