From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: by sourceware.org (Postfix, from userid 1129) id 781B23858D1E; Tue, 16 May 2023 23:35:09 +0000 (GMT) DKIM-Filter: OpenDKIM Filter v2.11.0 sourceware.org 781B23858D1E DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=sourceware.org; s=default; t=1684280109; bh=TI3z5XGnTZVYRi5Bl8pkWBdA8PN3U0zvVCM+9HDMRd8=; h=From:To:Subject:Date:From; b=Pm8XaqJjkPP+e83fyADM04jGZCBQJwtDOZTVCyOTcyFQg8UYVZ4wa03CeHRYN4Q16 G34l+MEMZ+ORd3ahsADAhrCNAvel6Ld6jl3JJ1O1zEcm8Lb0fFJm8CmzQluoxi46QM tZsEDJK/n27awwalfdwvgq/RCJL/h1kkXDOCgwn0= Content-Type: text/plain; charset="utf-8" MIME-Version: 1.0 Content-Transfer-Encoding: quoted-printable From: Joel Sherrill To: newlib-cvs@sourceware.org Subject: [newlib-cygwin] newlib: Add FreeBSD files for non LDBL_EQ_DBL support X-Act-Checkin: newlib-cygwin X-Git-Author: Jennifer Averett X-Git-Refname: refs/heads/master X-Git-Oldrev: 41fdb869f9984ca2f8600aaa980b4ed2eae2e980 X-Git-Newrev: c630a6a837fd581d741fb94602aed6a4a5ed9bf4 Message-Id: <20230516233509.781B23858D1E@sourceware.org> Date: Tue, 16 May 2023 23:35:09 +0000 (GMT) List-Id: https://sourceware.org/git/gitweb.cgi?p=3Dnewlib-cygwin.git;h=3Dc630a6a837f= d581d741fb94602aed6a4a5ed9bf4 commit c630a6a837fd581d741fb94602aed6a4a5ed9bf4 Author: Jennifer Averett Date: Fri May 5 13:39:16 2023 -0500 newlib: Add FreeBSD files for non LDBL_EQ_DBL support =20 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/end= ian.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 + * 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 PURP= OSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENT= IAL + * 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, STR= ICT + * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY W= AY + * 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 +#include +#include + +#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 litt= le + * endian to host byte order functions as detailed in byteorder(9). + */ +#if _BYTE_ORDER =3D=3D _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 !=3D _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 =3D=3D _LITTLE_ENDIAN */ + +/* Alignment-agnostic encode/decode bytestream to/from little/big endian. = */ + +static __inline uint16_t +be16dec(const void *pp) +{ + uint8_t const *p =3D (uint8_t const *)pp; + + return (((unsigned)p[0] << 8) | p[1]); +} + +static __inline uint32_t +be32dec(const void *pp) +{ + uint8_t const *p =3D (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 =3D (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 =3D (uint8_t const *)pp; + + return (((unsigned)p[1] << 8) | p[0]); +} + +static __inline uint32_t +le32dec(const void *pp) +{ + uint8_t const *p =3D (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 =3D (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 =3D (uint8_t *)pp; + + p[0] =3D (u >> 8) & 0xff; + p[1] =3D u & 0xff; +} + +static __inline void +be32enc(void *pp, uint32_t u) +{ + uint8_t *p =3D (uint8_t *)pp; + + p[0] =3D (u >> 24) & 0xff; + p[1] =3D (u >> 16) & 0xff; + p[2] =3D (u >> 8) & 0xff; + p[3] =3D u & 0xff; +} + +static __inline void +be64enc(void *pp, uint64_t u) +{ + uint8_t *p =3D (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 =3D (uint8_t *)pp; + + p[0] =3D u & 0xff; + p[1] =3D (u >> 8) & 0xff; +} + +static __inline void +le32enc(void *pp, uint32_t u) +{ + uint8_t *p =3D (uint8_t *)pp; + + p[0] =3D u & 0xff; + p[1] =3D (u >> 8) & 0xff; + p[2] =3D (u >> 16) & 0xff; + p[3] =3D (u >> 24) & 0xff; +} + +static __inline void +le64enc(void *pp, uint64_t u) +{ + uint8_t *p =3D (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/ma= chine/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 + * 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 PURP= OSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENT= IAL + * 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, STR= ICT + * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY W= AY + * 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] =3D (uint32_t)(u).bits.manl; \ + (a)[1] =3D (uint32_t)((u).bits.manl >> 32); \ + (a)[2] =3D (uint32_t)(u).bits.manh; \ + (a)[3] =3D (uint32_t)((u).bits.manh >> 32); \ +} while(0) diff --git a/newlib/libc/machine/i386/machine/_fpmath.h b/newlib/libc/machi= ne/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 + * 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 PURP= OSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENT= IAL + * 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, STR= ICT + * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY W= AY + * 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 &=3D ~LDBL_NBIT) + +#define LDBL_MANH_SIZE 32 +#define LDBL_MANL_SIZE 32 + +#define LDBL_TO_ARRAY32(u, a) do { \ + (a)[0] =3D (uint32_t)(u).bits.manl; \ + (a)[1] =3D (uint32_t)(u).bits.manh; \ +} while (0) diff --git a/newlib/libc/machine/x86_64/machine/_fpmath.h b/newlib/libc/mac= hine/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 + * 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 PURP= OSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENT= IAL + * 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, STR= ICT + * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY W= AY + * 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 &=3D ~LDBL_NBIT) + +#define LDBL_MANH_SIZE 32 +#define LDBL_MANL_SIZE 32 + +#define LDBL_TO_ARRAY32(u, a) do { \ + (a)[0] =3D (uint32_t)(u).bits.manl; \ + (a)[1] =3D (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 */ +/* + * =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D + * 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. + * =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D + * + */ + +#include +__FBSDID("$FreeBSD$"); + +/* + * See e_acosh.c for complete comments. + * + * Converted to long double by David Schultz and + * Bruce D. Evans. + */ + +#include +#ifdef __i386__ +#include +#endif + +#include "fpmath.h" +#include "math.h" +#include "math_private.h" + +/* EXP_LARGE is the threshold above which we use acosh(x) ~=3D log(2x). */ +#if LDBL_MANT_DIG =3D=3D 64 +#define EXP_LARGE 34 +#elif LDBL_MANT_DIG =3D=3D 113 +#define EXP_LARGE 58 +#else +#error "Unsupported long double format" +#endif + +#if LDBL_MAX_EXP !=3D 0x4000 +/* We also require the usual expsign encoding. */ +#error "Unsupported long double format" +#endif + +#define BIAS (LDBL_MAX_EXP - 1) + +static const double +one =3D 1.0; + +#if LDBL_MANT_DIG =3D=3D 64 +static const union IEEEl2bits +u_ln2 =3D LD80C(0xb17217f7d1cf79ac, -1, 6.93147180559945309417e-1L); +#define ln2 u_ln2.e +#elif LDBL_MANT_DIG =3D=3D 113 +static const long double +ln2 =3D 6.93147180559945309417232121458176568e-1L; /* 0x162e42fefa39ef357= 93c7673007e6.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 >=3D BIAS + EXP_LARGE) { /* x >=3D LARGE */ + if (hx >=3D 0x7fff) { /* x is inf, NaN or misnormal */ + RETURNI(x+x); + } else + RETURNI(logl(x)+ln2); /* acosh(huge)=3Dlog(2x), or misnormal */ + } else if (hx =3D=3D 0x3fff && x =3D=3D 1) { + RETURNI(0.0); /* acosh(1) =3D 0 */ + } else if (hx >=3D 0x4000) { /* LARGE > x >=3D 2, or misnormal */ + t=3Dx*x; + RETURNI(logl(2.0*x-one/(x+sqrtl(t-one)))); + } else { /* 1 +__FBSDID("$FreeBSD$"); + +/* + * See comments in e_acos.c. + * Converted to long double by David Schultz . + */ + +#include + +#include "invtrig.h" +#include "math.h" +#include "math_private.h" + +static const long double +one=3D 1.00000000000000000000e+00; + +#ifdef __i386__ +/* XXX Work around the fact that gcc truncates long double constants on i3= 86 */ +static volatile double +pi1 =3D 3.14159265358979311600e+00, /* 0x1.921fb54442d18p+1 */ +pi2 =3D 1.22514845490862001043e-16; /* 0x1.1a80000000000p-53 */ +#define pi ((long double)pi1 + pi2) +#else +static const long double +pi =3D 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 =3D x; + expsign =3D u.xbits.expsign; + expt =3D expsign & 0x7fff; + if(expt >=3D BIAS) { /* |x| >=3D 1 */ + if(expt=3D=3DBIAS && ((u.bits.manh&~LDBL_NBIT)|u.bits.manl)=3D=3D0) { + if (expsign>0) return 0.0; /* acos(1) =3D 0 */ + else return pi+2.0*pio2_lo; /* acos(-1)=3D pi */ + } + return (x-x)/(x-x); /* acos(|x|>1) is NaN */ + } + if(expt 0.5 */ + z =3D (one-x)*0.5; + s =3D sqrtl(z); + u.e =3D s; + u.bits.manl =3D 0; + df =3D u.e; + c =3D (z-df*df)/(s+df); + p =3D P(z); + q =3D Q(z); + r =3D p/q; + w =3D 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 */ +/* + * =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D + * 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. + * =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D + */ + +#include +__FBSDID("$FreeBSD$"); + +/* + * See comments in e_asin.c. + * Converted to long double by David Schultz . + */ + +#include + +#include "invtrig.h" +#include "math.h" +#include "math_private.h" + +static const long double +one =3D 1.00000000000000000000e+00, +huge =3D 1.000e+300; + +long double +asinl(long double x) +{ + union IEEEl2bits u; + long double t=3D0.0,w,p,q,c,r,s; + int16_t expsign, expt; + u.e =3D x; + expsign =3D u.xbits.expsign; + expt =3D expsign & 0x7fff; + if(expt >=3D BIAS) { /* |x|>=3D 1 */ + if(expt=3D=3DBIAS && ((u.bits.manh&~LDBL_NBIT)|u.bits.manl)=3D=3D0) + /* asin(1)=3D+-pi/2 with inexact */ + return x*pio2_hi+x*pio2_lo; + return (x-x)/(x-x); /* asin(|x|>1) is NaN */ + } else if (exptone) return x;/* return x with inexact if x!=3D0*/ + } + t =3D x*x; + p =3D P(t); + q =3D Q(t); + w =3D p/q; + return x+x*w; + } + /* 1> |x|>=3D 0.5 */ + w =3D one-fabsl(x); + t =3D w*0.5; + p =3D P(t); + q =3D Q(t); + s =3D sqrtl(t); + if(u.bits.manh>=3DTHRESH) { /* if |x| is close to 1 */ + w =3D p/q; + t =3D pio2_hi-(2.0*(s+s*w)-pio2_lo); + } else { + u.e =3D s; + u.bits.manl =3D 0; + w =3D u.e; + c =3D (t-w*w)/(s+w); + r =3D p/q; + p =3D 2.0*s*r-(pio2_lo-2.0*c); + q =3D pio4_hi-2.0*w; + t =3D 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 */ +/* + * =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D + * 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. + * =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D + * + */ + +#include +__FBSDID("$FreeBSD$"); + +/* + * See comments in e_atan2.c. + * Converted to long double by David Schultz . + */ + +#include + +#include "invtrig.h" +#include "math.h" +#include "math_private.h" + +static volatile long double +tiny =3D 1.0e-300; +static const long double +zero =3D 0.0; + +#ifdef __i386__ +/* XXX Work around the fact that gcc truncates long double constants on i3= 86 */ +static volatile double +pi1 =3D 3.14159265358979311600e+00, /* 0x1.921fb54442d18p+1 */ +pi2 =3D 1.22514845490862001043e-16; /* 0x1.1a80000000000p-53 */ +#define pi ((long double)pi1 + pi2) +#else +static const long double +pi =3D 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 =3D y; + expsigny =3D uy.xbits.expsign; + expty =3D expsigny & 0x7fff; + ux.e =3D x; + expsignx =3D ux.xbits.expsign; + exptx =3D expsignx & 0x7fff; + + if ((exptx=3D=3DBIAS+LDBL_MAX_EXP && + ((ux.bits.manh&~LDBL_NBIT)|ux.bits.manl)!=3D0) || /* x is NaN */ + (expty=3D=3DBIAS+LDBL_MAX_EXP && + ((uy.bits.manh&~LDBL_NBIT)|uy.bits.manl)!=3D0)) /* y is NaN */ + return nan_mix(x, y); + if (expsignx=3D=3DBIAS && ((ux.bits.manh&~LDBL_NBIT)|ux.bits.manl)=3D=3D0) + return atanl(y); /* x=3D1.0 */ + m =3D ((expsigny>>15)&1)|((expsignx>>14)&2); /* 2*sign(x)+sign(y) */ + + /* when y =3D 0 */ + if(expty=3D=3D0 && ((uy.bits.manh&~LDBL_NBIT)|uy.bits.manl)=3D=3D0) { + switch(m) { + case 0: + case 1: return y; /* atan(+-0,+anything)=3D+-0 */ + case 2: return pi+tiny;/* atan(+0,-anything) =3D pi */ + case 3: return -pi-tiny;/* atan(-0,-anything) =3D-pi */ + } + } + /* when x =3D 0 */ + if(exptx=3D=3D0 && ((ux.bits.manh&~LDBL_NBIT)|ux.bits.manl)=3D=3D0) + return (expsigny<0)? -pio2_hi-tiny: pio2_hi+tiny; + + /* when x is INF */ + if(exptx=3D=3DBIAS+LDBL_MAX_EXP) { + if(expty=3D=3DBIAS+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=3D=3DBIAS+LDBL_MAX_EXP) + return (expsigny<0)? -pio2_hi-tiny: pio2_hi+tiny; + + /* compute y/x */ + k =3D expty-exptx; + if(k > LDBL_MANT_DIG+2) { /* |y/x| huge */ + z=3Dpio2_hi+pio2_lo; + m&=3D1; + } + else if(expsignx<0&&k<-LDBL_MANT_DIG-2) z=3D0.0; /* |y/x| tiny, x<0 */ + else z=3Datanl(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 */ +/* + * =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D + * 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. + * =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D + * + */ + +#include +__FBSDID("$FreeBSD$"); + +/* + * See e_atanh.c for complete comments. + * + * Converted to long double by David Schultz and + * Bruce D. Evans. + */ + +#include +#ifdef __i386__ +#include +#endif + +#include "fpmath.h" +#include "math.h" +#include "math_private.h" + +/* EXP_TINY is the threshold below which we use atanh(x) ~=3D x. */ +#if LDBL_MANT_DIG =3D=3D 64 +#define EXP_TINY -34 +#elif LDBL_MANT_DIG =3D=3D 113 +#define EXP_TINY -58 +#else +#error "Unsupported long double format" +#endif + +#if LDBL_MAX_EXP !=3D 0x4000 +/* We also require the usual expsign encoding. */ +#error "Unsupported long double format" +#endif + +#define BIAS (LDBL_MAX_EXP - 1) + +static const double one =3D 1.0, huge =3D 1e300; +static const double zero =3D 0.0; + +long double +atanhl(long double x) +{ + long double t; + uint16_t hx, ix; + + ENTERI(); + GET_LDBL_EXPSIGN(hx, x); + ix =3D hx & 0x7fff; + if (ix >=3D 0x3fff) /* |x| >=3D 1, or NaN or misnormal */ + RETURNI(fabsl(x) =3D=3D 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 =3D x+x; + t =3D 0.5*log1pl(t+t*x/(one-x)); + } else + t =3D 0.5*log1pl((x+x)/(one-x)); + RETURNI((hx & 0x8000) =3D=3D 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 */ + +/* + * =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D + * 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. + * =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D + */ + +#include +__FBSDID("$FreeBSD$"); + +/* + * See e_cosh.c for complete comments. + * + * Converted to long double by Bruce D. Evans. + */ + +#include +#ifdef __i386__ +#include +#endif + +#include "fpmath.h" +#include "math.h" +#include "math_private.h" +#include "k_expl.h" + +#if LDBL_MAX_EXP !=3D 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 =3D 0x1p10000L, tiny =3D 0x1p-10000= L; +#if LDBL_MANT_DIG =3D=3D 64 +/* + * Domain [-1, 1], range ~[-1.8211e-21, 1.8211e-21]: + * |cosh(x) - c(x)| < 2**-68.8 + */ +static const union IEEEl2bits +C4u =3D LD80C(0xaaaaaaaaaaaaac78, -5, 4.16666666666666682297e-2L); +#define C4 C4u.e +static const double +C2 =3D 0.5, +C6 =3D 1.3888888888888616e-3, /* 0x16c16c16c16b99.0p-62 */ +C8 =3D 2.4801587301767953e-5, /* 0x1a01a01a027061.0p-68 */ +C10 =3D 2.7557319163300398e-7, /* 0x127e4fb6c9b55f.0p-74 */ +C12 =3D 2.0876768371393075e-9, /* 0x11eed99406a3f4.0p-81 */ +C14 =3D 1.1469537039374480e-11, /* 0x1938c67cd18c48.0p-89 */ +C16 =3D 4.8473490896852041e-14; /* 0x1b49c429701e45.0p-97 */ +#elif LDBL_MANT_DIG =3D=3D 113 +/* + * Domain [-1, 1], range ~[-2.3194e-37, 2.3194e-37]: + * |cosh(x) - c(x)| < 2**-121.69 + */ +static const long double +C4 =3D 4.16666666666666666666666666666666225e-2L, /* 0x1555555555555555= 555555555554e.0p-117L */ +C6 =3D 1.38888888888888888888888888889434831e-3L, /* 0x16c16c16c16c16c1= 6c16c16c1dd7a.0p-122L */ +C8 =3D 2.48015873015873015873015871870962089e-5L, /* 0x1a01a01a01a01a01= a01a017af2756.0p-128L */ +C10 =3D 2.75573192239858906525574318600800201e-7L, /* 0x127e4fb7789f5c72= ef01c8a040640.0p-134L */ +C12 =3D 2.08767569878680989791444691755468269e-9L, /* 0x11eed8eff8d897b5= 43d0679607399.0p-141L */ +C14=3D 1.14707455977297247387801189650495351e-11L, /* 0x193974a8c07c9d24= ae169a7fa9b54.0p-149L */ +C16 =3D 4.77947733238737883626416876486279985e-14L; /* 0x1ae7f3e733b814d= 4e1b90f5727fe4.0p-157L */ +static const double +C2 =3D 0.5, +C18 =3D 1.5619206968597871e-16, /* 0x16827863b9900b.0p-105 */ +C20 =3D 4.1103176218528049e-19, /* 0x1e542ba3d3c269.0p-114 */ +C22 =3D 8.8967926401641701e-22, /* 0x10ce399542a014.0p-122 */ +C24 =3D 1.6116681626523904e-24, /* 0x1f2c981d1f0cb7.0p-132 */ +C26 =3D 2.5022374732804632e-27; /* 0x18c7ecf8b2c4a0.0p-141 */ +#else +#error "Unsupported long double format" +#endif /* LDBL_MANT_DIG =3D=3D 64 */ + +/* log(2**16385 - 0.5) rounded up: */ +static const float +o_threshold =3D 1.13572168e4; /* 0xb174de.0p-10 */ + +long double +coshl(long double x) +{ + long double hi,lo,x2,x4; +#if LDBL_MANT_DIG =3D=3D 113 + double dx2; +#endif + uint16_t ix; + + GET_LDBL_EXPSIGN(ix,x); + ix &=3D 0x7fff; + + /* x is INF or NaN */ + if(ix>=3D0x7fff) return x*x; + + ENTERI(); + + /* |x| < 1, return 1 or c(x) */ + if(ix<0x3fff) { + if (ix 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 */ +/*- + * =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D + * 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. + * =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D + */ + +#include +__FBSDID("$FreeBSD$"); + +#include +#include + +#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 =3D 1.0, Zero[] =3D {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 =3D x; + uy.e =3D y; + sx =3D ux.bits.sign; + + /* purge off exception values */ + if((uy.bits.exp|uy.bits.manh|uy.bits.manl)=3D=3D0 || /* y=3D0 */ + (ux.bits.exp =3D=3D BIAS + LDBL_MAX_EXP) || /* or x not finite */ + (uy.bits.exp =3D=3D BIAS + LDBL_MAX_EXP && + ((uy.bits.manh&~LDBL_NBIT)|uy.bits.manl)!=3D0)) /* or y is NaN */ + return nan_mix_op(x, y, *)/nan_mix_op(x, y, *); + if(ux.bits.exp<=3Duy.bits.exp) { + if((ux.bits.exp>MANL_SHIFT); lx =3D lx+lx;} + else { + if ((hz|lz)=3D=3D0) /* return sign(x)*0 */ + return Zero[sx]; + hx =3D hz+hz+(lz>>MANL_SHIFT); lx =3D lz+lz; + } + } + hz=3Dhx-hy;lz=3Dlx-ly; if(lx=3D0) {hx=3Dhz;lx=3Dlz;} + + /* convert back to floating value and restore the sign */ + if((hx|lx)=3D=3D0) /* return sign(x)*0 */ + return Zero[sx]; + while(hx<(1ULL<>MANL_SHIFT); lx =3D lx+lx; + iy -=3D 1; + } + ux.bits.manh =3D hx; /* The mantissa is truncated here if needed. */ + ux.bits.manl =3D lx; + if (iy < LDBL_MIN_EXP) { + ux.bits.exp =3D iy + (BIAS + 512); + ux.e *=3D 0x1p-512; + } else { + ux.bits.exp =3D iy + BIAS; + } + x =3D 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 */ +/* + * =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D + * 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. + * =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D + */ + +#include +__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 + * 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 PURP= OSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENT= IAL + * 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, STR= ICT + * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY W= AY + * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF + * SUCH DAMAGE. + */ + +#include +__FBSDID("$FreeBSD$"); + +#include + +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 */ + +/* + * =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D + * 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. + * =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D + */ + +#include +__FBSDID("$FreeBSD$"); + +/* + * See e_sinh.c for complete comments. + * + * Converted to long double by Bruce D. Evans. + */ + +#include +#ifdef __i386__ +#include +#endif + +#include "fpmath.h" +#include "math.h" +#include "math_private.h" +#include "k_expl.h" + +#if LDBL_MAX_EXP !=3D 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 =3D 0x1p16383L; +#if LDBL_MANT_DIG =3D=3D 64 +/* + * Domain [-1, 1], range ~[-6.6749e-22, 6.6749e-22]: + * |sinh(x)/x - s(x)| < 2**-70.3 + */ +static const union IEEEl2bits +S3u =3D LD80C(0xaaaaaaaaaaaaaaaa, -3, 1.66666666666666666658e-1L); +#define S3 S3u.e +static const double +S5 =3D 8.3333333333333332e-3, /* 0x11111111111111.0p-59 */ +S7 =3D 1.9841269841270074e-4, /* 0x1a01a01a01a070.0p-65 */ +S9 =3D 2.7557319223873889e-6, /* 0x171de3a5565fe6.0p-71 */ +S11 =3D 2.5052108406704084e-8, /* 0x1ae6456857530f.0p-78 */ +S13 =3D 1.6059042748655297e-10, /* 0x161245fa910697.0p-85 */ +S15 =3D 7.6470006914396920e-13, /* 0x1ae7ce4eff2792.0p-93 */ +S17 =3D 2.8346142308424267e-15; /* 0x19882ce789ffc6.0p-101 */ +#elif LDBL_MANT_DIG =3D=3D 113 +/* + * Domain [-1, 1], range ~[-2.9673e-36, 2.9673e-36]: + * |sinh(x)/x - s(x)| < 2**-118.0 + */ +static const long double +S3 =3D 1.66666666666666666666666666666666033e-1L, /* 0x1555555555555555= 555555555553b.0p-115L */ +S5 =3D 8.33333333333333333333333333337643193e-3L, /* 0x1111111111111111= 11111111180f5.0p-119L */ +S7 =3D 1.98412698412698412698412697391263199e-4L, /* 0x1a01a01a01a01a01= a01a0176aad11.0p-125L */ +S9 =3D 2.75573192239858906525574406205464218e-6L, /* 0x171de3a556c7338f= aac243aaa9592.0p-131L */ +S11 =3D 2.50521083854417187749675637460977997e-8L, /* 0x1ae64567f544e38f= e59b3380d7413.0p-138L */ +S13 =3D 1.60590438368216146368737762431552702e-10L, /* 0x16124613a86d098= 059c7620850fc2.0p-145L */ +S15 =3D 7.64716373181980539786802470969096440e-13L, /* 0x1ae7f3e733b8141= 93af09ce723043.0p-153L */ +S17 =3D 2.81145725434775409870584280722701574e-15L; /* 0x1952c77030c3689= 8c3fd0b6dfc562.0p-161L */ +static const double +S19=3D 8.2206352435411005e-18, /* 0x12f49b4662b86d.0p-109 */ +S21=3D 1.9572943931418891e-20, /* 0x171b8f2fab9628.0p-118 */ +S23 =3D 3.8679983530666939e-23, /* 0x17617002b73afc.0p-127 */ +S25 =3D 6.5067867911512749e-26; /* 0x1423352626048a.0p-136 */ +#else +#error "Unsupported long double format" +#endif /* LDBL_MANT_DIG =3D=3D 64 */ + +/* log(2**16385 - 0.5) rounded up: */ +static const float +o_threshold =3D 1.13572168e4; /* 0xb174de.0p-10 */ + +long double +sinhl(long double x) +{ + long double hi,lo,x2,x4; +#if LDBL_MANT_DIG =3D=3D 113 + double dx2; +#endif + double s; + int16_t ix,jx; + + GET_LDBL_EXPSIGN(jx,x); + ix =3D jx&0x7fff; + + /* x is INF or NaN */ + if(ix>=3D0x7fff) return x+x; + + ENTERI(); + + s =3D 1; + if (jx<0) s =3D -1; + + /* |x| < 64, return x, s(x), or accurate s*(exp(|x|)/2-1/exp(|x|)/2) */ + if (ix<0x4005) { /* |x|<64 */ + if (ix1) RETURNI(x); /* sinh(tiny) =3D tiny with inexact */ + if (ix<0x3fff) { /* |x|<1 */ + x2 =3D x*x; +#if LDBL_MANT_DIG =3D=3D 64 + x4 =3D 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 =3D=3D 113 + dx2 =3D 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) <=3D 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 + * Copyright (c) 2002 David Schultz + * 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 PURP= OSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENT= IAL + * 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, STR= ICT + * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY W= AY + * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF + * SUCH DAMAGE. + * + * $FreeBSD$ + */ + +#ifndef _FPMATH_H_ +#define _FPMATH_H_ + +#include +#include "_fpmath.h" + +#ifndef _IEEE_WORD_ORDER +#define _IEEE_WORD_ORDER _BYTE_ORDER +#endif + +union IEEEf2bits { + float f; + struct { +#if _BYTE_ORDER =3D=3D _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 =3D=3D _LITTLE_ENDIAN +#if _IEEE_WORD_ORDER =3D=3D _LITTLE_ENDIAN + unsigned int manl :32; +#endif + unsigned int manh :20; + unsigned int exp :11; + unsigned int sign :1; +#if _IEEE_WORD_ORDER =3D=3D _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 @@ +/* + * =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D + * 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. + * =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D + */ + +/* + * from: @(#)fdlibm.h 5.1 93/09/24 + * $FreeBSD$ + */ + +#ifndef _MATH_PRIVATE_H_ +#define _MATH_PRIVATE_H_ + +#include +#include + +/* + * The original fdlibm code used statements like: + * n0 =3D ((*(int*)&one)>>29)^1; * index of high word * + * ix0 =3D *(n0+(int*)&x); * high word of x * + * ix1 =3D *((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 =3D=3D 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 =3D=3D 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 =3D=3D 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 =3D=3D 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 =3D (d); \ + (ix0) =3D ew_u.parts.msw; \ + (ix1) =3D 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 =3D (d); \ + (ix) =3D 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 =3D (d); \ + (i) =3D 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 =3D (d); \ + (i) =3D 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 =3D (ix0); \ + iw_u.parts.lsw =3D (ix1); \ + (d) =3D 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 =3D (ix); \ + (d) =3D 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 =3D (d); \ + sh_u.parts.msw =3D (v); \ + (d) =3D 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 =3D (d); \ + sl_u.parts.lsw =3D (v); \ + (d) =3D 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 =3D (d); \ + (i) =3D 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 =3D (i); \ + (d) =3D 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 =3D (d); \ + (ix0) =3D ew_u.xbits.expsign; \ + (ix1) =3D ew_u.xbits.man; \ +} while (0) + +/* + * Get expsign and mantissa as one 16 bit and two 64 bit ints from a 128 b= it + * long double. + */ + +#define EXTRACT_LDBL128_WORDS(ix0,ix1,ix2,d) \ +do { \ + union IEEEl2bits ew_u; \ + ew_u.e =3D (d); \ + (ix0) =3D ew_u.xbits.expsign; \ + (ix1) =3D ew_u.xbits.manh; \ + (ix2) =3D 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 =3D (d); \ + (i) =3D 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 =3D (ix0); \ + iw_u.xbits.man =3D (ix1); \ + (d) =3D 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 =3D (ix0); \ + iw_u.xbits.manh =3D (ix1); \ + iw_u.xbits.manl =3D (ix2); \ + (d) =3D 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 =3D (d); \ + se_u.xbits.expsign =3D (v); \ + (d) =3D se_u.e; \ +} while (0) + +#ifdef __i386__ +/* Long double constants are broken on i386. */ +#define LD80C(m, ex, v) { \ + .xbits.man =3D __CONCAT(m, ULL), \ + .xbits.expsign =3D (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 =3D (v), } +#endif + +#ifdef FLT_EVAL_METHOD +/* + * Attempt to get strict C99 semantics for assignment with non-C99 compile= rs. + */ +#if FLT_EVAL_METHOD =3D=3D 0 || __GNUC__ =3D=3D 0 +#define STRICT_ASSIGN(type, lval, rval) ((lval) =3D (rval)) +#else +#define STRICT_ASSIGN(type, lval, rval) do { \ + volatile type __lval; \ + \ + if (sizeof(type) >=3D sizeof(long double)) \ + (lval) =3D (rval); \ + else { \ + __lval =3D (rval); \ + (lval) =3D __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 =3D fpgetprec()) !=3D FP_PE) \ + fpsetprec(FP_PE) +#define RETURNI(x) do { \ + __retval =3D (x); \ + if (__oprec !=3D FP_PE) \ + fpsetprec(__oprec); \ + RETURNF(__retval); \ +} while (0) +#define ENTERV() \ + fp_prec_t __oprec; \ + \ + if ((__oprec =3D fpgetprec()) !=3D FP_PE) \ + fpsetprec(FP_PE) +#define RETURNV() do { \ + if (__oprec !=3D 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| >=3D |b| or + * a =3D=3D 0, but is slower. + */ +#define _2sum(a, b) do { \ + __typeof(a) __s, __w; \ + \ + __w =3D (a) + (b); \ + __s =3D __w - (a); \ + (b) =3D ((a) - (__w - __s)) + ((b) - __s); \ + (a) =3D __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 >=3D |b| or 0, b's type must be no larger than 'a's type, a= nd + * 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 =3D (a); \ + __ib =3D (b); \ + assert(__ia =3D=3D 0 || fabsl(__ia) >=3D fabsl(__ib)); \ + \ + __w =3D (a) + (b); \ + (b) =3D ((a) - __w) + (b); \ + (a) =3D __w; \ + \ + /* The next 2 assertions are weak if (a) is already long double. */ \ + assert((long double)__ia + __ib =3D=3D (long double)(a) + (b)); \ + __vw =3D __ia + __ib; \ + __r =3D __ia - __vw; \ + __r +=3D __ib; \ + assert(__vw =3D=3D (a) && __r =3D=3D (b)); \ +} while (0) +#else /* !DEBUG */ +#define _2sumF(a, b) do { \ + __typeof(a) __w; \ + \ + __w =3D (a) + (b); \ + (b) =3D ((a) - __w) + (b); \ + (a) =3D __w; \ +} while (0) +#endif /* DEBUG */ + +/* + * Set x +=3D 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 >=3D |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. Howeve= r, + * 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 =3D (c); \ + _2sumF(__tmp, (a)); \ + (b) +=3D (a); \ + (a) =3D __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 bo= th + * 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) =3D x; + IMAGPART(z) =3D y; + return (z.f); +} +#endif + +#ifndef CMPLX +static __inline double complex +CMPLX(double x, double y) +{ + double_complex z; + + REALPART(z) =3D x; + IMAGPART(z) =3D 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) =3D x; + IMAGPART(z) =3D 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 problem= s. + */ + +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) =3D=3D sizeof(float) && \ + sizeof(__float_t) =3D=3D sizeof(long double) ? irintf(x) : \ + sizeof(x) =3D=3D sizeof(double) && \ + sizeof(__double_t) =3D=3D sizeof(long double) ? irintd(x) : \ + sizeof(x) =3D=3D 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" : "=3Dm" (n) : "t" (x)); + return (n); +} + +static __inline int +irintd(double x) +{ + int n; + + __asm("fistl %0" : "=3Dm" (n) : "t" (x)); + return (n); +} +#endif + +#if defined(__amd64__) || defined(__i386__) +static __inline int +irintl(long double x) +{ + int n; + + __asm("fistl %0" : "=3Dm" (n) : "t" (x)); + return (n); +} +#endif + +#ifdef DEBUG +#if defined(__amd64__) || defined(__i386__) +#define breakpoint() asm("int $3") +#else +#include + +#define breakpoint() raise(SIGTRAP) +#endif +#endif + +/* Write a pari script to test things externally. */ +#ifdef DOPRINT +#include + +#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 ^=3D DOPRINT_SWIZZLE; \ + INSERT_LDBL80_WORDS(*xp, __hx, __lx); \ + printf("x =3D %.21Lg; ", (long double)*xp); \ +} while (0) +#define DOPRINT_END1(v) \ + printf("y =3D %.21Lg; z =3D 0; show(x, y, z);\n", (long double)(v)) +#define DOPRINT_END2(hi, lo) \ + printf("y =3D %.21Lg; z =3D %.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 ^=3D DOPRINT_SWIZZLE; \ + INSERT_WORDS(*xp, __hx, __lx); \ + printf("x =3D %.21Lg; ", (long double)*xp); \ +} while (0) +#define DOPRINT_END1(v) \ + printf("y =3D %.21Lg; z =3D 0; show(x, y, z);\n", (long double)(v)) +#define DOPRINT_END2(hi, lo) \ + printf("y =3D %.21Lg; z =3D %.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 ^=3D DOPRINT_SWIZZLE; \ + SET_FLOAT_WORD(*xp, __hx); \ + printf("x =3D %.21Lg; ", (long double)*xp); \ +} while (0) +#define DOPRINT_END1(v) \ + printf("y =3D %.21Lg; z =3D 0; show(x, y, z);\n", (long double)(v)) +#define DOPRINT_END2(hi, lo) \ + printf("y =3D %.21Lg; z =3D %.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 ^=3D DOPRINT_SWIZZLE; \ + __lx ^=3D DOPRINT_SWIZZLE_HIGH; \ + INSERT_LDBL128_WORDS(*xp, __hx, __lx, __llx); \ + printf("x =3D %.36Lg; ", (long double)*xp); \ +} while (0) +#define DOPRINT_END1(v) \ + printf("y =3D %.36Lg; z =3D 0; show(x, y, z);\n", (long double)(v)) +#define DOPRINT_END2(hi, lo) \ + printf("y =3D %.36Lg; z =3D %.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 =3D (x); \ + const __typeof (y) __y =3D (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 */ +/* + * =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D + * 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. + * =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D + */ + +#include +__FBSDID("$FreeBSD$"); + +/* + * See s_asinh.c for complete comments. + * + * Converted to long double by David Schultz and + * Bruce D. Evans. + */ + +#include +#ifdef __i386__ +#include +#endif + +#include "fpmath.h" +#include "math.h" +#include "math_private.h" + +/* EXP_LARGE is the threshold above which we use asinh(x) ~=3D log(2x). */ +/* EXP_TINY is the threshold below which we use asinh(x) ~=3D x. */ +#if LDBL_MANT_DIG =3D=3D 64 +#define EXP_LARGE 34 +#define EXP_TINY -34 +#elif LDBL_MANT_DIG =3D=3D 113 +#define EXP_LARGE 58 +#define EXP_TINY -58 +#else +#error "Unsupported long double format" +#endif + +#if LDBL_MAX_EXP !=3D 0x4000 +/* We also require the usual expsign encoding. */ +#error "Unsupported long double format" +#endif + +#define BIAS (LDBL_MAX_EXP - 1) + +static const double +one =3D 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ +huge=3D 1.00000000000000000000e+300; + +#if LDBL_MANT_DIG =3D=3D 64 +static const union IEEEl2bits +u_ln2 =3D LD80C(0xb17217f7d1cf79ac, -1, 6.93147180559945309417e-1L); +#define ln2 u_ln2.e +#elif LDBL_MANT_DIG =3D=3D 113 +static const long double +ln2 =3D 6.93147180559945309417232121458176568e-1L; /* 0x162e42fefa39ef357= 93c7673007e6.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 =3D hx & 0x7fff; + if (ix >=3D 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 >=3D BIAS + EXP_LARGE) { /* |x| >=3D LARGE, or misnormal */ + w =3D logl(fabsl(x))+ln2; + } else if (ix >=3D 0x4000) { /* LARGE > |x| >=3D 2.0, or misnormal */ + t =3D fabsl(x); + w =3D logl(2.0*t+one/(sqrtl(x*x+one)+t)); + } else { /* 2.0 > |x| >=3D TINY, or misnormal */ + t =3D x*x; + w =3Dlog1pl(fabsl(x)+t/(one+sqrtl(one+t))); + } + RETURNI((hx & 0x8000) =3D=3D 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 */ +/* + * =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D + * 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. + * =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D + */ + +#include +__FBSDID("$FreeBSD$"); + +/* + * See comments in s_atan.c. + * Converted to long double by David Schultz . + */ + +#include + +#include "invtrig.h" +#include "math.h" +#include "math_private.h" + +static const long double +one =3D 1.0, +huge =3D 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 =3D x; + expsign =3D u.xbits.expsign; + expt =3D expsign & 0x7fff; + if(expt >=3D ATAN_CONST) { /* if |x| is large, atan(x)~=3Dpi/2 */ + if(expt =3D=3D BIAS + LDBL_MAX_EXP && + ((u.bits.manh&~LDBL_NBIT)|u.bits.manl)!=3D0) + 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 =3D (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)~=3Dx */ + if(huge+x>one) return x; /* raise inexact */ + } + id =3D -1; + } else { + x =3D fabsl(x); + if (expman < (BIAS << 8) + 0x30) { /* |x| < 1.1875 */ + if (expman < ((BIAS - 1) << 8) + 0x60) { /* 7/16 <=3D|x|<11/16 */ + id =3D 0; x =3D (2.0*x-one)/(2.0+x); + } else { /* 11/16<=3D|x|< 19/16 */ + id =3D 1; x =3D (x-one)/(x+one); + } + } else { + if (expman < ((BIAS + 1) << 8) + 0x38) { /* |x| < 2.4375 */ + id =3D 2; x =3D (x-1.5)/(one+1.5*x); + } else { /* 2.4375 <=3D |x| < 2^ATAN_CONST */ + id =3D 3; x =3D -1.0/x; + } + }} + /* end of argument reduction */ + z =3D x*x; + w =3D z*z; + /* break sum aT[i]z**(i+1) into odd and even poly */ + s1 =3D z*T_even(w); + s2 =3D w*T_odd(w); + if (id<0) return x - x*(s1+s2); + else { + z =3D 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 @@ +/*- + * =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D + * 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. + * =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D + * + * 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 +__FBSDID("$FreeBSD$"); + +#include +#ifdef __i386__ +#include +#endif + +#include "fpmath.h" +#include "math.h" +#include "math_private.h" + +#define BIAS (LDBL_MAX_EXP - 1) + +static const unsigned + B1 =3D 709958130; /* B1 =3D (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 =3D x; + expsign =3D u.xbits.expsign; + k =3D expsign & 0x7fff; + + /* + * If x =3D +-Inf, then cbrt(x) =3D +-Inf. + * If x =3D NaN, then cbrt(x) =3D NaN. + */ + if (k =3D=3D BIAS + LDBL_MAX_EXP) + return (x + x); + + ENTERI(); + if (k =3D=3D 0) { + /* If x =3D +-0, then cbrt(x) =3D +-0. */ + if ((u.bits.manh | u.bits.manl) =3D=3D 0) + RETURNI(x); + /* Adjust subnormal numbers. */ + u.e *=3D 0x1.0p514; + k =3D u.bits.exp; + k -=3D BIAS + 514; + } else + k -=3D BIAS; + u.xbits.expsign =3D BIAS; + v.e =3D 1; + + x =3D u.e; + switch (k % 3) { + case 1: + case -2: + x =3D 2*x; + k--; + break; + case 2: + case -1: + x =3D 4*x; + k -=3D 2; + break; + } + v.xbits.expsign =3D (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 =3D x; + GET_FLOAT_WORD(hx, fx); + SET_FLOAT_WORD(ft, ((hx & 0x7fffffff) / 3 + B1)); + + /* ~16-bit estimate: */ + dx =3D x; + dt =3D ft; + dr =3D dt * dt * dt; + dt =3D dt * (dx + dx + dr) / (dx + dr + dr); + + /* ~47-bit estimate: */ + dr =3D dt * dt * dt; + dt =3D dt * (dx + dx + dr) / (dx + dr + dr); + +#if LDBL_MANT_DIG =3D=3D 64 + /* + * dt is cbrtl(x) to ~47 bits (after x has been reduced to 1 <=3D 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 =3D 0x1.0p32; + volatile double vd1 =3D 0x1.0p-31; + #define vd ((long double)vd2 + vd1) + + t =3D dt + vd - 0x1.0p32; +#elif LDBL_MANT_DIG =3D=3D 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 =3D 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=3Dt*t; /* t*t is exact */ + r=3Dx/s; /* error <=3D 0.5 ulps; |r| < |t| */ + w=3Dt+t; /* t+t is exact */ + r=3D(r-t)/(w+r); /* r-t is exact; w+r ~=3D 3*t */ + t=3Dt+t*r; /* error <=3D (0.5 + 0.5/3) * ulp */ + + t *=3D 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 @@ +/* + * =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D + * 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. + * =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D + * + * From: @(#)s_ceil.c 5.1 93/09/24 + */ + +#include +__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 +#include +#include + +#include "fpmath.h" + +#ifdef LDBL_IMPLICIT_NBIT +#define MANH_SIZE (LDBL_MANH_SIZE + 1) +#define INC_MANH(u, c) do { \ + uint64_t o =3D u.bits.manh; \ + u.bits.manh +=3D (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 =3D u.bits.manh; \ + u.bits.manh +=3D (c); \ + if (u.bits.manh < o) { \ + u.bits.exp++; \ + u.bits.manh |=3D 1llu << (LDBL_MANH_SIZE - 1); \ + } \ +} while (0) +#endif + +static const long double huge =3D 1.0e300; + +long double +ceill(long double x) +{ + union IEEEl2bits u =3D { .e =3D x }; + int e =3D u.bits.exp - LDBL_MAX_EXP + 1; + + if (e < MANH_SIZE - 1) { + if (e < 0) { /* raise inexact if x !=3D 0 */ + if (huge + x > 0.0) + if (u.bits.exp > 0 || + (u.bits.manh | u.bits.manl) !=3D 0) + u.e =3D u.bits.sign ? -0.0 : 1.0; + } else { + uint64_t m =3D ((1llu << MANH_SIZE) - 1) >> (e + 1); + if (((u.bits.manh & m) | u.bits.manl) =3D=3D 0) + return (x); /* x is integral */ + if (!u.bits.sign) { +#ifdef LDBL_IMPLICIT_NBIT + if (e =3D=3D 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 &=3D ~m; + u.bits.manl =3D 0; + } + } + } else if (e < LDBL_MANT_DIG - 1) { + uint64_t m =3D (uint64_t)-1 >> (64 - LDBL_MANT_DIG + e + 1); + if ((u.bits.manl & m) =3D=3D 0) + return (x); /* x is integral */ + if (!u.bits.sign) { + if (e =3D=3D MANH_SIZE - 1) + INC_MANH(u, 1); + else { + uint64_t o =3D u.bits.manl; + u.bits.manl +=3D 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 &=3D ~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 PURP= OSE + * ARE DISCLAIMED. IN NO EVENT SHALL AUTHOR OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENT= IAL + * 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, STR= ICT + * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY W= AY + * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF + * SUCH DAMAGE. + * + * $FreeBSD$ + */ + +#include + +#include "fpmath.h" + +long double +copysignl(long double x, long double y) +{ + union IEEEl2bits ux, uy; + + ux.e =3D x; + uy.e =3D y; + ux.bits.sign =3D 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 WARRANTI= ES + * 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 US= E, + * 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 +__FBSDID("$FreeBSD$"); + +/* + * Limited testing on pseudorandom numbers drawn within [-2e8:4e8] shows + * an accuracy of <=3D 0.7412 ULP. + */ + +#include +#ifdef __i386__ +#include +#endif + +#include "fpmath.h" +#include "math.h" +#include "math_private.h" +#if LDBL_MANT_DIG =3D=3D 64 +#include "../ld80/e_rem_pio2l.h" +static const union IEEEl2bits +pio4u =3D LD80C(0xc90fdaa22168c235, -00001, 7.85398163397448309628e-01L); +#define pio4 (pio4u.e) +#elif LDBL_MANT_DIG =3D=3D 113 +#include "../ld128/e_rem_pio2l.h" +long double pio4 =3D 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 =3D x; + z.bits.sign =3D 0; + + /* If x =3D +-0 or x is a subnormal number, then cos(x) =3D 1 */ + if (z.bits.exp =3D=3D 0) + return (1.0); + + /* If x =3D NaN or Inf, then cos(x) =3D NaN. */ + if (z.bits.exp =3D=3D 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 =3D __ieee754_rem_pio2l(x, y); + hi =3D y[0]; + lo =3D y[1]; + + switch (e0 & 3) { + case 0: + hi =3D __kernel_cosl(hi, lo); + break; + case 1: + hi =3D - __kernel_sinl(hi, lo, 1); + break; + case 2: + hi =3D - __kernel_cosl(hi, lo); + break; + case 3: + hi =3D __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=C3=B8rgrav + * 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 WARRANTI= ES + * 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 US= E, + * 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 + +#include "fpmath.h" + +long double +fabsl(long double x) +{ + union IEEEl2bits u; + + u.e =3D x; + u.bits.sign =3D 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 + * 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 PURP= OSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENT= IAL + * 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, STR= ICT + * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY W= AY + * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF + * SUCH DAMAGE. + */ + +#include +__FBSDID("$FreeBSD$"); + +#include + +#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 @@ +/* + * =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D + * 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. + * =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D + * + * From: @(#)s_floor.c 5.1 93/09/24 + */ + +#include +__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 +#include +#include + +#include "fpmath.h" + +#ifdef LDBL_IMPLICIT_NBIT +#define MANH_SIZE (LDBL_MANH_SIZE + 1) +#define INC_MANH(u, c) do { \ + uint64_t o =3D u.bits.manh; \ + u.bits.manh +=3D (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 =3D u.bits.manh; \ + u.bits.manh +=3D (c); \ + if (u.bits.manh < o) { \ + u.bits.exp++; \ + u.bits.manh |=3D 1llu << (LDBL_MANH_SIZE - 1); \ + } \ +} while (0) +#endif + +static const long double huge =3D 1.0e300; + +long double +floorl(long double x) +{ + union IEEEl2bits u =3D { .e =3D x }; + int e =3D u.bits.exp - LDBL_MAX_EXP + 1; + + if (e < MANH_SIZE - 1) { + if (e < 0) { /* raise inexact if x !=3D 0 */ + if (huge + x > 0.0) + if (u.bits.exp > 0 || + (u.bits.manh | u.bits.manl) !=3D 0) + u.e =3D u.bits.sign ? -1.0 : 0.0; + } else { + uint64_t m =3D ((1llu << MANH_SIZE) - 1) >> (e + 1); + if (((u.bits.manh & m) | u.bits.manl) =3D=3D 0) + return (x); /* x is integral */ + if (u.bits.sign) { +#ifdef LDBL_IMPLICIT_NBIT + if (e =3D=3D 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 &=3D ~m; + u.bits.manl =3D 0; + } + } + } else if (e < LDBL_MANT_DIG - 1) { + uint64_t m =3D (uint64_t)-1 >> (64 - LDBL_MANT_DIG + e + 1); + if ((u.bits.manl & m) =3D=3D 0) + return (x); /* x is integral */ + if (u.bits.sign) { + if (e =3D=3D MANH_SIZE - 1) + INC_MANH(u, 1); + else { + uint64_t o =3D u.bits.manl; + u.bits.manl +=3D 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 &=3D ~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 + * 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 PURP= OSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENT= IAL + * 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, STR= I[...] [diff truncated at 100000 bytes]