From fea2240d10b31df0708048e117338954d65683f2 Mon Sep 17 00:00:00 2001 From: das Date: Thu, 31 Jul 2008 22:41:26 +0000 Subject: Add implementations of acosl(), asinl(), atanl(), atan2l(), and cargl(). Reviewed by: bde sparc64 testing resources from: remko --- lib/msun/Makefile | 20 +++++---- lib/msun/Symbol.map | 5 +++ lib/msun/ld128/invtrig.c | 100 +++++++++++++++++++++++++++++++++++++++++ lib/msun/ld128/invtrig.h | 113 +++++++++++++++++++++++++++++++++++++++++++++++ lib/msun/ld80/invtrig.c | 82 ++++++++++++++++++++++++++++++++++ lib/msun/ld80/invtrig.h | 101 ++++++++++++++++++++++++++++++++++++++++++ lib/msun/man/acos.3 | 33 +++++++------- lib/msun/man/asin.3 | 33 +++++++------- lib/msun/man/atan.3 | 31 ++++++------- lib/msun/man/atan2.3 | 38 +++++++++++----- lib/msun/src/e_acos.c | 6 +++ lib/msun/src/e_acosl.c | 77 ++++++++++++++++++++++++++++++++ lib/msun/src/e_asin.c | 5 +++ lib/msun/src/e_asinl.c | 77 ++++++++++++++++++++++++++++++++ lib/msun/src/e_atan2.c | 6 +++ lib/msun/src/e_atan2l.c | 107 ++++++++++++++++++++++++++++++++++++++++++++ lib/msun/src/math.h | 6 +++ lib/msun/src/s_atan.c | 6 +++ lib/msun/src/s_atanl.c | 85 +++++++++++++++++++++++++++++++++++ lib/msun/src/s_cargl.c | 38 ++++++++++++++++ 20 files changed, 900 insertions(+), 69 deletions(-) create mode 100644 lib/msun/ld128/invtrig.c create mode 100644 lib/msun/ld128/invtrig.h create mode 100644 lib/msun/ld80/invtrig.c create mode 100644 lib/msun/ld80/invtrig.h create mode 100644 lib/msun/src/e_acosl.c create mode 100644 lib/msun/src/e_asinl.c create mode 100644 lib/msun/src/e_atan2l.c create mode 100644 lib/msun/src/s_atanl.c create mode 100644 lib/msun/src/s_cargl.c (limited to 'lib') diff --git a/lib/msun/Makefile b/lib/msun/Makefile index 404dd1c..11aabef 100644 --- a/lib/msun/Makefile +++ b/lib/msun/Makefile @@ -25,8 +25,10 @@ ARCH_SUBDIR= ${MACHINE_ARCH} # long double format .if ${LDBL_PREC} == 64 .PATH: ${.CURDIR}/ld80 +CFLAGS+= -I${.CURDIR}/ld80 .elif ${LDBL_PREC} == 113 .PATH: ${.CURDIR}/ld128 +CFLAGS+= -I${.CURDIR}/ld128 .endif .PATH: ${.CURDIR}/bsdsrc @@ -48,7 +50,7 @@ COMMON_SRCS= b_exp.c b_log.c b_tgamma.c \ e_sinh.c e_sinhf.c e_sqrt.c e_sqrtf.c fenv.c \ k_cos.c k_cosf.c k_rem_pio2.c k_sin.c k_sinf.c \ k_tan.c k_tanf.c \ - s_asinh.c s_asinhf.c s_atan.c s_atanf.c s_carg.c s_cargf.c \ + s_asinh.c s_asinhf.c s_atan.c s_atanf.c s_carg.c s_cargf.c s_cargl.c \ s_cbrt.c s_cbrtf.c s_ceil.c s_ceilf.c \ s_copysign.c s_copysignf.c s_cos.c s_cosf.c \ s_csqrt.c s_csqrtf.c s_erf.c s_erff.c \ @@ -82,9 +84,10 @@ SYMBOL_MAPS= ${SYM_MAPS} COMMON_SRCS+= s_copysignl.c s_fabsl.c s_llrintl.c s_lrintl.c s_modfl.c .if ${LDBL_PREC} != 53 # If long double != double use these; otherwise, we alias the double versions. -COMMON_SRCS+= e_fmodl.c e_hypotl.c e_remainderl.c e_sqrtl.c \ - k_cosl.c k_sinl.c k_tanl.c \ - s_ceill.c s_cosl.c s_csqrtl.c s_exp2l.c s_floorl.c s_fmal.c \ +COMMON_SRCS+= e_acosl.c e_asinl.c e_atan2l.c e_fmodl.c \ + e_hypotl.c e_remainderl.c e_sqrtl.c \ + invtrig.c k_cosl.c k_sinl.c k_tanl.c \ + s_atanl.c s_ceill.c s_cosl.c s_csqrtl.c s_exp2l.c s_floorl.c s_fmal.c \ s_frexpl.c s_logbl.c s_nanl.c s_nextafterl.c s_nexttoward.c \ s_remquol.c s_rintl.c s_scalbnl.c \ s_sinl.c s_tanl.c s_truncl.c w_cabsl.c @@ -117,13 +120,14 @@ MAN= acos.3 acosh.3 asin.3 asinh.3 atan.3 atan2.3 atanh.3 ceil.3 \ nextafter.3 remainder.3 rint.3 \ round.3 scalbn.3 signbit.3 sin.3 sinh.3 sqrt.3 tan.3 tanh.3 trunc.3 -MLINKS+=acos.3 acosf.3 +MLINKS+=acos.3 acosf.3 acos.3 acosl.3 MLINKS+=acosh.3 acoshf.3 -MLINKS+=asin.3 asinf.3 +MLINKS+=asin.3 asinf.3 asin.3 asinl.3 MLINKS+=asinh.3 asinhf.3 -MLINKS+=atan.3 atanf.3 +MLINKS+=atan.3 atanf.3 atan.3 atanl.3 MLINKS+=atanh.3 atanhf.3 -MLINKS+=atan2.3 atan2f.3 atan2.3 carg.3 atan2.3 cargf.3 +MLINKS+=atan2.3 atan2f.3 atan2.3 atan2l.3 \ + atan2.3 carg.3 atan2.3 cargf.3 atan2.3 cargl.3 MLINKS+=ceil.3 ceilf.3 ceil.3 ceill.3 MLINKS+=cimag.3 cimagf.3 cimag.3 cimagl.3 \ cimag.3 conj.3 cimag.3 conjf.3 cimag.3 conjl.3 \ diff --git a/lib/msun/Symbol.map b/lib/msun/Symbol.map index 526d7e3..e1e8c4a 100644 --- a/lib/msun/Symbol.map +++ b/lib/msun/Symbol.map @@ -209,4 +209,9 @@ FBSD_1.1 { remquol; remainderl; fmodl; + acosl; + asinl; + atan2l; + atanl; + cargl; }; diff --git a/lib/msun/ld128/invtrig.c b/lib/msun/ld128/invtrig.c new file mode 100644 index 0000000..4ceca8a --- /dev/null +++ b/lib/msun/ld128/invtrig.c @@ -0,0 +1,100 @@ +/*- + * 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 PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL + * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS + * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) + * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT + * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY + * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF + * SUCH DAMAGE. + */ + +#include +__FBSDID("$FreeBSD$"); + +#include "invtrig.h" + +/* + * asinl() and acosl() + */ +const long double +pS0 = 1.66666666666666666666666666666700314e-01L, +pS1 = -7.32816946414566252574527475428622708e-01L, +pS2 = 1.34215708714992334609030036562143589e+00L, +pS3 = -1.32483151677116409805070261790752040e+00L, +pS4 = 7.61206183613632558824485341162121989e-01L, +pS5 = -2.56165783329023486777386833928147375e-01L, +pS6 = 4.80718586374448793411019434585413855e-02L, +pS7 = -4.42523267167024279410230886239774718e-03L, +pS8 = 1.44551535183911458253205638280410064e-04L, +pS9 = -2.10558957916600254061591040482706179e-07L, +qS1 = -4.84690167848739751544716485245697428e+00L, +qS2 = 9.96619113536172610135016921140206980e+00L, +qS3 = -1.13177895428973036660836798461641458e+01L, +qS4 = 7.74004374389488266169304117714658761e+00L, +qS5 = -3.25871986053534084709023539900339905e+00L, +qS6 = 8.27830318881232209752469022352928864e-01L, +qS7 = -1.18768052702942805423330715206348004e-01L, +qS8 = 8.32600764660522313269101537926539470e-03L, +qS9 = -1.99407384882605586705979504567947007e-04L; + +/* + * atanl() + */ +const long double atanhi[] = { + 4.63647609000806116214256231461214397e-01L, + 7.85398163397448309615660845819875699e-01L, + 9.82793723247329067985710611014666038e-01L, + 1.57079632679489661923132169163975140e+00L, +}; + +const long double atanlo[] = { + 4.89509642257333492668618435220297706e-36L, + 2.16795253253094525619926100651083806e-35L, + -2.31288434538183565909319952098066272e-35L, + 4.33590506506189051239852201302167613e-35L, +}; + +const long double aT[] = { + 3.33333333333333333333333333333333125e-01L, + -1.99999999999999999999999999999180430e-01L, + 1.42857142857142857142857142125269827e-01L, + -1.11111111111111111111110834490810169e-01L, + 9.09090909090909090908522355708623681e-02L, + -7.69230769230769230696553844935357021e-02L, + 6.66666666666666660390096773046256096e-02L, + -5.88235294117646671706582985209643694e-02L, + 5.26315789473666478515847092020327506e-02L, + -4.76190476189855517021024424991436144e-02L, + 4.34782608678695085948531993458097026e-02L, + -3.99999999632663469330634215991142368e-02L, + 3.70370363987423702891250829918659723e-02L, + -3.44827496515048090726669907612335954e-02L, + 3.22579620681420149871973710852268528e-02L, + -3.03020767654269261041647570626778067e-02L, + 2.85641979882534783223403715930946138e-02L, + -2.69824879726738568189929461383741323e-02L, + 2.54194698498808542954187110873675769e-02L, + -2.35083879708189059926183138130183215e-02L, + 2.04832358998165364349957325067131428e-02L, + -1.54489555488544397858507248612362957e-02L, + 8.64492360989278761493037861575248038e-03L, + -2.58521121597609872727919154569765469e-03L, +}; + +const long double pi_lo = 8.67181013012378102479704402604335225e-35L; diff --git a/lib/msun/ld128/invtrig.h b/lib/msun/ld128/invtrig.h new file mode 100644 index 0000000..12f598b --- /dev/null +++ b/lib/msun/ld128/invtrig.h @@ -0,0 +1,113 @@ +/*- + * 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 PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL + * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS + * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) + * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT + * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY + * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF + * SUCH DAMAGE. + * + * $FreeBSD$ + */ + +#include + +#include "fpmath.h" + +#define BIAS (LDBL_MAX_EXP - 1) +#define MANH_SIZE (LDBL_MANH_SIZE + 1) + +/* Approximation thresholds. */ +#define ASIN_LINEAR (BIAS - 56) /* 2**-56 */ +#define ACOS_CONST (BIAS - 113) /* 2**-113 */ +#define ATAN_CONST (BIAS + 113) /* 2**113 */ +#define ATAN_LINEAR (BIAS - 56) /* 2**-56 */ + +/* 0.95 */ +#define THRESH ((0xe666666666666666ULL>>(64-(MANH_SIZE-1)))|LDBL_NBIT) + +/* Constants shared by the long double inverse trig functions. */ +#define pS0 _ItL_pS0 +#define pS1 _ItL_pS1 +#define pS2 _ItL_pS2 +#define pS3 _ItL_pS3 +#define pS4 _ItL_pS4 +#define pS5 _ItL_pS5 +#define pS6 _ItL_pS6 +#define pS7 _ItL_pS7 +#define pS8 _ItL_pS8 +#define pS9 _ItL_pS9 +#define qS1 _ItL_qS1 +#define qS2 _ItL_qS2 +#define qS3 _ItL_qS3 +#define qS4 _ItL_qS4 +#define qS5 _ItL_qS5 +#define qS6 _ItL_qS6 +#define qS7 _ItL_qS7 +#define qS8 _ItL_qS8 +#define qS9 _ItL_qS9 +#define atanhi _ItL_atanhi +#define atanlo _ItL_atanlo +#define aT _ItL_aT +#define pi_lo _ItL_pi_lo + +#define pio2_hi atanhi[3] +#define pio2_lo atanlo[3] +#define pio4_hi atanhi[1] + +/* Constants shared by the long double inverse trig functions. */ +extern const long double pS0, pS1, pS2, pS3, pS4, pS5, pS6, pS7, pS8, pS9; +extern const long double qS1, qS2, qS3, qS4, qS5, qS6, qS7, qS8, qS9; +extern const long double atanhi[], atanlo[], aT[]; +extern const long double pi_lo; + +static inline long double +P(long double x) +{ + + return (x * (pS0 + x * (pS1 + x * (pS2 + x * (pS3 + x * \ + (pS4 + x * (pS5 + x * (pS6 + x * (pS7 + x * (pS8 + x * \ + pS9)))))))))); +} + +static inline long double +Q(long double x) +{ + + return (1.0 + x * (qS1 + x * (qS2 + x * (qS3 + x * (qS4 + x * \ + (qS5 + x * (qS6 + x * (qS7 + x * (qS8 + x * qS9))))))))); +} + +static inline long double +T_even(long double x) +{ + + return (aT[0] + x * (aT[2] + x * (aT[4] + x * (aT[6] + x * \ + (aT[8] + x * (aT[10] + x * (aT[12] + x * (aT[14] + x * \ + (aT[16] + x * (aT[18] + x * (aT[20] + x * aT[22]))))))))))); +} + +static inline long double +T_odd(long double x) +{ + + return (aT[1] + x * (aT[3] + x * (aT[5] + x * (aT[7] + x * \ + (aT[9] + x * (aT[11] + x * (aT[13] + x * (aT[15] + x * \ + (aT[17] + x * (aT[19] + x * (aT[21] + x * aT[23]))))))))))); +} diff --git a/lib/msun/ld80/invtrig.c b/lib/msun/ld80/invtrig.c new file mode 100644 index 0000000..06d0092 --- /dev/null +++ b/lib/msun/ld80/invtrig.c @@ -0,0 +1,82 @@ +/*- + * 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 PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL + * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS + * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) + * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT + * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY + * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF + * SUCH DAMAGE. + */ + +#include +__FBSDID("$FreeBSD$"); + +#include "invtrig.h" + +/* + * asinl() and acosl() + */ +const long double +pS0 = 1.66666666666666666631e-01L, +pS1 = -4.16313987993683104320e-01L, +pS2 = 3.69068046323246813704e-01L, +pS3 = -1.36213932016738603108e-01L, +pS4 = 1.78324189708471965733e-02L, +pS5 = -2.19216428382605211588e-04L, +pS6 = -7.10526623669075243183e-06L, +qS1 = -2.94788392796209867269e+00L, +qS2 = 3.27309890266528636716e+00L, +qS3 = -1.68285799854822427013e+00L, +qS4 = 3.90699412641738801874e-01L, +qS5 = -3.14365703596053263322e-02L; + +/* + * atanl() + */ +const long double atanhi[] = { + 4.63647609000806116202e-01L, + 7.85398163397448309628e-01L, + 9.82793723247329067960e-01L, + 1.57079632679489661926e+00L, +}; + +const long double atanlo[] = { + 1.18469937025062860669e-20L, + -1.25413940316708300586e-20L, + 2.55232234165405176172e-20L, + -2.50827880633416601173e-20L, +}; + +const long double aT[] = { + 3.33333333333333333017e-01L, + -1.99999999999999632011e-01L, + 1.42857142857046531280e-01L, + -1.11111111100562372733e-01L, + 9.09090902935647302252e-02L, + -7.69230552476207730353e-02L, + 6.66661718042406260546e-02L, + -5.88158892835030888692e-02L, + 5.25499891539726639379e-02L, + -4.70119845393155721494e-02L, + 4.03539201366454414072e-02L, + -2.91303858419364158725e-02L, + 1.24822046299269234080e-02L, +}; + +const long double pi_lo = -5.01655761266833202345e-20L; diff --git a/lib/msun/ld80/invtrig.h b/lib/msun/ld80/invtrig.h new file mode 100644 index 0000000..069a408 --- /dev/null +++ b/lib/msun/ld80/invtrig.h @@ -0,0 +1,101 @@ +/*- + * 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 PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL + * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS + * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) + * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT + * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY + * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF + * SUCH DAMAGE. + * + * $FreeBSD$ + */ + +#include + +#include "fpmath.h" + +#define BIAS (LDBL_MAX_EXP - 1) +#define MANH_SIZE LDBL_MANH_SIZE + +/* Approximation thresholds. */ +#define ASIN_LINEAR (BIAS - 32) /* 2**-32 */ +#define ACOS_CONST (BIAS - 65) /* 2**-65 */ +#define ATAN_CONST (BIAS + 65) /* 2**65 */ +#define ATAN_LINEAR (BIAS - 32) /* 2**-32 */ + +/* 0.95 */ +#define THRESH ((0xe666666666666666ULL>>(64-(MANH_SIZE-1)))|LDBL_NBIT) + +/* Constants shared by the long double inverse trig functions. */ +#define pS0 _ItL_pS0 +#define pS1 _ItL_pS1 +#define pS2 _ItL_pS2 +#define pS3 _ItL_pS3 +#define pS4 _ItL_pS4 +#define pS5 _ItL_pS5 +#define pS6 _ItL_pS6 +#define qS1 _ItL_qS1 +#define qS2 _ItL_qS2 +#define qS3 _ItL_qS3 +#define qS4 _ItL_qS4 +#define qS5 _ItL_qS5 +#define atanhi _ItL_atanhi +#define atanlo _ItL_atanlo +#define aT _ItL_aT +#define pi_lo _ItL_pi_lo + +#define pio2_hi atanhi[3] +#define pio2_lo atanlo[3] +#define pio4_hi atanhi[1] + +extern const long double pS0, pS1, pS2, pS3, pS4, pS5, pS6; +extern const long double qS1, qS2, qS3, qS4, qS5; +extern const long double atanhi[], atanlo[], aT[]; +extern const long double pi_lo; + +static inline long double +P(long double x) +{ + + return (x * (pS0 + x * (pS1 + x * (pS2 + x * (pS3 + x * \ + (pS4 + x * (pS5 + x * pS6))))))); +} + +static inline long double +Q(long double x) +{ + + return (1.0 + x * (qS1 + x * (qS2 + x * (qS3 + x * (qS4 + x * qS5))))); +} + +static inline long double +T_even(long double x) +{ + + return (aT[0] + x * (aT[2] + x * (aT[4] + x * (aT[6] + x * \ + (aT[8] + x * (aT[10] + x * aT[12])))))); +} + +static inline long double +T_odd(long double x) +{ + + return (aT[1] + x * (aT[3] + x * (aT[5] + x * (aT[7] + x * \ + (aT[9] + x * aT[11]))))); +} diff --git a/lib/msun/man/acos.3 b/lib/msun/man/acos.3 index aa73699..b38a6e0 100644 --- a/lib/msun/man/acos.3 +++ b/lib/msun/man/acos.3 @@ -28,12 +28,13 @@ .\" from: @(#)acos.3 5.1 (Berkeley) 5/2/91 .\" $FreeBSD$ .\" -.Dd January 14, 2005 +.Dd July 31, 2008 .Dt ACOS 3 .Os .Sh NAME .Nm acos , -.Nm acosf +.Nm acosf , +.Nm acosl .Nd arc cosine functions .Sh LIBRARY .Lb libm @@ -43,23 +44,18 @@ .Fn acos "double x" .Ft float .Fn acosf "float x" +.Ft long double +.Fn acosl "long double x" .Sh DESCRIPTION The -.Fn acos -and the -.Fn acosf +.Fn acos , +.Fn acosf , +and +.Fn acosl functions compute the principal value of the arc cosine of .Fa x . -A domain error occurs for arguments not in the range -.Bq -1 , +1 . -For a discussion of error due to roundoff, see -.Xr math 3 . .Sh RETURN VALUES -The -.Fn acos -and the -.Fn acosf -functions return the arc cosine in the range +These functions return the arc cosine in the range .Bq 0 , \*(Pi radians. If: @@ -83,6 +79,9 @@ raises an invalid exception and returns an \*(Na. .Xr tanh 3 .Sh STANDARDS The -.Fn acos -function conforms to -.St -isoC . +.Fn acos , +.Fn acosf , +and +.Fn acosl +functions conform to +.St -isoC-99 . diff --git a/lib/msun/man/asin.3 b/lib/msun/man/asin.3 index f7e8404..d9cad73 100644 --- a/lib/msun/man/asin.3 +++ b/lib/msun/man/asin.3 @@ -28,12 +28,13 @@ .\" from: @(#)asin.3 5.1 (Berkeley) 5/2/91 .\" $FreeBSD$ .\" -.Dd January 14, 2005 +.Dd July 31, 2008 .Dt ASIN 3 .Os .Sh NAME .Nm asin , -.Nm asinf +.Nm asinf , +.Nm asinl .Nd arc sine functions .Sh LIBRARY .Lb libm @@ -43,23 +44,18 @@ .Fn asin "double x" .Ft float .Fn asinf "float x" +.Ft long double +.Fn asinl "long double x" .Sh DESCRIPTION The -.Fn asin -and the -.Fn asinf +.Fn asin , +.Fn asinf , +and +.Fn asinl functions compute the principal value of the arc sine of .Fa x . -A domain error occurs for arguments not in the range -.Bq -1 , +1 . -For a discussion of error due to roundoff, see -.Xr math 3 . .Sh RETURN VALUES -The -.Fn asin -and the -.Fn asinf -functions return the arc sine in the range +These functions return the arc sine in the range .Bk -words .Bq -\*(Pi/2 , +\*(Pi/2 .Ek @@ -85,6 +81,9 @@ raises an invalid exception and returns an \*(Na. .Xr tanh 3 .Sh STANDARDS The -.Fn asin -function conforms to -.St -isoC . +.Fn asin , +.Fn asinf , +and +.Fn asinl +functions conform to +.St -isoC-99 . diff --git a/lib/msun/man/atan.3 b/lib/msun/man/atan.3 index 2c2bd64..ef119f1 100644 --- a/lib/msun/man/atan.3 +++ b/lib/msun/man/atan.3 @@ -28,12 +28,13 @@ .\" from: @(#)atan.3 5.1 (Berkeley) 5/2/91 .\" $FreeBSD$ .\" -.Dd May 2, 1991 +.Dd July 31, 2008 .Dt ATAN 3 .Os .Sh NAME .Nm atan , -.Nm atanf +.Nm atanf , +.Nm atanl .Nd arc tangent functions of one variable .Sh LIBRARY .Lb libm @@ -43,21 +44,18 @@ .Fn atan "double x" .Ft float .Fn atanf "float x" +.Ft long double +.Fn atanl "long double x" .Sh DESCRIPTION The -.Fn atan -and the -.Fn atanf +.Fn atan , +.Fn atanf , +and +.Fn atanl functions compute the principal value of the arc tangent of .Fa x . -For a discussion of error due to roundoff, see -.Xr math 3 . .Sh RETURN VALUES -The -.Fn atan -and the -.Fn atanf -function returns the arc tangent in the range +These functions return the arc tangent in the range .Bk -words .Bq -\*(Pi/2 , +\*(Pi/2 .Ek @@ -75,6 +73,9 @@ radians. .Xr tanh 3 .Sh STANDARDS The -.Fn atan -function conforms to -.St -isoC . +.Fn atan , +.Fn atanf , +and +.Fn atanl +functions conform to +.St -isoC-99 . diff --git a/lib/msun/man/atan2.3 b/lib/msun/man/atan2.3 index dee65ae..06ac068 100644 --- a/lib/msun/man/atan2.3 +++ b/lib/msun/man/atan2.3 @@ -28,14 +28,16 @@ .\" from: @(#)atan2.3 5.1 (Berkeley) 5/2/91 .\" $FreeBSD$ .\" -.Dd December 12, 2007 +.Dd July 31, 2008 .Dt ATAN2 3 .Os .Sh NAME .Nm atan2 , .Nm atan2f , +.Nm atan2l , .Nm carg , -.Nm cargf +.Nm cargf , +.Nm cargl .Nd arc tangent and complex phase angle functions .Sh LIBRARY .Lb libm @@ -45,15 +47,21 @@ .Fn atan2 "double y" "double x" .Ft float .Fn atan2f "float y" "float x" +.Ft long double +.Fn atan2l "long double y" "long double x" +.In complex.h .Ft double .Fn carg "double complex z" .Ft float .Fn cargf "float complex z" +.Ft long double +.Fn cargl "long double complex z" .Sh DESCRIPTION The -.Fn atan2 -and the -.Fn atan2f +.Fn atan2 , +.Fn atan2f , +and +.Fn atan2l functions compute the principal value of the arc tangent of .Fa y/ Ns Ar x , using the signs of both arguments to determine the quadrant of @@ -66,9 +74,10 @@ the return value. .\} .Pp The -.Fn carg +.Fn carg , +.Fn cargf , and -.Fn cargf +.Fn cargl functions compute the complex argument (or phase angle) of .Fa z . The complex argument is the number \*(Th such that @@ -80,12 +89,15 @@ The call is equivalent to .Li atan2(cimag(z), creal(z)) , and similarly for -.Fn cargf . +.Fn cargf +and +.Fn cargl . .Sh RETURN VALUES The -.Fn atan2 -and the -.Fn atan2f +.Fn atan2 , +.Fn atan2f , +and +.Fn atan2l functions, if successful, return the arc tangent of .Fa y/ Ns Ar x @@ -210,8 +222,10 @@ r := \(sr(x\(**x+y\(**y);\0\0if r = 0 then x := copysign(1,x); The .Fn atan2 , .Fn atan2f , +.Fn atan2l , .Fn carg , +.Fn cargf , and -.Fn cargf +.Fn cargl functions conform to .St -isoC-99 . diff --git a/lib/msun/src/e_acos.c b/lib/msun/src/e_acos.c index 5b1b85a..1f6dca5 100644 --- a/lib/msun/src/e_acos.c +++ b/lib/msun/src/e_acos.c @@ -38,6 +38,8 @@ __FBSDID("$FreeBSD$"); * Function needed: sqrt */ +#include + #include "math.h" #include "math_private.h" @@ -103,3 +105,7 @@ __ieee754_acos(double x) return 2.0*(df+w); } } + +#if LDBL_MANT_DIG == 53 +__weak_reference(acos, acosl); +#endif diff --git a/lib/msun/src/e_acosl.c b/lib/msun/src/e_acosl.c new file mode 100644 index 0000000..4e9df52 --- /dev/null +++ b/lib/msun/src/e_acosl.c @@ -0,0 +1,77 @@ + +/* @(#)e_acos.c 1.3 95/01/18 */ +/* FreeBSD: head/lib/msun/src/e_acos.c 176451 2008-02-22 02:30:36Z das */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +#include +__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= 1.00000000000000000000e+00, +pi = 3.14159265358979323846264338327950280e+00L; + +long double +acosl(long double x) +{ + union IEEEl2bits u; + long double z,p,q,r,w,s,c,df; + int16_t expsign, expt; + u.e = x; + expsign = u.xbits.expsign; + expt = expsign & 0x7fff; + if(expt >= BIAS) { /* |x| >= 1 */ + if(expt==BIAS && ((u.bits.manh&~LDBL_NBIT)|u.bits.manl)==0) { + if (expsign>0) return 0.0; /* acos(1) = 0 */ + else return pi+2.0*pio2_lo; /* acos(-1)= pi */ + } + return (x-x)/(x-x); /* acos(|x|>1) is NaN */ + } + if(expt 0.5 */ + z = (one-x)*0.5; + s = sqrtl(z); + u.e = s; + u.bits.manl = 0; + df = u.e; + c = (z-df*df)/(s+df); + p = P(z); + q = Q(z); + r = p/q; + w = r*s+c; + return 2.0*(df+w); + } +} diff --git a/lib/msun/src/e_asin.c b/lib/msun/src/e_asin.c index 149a533..4458604 100644 --- a/lib/msun/src/e_asin.c +++ b/lib/msun/src/e_asin.c @@ -44,6 +44,7 @@ __FBSDID("$FreeBSD$"); * */ +#include #include "math.h" #include "math_private.h" @@ -110,3 +111,7 @@ __ieee754_asin(double x) } if(hx>0) return t; else return -t; } + +#if LDBL_MANT_DIG == 53 +__weak_reference(asin, asinl); +#endif diff --git a/lib/msun/src/e_asinl.c b/lib/msun/src/e_asinl.c new file mode 100644 index 0000000..ee7f5af --- /dev/null +++ b/lib/msun/src/e_asinl.c @@ -0,0 +1,77 @@ + +/* @(#)e_asin.c 1.3 95/01/18 */ +/* FreeBSD: head/lib/msun/src/e_asin.c 176451 2008-02-22 02:30:36Z das */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +#include +__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 = 1.00000000000000000000e+00, +huge = 1.000e+300; + +long double +asinl(long double x) +{ + union IEEEl2bits u; + long double t=0.0,w,p,q,c,r,s; + int16_t expsign, expt; + u.e = x; + expsign = u.xbits.expsign; + expt = expsign & 0x7fff; + if(expt >= BIAS) { /* |x|>= 1 */ + if(expt==BIAS && ((u.bits.manh&~LDBL_NBIT)|u.bits.manl)==0) + /* asin(1)=+-pi/2 with inexact */ + return x*pio2_hi+x*pio2_lo; + return (x-x)/(x-x); /* asin(|x|>1) is NaN */ + } else if (exptone) return x;/* return x with inexact if x!=0*/ + } else + t = x*x; + p = P(t); + q = Q(t); + w = p/q; + return x+x*w; + } + /* 1> |x|>= 0.5 */ + w = one-fabsl(x); + t = w*0.5; + p = P(t); + q = Q(t); + s = sqrtl(t); + if(u.bits.manh>=THRESH) { /* if |x| is close to 1 */ + w = p/q; + t = pio2_hi-(2.0*(s+s*w)-pio2_lo); + } else { + u.e = s; + u.bits.manl = 0; + w = u.e; + c = (t-w*w)/(s+w); + r = p/q; + p = 2.0*s*r-(pio2_lo-2.0*c); + q = pio4_hi-2.0*w; + t = pio4_hi-(p-q); + } + if(expsign>0) return t; else return -t; +} diff --git a/lib/msun/src/e_atan2.c b/lib/msun/src/e_atan2.c index ba0b618..399822a 100644 --- a/lib/msun/src/e_atan2.c +++ b/lib/msun/src/e_atan2.c @@ -42,6 +42,8 @@ __FBSDID("$FreeBSD$"); * to produce the hexadecimal values shown. */ +#include + #include "math.h" #include "math_private.h" @@ -123,3 +125,7 @@ __ieee754_atan2(double y, double x) return (z-pi_lo)-pi;/* atan(-,-) */ } } + +#if LDBL_MANT_DIG == 53 +__weak_reference(atan2, atan2l); +#endif diff --git a/lib/msun/src/e_atan2l.c b/lib/msun/src/e_atan2l.c new file mode 100644 index 0000000..8368bbf --- /dev/null +++ b/lib/msun/src/e_atan2l.c @@ -0,0 +1,107 @@ + +/* @(#)e_atan2.c 1.3 95/01/18 */ +/* FreeBSD: head/lib/msun/src/e_atan2.c 176451 2008-02-22 02:30:36Z das */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + * + */ + +#include +__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 = 1.0e-300; +static const long double +zero = 0.0, +pi = 3.14159265358979323846264338327950280e+00L; + +long double +atan2l(long double y, long double x) +{ + union IEEEl2bits ux, uy; + long double z; + int32_t k,m; + int16_t exptx, expsignx, expty, expsigny; + + uy.e = y; + expsigny = uy.xbits.expsign; + expty = expsigny & 0x7fff; + ux.e = x; + expsignx = ux.xbits.expsign; + exptx = expsignx & 0x7fff; + + if ((exptx==BIAS+LDBL_MAX_EXP && + ((ux.bits.manh&~LDBL_NBIT)|ux.bits.manl)!=0) || /* x is NaN */ + (expty==BIAS+LDBL_MAX_EXP && + ((uy.bits.manh&~LDBL_NBIT)|uy.bits.manl)!=0)) /* y is NaN */ + return x+y; + if (expsignx==BIAS && ((ux.bits.manh&~LDBL_NBIT)|ux.bits.manl)==0) + return atanl(y); /* x=1.0 */ + m = ((expsigny>>15)&1)|((expsignx>>14)&2); /* 2*sign(x)+sign(y) */ + + /* when y = 0 */ + if(expty==0 && ((uy.bits.manh&~LDBL_NBIT)|uy.bits.manl)==0) { + switch(m) { + case 0: + case 1: return y; /* atan(+-0,+anything)=+-0 */ + case 2: return pi+tiny;/* atan(+0,-anything) = pi */ + case 3: return -pi-tiny;/* atan(-0,-anything) =-pi */ + } + } + /* when x = 0 */ + if(exptx==0 && ((ux.bits.manh&~LDBL_NBIT)|ux.bits.manl)==0) + return (expsigny<0)? -pio2_hi-tiny: pio2_hi+tiny; + + /* when x is INF */ + if(exptx==BIAS+LDBL_MAX_EXP) { + if(expty==BIAS+LDBL_MAX_EXP) { + switch(m) { + case 0: return pio2_hi*0.5+tiny;/* atan(+INF,+INF) */ + case 1: return -pio2_hi*0.5-tiny;/* atan(-INF,+INF) */ + case 2: return 1.5*pio2_hi+tiny;/*atan(+INF,-INF)*/ + case 3: return -1.5*pio2_hi-tiny;/*atan(-INF,-INF)*/ + } + } else { + switch(m) { + case 0: return zero ; /* atan(+...,+INF) */ + case 1: return -zero ; /* atan(-...,+INF) */ + case 2: return pi+tiny ; /* atan(+...,-INF) */ + case 3: return -pi-tiny ; /* atan(-...,-INF) */ + } + } + } + /* when y is INF */ + if(expty==BIAS+LDBL_MAX_EXP) + return (expsigny<0)? -pio2_hi-tiny: pio2_hi+tiny; + + /* compute y/x */ + k = expty-exptx; + if(k > LDBL_MANT_DIG+2) z=pio2_hi+pio2_lo; /* |y/x| huge */ + else if(expsignx<0&&k<-LDBL_MANT_DIG-2) z=0.0; /* |y/x| tiny, x<0 */ + else z=atanl(fabsl(y/x)); /* safe to do y/x */ + switch (m) { + case 0: return z ; /* atan(+,+) */ + case 1: return -z ; /* atan(-,+) */ + case 2: return pi-(z-pi_lo);/* atan(+,-) */ + default: /* case 3 */ + return (z-pi_lo)-pi;/* atan(-,-) */ + } +} diff --git a/lib/msun/src/math.h b/lib/msun/src/math.h index 97d1a58..374f95f 100644 --- a/lib/msun/src/math.h +++ b/lib/msun/src/math.h @@ -394,12 +394,18 @@ float significandf(float); #if __ISO_C_VISIBLE >= 1999 #if 0 long double acoshl(long double); +#endif long double acosl(long double); +#if 0 long double asinhl(long double); +#endif long double asinl(long double); long double atan2l(long double, long double); +#if 0 long double atanhl(long double); +#endif long double atanl(long double); +#if 0 long double cbrtl(long double); #endif long double ceill(long double); diff --git a/lib/msun/src/s_atan.c b/lib/msun/src/s_atan.c index 01e7bea..24daed5 100644 --- a/lib/msun/src/s_atan.c +++ b/lib/msun/src/s_atan.c @@ -33,6 +33,8 @@ __FBSDID("$FreeBSD$"); * to produce the hexadecimal values shown. */ +#include + #include "math.h" #include "math_private.h" @@ -116,3 +118,7 @@ atan(double x) return (hx<0)? -z:z; } } + +#if LDBL_MANT_DIG == 53 +__weak_reference(atan, atanl); +#endif diff --git a/lib/msun/src/s_atanl.c b/lib/msun/src/s_atanl.c new file mode 100644 index 0000000..ff29c3c --- /dev/null +++ b/lib/msun/src/s_atanl.c @@ -0,0 +1,85 @@ +/* @(#)s_atan.c 5.1 93/09/24 */ +/* FreeBSD: head/lib/msun/src/s_atan.c 176451 2008-02-22 02:30:36Z das */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +#include +__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 = 1.0, +huge = 1.0e300; + +long double +atanl(long double x) +{ + union IEEEl2bits u; + long double w,s1,s2,z; + int id; + int16_t expsign, expt; + int32_t expman; + + u.e = x; + expsign = u.xbits.expsign; + expt = expsign & 0x7fff; + if(expt >= ATAN_CONST) { /* if |x| is large, atan(x)~=pi/2 */ + if(expt == BIAS + LDBL_MAX_EXP && + ((u.bits.manh&~LDBL_NBIT)|u.bits.manl)!=0) + return x+x; /* NaN */ + if(expsign>0) return atanhi[3]+atanlo[3]; + else return -atanhi[3]-atanlo[3]; + } + /* Extract the exponent and the first few bits of the mantissa. */ + /* XXX There should be a more convenient way to do this. */ + expman = (expt << 8) | ((u.bits.manh >> (MANH_SIZE - 9)) & 0xff); + if (expman < ((BIAS - 2) << 8) + 0xc0) { /* |x| < 0.4375 */ + if (expt < ATAN_LINEAR) { /* if |x| is small, atanl(x)~=x */ + if(huge+x>one) return x; /* raise inexact */ + } + id = -1; + } else { + x = fabsl(x); + if (expman < (BIAS << 8) + 0x30) { /* |x| < 1.1875 */ + if (expman < ((BIAS - 1) << 8) + 0x60) { /* 7/16 <=|x|<11/16 */ + id = 0; x = (2.0*x-one)/(2.0+x); + } else { /* 11/16<=|x|< 19/16 */ + id = 1; x = (x-one)/(x+one); + } + } else { + if (expman < ((BIAS + 1) << 8) + 0x38) { /* |x| < 2.4375 */ + id = 2; x = (x-1.5)/(one+1.5*x); + } else { /* 2.4375 <= |x| < 2^ATAN_CONST */ + id = 3; x = -1.0/x; + } + }} + /* end of argument reduction */ + z = x*x; + w = z*z; + /* break sum aT[i]z**(i+1) into odd and even poly */ + s1 = z*T_even(w); + s2 = w*T_odd(w); + if (id<0) return x - x*(s1+s2); + else { + z = atanhi[id] - ((x*(s1+s2) - atanlo[id]) - x); + return (expsign<0)? -z:z; + } +} diff --git a/lib/msun/src/s_cargl.c b/lib/msun/src/s_cargl.c new file mode 100644 index 0000000..0555083 --- /dev/null +++ b/lib/msun/src/s_cargl.c @@ -0,0 +1,38 @@ +/*- + * Copyright (c) 2005-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 PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL + * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS + * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) + * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT + * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY + * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF + * SUCH DAMAGE. + */ + +#include +__FBSDID("$FreeBSD$"); + +#include +#include + +long double +cargl(long double complex z) +{ + + return (atan2l(cimagl(z), creall(z))); +} -- cgit v1.1