public inbox for newlib-cvs@sourceware.org
help / color / mirror / Atom feed
* [newlib-cygwin] newlib: Add FreeBSD files for non LDBL_EQ_DBL support
@ 2023-05-16 23:35 Joel Sherrill
0 siblings, 0 replies; only message in thread
From: Joel Sherrill @ 2023-05-16 23:35 UTC (permalink / raw)
To: newlib-cvs
https://sourceware.org/git/gitweb.cgi?p=newlib-cygwin.git;h=c630a6a837fd581d741fb94602aed6a4a5ed9bf4
commit c630a6a837fd581d741fb94602aed6a4a5ed9bf4
Author: Jennifer Averett <jennifer.averett@oarcorp.com>
Date: Fri May 5 13:39:16 2023 -0500
newlib: Add FreeBSD files for non LDBL_EQ_DBL support
FreeBSD files to add long double support for i386,
aarch64 and x86_64.
Diff:
---
newlib/libc/include/sys/endian.h | 207 ++++++
newlib/libc/machine/aarch64/machine/_fpmath.h | 58 ++
newlib/libc/machine/i386/machine/_fpmath.h | 56 ++
newlib/libc/machine/x86_64/machine/_fpmath.h | 57 ++
newlib/libm/ld/e_acoshl.c | 89 +++
newlib/libm/ld/e_acosl.c | 87 +++
newlib/libm/ld/e_asinl.c | 77 +++
newlib/libm/ld/e_atan2l.c | 120 ++++
newlib/libm/ld/e_atanhl.c | 74 +++
newlib/libm/ld/e_coshl.c | 132 ++++
newlib/libm/ld/e_fmodl.c | 149 +++++
newlib/libm/ld/e_lgammal.c | 25 +
newlib/libm/ld/e_remainderl.c | 40 ++
newlib/libm/ld/e_sinhl.c | 134 ++++
newlib/libm/ld/fpmath.h | 82 +++
newlib/libm/ld/math_private.h | 924 ++++++++++++++++++++++++++
newlib/libm/ld/s_asinhl.c | 91 +++
newlib/libm/ld/s_atanl.c | 85 +++
newlib/libm/ld/s_cbrtl.c | 143 ++++
newlib/libm/ld/s_ceill.c | 101 +++
newlib/libm/ld/s_copysignl.c | 44 ++
newlib/libm/ld/s_cosl.c | 102 +++
newlib/libm/ld/s_fabsl.c | 45 ++
newlib/libm/ld/s_fdim.c | 48 ++
newlib/libm/ld/s_floorl.c | 101 +++
newlib/libm/ld/s_fmal.c | 274 ++++++++
newlib/libm/ld/s_fmaxl.c | 57 ++
newlib/libm/ld/s_fminl.c | 57 ++
newlib/libm/ld/s_frexpl.c | 64 ++
newlib/libm/ld/s_ilogbl.c | 53 ++
newlib/libm/ld/s_llrintl.c | 9 +
newlib/libm/ld/s_llroundl.c | 11 +
newlib/libm/ld/s_logbl.c | 54 ++
newlib/libm/ld/s_lrint.c | 60 ++
newlib/libm/ld/s_lrintl.c | 9 +
newlib/libm/ld/s_lround.c | 70 ++
newlib/libm/ld/s_lroundl.c | 11 +
newlib/libm/ld/s_modfl.c | 103 +++
newlib/libm/ld/s_nearbyint.c | 61 ++
newlib/libm/ld/s_nextafterl.c | 80 +++
newlib/libm/ld/s_nexttoward.c | 72 ++
newlib/libm/ld/s_nexttowardf.c | 59 ++
newlib/libm/ld/s_remquol.c | 173 +++++
newlib/libm/ld/s_rintl.c | 92 +++
newlib/libm/ld/s_roundl.c | 64 ++
newlib/libm/ld/s_scalbln.c | 56 ++
newlib/libm/ld/s_scalbnl.c | 48 ++
newlib/libm/ld/s_sinl.c | 95 +++
newlib/libm/ld/s_tanhl.c | 174 +++++
newlib/libm/ld/s_tanl.c | 97 +++
newlib/libm/ld/s_truncl.c | 68 ++
newlib/libm/ld128/b_tgammal.c | 57 ++
newlib/libm/ld128/e_lgammal_r.c | 330 +++++++++
newlib/libm/ld128/e_powl.c | 443 ++++++++++++
newlib/libm/ld128/e_rem_pio2l.h | 135 ++++
newlib/libm/ld128/invtrig.c | 102 +++
newlib/libm/ld128/invtrig.h | 115 ++++
newlib/libm/ld128/k_cosl.c | 59 ++
newlib/libm/ld128/k_expl.h | 324 +++++++++
newlib/libm/ld128/k_sinl.c | 59 ++
newlib/libm/ld128/s_erfl.c | 329 +++++++++
newlib/libm/ld128/s_exp2l.c | 429 ++++++++++++
newlib/libm/ld128/s_expl.c | 326 +++++++++
newlib/libm/ld128/s_logl.c | 740 +++++++++++++++++++++
newlib/libm/ld80/b_expl.c | 113 ++++
newlib/libm/ld80/b_logl.c | 375 +++++++++++
newlib/libm/ld80/b_tgammal.c | 419 ++++++++++++
newlib/libm/ld80/e_lgammal_r.c | 358 ++++++++++
newlib/libm/ld80/e_powl.c | 662 ++++++++++++++++++
newlib/libm/ld80/e_rem_pio2l.h | 143 ++++
newlib/libm/ld80/invtrig.c | 84 +++
newlib/libm/ld80/invtrig.h | 116 ++++
newlib/libm/ld80/k_cosl.c | 78 +++
newlib/libm/ld80/k_cospil.h | 42 ++
newlib/libm/ld80/k_expl.h | 301 +++++++++
newlib/libm/ld80/k_sinl.c | 62 ++
newlib/libm/ld80/k_sinpil.h | 42 ++
newlib/libm/ld80/s_cospil.c | 129 ++++
newlib/libm/ld80/s_erfl.c | 337 ++++++++++
newlib/libm/ld80/s_exp2l.c | 290 ++++++++
newlib/libm/ld80/s_expl.c | 279 ++++++++
newlib/libm/ld80/s_logl.c | 722 ++++++++++++++++++++
newlib/libm/ld80/s_sinpil.c | 140 ++++
83 files changed, 13182 insertions(+)
diff --git a/newlib/libc/include/sys/endian.h b/newlib/libc/include/sys/endian.h
new file mode 100644
index 000000000..3685322b6
--- /dev/null
+++ b/newlib/libc/include/sys/endian.h
@@ -0,0 +1,207 @@
+/*-
+ * Copyright (c) 2002 Thomas Moestl <tmm@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.
+ *
+ * $FreeBSD: head/sys/sys/endian.h 208331 2010-05-20 06:16:13Z phk $
+ */
+
+#ifndef _SYS_ENDIAN_H_
+#define _SYS_ENDIAN_H_
+
+#include <sys/cdefs.h>
+#include <sys/_types.h>
+#include <machine/endian.h>
+
+#ifndef _UINT8_T_DECLARED
+typedef __uint8_t uint8_t;
+#define _UINT8_T_DECLARED
+#endif
+
+#ifndef _UINT16_T_DECLARED
+typedef __uint16_t uint16_t;
+#define _UINT16_T_DECLARED
+#endif
+
+#ifndef _UINT32_T_DECLARED
+typedef __uint32_t uint32_t;
+#define _UINT32_T_DECLARED
+#endif
+
+#ifndef _UINT64_T_DECLARED
+typedef __uint64_t uint64_t;
+#define _UINT64_T_DECLARED
+#endif
+
+/*
+ * General byte order swapping functions.
+ */
+#define bswap16(x) __bswap16(x)
+#define bswap32(x) __bswap32(x)
+#define bswap64(x) __bswap64(x)
+
+/*
+ * Host to big endian, host to little endian, big endian to host, and little
+ * endian to host byte order functions as detailed in byteorder(9).
+ */
+#if _BYTE_ORDER == _LITTLE_ENDIAN
+#define htobe16(x) bswap16((x))
+#define htobe32(x) bswap32((x))
+#define htobe64(x) bswap64((x))
+#define htole16(x) ((uint16_t)(x))
+#define htole32(x) ((uint32_t)(x))
+#define htole64(x) ((uint64_t)(x))
+
+#define be16toh(x) bswap16((x))
+#define be32toh(x) bswap32((x))
+#define be64toh(x) bswap64((x))
+#define le16toh(x) ((uint16_t)(x))
+#define le32toh(x) ((uint32_t)(x))
+#define le64toh(x) ((uint64_t)(x))
+#else /* _BYTE_ORDER != _LITTLE_ENDIAN */
+#define htobe16(x) ((uint16_t)(x))
+#define htobe32(x) ((uint32_t)(x))
+#define htobe64(x) ((uint64_t)(x))
+#define htole16(x) bswap16((x))
+#define htole32(x) bswap32((x))
+#define htole64(x) bswap64((x))
+
+#define be16toh(x) ((uint16_t)(x))
+#define be32toh(x) ((uint32_t)(x))
+#define be64toh(x) ((uint64_t)(x))
+#define le16toh(x) bswap16((x))
+#define le32toh(x) bswap32((x))
+#define le64toh(x) bswap64((x))
+#endif /* _BYTE_ORDER == _LITTLE_ENDIAN */
+
+/* Alignment-agnostic encode/decode bytestream to/from little/big endian. */
+
+static __inline uint16_t
+be16dec(const void *pp)
+{
+ uint8_t const *p = (uint8_t const *)pp;
+
+ return (((unsigned)p[0] << 8) | p[1]);
+}
+
+static __inline uint32_t
+be32dec(const void *pp)
+{
+ uint8_t const *p = (uint8_t const *)pp;
+
+ return (((uint32_t)p[0] << 24) | ((uint32_t)p[1] << 16) |
+ ((uint32_t)p[2] << 8) | p[3]);
+}
+
+static __inline uint64_t
+be64dec(const void *pp)
+{
+ uint8_t const *p = (uint8_t const *)pp;
+
+ return (((uint64_t)be32dec(p) << 32) | be32dec(p + 4));
+}
+
+static __inline uint16_t
+le16dec(const void *pp)
+{
+ uint8_t const *p = (uint8_t const *)pp;
+
+ return (((unsigned)p[1] << 8) | p[0]);
+}
+
+static __inline uint32_t
+le32dec(const void *pp)
+{
+ uint8_t const *p = (uint8_t const *)pp;
+
+ return (((uint32_t)p[3] << 24) | ((uint32_t)p[2] << 16) |
+ ((uint32_t)p[1] << 8) | p[0]);
+}
+
+static __inline uint64_t
+le64dec(const void *pp)
+{
+ uint8_t const *p = (uint8_t const *)pp;
+
+ return (((uint64_t)le32dec(p + 4) << 32) | le32dec(p));
+}
+
+static __inline void
+be16enc(void *pp, uint16_t u)
+{
+ uint8_t *p = (uint8_t *)pp;
+
+ p[0] = (u >> 8) & 0xff;
+ p[1] = u & 0xff;
+}
+
+static __inline void
+be32enc(void *pp, uint32_t u)
+{
+ uint8_t *p = (uint8_t *)pp;
+
+ p[0] = (u >> 24) & 0xff;
+ p[1] = (u >> 16) & 0xff;
+ p[2] = (u >> 8) & 0xff;
+ p[3] = u & 0xff;
+}
+
+static __inline void
+be64enc(void *pp, uint64_t u)
+{
+ uint8_t *p = (uint8_t *)pp;
+
+ be32enc(p, (uint32_t)(u >> 32));
+ be32enc(p + 4, (uint32_t)(u & 0xffffffffU));
+}
+
+static __inline void
+le16enc(void *pp, uint16_t u)
+{
+ uint8_t *p = (uint8_t *)pp;
+
+ p[0] = u & 0xff;
+ p[1] = (u >> 8) & 0xff;
+}
+
+static __inline void
+le32enc(void *pp, uint32_t u)
+{
+ uint8_t *p = (uint8_t *)pp;
+
+ p[0] = u & 0xff;
+ p[1] = (u >> 8) & 0xff;
+ p[2] = (u >> 16) & 0xff;
+ p[3] = (u >> 24) & 0xff;
+}
+
+static __inline void
+le64enc(void *pp, uint64_t u)
+{
+ uint8_t *p = (uint8_t *)pp;
+
+ le32enc(p, (uint32_t)(u & 0xffffffffU));
+ le32enc(p + 4, (uint32_t)(u >> 32));
+}
+
+#endif /* _SYS_ENDIAN_H_ */
diff --git a/newlib/libc/machine/aarch64/machine/_fpmath.h b/newlib/libc/machine/aarch64/machine/_fpmath.h
new file mode 100644
index 000000000..71d0a7152
--- /dev/null
+++ b/newlib/libc/machine/aarch64/machine/_fpmath.h
@@ -0,0 +1,58 @@
+/*-
+ * Copyright (c) 2002, 2003 David Schultz <das@FreeBSD.ORG>
+ * Copyright (2) 2014 The FreeBSD Foundation
+ * 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.
+ *
+ * $FreeBSD$
+ */
+
+union IEEEl2bits {
+ long double e;
+ struct {
+ unsigned long manl :64;
+ unsigned long manh :48;
+ unsigned int exp :15;
+ unsigned int sign :1;
+ } bits;
+ /* TODO andrew: Check the packing here */
+ struct {
+ unsigned long manl :64;
+ unsigned long manh :48;
+ unsigned int expsign :16;
+ } xbits;
+};
+
+#define LDBL_NBIT 0
+#define LDBL_IMPLICIT_NBIT
+#define mask_nbit_l(u) ((void)0)
+
+#define LDBL_MANH_SIZE 48
+#define LDBL_MANL_SIZE 64
+
+#define LDBL_TO_ARRAY32(u, a) do { \
+ (a)[0] = (uint32_t)(u).bits.manl; \
+ (a)[1] = (uint32_t)((u).bits.manl >> 32); \
+ (a)[2] = (uint32_t)(u).bits.manh; \
+ (a)[3] = (uint32_t)((u).bits.manh >> 32); \
+} while(0)
diff --git a/newlib/libc/machine/i386/machine/_fpmath.h b/newlib/libc/machine/i386/machine/_fpmath.h
new file mode 100644
index 000000000..874439c34
--- /dev/null
+++ b/newlib/libc/machine/i386/machine/_fpmath.h
@@ -0,0 +1,56 @@
+/*-
+ * SPDX-License-Identifier: BSD-2-Clause-FreeBSD
+ *
+ * Copyright (c) 2002, 2003 David Schultz <das@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.
+ *
+ * $FreeBSD$
+ */
+
+union IEEEl2bits {
+ long double e;
+ struct {
+ unsigned int manl :32;
+ unsigned int manh :32;
+ unsigned int exp :15;
+ unsigned int sign :1;
+ unsigned int junk :16;
+ } bits;
+ struct {
+ unsigned long long man :64;
+ unsigned int expsign :16;
+ unsigned int junk :16;
+ } xbits;
+};
+
+#define LDBL_NBIT 0x80000000
+#define mask_nbit_l(u) ((u).bits.manh &= ~LDBL_NBIT)
+
+#define LDBL_MANH_SIZE 32
+#define LDBL_MANL_SIZE 32
+
+#define LDBL_TO_ARRAY32(u, a) do { \
+ (a)[0] = (uint32_t)(u).bits.manl; \
+ (a)[1] = (uint32_t)(u).bits.manh; \
+} while (0)
diff --git a/newlib/libc/machine/x86_64/machine/_fpmath.h b/newlib/libc/machine/x86_64/machine/_fpmath.h
new file mode 100644
index 000000000..8be7b7dba
--- /dev/null
+++ b/newlib/libc/machine/x86_64/machine/_fpmath.h
@@ -0,0 +1,57 @@
+/*-
+ * SPDX-License-Identifier: BSD-2-Clause-FreeBSD
+ *
+ * Copyright (c) 2002, 2003 David Schultz <das@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.
+ *
+ * $FreeBSD$
+ */
+
+union IEEEl2bits {
+ long double e;
+ struct {
+ unsigned int manl :32;
+ unsigned int manh :32;
+ unsigned int exp :15;
+ unsigned int sign :1;
+ unsigned int junkl :16;
+ unsigned int junkh :32;
+ } bits;
+ struct {
+ unsigned long man :64;
+ unsigned int expsign :16;
+ unsigned long junk :48;
+ } xbits;
+};
+
+#define LDBL_NBIT 0x80000000
+#define mask_nbit_l(u) ((u).bits.manh &= ~LDBL_NBIT)
+
+#define LDBL_MANH_SIZE 32
+#define LDBL_MANL_SIZE 32
+
+#define LDBL_TO_ARRAY32(u, a) do { \
+ (a)[0] = (uint32_t)(u).bits.manl; \
+ (a)[1] = (uint32_t)(u).bits.manh; \
+} while (0)
diff --git a/newlib/libm/ld/e_acoshl.c b/newlib/libm/ld/e_acoshl.c
new file mode 100644
index 000000000..59faeb0fb
--- /dev/null
+++ b/newlib/libm/ld/e_acoshl.c
@@ -0,0 +1,89 @@
+/* from: FreeBSD: head/lib/msun/src/e_acosh.c 176451 2008-02-22 02:30:36Z das */
+
+/* @(#)e_acosh.c 1.3 95/01/18 */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ *
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/*
+ * See e_acosh.c for complete comments.
+ *
+ * Converted to long double by David Schultz <das@FreeBSD.ORG> and
+ * Bruce D. Evans.
+ */
+
+#include <float.h>
+#ifdef __i386__
+#include <ieeefp.h>
+#endif
+
+#include "fpmath.h"
+#include "math.h"
+#include "math_private.h"
+
+/* EXP_LARGE is the threshold above which we use acosh(x) ~= log(2x). */
+#if LDBL_MANT_DIG == 64
+#define EXP_LARGE 34
+#elif LDBL_MANT_DIG == 113
+#define EXP_LARGE 58
+#else
+#error "Unsupported long double format"
+#endif
+
+#if LDBL_MAX_EXP != 0x4000
+/* We also require the usual expsign encoding. */
+#error "Unsupported long double format"
+#endif
+
+#define BIAS (LDBL_MAX_EXP - 1)
+
+static const double
+one = 1.0;
+
+#if LDBL_MANT_DIG == 64
+static const union IEEEl2bits
+u_ln2 = LD80C(0xb17217f7d1cf79ac, -1, 6.93147180559945309417e-1L);
+#define ln2 u_ln2.e
+#elif LDBL_MANT_DIG == 113
+static const long double
+ln2 = 6.93147180559945309417232121458176568e-1L; /* 0x162e42fefa39ef35793c7673007e6.0p-113 */
+#else
+#error "Unsupported long double format"
+#endif
+
+long double
+acoshl(long double x)
+{
+ long double t;
+ int16_t hx;
+
+ ENTERI();
+ GET_LDBL_EXPSIGN(hx, x);
+ if (hx < 0x3fff) { /* x < 1, or misnormal */
+ RETURNI((x-x)/(x-x));
+ } else if (hx >= BIAS + EXP_LARGE) { /* x >= LARGE */
+ if (hx >= 0x7fff) { /* x is inf, NaN or misnormal */
+ RETURNI(x+x);
+ } else
+ RETURNI(logl(x)+ln2); /* acosh(huge)=log(2x), or misnormal */
+ } else if (hx == 0x3fff && x == 1) {
+ RETURNI(0.0); /* acosh(1) = 0 */
+ } else if (hx >= 0x4000) { /* LARGE > x >= 2, or misnormal */
+ t=x*x;
+ RETURNI(logl(2.0*x-one/(x+sqrtl(t-one))));
+ } else { /* 1<x<2 */
+ t = x-one;
+ RETURNI(log1pl(t+sqrtl(2.0*t+t*t)));
+ }
+}
diff --git a/newlib/libm/ld/e_acosl.c b/newlib/libm/ld/e_acosl.c
new file mode 100644
index 000000000..3e82e98ba
--- /dev/null
+++ b/newlib/libm/ld/e_acosl.c
@@ -0,0 +1,87 @@
+
+/* @(#)e_acos.c 1.3 95/01/18 */
+/* FreeBSD: head/lib/msun/src/e_acos.c 176451 2008-02-22 02:30:36Z das */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/*
+ * See comments in e_acos.c.
+ * Converted to long double by David Schultz <das@FreeBSD.ORG>.
+ */
+
+#include <float.h>
+
+#include "invtrig.h"
+#include "math.h"
+#include "math_private.h"
+
+static const long double
+one= 1.00000000000000000000e+00;
+
+#ifdef __i386__
+/* XXX Work around the fact that gcc truncates long double constants on i386 */
+static volatile double
+pi1 = 3.14159265358979311600e+00, /* 0x1.921fb54442d18p+1 */
+pi2 = 1.22514845490862001043e-16; /* 0x1.1a80000000000p-53 */
+#define pi ((long double)pi1 + pi2)
+#else
+static const long double
+pi = 3.14159265358979323846264338327950280e+00L;
+#endif
+
+long double
+acosl(long double x)
+{
+ union IEEEl2bits u;
+ long double z,p,q,r,w,s,c,df;
+ int16_t expsign, expt;
+ u.e = x;
+ expsign = u.xbits.expsign;
+ expt = expsign & 0x7fff;
+ if(expt >= BIAS) { /* |x| >= 1 */
+ if(expt==BIAS && ((u.bits.manh&~LDBL_NBIT)|u.bits.manl)==0) {
+ if (expsign>0) return 0.0; /* acos(1) = 0 */
+ else return pi+2.0*pio2_lo; /* acos(-1)= pi */
+ }
+ return (x-x)/(x-x); /* acos(|x|>1) is NaN */
+ }
+ if(expt<BIAS-1) { /* |x| < 0.5 */
+ if(expt<ACOS_CONST) return pio2_hi+pio2_lo;/*x tiny: acosl=pi/2*/
+ z = x*x;
+ p = P(z);
+ q = Q(z);
+ r = p/q;
+ return pio2_hi - (x - (pio2_lo-x*r));
+ } else if (expsign<0) { /* x < -0.5 */
+ z = (one+x)*0.5;
+ p = P(z);
+ q = Q(z);
+ s = sqrtl(z);
+ r = p/q;
+ w = r*s-pio2_lo;
+ return pi - 2.0*(s+w);
+ } else { /* x > 0.5 */
+ z = (one-x)*0.5;
+ s = sqrtl(z);
+ u.e = s;
+ u.bits.manl = 0;
+ df = u.e;
+ c = (z-df*df)/(s+df);
+ p = P(z);
+ q = Q(z);
+ r = p/q;
+ w = r*s+c;
+ return 2.0*(df+w);
+ }
+}
diff --git a/newlib/libm/ld/e_asinl.c b/newlib/libm/ld/e_asinl.c
new file mode 100644
index 000000000..066a591f4
--- /dev/null
+++ b/newlib/libm/ld/e_asinl.c
@@ -0,0 +1,77 @@
+
+/* @(#)e_asin.c 1.3 95/01/18 */
+/* FreeBSD: head/lib/msun/src/e_asin.c 176451 2008-02-22 02:30:36Z das */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/*
+ * See comments in e_asin.c.
+ * Converted to long double by David Schultz <das@FreeBSD.ORG>.
+ */
+
+#include <float.h>
+
+#include "invtrig.h"
+#include "math.h"
+#include "math_private.h"
+
+static const long double
+one = 1.00000000000000000000e+00,
+huge = 1.000e+300;
+
+long double
+asinl(long double x)
+{
+ union IEEEl2bits u;
+ long double t=0.0,w,p,q,c,r,s;
+ int16_t expsign, expt;
+ u.e = x;
+ expsign = u.xbits.expsign;
+ expt = expsign & 0x7fff;
+ if(expt >= BIAS) { /* |x|>= 1 */
+ if(expt==BIAS && ((u.bits.manh&~LDBL_NBIT)|u.bits.manl)==0)
+ /* asin(1)=+-pi/2 with inexact */
+ return x*pio2_hi+x*pio2_lo;
+ return (x-x)/(x-x); /* asin(|x|>1) is NaN */
+ } else if (expt<BIAS-1) { /* |x|<0.5 */
+ if(expt<ASIN_LINEAR) { /* if |x| is small, asinl(x)=x */
+ if(huge+x>one) return x;/* return x with inexact if x!=0*/
+ }
+ t = x*x;
+ p = P(t);
+ q = Q(t);
+ w = p/q;
+ return x+x*w;
+ }
+ /* 1> |x|>= 0.5 */
+ w = one-fabsl(x);
+ t = w*0.5;
+ p = P(t);
+ q = Q(t);
+ s = sqrtl(t);
+ if(u.bits.manh>=THRESH) { /* if |x| is close to 1 */
+ w = p/q;
+ t = pio2_hi-(2.0*(s+s*w)-pio2_lo);
+ } else {
+ u.e = s;
+ u.bits.manl = 0;
+ w = u.e;
+ c = (t-w*w)/(s+w);
+ r = p/q;
+ p = 2.0*s*r-(pio2_lo-2.0*c);
+ q = pio4_hi-2.0*w;
+ t = pio4_hi-(p-q);
+ }
+ if(expsign>0) return t; else return -t;
+}
diff --git a/newlib/libm/ld/e_atan2l.c b/newlib/libm/ld/e_atan2l.c
new file mode 100644
index 000000000..619bfb595
--- /dev/null
+++ b/newlib/libm/ld/e_atan2l.c
@@ -0,0 +1,120 @@
+
+/* @(#)e_atan2.c 1.3 95/01/18 */
+/* FreeBSD: head/lib/msun/src/e_atan2.c 176451 2008-02-22 02:30:36Z das */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ *
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/*
+ * See comments in e_atan2.c.
+ * Converted to long double by David Schultz <das@FreeBSD.ORG>.
+ */
+
+#include <float.h>
+
+#include "invtrig.h"
+#include "math.h"
+#include "math_private.h"
+
+static volatile long double
+tiny = 1.0e-300;
+static const long double
+zero = 0.0;
+
+#ifdef __i386__
+/* XXX Work around the fact that gcc truncates long double constants on i386 */
+static volatile double
+pi1 = 3.14159265358979311600e+00, /* 0x1.921fb54442d18p+1 */
+pi2 = 1.22514845490862001043e-16; /* 0x1.1a80000000000p-53 */
+#define pi ((long double)pi1 + pi2)
+#else
+static const long double
+pi = 3.14159265358979323846264338327950280e+00L;
+#endif
+
+long double
+atan2l(long double y, long double x)
+{
+ union IEEEl2bits ux, uy;
+ long double z;
+ int32_t k,m;
+ int16_t exptx, expsignx, expty, expsigny;
+
+ uy.e = y;
+ expsigny = uy.xbits.expsign;
+ expty = expsigny & 0x7fff;
+ ux.e = x;
+ expsignx = ux.xbits.expsign;
+ exptx = expsignx & 0x7fff;
+
+ if ((exptx==BIAS+LDBL_MAX_EXP &&
+ ((ux.bits.manh&~LDBL_NBIT)|ux.bits.manl)!=0) || /* x is NaN */
+ (expty==BIAS+LDBL_MAX_EXP &&
+ ((uy.bits.manh&~LDBL_NBIT)|uy.bits.manl)!=0)) /* y is NaN */
+ return nan_mix(x, y);
+ if (expsignx==BIAS && ((ux.bits.manh&~LDBL_NBIT)|ux.bits.manl)==0)
+ return atanl(y); /* x=1.0 */
+ m = ((expsigny>>15)&1)|((expsignx>>14)&2); /* 2*sign(x)+sign(y) */
+
+ /* when y = 0 */
+ if(expty==0 && ((uy.bits.manh&~LDBL_NBIT)|uy.bits.manl)==0) {
+ switch(m) {
+ case 0:
+ case 1: return y; /* atan(+-0,+anything)=+-0 */
+ case 2: return pi+tiny;/* atan(+0,-anything) = pi */
+ case 3: return -pi-tiny;/* atan(-0,-anything) =-pi */
+ }
+ }
+ /* when x = 0 */
+ if(exptx==0 && ((ux.bits.manh&~LDBL_NBIT)|ux.bits.manl)==0)
+ return (expsigny<0)? -pio2_hi-tiny: pio2_hi+tiny;
+
+ /* when x is INF */
+ if(exptx==BIAS+LDBL_MAX_EXP) {
+ if(expty==BIAS+LDBL_MAX_EXP) {
+ switch(m) {
+ case 0: return pio2_hi*0.5+tiny;/* atan(+INF,+INF) */
+ case 1: return -pio2_hi*0.5-tiny;/* atan(-INF,+INF) */
+ case 2: return 1.5*pio2_hi+tiny;/*atan(+INF,-INF)*/
+ case 3: return -1.5*pio2_hi-tiny;/*atan(-INF,-INF)*/
+ }
+ } else {
+ switch(m) {
+ case 0: return zero ; /* atan(+...,+INF) */
+ case 1: return -zero ; /* atan(-...,+INF) */
+ case 2: return pi+tiny ; /* atan(+...,-INF) */
+ case 3: return -pi-tiny ; /* atan(-...,-INF) */
+ }
+ }
+ }
+ /* when y is INF */
+ if(expty==BIAS+LDBL_MAX_EXP)
+ return (expsigny<0)? -pio2_hi-tiny: pio2_hi+tiny;
+
+ /* compute y/x */
+ k = expty-exptx;
+ if(k > LDBL_MANT_DIG+2) { /* |y/x| huge */
+ z=pio2_hi+pio2_lo;
+ m&=1;
+ }
+ else if(expsignx<0&&k<-LDBL_MANT_DIG-2) z=0.0; /* |y/x| tiny, x<0 */
+ else z=atanl(fabsl(y/x)); /* safe to do y/x */
+ switch (m) {
+ case 0: return z ; /* atan(+,+) */
+ case 1: return -z ; /* atan(-,+) */
+ case 2: return pi-(z-pi_lo);/* atan(+,-) */
+ default: /* case 3 */
+ return (z-pi_lo)-pi;/* atan(-,-) */
+ }
+}
diff --git a/newlib/libm/ld/e_atanhl.c b/newlib/libm/ld/e_atanhl.c
new file mode 100644
index 000000000..a888426d0
--- /dev/null
+++ b/newlib/libm/ld/e_atanhl.c
@@ -0,0 +1,74 @@
+/* from: FreeBSD: head/lib/msun/src/e_atanh.c 176451 2008-02-22 02:30:36Z das */
+
+/* @(#)e_atanh.c 1.3 95/01/18 */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ *
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/*
+ * See e_atanh.c for complete comments.
+ *
+ * Converted to long double by David Schultz <das@FreeBSD.ORG> and
+ * Bruce D. Evans.
+ */
+
+#include <float.h>
+#ifdef __i386__
+#include <ieeefp.h>
+#endif
+
+#include "fpmath.h"
+#include "math.h"
+#include "math_private.h"
+
+/* EXP_TINY is the threshold below which we use atanh(x) ~= x. */
+#if LDBL_MANT_DIG == 64
+#define EXP_TINY -34
+#elif LDBL_MANT_DIG == 113
+#define EXP_TINY -58
+#else
+#error "Unsupported long double format"
+#endif
+
+#if LDBL_MAX_EXP != 0x4000
+/* We also require the usual expsign encoding. */
+#error "Unsupported long double format"
+#endif
+
+#define BIAS (LDBL_MAX_EXP - 1)
+
+static const double one = 1.0, huge = 1e300;
+static const double zero = 0.0;
+
+long double
+atanhl(long double x)
+{
+ long double t;
+ uint16_t hx, ix;
+
+ ENTERI();
+ GET_LDBL_EXPSIGN(hx, x);
+ ix = hx & 0x7fff;
+ if (ix >= 0x3fff) /* |x| >= 1, or NaN or misnormal */
+ RETURNI(fabsl(x) == 1 ? x / zero : (x - x) / (x - x));
+ if (ix < BIAS + EXP_TINY && (huge + x) > zero)
+ RETURNI(x); /* x is tiny */
+ SET_LDBL_EXPSIGN(x, ix);
+ if (ix < 0x3ffe) { /* |x| < 0.5, or misnormal */
+ t = x+x;
+ t = 0.5*log1pl(t+t*x/(one-x));
+ } else
+ t = 0.5*log1pl((x+x)/(one-x));
+ RETURNI((hx & 0x8000) == 0 ? t : -t);
+}
diff --git a/newlib/libm/ld/e_coshl.c b/newlib/libm/ld/e_coshl.c
new file mode 100644
index 000000000..4e3b28311
--- /dev/null
+++ b/newlib/libm/ld/e_coshl.c
@@ -0,0 +1,132 @@
+/* from: FreeBSD: head/lib/msun/src/e_coshl.c XXX */
+
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/*
+ * See e_cosh.c for complete comments.
+ *
+ * Converted to long double by Bruce D. Evans.
+ */
+
+#include <float.h>
+#ifdef __i386__
+#include <ieeefp.h>
+#endif
+
+#include "fpmath.h"
+#include "math.h"
+#include "math_private.h"
+#include "k_expl.h"
+
+#if LDBL_MAX_EXP != 0x4000
+/* We also require the usual expsign encoding. */
+#error "Unsupported long double format"
+#endif
+
+#define BIAS (LDBL_MAX_EXP - 1)
+
+static const volatile long double huge = 0x1p10000L, tiny = 0x1p-10000L;
+#if LDBL_MANT_DIG == 64
+/*
+ * Domain [-1, 1], range ~[-1.8211e-21, 1.8211e-21]:
+ * |cosh(x) - c(x)| < 2**-68.8
+ */
+static const union IEEEl2bits
+C4u = LD80C(0xaaaaaaaaaaaaac78, -5, 4.16666666666666682297e-2L);
+#define C4 C4u.e
+static const double
+C2 = 0.5,
+C6 = 1.3888888888888616e-3, /* 0x16c16c16c16b99.0p-62 */
+C8 = 2.4801587301767953e-5, /* 0x1a01a01a027061.0p-68 */
+C10 = 2.7557319163300398e-7, /* 0x127e4fb6c9b55f.0p-74 */
+C12 = 2.0876768371393075e-9, /* 0x11eed99406a3f4.0p-81 */
+C14 = 1.1469537039374480e-11, /* 0x1938c67cd18c48.0p-89 */
+C16 = 4.8473490896852041e-14; /* 0x1b49c429701e45.0p-97 */
+#elif LDBL_MANT_DIG == 113
+/*
+ * Domain [-1, 1], range ~[-2.3194e-37, 2.3194e-37]:
+ * |cosh(x) - c(x)| < 2**-121.69
+ */
+static const long double
+C4 = 4.16666666666666666666666666666666225e-2L, /* 0x1555555555555555555555555554e.0p-117L */
+C6 = 1.38888888888888888888888888889434831e-3L, /* 0x16c16c16c16c16c16c16c16c1dd7a.0p-122L */
+C8 = 2.48015873015873015873015871870962089e-5L, /* 0x1a01a01a01a01a01a01a017af2756.0p-128L */
+C10 = 2.75573192239858906525574318600800201e-7L, /* 0x127e4fb7789f5c72ef01c8a040640.0p-134L */
+C12 = 2.08767569878680989791444691755468269e-9L, /* 0x11eed8eff8d897b543d0679607399.0p-141L */
+C14= 1.14707455977297247387801189650495351e-11L, /* 0x193974a8c07c9d24ae169a7fa9b54.0p-149L */
+C16 = 4.77947733238737883626416876486279985e-14L; /* 0x1ae7f3e733b814d4e1b90f5727fe4.0p-157L */
+static const double
+C2 = 0.5,
+C18 = 1.5619206968597871e-16, /* 0x16827863b9900b.0p-105 */
+C20 = 4.1103176218528049e-19, /* 0x1e542ba3d3c269.0p-114 */
+C22 = 8.8967926401641701e-22, /* 0x10ce399542a014.0p-122 */
+C24 = 1.6116681626523904e-24, /* 0x1f2c981d1f0cb7.0p-132 */
+C26 = 2.5022374732804632e-27; /* 0x18c7ecf8b2c4a0.0p-141 */
+#else
+#error "Unsupported long double format"
+#endif /* LDBL_MANT_DIG == 64 */
+
+/* log(2**16385 - 0.5) rounded up: */
+static const float
+o_threshold = 1.13572168e4; /* 0xb174de.0p-10 */
+
+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);
+ ix &= 0x7fff;
+
+ /* x is INF or NaN */
+ if(ix>=0x7fff) return x*x;
+
+ ENTERI();
+
+ /* |x| < 1, return 1 or c(x) */
+ if(ix<0x3fff) {
+ if (ix<BIAS-(LDBL_MANT_DIG+1)/2) /* |x| < TINY */
+ RETURNI(1+tiny); /* cosh(tiny) = 1(+) with inexact */
+ x2 = x*x;
+#if LDBL_MANT_DIG == 64
+ x4 = x2*x2;
+ RETURNI(((C16*x2 + C14)*x4 + (C12*x2 + C10))*(x4*x4*x2) +
+ ((C8*x2 + C6)*x2 + C4)*x4 + C2*x2 + 1);
+#elif LDBL_MANT_DIG == 113
+ dx2 = x2;
+ RETURNI((((((((((((C26*dx2 + C24)*dx2 + C22)*dx2 +
+ C20)*x2 + C18)*x2 +
+ C16)*x2 + C14)*x2 + C12)*x2 + C10)*x2 + C8)*x2 + C6)*x2 +
+ C4)*(x2*x2) + C2*x2 + 1);
+#endif
+ }
+
+ /* |x| in [1, 64), return accurate exp(|x|)/2+1/exp(|x|)/2 */
+ if (ix < 0x4005) {
+ k_hexpl(fabsl(x), &hi, &lo);
+ RETURNI(lo + 0.25/(hi + lo) + hi);
+ }
+
+ /* |x| in [64, o_threshold], return correctly-overflowing exp(|x|)/2 */
+ if (fabsl(x) <= o_threshold)
+ RETURNI(hexpl(fabsl(x)));
+
+ /* |x| > o_threshold, cosh(x) overflow */
+ RETURNI(huge*huge);
+}
diff --git a/newlib/libm/ld/e_fmodl.c b/newlib/libm/ld/e_fmodl.c
new file mode 100644
index 000000000..c7689b9ad
--- /dev/null
+++ b/newlib/libm/ld/e_fmodl.c
@@ -0,0 +1,149 @@
+/* @(#)e_fmod.c 1.3 95/01/18 */
+/*-
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+#include <float.h>
+#include <stdint.h>
+
+#include "fpmath.h"
+#include "math.h"
+#include "math_private.h"
+
+#define BIAS (LDBL_MAX_EXP - 1)
+
+#if LDBL_MANL_SIZE > 32
+typedef uint64_t manl_t;
+#else
+typedef uint32_t manl_t;
+#endif
+
+#if LDBL_MANH_SIZE > 32
+typedef uint64_t manh_t;
+#else
+typedef uint32_t manh_t;
+#endif
+
+/*
+ * These macros add and remove an explicit integer bit in front of the
+ * fractional mantissa, if the architecture doesn't have such a bit by
+ * default already.
+ */
+#ifdef LDBL_IMPLICIT_NBIT
+#define SET_NBIT(hx) ((hx) | (1ULL << LDBL_MANH_SIZE))
+#define HFRAC_BITS LDBL_MANH_SIZE
+#else
+#define SET_NBIT(hx) (hx)
+#define HFRAC_BITS (LDBL_MANH_SIZE - 1)
+#endif
+
+#define MANL_SHIFT (LDBL_MANL_SIZE - 1)
+
+static const long double one = 1.0, Zero[] = {0.0, -0.0,};
+
+/*
+ * fmodl(x,y)
+ * Return x mod y in exact arithmetic
+ * Method: shift and subtract
+ *
+ * Assumptions:
+ * - The low part of the mantissa fits in a manl_t exactly.
+ * - The high part of the mantissa fits in an int64_t with enough room
+ * for an explicit integer bit in front of the fractional bits.
+ */
+long double
+fmodl(long double x, long double y)
+{
+ union IEEEl2bits ux, uy;
+ int64_t hx,hz; /* We need a carry bit even if LDBL_MANH_SIZE is 32. */
+ manh_t hy;
+ manl_t lx,ly,lz;
+ int ix,iy,n,sx;
+
+ ux.e = x;
+ uy.e = y;
+ sx = ux.bits.sign;
+
+ /* purge off exception values */
+ if((uy.bits.exp|uy.bits.manh|uy.bits.manl)==0 || /* y=0 */
+ (ux.bits.exp == BIAS + LDBL_MAX_EXP) || /* or x not finite */
+ (uy.bits.exp == BIAS + LDBL_MAX_EXP &&
+ ((uy.bits.manh&~LDBL_NBIT)|uy.bits.manl)!=0)) /* or y is NaN */
+ return nan_mix_op(x, y, *)/nan_mix_op(x, y, *);
+ if(ux.bits.exp<=uy.bits.exp) {
+ if((ux.bits.exp<uy.bits.exp) ||
+ (ux.bits.manh<=uy.bits.manh &&
+ (ux.bits.manh<uy.bits.manh ||
+ ux.bits.manl<uy.bits.manl))) {
+ return x; /* |x|<|y| return x or x-y */
+ }
+ if(ux.bits.manh==uy.bits.manh && ux.bits.manl==uy.bits.manl) {
+ return Zero[sx]; /* |x|=|y| return x*0*/
+ }
+ }
+
+ /* determine ix = ilogb(x) */
+ if(ux.bits.exp == 0) { /* subnormal x */
+ ux.e *= 0x1.0p512;
+ ix = ux.bits.exp - (BIAS + 512);
+ } else {
+ ix = ux.bits.exp - BIAS;
+ }
+
+ /* determine iy = ilogb(y) */
+ if(uy.bits.exp == 0) { /* subnormal y */
+ uy.e *= 0x1.0p512;
+ iy = uy.bits.exp - (BIAS + 512);
+ } else {
+ iy = uy.bits.exp - BIAS;
+ }
+
+ /* set up {hx,lx}, {hy,ly} and align y to x */
+ hx = SET_NBIT(ux.bits.manh);
+ hy = SET_NBIT(uy.bits.manh);
+ lx = ux.bits.manl;
+ ly = uy.bits.manl;
+
+ /* fix point fmod */
+ n = ix - iy;
+
+ while(n--) {
+ hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
+ if(hz<0){hx = hx+hx+(lx>>MANL_SHIFT); lx = lx+lx;}
+ else {
+ if ((hz|lz)==0) /* return sign(x)*0 */
+ return Zero[sx];
+ hx = hz+hz+(lz>>MANL_SHIFT); lx = lz+lz;
+ }
+ }
+ hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
+ if(hz>=0) {hx=hz;lx=lz;}
+
+ /* convert back to floating value and restore the sign */
+ if((hx|lx)==0) /* return sign(x)*0 */
+ return Zero[sx];
+ while(hx<(1ULL<<HFRAC_BITS)) { /* normalize x */
+ hx = hx+hx+(lx>>MANL_SHIFT); lx = lx+lx;
+ iy -= 1;
+ }
+ ux.bits.manh = hx; /* The mantissa is truncated here if needed. */
+ ux.bits.manl = lx;
+ if (iy < LDBL_MIN_EXP) {
+ ux.bits.exp = iy + (BIAS + 512);
+ ux.e *= 0x1p-512;
+ } else {
+ ux.bits.exp = iy + BIAS;
+ }
+ x = ux.e * one; /* create necessary signal */
+ return x; /* exact output */
+}
diff --git a/newlib/libm/ld/e_lgammal.c b/newlib/libm/ld/e_lgammal.c
new file mode 100644
index 000000000..ebc2fc78c
--- /dev/null
+++ b/newlib/libm/ld/e_lgammal.c
@@ -0,0 +1,25 @@
+/* @(#)e_lgamma.c 1.3 95/01/18 */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+#include "math.h"
+#include "math_private.h"
+
+extern int signgam;
+
+long double
+lgammal(long double x)
+{
+ return lgammal_r(x,&signgam);
+}
diff --git a/newlib/libm/ld/e_remainderl.c b/newlib/libm/ld/e_remainderl.c
new file mode 100644
index 000000000..4a6786309
--- /dev/null
+++ b/newlib/libm/ld/e_remainderl.c
@@ -0,0 +1,40 @@
+/*-
+ * SPDX-License-Identifier: BSD-2-Clause-FreeBSD
+ *
+ * Copyright (c) 2008 David Schultz <das@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.
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+#include <math.h>
+
+long double
+remainderl(long double x, long double y)
+{
+ int quo;
+
+ return (remquol(x, y, &quo));
+}
diff --git a/newlib/libm/ld/e_sinhl.c b/newlib/libm/ld/e_sinhl.c
new file mode 100644
index 000000000..38d3df195
--- /dev/null
+++ b/newlib/libm/ld/e_sinhl.c
@@ -0,0 +1,134 @@
+/* from: FreeBSD: head/lib/msun/src/e_sinhl.c XXX */
+
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/*
+ * See e_sinh.c for complete comments.
+ *
+ * Converted to long double by Bruce D. Evans.
+ */
+
+#include <float.h>
+#ifdef __i386__
+#include <ieeefp.h>
+#endif
+
+#include "fpmath.h"
+#include "math.h"
+#include "math_private.h"
+#include "k_expl.h"
+
+#if LDBL_MAX_EXP != 0x4000
+/* We also require the usual expsign encoding. */
+#error "Unsupported long double format"
+#endif
+
+#define BIAS (LDBL_MAX_EXP - 1)
+
+static const long double shuge = 0x1p16383L;
+#if LDBL_MANT_DIG == 64
+/*
+ * Domain [-1, 1], range ~[-6.6749e-22, 6.6749e-22]:
+ * |sinh(x)/x - s(x)| < 2**-70.3
+ */
+static const union IEEEl2bits
+S3u = LD80C(0xaaaaaaaaaaaaaaaa, -3, 1.66666666666666666658e-1L);
+#define S3 S3u.e
+static const double
+S5 = 8.3333333333333332e-3, /* 0x11111111111111.0p-59 */
+S7 = 1.9841269841270074e-4, /* 0x1a01a01a01a070.0p-65 */
+S9 = 2.7557319223873889e-6, /* 0x171de3a5565fe6.0p-71 */
+S11 = 2.5052108406704084e-8, /* 0x1ae6456857530f.0p-78 */
+S13 = 1.6059042748655297e-10, /* 0x161245fa910697.0p-85 */
+S15 = 7.6470006914396920e-13, /* 0x1ae7ce4eff2792.0p-93 */
+S17 = 2.8346142308424267e-15; /* 0x19882ce789ffc6.0p-101 */
+#elif LDBL_MANT_DIG == 113
+/*
+ * Domain [-1, 1], range ~[-2.9673e-36, 2.9673e-36]:
+ * |sinh(x)/x - s(x)| < 2**-118.0
+ */
+static const long double
+S3 = 1.66666666666666666666666666666666033e-1L, /* 0x1555555555555555555555555553b.0p-115L */
+S5 = 8.33333333333333333333333333337643193e-3L, /* 0x111111111111111111111111180f5.0p-119L */
+S7 = 1.98412698412698412698412697391263199e-4L, /* 0x1a01a01a01a01a01a01a0176aad11.0p-125L */
+S9 = 2.75573192239858906525574406205464218e-6L, /* 0x171de3a556c7338faac243aaa9592.0p-131L */
+S11 = 2.50521083854417187749675637460977997e-8L, /* 0x1ae64567f544e38fe59b3380d7413.0p-138L */
+S13 = 1.60590438368216146368737762431552702e-10L, /* 0x16124613a86d098059c7620850fc2.0p-145L */
+S15 = 7.64716373181980539786802470969096440e-13L, /* 0x1ae7f3e733b814193af09ce723043.0p-153L */
+S17 = 2.81145725434775409870584280722701574e-15L; /* 0x1952c77030c36898c3fd0b6dfc562.0p-161L */
+static const double
+S19= 8.2206352435411005e-18, /* 0x12f49b4662b86d.0p-109 */
+S21= 1.9572943931418891e-20, /* 0x171b8f2fab9628.0p-118 */
+S23 = 3.8679983530666939e-23, /* 0x17617002b73afc.0p-127 */
+S25 = 6.5067867911512749e-26; /* 0x1423352626048a.0p-136 */
+#else
+#error "Unsupported long double format"
+#endif /* LDBL_MANT_DIG == 64 */
+
+/* log(2**16385 - 0.5) rounded up: */
+static const float
+o_threshold = 1.13572168e4; /* 0xb174de.0p-10 */
+
+long double
+sinhl(long double x)
+{
+ long double hi,lo,x2,x4;
+#if LDBL_MANT_DIG == 113
+ double dx2;
+#endif
+ double s;
+ int16_t ix,jx;
+
+ GET_LDBL_EXPSIGN(jx,x);
+ ix = jx&0x7fff;
+
+ /* x is INF or NaN */
+ if(ix>=0x7fff) return x+x;
+
+ ENTERI();
+
+ s = 1;
+ if (jx<0) s = -1;
+
+ /* |x| < 64, return x, s(x), or accurate s*(exp(|x|)/2-1/exp(|x|)/2) */
+ if (ix<0x4005) { /* |x|<64 */
+ if (ix<BIAS-(LDBL_MANT_DIG+1)/2) /* |x|<TINY */
+ if(shuge+x>1) RETURNI(x); /* sinh(tiny) = tiny with inexact */
+ if (ix<0x3fff) { /* |x|<1 */
+ x2 = x*x;
+#if LDBL_MANT_DIG == 64
+ x4 = x2*x2;
+ RETURNI(((S17*x2 + S15)*x4 + (S13*x2 + S11))*(x2*x*x4*x4) +
+ ((S9*x2 + S7)*x2 + S5)*(x2*x*x2) + S3*(x2*x) + x);
+#elif LDBL_MANT_DIG == 113
+ dx2 = x2;
+ RETURNI(((((((((((S25*dx2 + S23)*dx2 +
+ S21)*x2 + S19)*x2 +
+ S17)*x2 + S15)*x2 + S13)*x2 + S11)*x2 + S9)*x2 + S7)*x2 +
+ S5)* (x2*x*x2) +
+ S3*(x2*x) + x);
+#endif
+ }
+ k_hexpl(fabsl(x), &hi, &lo);
+ RETURNI(s*(lo - 0.25/(hi + lo) + hi));
+ }
+
+ /* |x| in [64, o_threshold], return correctly-overflowing s*exp(|x|)/2 */
+ if (fabsl(x) <= o_threshold)
+ RETURNI(s*hexpl(fabsl(x)));
+
+ /* |x| > o_threshold, sinh(x) overflow */
+ return x*shuge;
+}
diff --git a/newlib/libm/ld/fpmath.h b/newlib/libm/ld/fpmath.h
new file mode 100644
index 000000000..ce935eb40
--- /dev/null
+++ b/newlib/libm/ld/fpmath.h
@@ -0,0 +1,82 @@
+/*-
+ * SPDX-License-Identifier: BSD-2-Clause-FreeBSD
+ *
+ * Copyright (c) 2003 Mike Barcroft <mike@FreeBSD.org>
+ * Copyright (c) 2002 David Schultz <das@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.
+ *
+ * $FreeBSD$
+ */
+
+#ifndef _FPMATH_H_
+#define _FPMATH_H_
+
+#include <sys/endian.h>
+#include "_fpmath.h"
+
+#ifndef _IEEE_WORD_ORDER
+#define _IEEE_WORD_ORDER _BYTE_ORDER
+#endif
+
+union IEEEf2bits {
+ float f;
+ struct {
+#if _BYTE_ORDER == _LITTLE_ENDIAN
+ unsigned int man :23;
+ unsigned int exp :8;
+ unsigned int sign :1;
+#else /* _BIG_ENDIAN */
+ unsigned int sign :1;
+ unsigned int exp :8;
+ unsigned int man :23;
+#endif
+ } bits;
+};
+
+#define DBL_MANH_SIZE 20
+#define DBL_MANL_SIZE 32
+
+union IEEEd2bits {
+ double d;
+ struct {
+#if _BYTE_ORDER == _LITTLE_ENDIAN
+#if _IEEE_WORD_ORDER == _LITTLE_ENDIAN
+ unsigned int manl :32;
+#endif
+ unsigned int manh :20;
+ unsigned int exp :11;
+ unsigned int sign :1;
+#if _IEEE_WORD_ORDER == _BIG_ENDIAN
+ unsigned int manl :32;
+#endif
+#else /* _BIG_ENDIAN */
+ unsigned int sign :1;
+ unsigned int exp :11;
+ unsigned int manh :20;
+ unsigned int manl :32;
+#endif
+ } bits;
+};
+
+#endif /* !_FPMATH_H */
diff --git a/newlib/libm/ld/math_private.h b/newlib/libm/ld/math_private.h
new file mode 100644
index 000000000..c6d06fb14
--- /dev/null
+++ b/newlib/libm/ld/math_private.h
@@ -0,0 +1,924 @@
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+
+/*
+ * from: @(#)fdlibm.h 5.1 93/09/24
+ * $FreeBSD$
+ */
+
+#ifndef _MATH_PRIVATE_H_
+#define _MATH_PRIVATE_H_
+
+#include <sys/types.h>
+#include <machine/endian.h>
+
+/*
+ * The original fdlibm code used statements like:
+ * n0 = ((*(int*)&one)>>29)^1; * index of high word *
+ * ix0 = *(n0+(int*)&x); * high word of x *
+ * ix1 = *((1-n0)+(int*)&x); * low word of x *
+ * to dig two 32 bit words out of the 64 bit IEEE floating point
+ * value. That is non-ANSI, and, moreover, the gcc instruction
+ * scheduler gets it wrong. We instead use the following macros.
+ * Unlike the original code, we determine the endianness at compile
+ * time, not at run time; I don't see much benefit to selecting
+ * endianness at run time.
+ */
+
+/*
+ * A union which permits us to convert between a double and two 32 bit
+ * ints.
+ */
+
+#ifdef __arm__
+#if defined(__VFP_FP__) || defined(__ARM_EABI__)
+#define IEEE_WORD_ORDER BYTE_ORDER
+#else
+#define IEEE_WORD_ORDER BIG_ENDIAN
+#endif
+#else /* __arm__ */
+#define IEEE_WORD_ORDER BYTE_ORDER
+#endif
+
+/* A union which permits us to convert between a long double and
+ four 32 bit ints. */
+
+#if IEEE_WORD_ORDER == BIG_ENDIAN
+
+typedef union
+{
+ long double value;
+ struct {
+ u_int32_t mswhi;
+ u_int32_t mswlo;
+ u_int32_t lswhi;
+ u_int32_t lswlo;
+ } parts32;
+ struct {
+ u_int64_t msw;
+ u_int64_t lsw;
+ } parts64;
+} ieee_quad_shape_type;
+
+#endif
+
+#if IEEE_WORD_ORDER == LITTLE_ENDIAN
+
+typedef union
+{
+ long double value;
+ struct {
+ u_int32_t lswlo;
+ u_int32_t lswhi;
+ u_int32_t mswlo;
+ u_int32_t mswhi;
+ } parts32;
+ struct {
+ u_int64_t lsw;
+ u_int64_t msw;
+ } parts64;
+} ieee_quad_shape_type;
+
+#endif
+
+#if IEEE_WORD_ORDER == BIG_ENDIAN
+
+typedef union
+{
+ double value;
+ struct
+ {
+ u_int32_t msw;
+ u_int32_t lsw;
+ } parts;
+ struct
+ {
+ u_int64_t w;
+ } xparts;
+} ieee_double_shape_type;
+
+#endif
+
+#if IEEE_WORD_ORDER == LITTLE_ENDIAN
+
+typedef union
+{
+ double value;
+ struct
+ {
+ u_int32_t lsw;
+ u_int32_t msw;
+ } parts;
+ struct
+ {
+ u_int64_t w;
+ } xparts;
+} ieee_double_shape_type;
+
+#endif
+
+/* Get two 32 bit ints from a double. */
+
+#define EXTRACT_WORDS(ix0,ix1,d) \
+do { \
+ ieee_double_shape_type ew_u; \
+ ew_u.value = (d); \
+ (ix0) = ew_u.parts.msw; \
+ (ix1) = ew_u.parts.lsw; \
+} while (0)
+
+/* Get a 64-bit int from a double. */
+#define EXTRACT_WORD64(ix,d) \
+do { \
+ ieee_double_shape_type ew_u; \
+ ew_u.value = (d); \
+ (ix) = ew_u.xparts.w; \
+} while (0)
+
+/* Get the more significant 32 bit int from a double. */
+
+#define GET_HIGH_WORD(i,d) \
+do { \
+ ieee_double_shape_type gh_u; \
+ gh_u.value = (d); \
+ (i) = gh_u.parts.msw; \
+} while (0)
+
+/* Get the less significant 32 bit int from a double. */
+
+#define GET_LOW_WORD(i,d) \
+do { \
+ ieee_double_shape_type gl_u; \
+ gl_u.value = (d); \
+ (i) = gl_u.parts.lsw; \
+} while (0)
+
+/* Set a double from two 32 bit ints. */
+
+#define INSERT_WORDS(d,ix0,ix1) \
+do { \
+ ieee_double_shape_type iw_u; \
+ iw_u.parts.msw = (ix0); \
+ iw_u.parts.lsw = (ix1); \
+ (d) = iw_u.value; \
+} while (0)
+
+/* Set a double from a 64-bit int. */
+#define INSERT_WORD64(d,ix) \
+do { \
+ ieee_double_shape_type iw_u; \
+ iw_u.xparts.w = (ix); \
+ (d) = iw_u.value; \
+} while (0)
+
+/* Set the more significant 32 bits of a double from an int. */
+
+#define SET_HIGH_WORD(d,v) \
+do { \
+ ieee_double_shape_type sh_u; \
+ sh_u.value = (d); \
+ sh_u.parts.msw = (v); \
+ (d) = sh_u.value; \
+} while (0)
+
+/* Set the less significant 32 bits of a double from an int. */
+
+#define SET_LOW_WORD(d,v) \
+do { \
+ ieee_double_shape_type sl_u; \
+ sl_u.value = (d); \
+ sl_u.parts.lsw = (v); \
+ (d) = sl_u.value; \
+} while (0)
+
+/*
+ * A union which permits us to convert between a float and a 32 bit
+ * int.
+ */
+
+typedef union
+{
+ float value;
+ /* FIXME: Assumes 32 bit int. */
+ unsigned int word;
+} ieee_float_shape_type;
+
+/* Get a 32 bit int from a float. */
+
+#define GET_FLOAT_WORD(i,d) \
+do { \
+ ieee_float_shape_type gf_u; \
+ gf_u.value = (d); \
+ (i) = gf_u.word; \
+} while (0)
+
+/* Set a float from a 32 bit int. */
+
+#define SET_FLOAT_WORD(d,i) \
+do { \
+ ieee_float_shape_type sf_u; \
+ sf_u.word = (i); \
+ (d) = sf_u.value; \
+} while (0)
+
+/*
+ * Get expsign and mantissa as 16 bit and 64 bit ints from an 80 bit long
+ * double.
+ */
+
+#define EXTRACT_LDBL80_WORDS(ix0,ix1,d) \
+do { \
+ union IEEEl2bits ew_u; \
+ ew_u.e = (d); \
+ (ix0) = ew_u.xbits.expsign; \
+ (ix1) = ew_u.xbits.man; \
+} while (0)
+
+/*
+ * Get expsign and mantissa as one 16 bit and two 64 bit ints from a 128 bit
+ * long double.
+ */
+
+#define EXTRACT_LDBL128_WORDS(ix0,ix1,ix2,d) \
+do { \
+ union IEEEl2bits ew_u; \
+ ew_u.e = (d); \
+ (ix0) = ew_u.xbits.expsign; \
+ (ix1) = ew_u.xbits.manh; \
+ (ix2) = ew_u.xbits.manl; \
+} while (0)
+
+/* Get expsign as a 16 bit int from a long double. */
+
+#define GET_LDBL_EXPSIGN(i,d) \
+do { \
+ union IEEEl2bits ge_u; \
+ ge_u.e = (d); \
+ (i) = ge_u.xbits.expsign; \
+} while (0)
+
+/*
+ * Set an 80 bit long double from a 16 bit int expsign and a 64 bit int
+ * mantissa.
+ */
+
+#define INSERT_LDBL80_WORDS(d,ix0,ix1) \
+do { \
+ union IEEEl2bits iw_u; \
+ iw_u.xbits.expsign = (ix0); \
+ iw_u.xbits.man = (ix1); \
+ (d) = iw_u.e; \
+} while (0)
+
+/*
+ * Set a 128 bit long double from a 16 bit int expsign and two 64 bit ints
+ * comprising the mantissa.
+ */
+
+#define INSERT_LDBL128_WORDS(d,ix0,ix1,ix2) \
+do { \
+ union IEEEl2bits iw_u; \
+ iw_u.xbits.expsign = (ix0); \
+ iw_u.xbits.manh = (ix1); \
+ iw_u.xbits.manl = (ix2); \
+ (d) = iw_u.e; \
+} while (0)
+
+/* Set expsign of a long double from a 16 bit int. */
+
+#define SET_LDBL_EXPSIGN(d,v) \
+do { \
+ union IEEEl2bits se_u; \
+ se_u.e = (d); \
+ se_u.xbits.expsign = (v); \
+ (d) = se_u.e; \
+} while (0)
+
+#ifdef __i386__
+/* Long double constants are broken on i386. */
+#define LD80C(m, ex, v) { \
+ .xbits.man = __CONCAT(m, ULL), \
+ .xbits.expsign = (0x3fff + (ex)) | ((v) < 0 ? 0x8000 : 0), \
+}
+#else
+/* The above works on non-i386 too, but we use this to check v. */
+#define LD80C(m, ex, v) { .e = (v), }
+#endif
+
+#ifdef FLT_EVAL_METHOD
+/*
+ * Attempt to get strict C99 semantics for assignment with non-C99 compilers.
+ */
+#if FLT_EVAL_METHOD == 0 || __GNUC__ == 0
+#define STRICT_ASSIGN(type, lval, rval) ((lval) = (rval))
+#else
+#define STRICT_ASSIGN(type, lval, rval) do { \
+ volatile type __lval; \
+ \
+ if (sizeof(type) >= sizeof(long double)) \
+ (lval) = (rval); \
+ else { \
+ __lval = (rval); \
+ (lval) = __lval; \
+ } \
+} while (0)
+#endif
+#endif /* FLT_EVAL_METHOD */
+
+/* Support switching the mode to FP_PE if necessary. */
+#if defined(__i386__) && !defined(NO_FPSETPREC)
+#define ENTERI() ENTERIT(long double)
+#define ENTERIT(returntype) \
+ returntype __retval; \
+ fp_prec_t __oprec; \
+ \
+ if ((__oprec = fpgetprec()) != FP_PE) \
+ fpsetprec(FP_PE)
+#define RETURNI(x) do { \
+ __retval = (x); \
+ if (__oprec != FP_PE) \
+ 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()
+#define ENTERIT(x)
+#define RETURNI(x) RETURNF(x)
+#define ENTERV()
+#define RETURNV() return
+#endif
+
+/* Default return statement if hack*_t() is not used. */
+#define RETURNF(v) return (v)
+
+/*
+ * 2sum gives the same result as 2sumF without requiring |a| >= |b| or
+ * a == 0, but is slower.
+ */
+#define _2sum(a, b) do { \
+ __typeof(a) __s, __w; \
+ \
+ __w = (a) + (b); \
+ __s = __w - (a); \
+ (b) = ((a) - (__w - __s)) + ((b) - __s); \
+ (a) = __w; \
+} while (0)
+
+/*
+ * 2sumF algorithm.
+ *
+ * "Normalize" the terms in the infinite-precision expression a + b for
+ * the sum of 2 floating point values so that b is as small as possible
+ * relative to 'a'. (The resulting 'a' is the value of the expression in
+ * the same precision as 'a' and the resulting b is the rounding error.)
+ * |a| must be >= |b| or 0, b's type must be no larger than 'a's type, and
+ * exponent overflow or underflow must not occur. This uses a Theorem of
+ * Dekker (1971). See Knuth (1981) 4.2.2 Theorem C. The name "TwoSum"
+ * is apparently due to Skewchuk (1997).
+ *
+ * For this to always work, assignment of a + b to 'a' must not retain any
+ * extra precision in a + b. This is required by C standards but broken
+ * in many compilers. The brokenness cannot be worked around using
+ * STRICT_ASSIGN() like we do elsewhere, since the efficiency of this
+ * algorithm would be destroyed by non-null strict assignments. (The
+ * compilers are correct to be broken -- the efficiency of all floating
+ * point code calculations would be destroyed similarly if they forced the
+ * conversions.)
+ *
+ * Fortunately, a case that works well can usually be arranged by building
+ * any extra precision into the type of 'a' -- 'a' should have type float_t,
+ * double_t or long double. b's type should be no larger than 'a's type.
+ * Callers should use these types with scopes as large as possible, to
+ * reduce their own extra-precision and efficiciency problems. In
+ * particular, they shouldn't convert back and forth just to call here.
+ */
+#ifdef DEBUG
+#define _2sumF(a, b) do { \
+ __typeof(a) __w; \
+ volatile __typeof(a) __ia, __ib, __r, __vw; \
+ \
+ __ia = (a); \
+ __ib = (b); \
+ assert(__ia == 0 || fabsl(__ia) >= fabsl(__ib)); \
+ \
+ __w = (a) + (b); \
+ (b) = ((a) - __w) + (b); \
+ (a) = __w; \
+ \
+ /* The next 2 assertions are weak if (a) is already long double. */ \
+ assert((long double)__ia + __ib == (long double)(a) + (b)); \
+ __vw = __ia + __ib; \
+ __r = __ia - __vw; \
+ __r += __ib; \
+ assert(__vw == (a) && __r == (b)); \
+} while (0)
+#else /* !DEBUG */
+#define _2sumF(a, b) do { \
+ __typeof(a) __w; \
+ \
+ __w = (a) + (b); \
+ (b) = ((a) - __w) + (b); \
+ (a) = __w; \
+} while (0)
+#endif /* DEBUG */
+
+/*
+ * Set x += c, where x is represented in extra precision as a + b.
+ * x must be sufficiently normalized and sufficiently larger than c,
+ * and the result is then sufficiently normalized.
+ *
+ * The details of ordering are that |a| must be >= |c| (so that (a, c)
+ * can be normalized without extra work to swap 'a' with c). The details of
+ * the normalization are that b must be small relative to the normalized 'a'.
+ * Normalization of (a, c) makes the normalized c tiny relative to the
+ * normalized a, so b remains small relative to 'a' in the result. However,
+ * b need not ever be tiny relative to 'a'. For example, b might be about
+ * 2**20 times smaller than 'a' to give about 20 extra bits of precision.
+ * That is usually enough, and adding c (which by normalization is about
+ * 2**53 times smaller than a) cannot change b significantly. However,
+ * cancellation of 'a' with c in normalization of (a, c) may reduce 'a'
+ * significantly relative to b. The caller must ensure that significant
+ * cancellation doesn't occur, either by having c of the same sign as 'a',
+ * or by having |c| a few percent smaller than |a|. Pre-normalization of
+ * (a, b) may help.
+ *
+ * This is a variant of an algorithm of Kahan (see Knuth (1981) 4.2.2
+ * exercise 19). We gain considerable efficiency by requiring the terms to
+ * be sufficiently normalized and sufficiently increasing.
+ */
+#define _3sumF(a, b, c) do { \
+ __typeof(a) __tmp; \
+ \
+ __tmp = (c); \
+ _2sumF(__tmp, (a)); \
+ (b) += (a); \
+ (a) = __tmp; \
+} while (0)
+
+/*
+ * Common routine to process the arguments to nan(), nanf(), and nanl().
+ */
+void _scan_nan(uint32_t *__words, int __num_words, const char *__s);
+
+/*
+ * Mix 0, 1 or 2 NaNs. First add 0 to each arg. This normally just turns
+ * signaling NaNs into quiet NaNs by setting a quiet bit. We do this
+ * because we want to never return a signaling NaN, and also because we
+ * don't want the quiet bit to affect the result. Then mix the converted
+ * args using the specified operation.
+ *
+ * When one arg is NaN, the result is typically that arg quieted. When both
+ * args are NaNs, the result is typically the quietening of the arg whose
+ * mantissa is largest after quietening. When neither arg is NaN, the
+ * result may be NaN because it is indeterminate, or finite for subsequent
+ * construction of a NaN as the indeterminate 0.0L/0.0L.
+ *
+ * Technical complications: the result in bits after rounding to the final
+ * precision might depend on the runtime precision and/or on compiler
+ * optimizations, especially when different register sets are used for
+ * different precisions. Try to make the result not depend on at least the
+ * runtime precision by always doing the main mixing step in long double
+ * precision. Try to reduce dependencies on optimizations by adding the
+ * the 0's in different precisions (unless everything is in long double
+ * precision).
+ */
+#define nan_mix(x, y) (nan_mix_op((x), (y), +))
+#define nan_mix_op(x, y, op) (((x) + 0.0L) op ((y) + 0))
+
+#ifdef _COMPLEX_H
+
+/*
+ * C99 specifies that complex numbers have the same representation as
+ * an array of two elements, where the first element is the real part
+ * and the second element is the imaginary part.
+ */
+typedef union {
+ float complex f;
+ float a[2];
+} float_complex;
+typedef union {
+ double complex f;
+ double a[2];
+} double_complex;
+typedef union {
+ long double complex f;
+ long double a[2];
+} long_double_complex;
+#define REALPART(z) ((z).a[0])
+#define IMAGPART(z) ((z).a[1])
+
+/*
+ * Inline functions that can be used to construct complex values.
+ *
+ * The C99 standard intends x+I*y to be used for this, but x+I*y is
+ * currently unusable in general since gcc introduces many overflow,
+ * underflow, sign and efficiency bugs by rewriting I*y as
+ * (0.0+I)*(y+0.0*I) and laboriously computing the full complex product.
+ * In particular, I*Inf is corrupted to NaN+I*Inf, and I*-0 is corrupted
+ * to -0.0+I*0.0.
+ *
+ * The C11 standard introduced the macros CMPLX(), CMPLXF() and CMPLXL()
+ * to construct complex values. Compilers that conform to the C99
+ * standard require the following functions to avoid the above issues.
+ */
+
+#ifndef CMPLXF
+static __inline float complex
+CMPLXF(float x, float y)
+{
+ float_complex z;
+
+ REALPART(z) = x;
+ IMAGPART(z) = y;
+ return (z.f);
+}
+#endif
+
+#ifndef CMPLX
+static __inline double complex
+CMPLX(double x, double y)
+{
+ double_complex z;
+
+ REALPART(z) = x;
+ IMAGPART(z) = y;
+ return (z.f);
+}
+#endif
+
+#ifndef CMPLXL
+static __inline long double complex
+CMPLXL(long double x, long double y)
+{
+ long_double_complex z;
+
+ REALPART(z) = x;
+ IMAGPART(z) = y;
+ return (z.f);
+}
+#endif
+
+#endif /* _COMPLEX_H */
+
+/*
+ * The rnint() family rounds to the nearest integer for a restricted range
+ * range of args (up to about 2**MANT_DIG). We assume that the current
+ * rounding mode is FE_TONEAREST so that this can be done efficiently.
+ * Extra precision causes more problems in practice, and we only centralize
+ * this here to reduce those problems, and have not solved the efficiency
+ * problems. The exp2() family uses a more delicate version of this that
+ * requires extracting bits from the intermediate value, so it is not
+ * centralized here and should copy any solution of the efficiency problems.
+ */
+
+static inline double
+rnint(__double_t x)
+{
+ /*
+ * This casts to double to kill any extra precision. This depends
+ * on the cast being applied to a double_t to avoid compiler bugs
+ * (this is a cleaner version of STRICT_ASSIGN()). This is
+ * inefficient if there actually is extra precision, but is hard
+ * to improve on. We use double_t in the API to minimise conversions
+ * for just calling here. Note that we cannot easily change the
+ * magic number to the one that works directly with double_t, since
+ * the rounding precision is variable at runtime on x86 so the
+ * magic number would need to be variable. Assuming that the
+ * rounding precision is always the default is too fragile. This
+ * and many other complications will move when the default is
+ * changed to FP_PE.
+ */
+ return ((double)(x + 0x1.8p52) - 0x1.8p52);
+}
+
+static inline float
+rnintf(__float_t x)
+{
+ /*
+ * As for rnint(), except we could just call that to handle the
+ * extra precision case, usually without losing efficiency.
+ */
+ return ((float)(x + 0x1.8p23F) - 0x1.8p23F);
+}
+
+#ifdef LDBL_MANT_DIG
+/*
+ * The complications for extra precision are smaller for rnintl() since it
+ * can safely assume that the rounding precision has been increased from
+ * its default to FP_PE on x86. We don't exploit that here to get small
+ * optimizations from limiting the rangle to double. We just need it for
+ * the magic number to work with long doubles. ld128 callers should use
+ * rnint() instead of this if possible. ld80 callers should prefer
+ * rnintl() since for amd64 this avoids swapping the register set, while
+ * for i386 it makes no difference (assuming FP_PE), and for other arches
+ * it makes little difference.
+ */
+static inline long double
+rnintl(long double x)
+{
+ return (x + __CONCAT(0x1.8p, LDBL_MANT_DIG) / 2 -
+ __CONCAT(0x1.8p, LDBL_MANT_DIG) / 2);
+}
+#endif /* LDBL_MANT_DIG */
+
+/*
+ * irint() and i64rint() give the same result as casting to their integer
+ * return type provided their arg is a floating point integer. They can
+ * sometimes be more efficient because no rounding is required.
+ */
+#if defined(amd64) || defined(__i386__)
+#define irint(x) \
+ (sizeof(x) == sizeof(float) && \
+ sizeof(__float_t) == sizeof(long double) ? irintf(x) : \
+ sizeof(x) == sizeof(double) && \
+ sizeof(__double_t) == sizeof(long double) ? irintd(x) : \
+ sizeof(x) == sizeof(long double) ? irintl(x) : (int)(x))
+#else
+#define irint(x) ((int)(x))
+#endif
+
+#define i64rint(x) ((int64_t)(x)) /* only needed for ld128 so not opt. */
+
+#if defined(__i386__)
+static __inline int
+irintf(float x)
+{
+ int n;
+
+ __asm("fistl %0" : "=m" (n) : "t" (x));
+ return (n);
+}
+
+static __inline int
+irintd(double x)
+{
+ int n;
+
+ __asm("fistl %0" : "=m" (n) : "t" (x));
+ return (n);
+}
+#endif
+
+#if defined(__amd64__) || defined(__i386__)
+static __inline int
+irintl(long double x)
+{
+ int n;
+
+ __asm("fistl %0" : "=m" (n) : "t" (x));
+ return (n);
+}
+#endif
+
+#ifdef DEBUG
+#if defined(__amd64__) || defined(__i386__)
+#define breakpoint() asm("int $3")
+#else
+#include <signal.h>
+
+#define breakpoint() raise(SIGTRAP)
+#endif
+#endif
+
+/* Write a pari script to test things externally. */
+#ifdef DOPRINT
+#include <stdio.h>
+
+#ifndef DOPRINT_SWIZZLE
+#define DOPRINT_SWIZZLE 0
+#endif
+
+#ifdef DOPRINT_LD80
+
+#define DOPRINT_START(xp) do { \
+ uint64_t __lx; \
+ uint16_t __hx; \
+ \
+ /* Hack to give more-problematic args. */ \
+ EXTRACT_LDBL80_WORDS(__hx, __lx, *xp); \
+ __lx ^= DOPRINT_SWIZZLE; \
+ INSERT_LDBL80_WORDS(*xp, __hx, __lx); \
+ printf("x = %.21Lg; ", (long double)*xp); \
+} while (0)
+#define DOPRINT_END1(v) \
+ printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v))
+#define DOPRINT_END2(hi, lo) \
+ printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n", \
+ (long double)(hi), (long double)(lo))
+
+#elif defined(DOPRINT_D64)
+
+#define DOPRINT_START(xp) do { \
+ uint32_t __hx, __lx; \
+ \
+ EXTRACT_WORDS(__hx, __lx, *xp); \
+ __lx ^= DOPRINT_SWIZZLE; \
+ INSERT_WORDS(*xp, __hx, __lx); \
+ printf("x = %.21Lg; ", (long double)*xp); \
+} while (0)
+#define DOPRINT_END1(v) \
+ printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v))
+#define DOPRINT_END2(hi, lo) \
+ printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n", \
+ (long double)(hi), (long double)(lo))
+
+#elif defined(DOPRINT_F32)
+
+#define DOPRINT_START(xp) do { \
+ uint32_t __hx; \
+ \
+ GET_FLOAT_WORD(__hx, *xp); \
+ __hx ^= DOPRINT_SWIZZLE; \
+ SET_FLOAT_WORD(*xp, __hx); \
+ printf("x = %.21Lg; ", (long double)*xp); \
+} while (0)
+#define DOPRINT_END1(v) \
+ printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v))
+#define DOPRINT_END2(hi, lo) \
+ printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n", \
+ (long double)(hi), (long double)(lo))
+
+#else /* !DOPRINT_LD80 && !DOPRINT_D64 (LD128 only) */
+
+#ifndef DOPRINT_SWIZZLE_HIGH
+#define DOPRINT_SWIZZLE_HIGH 0
+#endif
+
+#define DOPRINT_START(xp) do { \
+ uint64_t __lx, __llx; \
+ uint16_t __hx; \
+ \
+ EXTRACT_LDBL128_WORDS(__hx, __lx, __llx, *xp); \
+ __llx ^= DOPRINT_SWIZZLE; \
+ __lx ^= DOPRINT_SWIZZLE_HIGH; \
+ INSERT_LDBL128_WORDS(*xp, __hx, __lx, __llx); \
+ printf("x = %.36Lg; ", (long double)*xp); \
+} while (0)
+#define DOPRINT_END1(v) \
+ printf("y = %.36Lg; z = 0; show(x, y, z);\n", (long double)(v))
+#define DOPRINT_END2(hi, lo) \
+ printf("y = %.36Lg; z = %.36Lg; show(x, y, z);\n", \
+ (long double)(hi), (long double)(lo))
+
+#endif /* DOPRINT_LD80 */
+
+#else /* !DOPRINT */
+#define DOPRINT_START(xp)
+#define DOPRINT_END1(v)
+#define DOPRINT_END2(hi, lo)
+#endif /* DOPRINT */
+
+#define RETURNP(x) do { \
+ DOPRINT_END1(x); \
+ RETURNF(x); \
+} while (0)
+#define RETURNPI(x) do { \
+ DOPRINT_END1(x); \
+ RETURNI(x); \
+} while (0)
+#define RETURN2P(x, y) do { \
+ DOPRINT_END2((x), (y)); \
+ RETURNF((x) + (y)); \
+} while (0)
+#define RETURN2PI(x, y) do { \
+ DOPRINT_END2((x), (y)); \
+ RETURNI((x) + (y)); \
+} while (0)
+#ifdef STRUCT_RETURN
+#define RETURNSP(rp) do { \
+ if (!(rp)->lo_set) \
+ RETURNP((rp)->hi); \
+ RETURN2P((rp)->hi, (rp)->lo); \
+} while (0)
+#define RETURNSPI(rp) do { \
+ if (!(rp)->lo_set) \
+ RETURNPI((rp)->hi); \
+ RETURN2PI((rp)->hi, (rp)->lo); \
+} while (0)
+#endif
+#define SUM2P(x, y) ({ \
+ const __typeof (x) __x = (x); \
+ const __typeof (y) __y = (y); \
+ \
+ DOPRINT_END2(__x, __y); \
+ __x + __y; \
+})
+
+/*
+ * ieee style elementary functions
+ *
+ * We rename functions here to improve other sources' diffability
+ * against fdlibm.
+ */
+#define __ieee754_sqrt sqrt
+#define __ieee754_acos acos
+#define __ieee754_acosh acosh
+#define __ieee754_log log
+#define __ieee754_log2 log2
+#define __ieee754_atanh atanh
+#define __ieee754_asin asin
+#define __ieee754_atan2 atan2
+#define __ieee754_exp exp
+#define __ieee754_cosh cosh
+#define __ieee754_fmod fmod
+#define __ieee754_pow pow
+#define __ieee754_lgamma lgamma
+#define __ieee754_gamma gamma
+#define __ieee754_lgamma_r lgamma_r
+#define __ieee754_gamma_r gamma_r
+#define __ieee754_log10 log10
+#define __ieee754_sinh sinh
+#define __ieee754_hypot hypot
+#define __ieee754_j0 j0
+#define __ieee754_j1 j1
+#define __ieee754_y0 y0
+#define __ieee754_y1 y1
+#define __ieee754_jn jn
+#define __ieee754_yn yn
+#define __ieee754_remainder remainder
+#define __ieee754_scalb scalb
+#define __ieee754_sqrtf sqrtf
+#define __ieee754_acosf acosf
+#define __ieee754_acoshf acoshf
+#define __ieee754_logf logf
+#define __ieee754_atanhf atanhf
+#define __ieee754_asinf asinf
+#define __ieee754_atan2f atan2f
+#define __ieee754_expf expf
+#define __ieee754_coshf coshf
+#define __ieee754_fmodf fmodf
+#define __ieee754_powf powf
+#define __ieee754_lgammaf lgammaf
+#define __ieee754_gammaf gammaf
+#define __ieee754_lgammaf_r lgammaf_r
+#define __ieee754_gammaf_r gammaf_r
+#define __ieee754_log10f log10f
+#define __ieee754_log2f log2f
+#define __ieee754_sinhf sinhf
+#define __ieee754_hypotf hypotf
+#define __ieee754_j0f j0f
+#define __ieee754_j1f j1f
+#define __ieee754_y0f y0f
+#define __ieee754_y1f y1f
+#define __ieee754_jnf jnf
+#define __ieee754_ynf ynf
+#define __ieee754_remainderf remainderf
+#define __ieee754_scalbf scalbf
+
+/* fdlibm kernel function */
+int __kernel_rem_pio2(double*,double*,int,int,int);
+
+/* double precision kernel functions */
+#ifndef INLINE_REM_PIO2
+int __ieee754_rem_pio2(double,double*);
+#endif
+double __kernel_sin(double,double,int);
+double __kernel_cos(double,double);
+double __kernel_tan(double,double,int);
+double __ldexp_exp(double,int);
+#ifdef _COMPLEX_H
+double complex __ldexp_cexp(double complex,int);
+#endif
+
+/* float precision kernel functions */
+#ifndef INLINE_REM_PIO2F
+int __ieee754_rem_pio2f(float,double*);
+#endif
+#ifndef INLINE_KERNEL_SINDF
+float __kernel_sindf(double);
+#endif
+#ifndef INLINE_KERNEL_COSDF
+float __kernel_cosdf(double);
+#endif
+#ifndef INLINE_KERNEL_TANDF
+float __kernel_tandf(double,int);
+#endif
+float __ldexp_expf(float,int);
+#ifdef _COMPLEX_H
+float complex __ldexp_cexpf(float complex,int);
+#endif
+
+/* long double precision kernel functions */
+long double __kernel_sinl(long double, long double, int);
+long double __kernel_cosl(long double, long double);
+long double __kernel_tanl(long double, long double, int);
+
+#endif /* !_MATH_PRIVATE_H_ */
diff --git a/newlib/libm/ld/s_asinhl.c b/newlib/libm/ld/s_asinhl.c
new file mode 100644
index 000000000..ba28f599c
--- /dev/null
+++ b/newlib/libm/ld/s_asinhl.c
@@ -0,0 +1,91 @@
+/* from: FreeBSD: head/lib/msun/src/e_acosh.c 176451 2008-02-22 02:30:36Z das */
+
+/* @(#)s_asinh.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/*
+ * See s_asinh.c for complete comments.
+ *
+ * Converted to long double by David Schultz <das@FreeBSD.ORG> and
+ * Bruce D. Evans.
+ */
+
+#include <float.h>
+#ifdef __i386__
+#include <ieeefp.h>
+#endif
+
+#include "fpmath.h"
+#include "math.h"
+#include "math_private.h"
+
+/* EXP_LARGE is the threshold above which we use asinh(x) ~= log(2x). */
+/* EXP_TINY is the threshold below which we use asinh(x) ~= x. */
+#if LDBL_MANT_DIG == 64
+#define EXP_LARGE 34
+#define EXP_TINY -34
+#elif LDBL_MANT_DIG == 113
+#define EXP_LARGE 58
+#define EXP_TINY -58
+#else
+#error "Unsupported long double format"
+#endif
+
+#if LDBL_MAX_EXP != 0x4000
+/* We also require the usual expsign encoding. */
+#error "Unsupported long double format"
+#endif
+
+#define BIAS (LDBL_MAX_EXP - 1)
+
+static const double
+one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
+huge= 1.00000000000000000000e+300;
+
+#if LDBL_MANT_DIG == 64
+static const union IEEEl2bits
+u_ln2 = LD80C(0xb17217f7d1cf79ac, -1, 6.93147180559945309417e-1L);
+#define ln2 u_ln2.e
+#elif LDBL_MANT_DIG == 113
+static const long double
+ln2 = 6.93147180559945309417232121458176568e-1L; /* 0x162e42fefa39ef35793c7673007e6.0p-113 */
+#else
+#error "Unsupported long double format"
+#endif
+
+long double
+asinhl(long double x)
+{
+ long double t, w;
+ uint16_t hx, ix;
+
+ ENTERI();
+ GET_LDBL_EXPSIGN(hx, x);
+ ix = hx & 0x7fff;
+ if (ix >= 0x7fff) RETURNI(x+x); /* x is inf, NaN or misnormal */
+ if (ix < BIAS + EXP_TINY) { /* |x| < TINY, or misnormal */
+ if (huge + x > one) RETURNI(x); /* return x inexact except 0 */
+ }
+ if (ix >= BIAS + EXP_LARGE) { /* |x| >= LARGE, or misnormal */
+ w = logl(fabsl(x))+ln2;
+ } else if (ix >= 0x4000) { /* LARGE > |x| >= 2.0, or misnormal */
+ t = fabsl(x);
+ w = logl(2.0*t+one/(sqrtl(x*x+one)+t));
+ } else { /* 2.0 > |x| >= TINY, or misnormal */
+ t = x*x;
+ w =log1pl(fabsl(x)+t/(one+sqrtl(one+t)));
+ }
+ RETURNI((hx & 0x8000) == 0 ? w : -w);
+}
diff --git a/newlib/libm/ld/s_atanl.c b/newlib/libm/ld/s_atanl.c
new file mode 100644
index 000000000..ff29c3ce8
--- /dev/null
+++ b/newlib/libm/ld/s_atanl.c
@@ -0,0 +1,85 @@
+/* @(#)s_atan.c 5.1 93/09/24 */
+/* FreeBSD: head/lib/msun/src/s_atan.c 176451 2008-02-22 02:30:36Z das */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/*
+ * See comments in s_atan.c.
+ * Converted to long double by David Schultz <das@FreeBSD.ORG>.
+ */
+
+#include <float.h>
+
+#include "invtrig.h"
+#include "math.h"
+#include "math_private.h"
+
+static const long double
+one = 1.0,
+huge = 1.0e300;
+
+long double
+atanl(long double x)
+{
+ union IEEEl2bits u;
+ long double w,s1,s2,z;
+ int id;
+ int16_t expsign, expt;
+ int32_t expman;
+
+ u.e = x;
+ expsign = u.xbits.expsign;
+ expt = expsign & 0x7fff;
+ if(expt >= ATAN_CONST) { /* if |x| is large, atan(x)~=pi/2 */
+ if(expt == BIAS + LDBL_MAX_EXP &&
+ ((u.bits.manh&~LDBL_NBIT)|u.bits.manl)!=0)
+ return x+x; /* NaN */
+ if(expsign>0) return atanhi[3]+atanlo[3];
+ else return -atanhi[3]-atanlo[3];
+ }
+ /* Extract the exponent and the first few bits of the mantissa. */
+ /* XXX There should be a more convenient way to do this. */
+ expman = (expt << 8) | ((u.bits.manh >> (MANH_SIZE - 9)) & 0xff);
+ if (expman < ((BIAS - 2) << 8) + 0xc0) { /* |x| < 0.4375 */
+ if (expt < ATAN_LINEAR) { /* if |x| is small, atanl(x)~=x */
+ if(huge+x>one) return x; /* raise inexact */
+ }
+ id = -1;
+ } else {
+ x = fabsl(x);
+ if (expman < (BIAS << 8) + 0x30) { /* |x| < 1.1875 */
+ if (expman < ((BIAS - 1) << 8) + 0x60) { /* 7/16 <=|x|<11/16 */
+ id = 0; x = (2.0*x-one)/(2.0+x);
+ } else { /* 11/16<=|x|< 19/16 */
+ id = 1; x = (x-one)/(x+one);
+ }
+ } else {
+ if (expman < ((BIAS + 1) << 8) + 0x38) { /* |x| < 2.4375 */
+ id = 2; x = (x-1.5)/(one+1.5*x);
+ } else { /* 2.4375 <= |x| < 2^ATAN_CONST */
+ id = 3; x = -1.0/x;
+ }
+ }}
+ /* end of argument reduction */
+ z = x*x;
+ w = z*z;
+ /* break sum aT[i]z**(i+1) into odd and even poly */
+ s1 = z*T_even(w);
+ s2 = w*T_odd(w);
+ if (id<0) return x - x*(s1+s2);
+ else {
+ z = atanhi[id] - ((x*(s1+s2) - atanlo[id]) - x);
+ return (expsign<0)? -z:z;
+ }
+}
diff --git a/newlib/libm/ld/s_cbrtl.c b/newlib/libm/ld/s_cbrtl.c
new file mode 100644
index 000000000..5c4a98a87
--- /dev/null
+++ b/newlib/libm/ld/s_cbrtl.c
@@ -0,0 +1,143 @@
+/*-
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ * Copyright (c) 2009-2011, Bruce D. Evans, Steven G. Kargl, David Schultz.
+ *
+ * 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.
+ * ====================================================
+ *
+ * The argument reduction and testing for exceptional cases was
+ * written by Steven G. Kargl with input from Bruce D. Evans
+ * and David A. Schultz.
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+#include <float.h>
+#ifdef __i386__
+#include <ieeefp.h>
+#endif
+
+#include "fpmath.h"
+#include "math.h"
+#include "math_private.h"
+
+#define BIAS (LDBL_MAX_EXP - 1)
+
+static const unsigned
+ B1 = 709958130; /* B1 = (127-127.0/3-0.03306235651)*2**23 */
+
+long double
+cbrtl(long double x)
+{
+ union IEEEl2bits u, v;
+ long double r, s, t, w;
+ double dr, dt, dx;
+ float ft, fx;
+ uint32_t hx;
+ uint16_t expsign;
+ int k;
+
+ u.e = x;
+ expsign = u.xbits.expsign;
+ k = expsign & 0x7fff;
+
+ /*
+ * If x = +-Inf, then cbrt(x) = +-Inf.
+ * If x = NaN, then cbrt(x) = NaN.
+ */
+ if (k == BIAS + LDBL_MAX_EXP)
+ return (x + x);
+
+ ENTERI();
+ if (k == 0) {
+ /* If x = +-0, then cbrt(x) = +-0. */
+ if ((u.bits.manh | u.bits.manl) == 0)
+ RETURNI(x);
+ /* Adjust subnormal numbers. */
+ u.e *= 0x1.0p514;
+ k = u.bits.exp;
+ k -= BIAS + 514;
+ } else
+ k -= BIAS;
+ u.xbits.expsign = BIAS;
+ v.e = 1;
+
+ x = u.e;
+ switch (k % 3) {
+ case 1:
+ case -2:
+ x = 2*x;
+ k--;
+ break;
+ case 2:
+ case -1:
+ x = 4*x;
+ k -= 2;
+ break;
+ }
+ v.xbits.expsign = (expsign & 0x8000) | (BIAS + k / 3);
+
+ /*
+ * The following is the guts of s_cbrtf, with the handling of
+ * special values removed and extra care for accuracy not taken,
+ * but with most of the extra accuracy not discarded.
+ */
+
+ /* ~5-bit estimate: */
+ fx = x;
+ GET_FLOAT_WORD(hx, fx);
+ SET_FLOAT_WORD(ft, ((hx & 0x7fffffff) / 3 + B1));
+
+ /* ~16-bit estimate: */
+ dx = x;
+ dt = ft;
+ dr = dt * dt * dt;
+ dt = dt * (dx + dx + dr) / (dx + dr + dr);
+
+ /* ~47-bit estimate: */
+ dr = dt * dt * dt;
+ dt = dt * (dx + dx + dr) / (dx + dr + dr);
+
+#if LDBL_MANT_DIG == 64
+ /*
+ * dt is cbrtl(x) to ~47 bits (after x has been reduced to 1 <= x < 8).
+ * Round it away from zero to 32 bits (32 so that t*t is exact, and
+ * away from zero for technical reasons).
+ */
+ volatile double vd2 = 0x1.0p32;
+ volatile double vd1 = 0x1.0p-31;
+ #define vd ((long double)vd2 + vd1)
+
+ t = dt + vd - 0x1.0p32;
+#elif LDBL_MANT_DIG == 113
+ /*
+ * Round dt away from zero to 47 bits. Since we don't trust the 47,
+ * add 2 47-bit ulps instead of 1 to round up. Rounding is slow and
+ * might be avoidable in this case, since on most machines dt will
+ * have been evaluated in 53-bit precision and the technical reasons
+ * for rounding up might not apply to either case in cbrtl() since
+ * dt is much more accurate than needed.
+ */
+ t = dt + 0x2.0p-46 + 0x1.0p60L - 0x1.0p60;
+#else
+#error "Unsupported long double format"
+#endif
+
+ /*
+ * Final step Newton iteration to 64 or 113 bits with
+ * error < 0.667 ulps
+ */
+ s=t*t; /* t*t is exact */
+ r=x/s; /* error <= 0.5 ulps; |r| < |t| */
+ w=t+t; /* t+t is exact */
+ r=(r-t)/(w+r); /* r-t is exact; w+r ~= 3*t */
+ t=t+t*r; /* error <= (0.5 + 0.5/3) * ulp */
+
+ t *= v.e;
+ RETURNI(t);
+}
diff --git a/newlib/libm/ld/s_ceill.c b/newlib/libm/ld/s_ceill.c
new file mode 100644
index 000000000..2d1045fe6
--- /dev/null
+++ b/newlib/libm/ld/s_ceill.c
@@ -0,0 +1,101 @@
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ *
+ * From: @(#)s_ceil.c 5.1 93/09/24
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/*
+ * ceill(x)
+ * Return x rounded toward -inf to integral value
+ * Method:
+ * Bit twiddling.
+ * Exception:
+ * Inexact flag raised if x not equal to ceill(x).
+ */
+
+#include <float.h>
+#include <math.h>
+#include <stdint.h>
+
+#include "fpmath.h"
+
+#ifdef LDBL_IMPLICIT_NBIT
+#define MANH_SIZE (LDBL_MANH_SIZE + 1)
+#define INC_MANH(u, c) do { \
+ uint64_t o = u.bits.manh; \
+ u.bits.manh += (c); \
+ if (u.bits.manh < o) \
+ u.bits.exp++; \
+} while (0)
+#else
+#define MANH_SIZE LDBL_MANH_SIZE
+#define INC_MANH(u, c) do { \
+ uint64_t o = u.bits.manh; \
+ u.bits.manh += (c); \
+ if (u.bits.manh < o) { \
+ u.bits.exp++; \
+ u.bits.manh |= 1llu << (LDBL_MANH_SIZE - 1); \
+ } \
+} while (0)
+#endif
+
+static const long double huge = 1.0e300;
+
+long double
+ceill(long double x)
+{
+ union IEEEl2bits u = { .e = x };
+ int e = u.bits.exp - LDBL_MAX_EXP + 1;
+
+ if (e < MANH_SIZE - 1) {
+ if (e < 0) { /* raise inexact if x != 0 */
+ if (huge + x > 0.0)
+ if (u.bits.exp > 0 ||
+ (u.bits.manh | u.bits.manl) != 0)
+ u.e = u.bits.sign ? -0.0 : 1.0;
+ } else {
+ uint64_t m = ((1llu << MANH_SIZE) - 1) >> (e + 1);
+ if (((u.bits.manh & m) | u.bits.manl) == 0)
+ return (x); /* x is integral */
+ if (!u.bits.sign) {
+#ifdef LDBL_IMPLICIT_NBIT
+ if (e == 0)
+ u.bits.exp++;
+ else
+#endif
+ INC_MANH(u, 1llu << (MANH_SIZE - e - 1));
+ }
+ if (huge + x > 0.0) { /* raise inexact flag */
+ u.bits.manh &= ~m;
+ u.bits.manl = 0;
+ }
+ }
+ } else if (e < LDBL_MANT_DIG - 1) {
+ uint64_t m = (uint64_t)-1 >> (64 - LDBL_MANT_DIG + e + 1);
+ if ((u.bits.manl & m) == 0)
+ return (x); /* x is integral */
+ if (!u.bits.sign) {
+ if (e == MANH_SIZE - 1)
+ INC_MANH(u, 1);
+ else {
+ uint64_t o = u.bits.manl;
+ u.bits.manl += 1llu << (LDBL_MANT_DIG - e - 1);
+ if (u.bits.manl < o) /* got a carry */
+ INC_MANH(u, 1);
+ }
+ }
+ if (huge + x > 0.0) /* raise inexact flag */
+ u.bits.manl &= ~m;
+ }
+ return (u.e);
+}
diff --git a/newlib/libm/ld/s_copysignl.c b/newlib/libm/ld/s_copysignl.c
new file mode 100644
index 000000000..bd6744705
--- /dev/null
+++ b/newlib/libm/ld/s_copysignl.c
@@ -0,0 +1,44 @@
+/*-
+ * SPDX-License-Identifier: BSD-2-Clause-FreeBSD
+ *
+ * Copyright (c) 2004 Stefan Farfeleder
+ * 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 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 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.
+ *
+ * $FreeBSD$
+ */
+
+#include <math.h>
+
+#include "fpmath.h"
+
+long double
+copysignl(long double x, long double y)
+{
+ union IEEEl2bits ux, uy;
+
+ ux.e = x;
+ uy.e = y;
+ ux.bits.sign = uy.bits.sign;
+ return (ux.e);
+}
diff --git a/newlib/libm/ld/s_cosl.c b/newlib/libm/ld/s_cosl.c
new file mode 100644
index 000000000..061de295a
--- /dev/null
+++ b/newlib/libm/ld/s_cosl.c
@@ -0,0 +1,102 @@
+/*-
+ * SPDX-License-Identifier: BSD-2-Clause-FreeBSD
+ *
+ * Copyright (c) 2007 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.
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/*
+ * Limited testing on pseudorandom numbers drawn within [-2e8:4e8] shows
+ * an accuracy of <= 0.7412 ULP.
+ */
+
+#include <float.h>
+#ifdef __i386__
+#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
+
+long double
+cosl(long double x)
+{
+ union IEEEl2bits z;
+ int e0;
+ long double y[2];
+ long double hi, lo;
+
+ z.e = x;
+ z.bits.sign = 0;
+
+ /* If x = +-0 or x is a subnormal number, then cos(x) = 1 */
+ if (z.bits.exp == 0)
+ return (1.0);
+
+ /* If x = NaN or Inf, then cos(x) = NaN. */
+ if (z.bits.exp == 32767)
+ return ((x - x) / (x - x));
+
+ ENTERI();
+
+ /* Optimize the case where x is already within range. */
+ if (z.e < pio4)
+ RETURNI(__kernel_cosl(z.e, 0));
+
+ e0 = __ieee754_rem_pio2l(x, y);
+ hi = y[0];
+ lo = y[1];
+
+ switch (e0 & 3) {
+ case 0:
+ hi = __kernel_cosl(hi, lo);
+ break;
+ case 1:
+ hi = - __kernel_sinl(hi, lo, 1);
+ break;
+ case 2:
+ hi = - __kernel_cosl(hi, lo);
+ break;
+ case 3:
+ hi = __kernel_sinl(hi, lo, 1);
+ break;
+ }
+
+ RETURNI(hi);
+}
diff --git a/newlib/libm/ld/s_fabsl.c b/newlib/libm/ld/s_fabsl.c
new file mode 100644
index 000000000..5076d8a9b
--- /dev/null
+++ b/newlib/libm/ld/s_fabsl.c
@@ -0,0 +1,45 @@
+/*-
+ * SPDX-License-Identifier: BSD-3-Clause
+ *
+ * Copyright (c) 2003 Dag-Erling Smørgrav
+ * 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
+ * in this position and unchanged.
+ * 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.
+ * 3. The name of the author may not be used to endorse or promote products
+ * derived from this software without specific prior written permission.
+ *
+ * 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.
+ *
+ * $FreeBSD$
+ */
+
+#include <math.h>
+
+#include "fpmath.h"
+
+long double
+fabsl(long double x)
+{
+ union IEEEl2bits u;
+
+ u.e = x;
+ u.bits.sign = 0;
+ return (u.e);
+}
diff --git a/newlib/libm/ld/s_fdim.c b/newlib/libm/ld/s_fdim.c
new file mode 100644
index 000000000..c40c3e9d3
--- /dev/null
+++ b/newlib/libm/ld/s_fdim.c
@@ -0,0 +1,48 @@
+/*-
+ * SPDX-License-Identifier: BSD-2-Clause-FreeBSD
+ *
+ * Copyright (c) 2004 David Schultz <das@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.
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+#include <math.h>
+
+#define DECL(type, fn) \
+type \
+fn(type x, type y) \
+{ \
+ \
+ if (isnan(x)) \
+ return (x); \
+ if (isnan(y)) \
+ return (y); \
+ return (x > y ? x - y : 0.0); \
+}
+
+DECL(double, fdim)
+DECL(float, fdimf)
+DECL(long double, fdiml)
diff --git a/newlib/libm/ld/s_floorl.c b/newlib/libm/ld/s_floorl.c
new file mode 100644
index 000000000..6cec3e781
--- /dev/null
+++ b/newlib/libm/ld/s_floorl.c
@@ -0,0 +1,101 @@
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ *
+ * From: @(#)s_floor.c 5.1 93/09/24
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/*
+ * floorl(x)
+ * Return x rounded toward -inf to integral value
+ * Method:
+ * Bit twiddling.
+ * Exception:
+ * Inexact flag raised if x not equal to floorl(x).
+ */
+
+#include <float.h>
+#include <math.h>
+#include <stdint.h>
+
+#include "fpmath.h"
+
+#ifdef LDBL_IMPLICIT_NBIT
+#define MANH_SIZE (LDBL_MANH_SIZE + 1)
+#define INC_MANH(u, c) do { \
+ uint64_t o = u.bits.manh; \
+ u.bits.manh += (c); \
+ if (u.bits.manh < o) \
+ u.bits.exp++; \
+} while (0)
+#else
+#define MANH_SIZE LDBL_MANH_SIZE
+#define INC_MANH(u, c) do { \
+ uint64_t o = u.bits.manh; \
+ u.bits.manh += (c); \
+ if (u.bits.manh < o) { \
+ u.bits.exp++; \
+ u.bits.manh |= 1llu << (LDBL_MANH_SIZE - 1); \
+ } \
+} while (0)
+#endif
+
+static const long double huge = 1.0e300;
+
+long double
+floorl(long double x)
+{
+ union IEEEl2bits u = { .e = x };
+ int e = u.bits.exp - LDBL_MAX_EXP + 1;
+
+ if (e < MANH_SIZE - 1) {
+ if (e < 0) { /* raise inexact if x != 0 */
+ if (huge + x > 0.0)
+ if (u.bits.exp > 0 ||
+ (u.bits.manh | u.bits.manl) != 0)
+ u.e = u.bits.sign ? -1.0 : 0.0;
+ } else {
+ uint64_t m = ((1llu << MANH_SIZE) - 1) >> (e + 1);
+ if (((u.bits.manh & m) | u.bits.manl) == 0)
+ return (x); /* x is integral */
+ if (u.bits.sign) {
+#ifdef LDBL_IMPLICIT_NBIT
+ if (e == 0)
+ u.bits.exp++;
+ else
+#endif
+ INC_MANH(u, 1llu << (MANH_SIZE - e - 1));
+ }
+ if (huge + x > 0.0) { /* raise inexact flag */
+ u.bits.manh &= ~m;
+ u.bits.manl = 0;
+ }
+ }
+ } else if (e < LDBL_MANT_DIG - 1) {
+ uint64_t m = (uint64_t)-1 >> (64 - LDBL_MANT_DIG + e + 1);
+ if ((u.bits.manl & m) == 0)
+ return (x); /* x is integral */
+ if (u.bits.sign) {
+ if (e == MANH_SIZE - 1)
+ INC_MANH(u, 1);
+ else {
+ uint64_t o = u.bits.manl;
+ u.bits.manl += 1llu << (LDBL_MANT_DIG - e - 1);
+ if (u.bits.manl < o) /* got a carry */
+ INC_MANH(u, 1);
+ }
+ }
+ if (huge + x > 0.0) /* raise inexact flag */
+ u.bits.manl &= ~m;
+ }
+ return (u.e);
+}
diff --git a/newlib/libm/ld/s_fmal.c b/newlib/libm/ld/s_fmal.c
new file mode 100644
index 000000000..a379346c1
--- /dev/null
+++ b/newlib/libm/ld/s_fmal.c
@@ -0,0 +1,274 @@
+/*-
+ * SPDX-License-Identifier: BSD-2-Clause-FreeBSD
+ *
+ * Copyright (c) 2005-2011 David Schultz <das@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, STRI[...]
[diff truncated at 100000 bytes]
^ permalink raw reply [flat|nested] only message in thread
only message in thread, other threads:[~2023-05-16 23:35 UTC | newest]
Thread overview: (only message) (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
2023-05-16 23:35 [newlib-cygwin] newlib: Add FreeBSD files for non LDBL_EQ_DBL support Joel Sherrill
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).