libm: Sync with FreeBSD (gains 6 long double functions)
authorJohn Marino <draco@marino.st>
Sat, 29 Nov 2014 21:28:11 +0000 (22:28 +0100)
committerJohn Marino <draco@marino.st>
Sat, 29 Nov 2014 23:00:28 +0000 (00:00 +0100)
The following functions have been implemented:
  - coshl
  - erfcl
  - erfl
  - lgammal
  - sinhl
  - tanhl

Before these were approximated with the double versions using the
"imprecise" macros.  I've left the old ones in place (unlike FreeBSD)
but with symbol versioning so that libraries built with the earlier
versions can link to them.  In other words, there are two versions of
these 6 functions, Df306.1 and DF402.0.

55 files changed:
lib/libc/Versions.def
lib/libm/Makefile
lib/libm/Symbol.map
lib/libm/ld80/e_lgammal_r.c [new file with mode: 0644]
lib/libm/ld80/k_expl.h [copied from lib/libm/ld80/s_expl.c with 57% similarity]
lib/libm/ld80/s_erfl.c [new file with mode: 0644]
lib/libm/ld80/s_expl.c
lib/libm/man/cosh.3
lib/libm/man/erf.3
lib/libm/man/lgamma.3
lib/libm/man/sinh.3
lib/libm/man/tanh.3
lib/libm/src/catrig.c
lib/libm/src/catrigf.c
lib/libm/src/e_cosh.c
lib/libm/src/e_coshl.c [new file with mode: 0644]
lib/libm/src/e_gamma.c
lib/libm/src/e_lgamma.c
lib/libm/src/e_lgamma_r.c
lib/libm/src/e_lgammaf_r.c
lib/libm/src/e_lgammal.c [copied from lib/libm/src/e_lgamma.c with 56% similarity]
lib/libm/src/e_pow.c
lib/libm/src/e_remainder.c [new file with mode: 0644]
lib/libm/src/e_remainderf.c [new file with mode: 0644]
lib/libm/src/e_remainderl.c [copied from lib/libm/src/imprecise.c with 56% similarity]
lib/libm/src/e_sinh.c
lib/libm/src/e_sinhl.c [new file with mode: 0644]
lib/libm/src/e_sqrt.c [new file with mode: 0644]
lib/libm/src/e_sqrtf.c [new file with mode: 0644]
lib/libm/src/e_sqrtl.c [new file with mode: 0644]
lib/libm/src/fenv-softfloat.h [new file with mode: 0644]
lib/libm/src/imprecise.c
lib/libm/src/math.h
lib/libm/src/s_erf.c
lib/libm/src/s_erff.c
lib/libm/src/s_llrint.c [new file with mode: 0644]
lib/libm/src/s_llrintf.c [new file with mode: 0644]
lib/libm/src/s_llrintl.c [new file with mode: 0644]
lib/libm/src/s_logbl.c [new file with mode: 0644]
lib/libm/src/s_lrint.c [copied from lib/libm/src/imprecise.c with 56% similarity]
lib/libm/src/s_lrintf.c [new file with mode: 0644]
lib/libm/src/s_lrintl.c [new file with mode: 0644]
lib/libm/src/s_remquo.c [new file with mode: 0644]
lib/libm/src/s_remquof.c [new file with mode: 0644]
lib/libm/src/s_remquol.c [new file with mode: 0644]
lib/libm/src/s_rintl.c [new file with mode: 0644]
lib/libm/src/s_round.c
lib/libm/src/s_roundf.c
lib/libm/src/s_roundl.c
lib/libm/src/s_scalbn.c [new file with mode: 0644]
lib/libm/src/s_scalbnf.c [new file with mode: 0644]
lib/libm/src/s_scalbnl.c [new file with mode: 0644]
lib/libm/src/s_tanh.c
lib/libm/src/s_tanhf.c
lib/libm/src/s_tanhl.c [new file with mode: 0644]

index 0c9cd87..13301f7 100644 (file)
@@ -13,6 +13,10 @@ DF306.0 {
 DF306.1 {
 } DF306.0;
 
+# version for DragonFly 4.1
+DF402.0 {
+} DF306.1;
+
 # This is our private namespace.  Any global interfaces that are
 # strictly for use only by other DragonFly applications and libraries
 # are listed here.  We use a separate namespace so we can write
@@ -20,4 +24,4 @@ DF306.1 {
 #
 # Please do NOT increment the version of this namespace.
 DFprivate_1.0 {
-} DF306.1;
+} DF402.0;
index 4af61d0..3b25ca7 100644 (file)
@@ -35,8 +35,8 @@ COMMON_SRCS= b_exp.c b_log.c b_tgamma.c \
        e_jn.c e_jnf.c e_lgamma.c e_lgamma_r.c e_lgammaf.c e_lgammaf_r.c \
        e_log.c e_log10.c e_log10f.c e_log2.c e_log2f.c e_logf.c \
        e_pow.c e_powf.c e_rem_pio2.c \
-       e_rem_pio2f.c e_scalb.c e_scalbf.c \
-       e_sinh.c e_sinhf.c fenv.c \
+       e_rem_pio2f.c e_remainder.c e_remainderf.c e_scalb.c e_scalbf.c \
+       e_sinh.c e_sinhf.c e_sqrt.c e_sqrtf.c fenv.c \
        imprecise.c \
        k_cos.c k_cosf.c k_exp.c k_expf.c k_rem_pio2.c k_sin.c k_sinf.c \
        k_tan.c k_tanf.c \
@@ -50,13 +50,13 @@ COMMON_SRCS= b_exp.c b_log.c b_tgamma.c \
        s_fmax.c s_fmaxf.c s_fmaxl.c s_fmin.c \
        s_fminf.c s_fminl.c s_frexp.c s_frexpf.c s_ilogb.c s_ilogbf.c \
        s_ilogbl.c s_isfinite.c s_isnan.c s_isnormal.c \
-       s_llround.c s_llroundf.c s_llroundl.c \
-       s_log1p.c s_log1pf.c s_logb.c s_logbf.c \
+       s_llrint.c s_llrintf.c s_llround.c s_llroundf.c s_llroundl.c \
+       s_log1p.c s_log1pf.c s_logb.c s_logbf.c s_lrint.c s_lrintf.c \
        s_lround.c s_lroundf.c s_lroundl.c s_modff.c \
        s_nan.c s_nearbyint.c s_nextafter.c s_nextafterf.c \
-       s_nexttowardf.c \
-       s_rint.c s_rintf.c s_round.c s_roundf.c s_roundl.c \
-       s_scalbln.c s_signbit.c \
+       s_nexttowardf.c s_remquo.c s_remquof.c \
+       s_rint.c s_rintf.c s_round.c s_roundf.c \
+       s_scalbln.c s_scalbn.c s_scalbnf.c s_signbit.c \
        s_signgam.c s_significand.c s_significandf.c s_sin.c s_sinf.c \
        s_tan.c s_tanf.c s_tanh.c s_tanhf.c s_tgammaf.c s_trunc.c s_truncf.c \
        w_cabs.c w_cabsf.c w_drem.c w_dremf.c
@@ -70,16 +70,18 @@ VERSION_DEF=        ${LIBCDIR}/Versions.def
 SYMBOL_MAPS=   ${SYM_MAPS}
 
 # C99 long double functions
-COMMON_SRCS+=  s_copysignl.c s_fabsl.c s_modfl.c
+COMMON_SRCS+=  s_copysignl.c s_fabsl.c s_llrintl.c s_lrintl.c s_modfl.c
 # If long double != double use these; otherwise, we alias the double versions.
 COMMON_SRCS+=  e_acoshl.c e_acosl.c e_asinl.c e_atan2l.c e_atanhl.c \
-       e_fmodl.c e_hypotl.c \
+       e_coshl.c e_fmodl.c e_hypotl.c \
+       e_lgammal.c e_lgammal_r.c \
+       e_remainderl.c e_sinhl.c e_sqrtl.c \
        invtrig.c k_cosl.c k_sinl.c k_tanl.c \
        s_asinhl.c s_atanl.c s_cbrtl.c s_ceill.c s_cosl.c s_cprojl.c \
-       s_csqrtl.c s_exp2l.c s_expl.c s_floorl.c s_fmal.c \
-       s_frexpl.c s_logl.c s_nanl.c s_nextafterl.c \
-       s_nexttoward.c \
-       s_sinl.c s_tanl.c s_truncl.c w_cabsl.c
+       s_csqrtl.c s_erfl.c s_exp2l.c s_expl.c s_floorl.c s_fmal.c \
+       s_frexpl.c s_logbl.c s_logl.c s_nanl.c s_nextafterl.c \
+       s_nexttoward.c s_remquol.c s_rintl.c s_roundl.c s_scalbnl.c \
+       s_sinl.c s_tanhl.c s_tanl.c s_truncl.c w_cabsl.c
 
 # C99 complex functions
 COMMON_SRCS+=  catrig.c catrigf.c \
@@ -104,13 +106,15 @@ SRCS=     ${COMMON_SRCS} ${ARCH_SRCS}
 INCS+= complex.h fenv.h math.h
 
 MAN=   acos.3 acosh.3 asin.3 asinh.3 atan.3 atan2.3 atanh.3 \
-       ceil.3 cacos.3 ccos.3 ccosh.3 cexp.3 cimag.3 complex.3 copysign.3 \
-       cos.3 cosh.3 csqrt.3 erf.3 exp.3 fabs.3 fdim.3 \
-       feclearexcept.3 feenableexcept.3 fegetenv.3 fegetround.3 fenv.3 \
-       floor.3 fma.3 fmax.3 fmod.3 hypot.3 ieee.3 ieee_test.3 ilogb.3 \
-       j0.3 lgamma.3 log.3 lrint.3 lround.3 math.3 nan.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
+       ceil.3 cacos.3 ccos.3 ccosh.3 cexp.3 \
+       cimag.3 copysign.3 cos.3 cosh.3 csqrt.3 erf.3 exp.3 fabs.3 fdim.3 \
+       feclearexcept.3 feenableexcept.3 fegetenv.3 \
+       fegetround.3 fenv.3 floor.3 \
+       fma.3 fmax.3 fmod.3 hypot.3 ieee.3 ieee_test.3 ilogb.3 j0.3 \
+       lgamma.3 log.3 lrint.3 lround.3 math.3 nan.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 \
+       complex.3
 
 MLINKS+=acos.3 acosf.3 \
        acos.3 acosl.3
@@ -122,13 +126,13 @@ MLINKS+=asinh.3 asinhf.3 \
        asinh.3 asinhl.3
 MLINKS+=atan.3 atanf.3 \
        atan.3 atanl.3
+MLINKS+=atanh.3 atanhf.3 \
+       atanh.3 atanhl.3
 MLINKS+=atan2.3 atan2f.3 \
        atan2.3 atan2l.3 \
        atan2.3 carg.3 \
        atan2.3 cargf.3 \
        atan2.3 cargl.3
-MLINKS+=atanh.3 atanhf.3 \
-       atanh.3 atanhl.3
 MLINKS+=cacos.3 cacosf.3 \
        cacos.3 cacosh.3 \
        cacos.3 cacoshf.3 \
@@ -168,22 +172,25 @@ MLINKS+=copysign.3 copysignf.3 \
        copysign.3 copysignl.3
 MLINKS+=cos.3 cosf.3 \
        cos.3 cosl.3
-MLINKS+=cosh.3 coshf.3
+MLINKS+=cosh.3 coshf.3 \
+       cosh.3 coshl.3
 MLINKS+=csqrt.3 csqrtf.3 \
        csqrt.3 csqrtl.3
-MLINKS+=erf.3 erff.3 \
-       erf.3 erfc.3 \
-       erf.3 erfcf.3
-MLINKS+=exp.3 expf.3 \
-       exp.3 expl.3 \
-       exp.3 exp2.3 \
-       exp.3 exp2f.3 \
-       exp.3 exp2l.3 \
-       exp.3 expm1.3 \
+MLINKS+=erf.3 erfc.3 \
+       erf.3 erff.3 \
+       erf.3 erfcf.3 \
+       erf.3 erfl.3 \
+       erf.3 erfcl.3
+MLINKS+=exp.3 expm1.3 \
        exp.3 expm1f.3 \
        exp.3 expm1l.3 \
        exp.3 pow.3 \
-       exp.3 powf.3
+       exp.3 powf.3 \
+       exp.3 exp2.3 \
+       exp.3 exp2f.3 \
+       exp.3 exp2l.3 \
+       exp.3 expf.3 \
+       exp.3 expl.3
 MLINKS+=fabs.3 fabsf.3 \
        fabs.3 fabsl.3
 MLINKS+=fdim.3 fdimf.3 \
@@ -214,47 +221,43 @@ MLINKS+=hypot.3 cabs.3 \
        hypot.3 cabsl.3 \
        hypot.3 hypotf.3 \
        hypot.3 hypotl.3
-MLINKS+=ieee.3 libm.3
 MLINKS+=ieee_test.3 scalb.3 \
-       ieee_test.3 scalbf.3 \
-       ieee_test.3 significand.3 \
+       ieee_test.3 scalbf.3
+MLINKS+=ieee_test.3 significand.3 \
        ieee_test.3 significandf.3
 MLINKS+=ilogb.3 ilogbf.3 \
        ilogb.3 ilogbl.3 \
        ilogb.3 logb.3 \
        ilogb.3 logbf.3 \
        ilogb.3 logbl.3
-MLINKS+=j0.3 j0f.3 \
-       j0.3 j1.3 \
-       j0.3 j1f.3 \
+MLINKS+=j0.3 j1.3 \
        j0.3 jn.3 \
-       j0.3 jnf.3 \
        j0.3 y0.3 \
-       j0.3 y0f.3 \
        j0.3 y1.3 \
        j0.3 y1f.3 \
-       j0.3 yn.3 \
+       j0.3 yn.3
+MLINKS+=j0.3 j0f.3 \
+       j0.3 j1f.3 \
+       j0.3 jnf.3 \
+       j0.3 y0f.3 \
        j0.3 ynf.3
-MLINKS+=lgamma.3 lgamma_r.3 \
-       lgamma.3 lgammaf.3 \
-       lgamma.3 lgammaf_r.3 \
-       lgamma.3 gamma.3 \
-       lgamma.3 gamma_r.3 \
+MLINKS+=lgamma.3 gamma.3 \
        lgamma.3 gammaf.3 \
-       lgamma.3 gammaf_r.3 \
+       lgamma.3 lgammaf.3 \
+       lgamma.3 lgammal.3 \
        lgamma.3 tgamma.3 \
        lgamma.3 tgammaf.3
-MLINKS+=log.3 logf.3 \
-       log.3 logl.3 \
-       log.3 log10.3 \
+MLINKS+=log.3 log10.3 \
        log.3 log10f.3 \
        log.3 log10l.3 \
-       log.3 log2.3 \
-       log.3 log2f.3 \
-       log.3 log2l.3 \
        log.3 log1p.3 \
        log.3 log1pf.3 \
-       log.3 log1pl.3
+       log.3 log1pl.3 \
+       log.3 logf.3 \
+       log.3 logl.3 \
+       log.3 log2.3 \
+       log.3 log2f.3 \
+       log.3 log2l.3
 MLINKS+=lrint.3 llrint.3 \
        lrint.3 llrintf.3 \
        lrint.3 llrintl.3 \
@@ -268,10 +271,10 @@ MLINKS+=lround.3 llround.3 \
 MLINKS+=nan.3 nanf.3 \
        nan.3 nanl.3
 MLINKS+=nextafter.3 nextafterf.3 \
-       nextafter.3 nextafterl.3 \
-       nextafter.3 nexttoward.3 \
-       nextafter.3 nexttowardf.3 \
-       nextafter.3 nexttowardl.3
+       nextafter.3 nextafterl.3
+MLINKS+=nextafter.3 nexttoward.3 \
+       nextafter.3 nexttowardf.3
+MLINKS+=nextafter.3 nexttowardl.3
 MLINKS+=remainder.3 remainderf.3 \
        remainder.3 remainderl.3 \
        remainder.3 remquo.3 \
@@ -286,12 +289,13 @@ MLINKS+=round.3 roundf.3 \
        round.3 roundl.3
 MLINKS+=scalbn.3 scalbln.3 \
        scalbn.3 scalblnf.3 \
-       scalbn.3 scalblnl.3 \
-       scalbn.3 scalbnf.3 \
+       scalbn.3 scalblnl.3
+MLINKS+=scalbn.3 scalbnf.3 \
        scalbn.3 scalbnl.3
 MLINKS+=sin.3 sinf.3 \
        sin.3 sinl.3
-MLINKS+=sinh.3 sinhf.3
+MLINKS+=sinh.3 sinhf.3 \
+       sinh.3 sinhl.3
 MLINKS+=sqrt.3 cbrt.3 \
        sqrt.3 cbrtf.3 \
        sqrt.3 cbrtl.3 \
@@ -299,7 +303,8 @@ MLINKS+=sqrt.3 cbrt.3 \
        sqrt.3 sqrtl.3
 MLINKS+=tan.3 tanf.3 \
        tan.3 tanl.3
-MLINKS+=tanh.3 tanhf.3
+MLINKS+=tanh.3 tanhf.3 \
+       tanh.3 tanhl.3
 MLINKS+=trunc.3 truncf.3 \
        trunc.3 truncl.3
 
index 6418f90..58b312c 100644 (file)
@@ -264,12 +264,28 @@ DF306.0 {
  */
 
 DF306.1 {
+       /*
+        * These old versions are implemented with __asm__ .symver
+        * coshl;
+        * erfcl;
+        * erfl;
+        * lgammal;
+        * sinhl;
+        * tanhl;
+        */
+       powl;
+       tgammal;
+};
+
+/*
+ * Implemented in 4.1-dev branch
+ */
+
+DF402.0 {
        coshl;
        erfcl;
        erfl;
        lgammal;
-       powl;
        sinhl;
        tanhl;
-       tgammal;
 };
diff --git a/lib/libm/ld80/e_lgammal_r.c b/lib/libm/ld80/e_lgammal_r.c
new file mode 100644 (file)
index 0000000..1a053d8
--- /dev/null
@@ -0,0 +1,356 @@
+/* @(#)e_lgamma_r.c 1.3 95/01/18 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+
+/*
+ * See e_lgamma_r.c for complete comments.
+ *
+ * Converted to long double by Steven G. Kargl.
+ */
+
+#ifdef __i386__
+#include <ieeefp.h>
+#endif
+
+#include "fpmath.h"
+#include "math.h"
+#include "math_private.h"
+
+static const volatile double vzero = 0;
+
+static const double
+zero=  0,
+half=  0.5,
+one =  1;
+
+static const union IEEEl2bits
+piu = LD80C(0xc90fdaa22168c235, 1,  3.14159265358979323851e+00L);
+#define        pi      (piu.e)
+/*
+ * Domain y in [0x1p-70, 0.27], range ~[-4.5264e-22, 4.5264e-22]:
+ * |(lgamma(2 - y) + y / 2) / y - a(y)| < 2**-70.9
+ */
+static const union IEEEl2bits
+a0u = LD80C(0x9e233f1bed863d26, -4,  7.72156649015328606028e-02L),
+a1u = LD80C(0xa51a6625307d3249, -2,  3.22467033424113218889e-01L),
+a2u = LD80C(0x89f000d2abafda8c, -4,  6.73523010531979398946e-02L),
+a3u = LD80C(0xa8991563eca75f26, -6,  2.05808084277991211934e-02L),
+a4u = LD80C(0xf2027e10634ce6b6, -8,  7.38555102796070454026e-03L),
+a5u = LD80C(0xbd6eb76dd22187f4, -9,  2.89051035162703932972e-03L),
+a6u = LD80C(0x9c562ab05e0458ed, -10,  1.19275351624639999297e-03L),
+a7u = LD80C(0x859baed93ee48e46, -11,  5.09674593842117925320e-04L),
+a8u = LD80C(0xe9f28a4432949af2, -13,  2.23109648015769155122e-04L),
+a9u = LD80C(0xd12ad0d9b93c6bb0, -14,  9.97387167479808509830e-05L),
+a10u= LD80C(0xb7522643c78a219b, -15,  4.37071076331030136818e-05L),
+a11u= LD80C(0xca024dcdece2cb79, -16,  2.40813493372040143061e-05L),
+a12u= LD80C(0xbb90fb6968ebdbf9, -19,  2.79495621083634031729e-06L),
+a13u= LD80C(0xba1c9ffeeae07b37, -17,  1.10931287015513924136e-05L);
+#define        a0      (a0u.e)
+#define        a1      (a1u.e)
+#define        a2      (a2u.e)
+#define        a3      (a3u.e)
+#define        a4      (a4u.e)
+#define        a5      (a5u.e)
+#define        a6      (a6u.e)
+#define        a7      (a7u.e)
+#define        a8      (a8u.e)
+#define        a9      (a9u.e)
+#define        a10     (a10u.e)
+#define        a11     (a11u.e)
+#define        a12     (a12u.e)
+#define        a13     (a13u.e)
+/*
+ * Domain x in [tc-0.24, tc+0.28], range ~[-6.1205e-22, 6.1205e-22]:
+ * |(lgamma(x) - tf) -  t(x - tc)| < 2**-70.5
+ */
+static const union IEEEl2bits
+tcu  = LD80C(0xbb16c31ab5f1fb71, 0,  1.46163214496836234128e+00L),
+tfu  = LD80C(0xf8cdcde61c520e0f, -4, -1.21486290535849608093e-01L),
+ttu  = LD80C(0xd46ee54b27d4de99, -69, -2.81152980996018785880e-21L),
+t0u  = LD80C(0x80b9406556a62a6b, -68,  3.40728634996055147231e-21L),
+t1u  = LD80C(0xc7e9c6f6df3f8c39, -67, -1.05833162742737073665e-20L),
+t2u  = LD80C(0xf7b95e4771c55d51, -2,  4.83836122723810583532e-01L),
+t3u  = LD80C(0x97213c6e35e119ff, -3, -1.47587722994530691476e-01L),
+t4u  = LD80C(0x845a14a6a81dc94b, -4,  6.46249402389135358063e-02L),
+t5u  = LD80C(0x864d46fa89997796, -5, -3.27885410884846056084e-02L),
+t6u  = LD80C(0x93373cbd00297438, -6,  1.79706751150707171293e-02L),
+t7u  = LD80C(0xa8fcfca7eddc8d1d, -7, -1.03142230361450732547e-02L),
+t8u  = LD80C(0xc7e7015ff4bc45af, -8,  6.10053603296546099193e-03L),
+t9u  = LD80C(0xf178d2247adc5093, -9, -3.68456964904901200152e-03L),
+t10u = LD80C(0x94188d58f12e5e9f, -9,  2.25976420273774583089e-03L),
+t11u = LD80C(0xb7cbaef14e1406f1, -10, -1.40224943666225639823e-03L),
+t12u = LD80C(0xe63a671e6704ea4d, -11,  8.78250640744776944887e-04L),
+t13u = LD80C(0x914b6c9cae61783e, -11, -5.54255012657716808811e-04L),
+t14u = LD80C(0xb858f5bdb79276fe, -12,  3.51614951536825927370e-04L),
+t15u = LD80C(0xea73e744c34b9591, -13, -2.23591563824520112236e-04L),
+t16u = LD80C(0x99aeabb0d67ba835, -13,  1.46562869351659194136e-04L),
+t17u = LD80C(0xd7c6938325db2024, -14, -1.02889866046435680588e-04L),
+t18u = LD80C(0xe24cb1e3b0474775, -15,  5.39540265505221957652e-05L);
+#define        tc      (tcu.e)
+#define        tf      (tfu.e)
+#define        tt      (ttu.e)
+#define        t0      (t0u.e)
+#define        t1      (t1u.e)
+#define        t2      (t2u.e)
+#define        t3      (t3u.e)
+#define        t4      (t4u.e)
+#define        t5      (t5u.e)
+#define        t6      (t6u.e)
+#define        t7      (t7u.e)
+#define        t8      (t8u.e)
+#define        t9      (t9u.e)
+#define        t10     (t10u.e)
+#define        t11     (t11u.e)
+#define        t12     (t12u.e)
+#define        t13     (t13u.e)
+#define        t14     (t14u.e)
+#define        t15     (t15u.e)
+#define        t16     (t16u.e)
+#define        t17     (t17u.e)
+#define        t18     (t18u.e)
+/*
+ * Domain y in [-0.1, 0.232], range ~[-8.1938e-22, 8.3815e-22]:
+ * |(lgamma(1 + y) + 0.5 * y) / y - u(y) / v(y)| < 2**-71.2
+ */
+static const union IEEEl2bits
+u0u = LD80C(0x9e233f1bed863d27, -4, -7.72156649015328606095e-02L),
+u1u = LD80C(0x98280ee45e4ddd3d, -1,  5.94361239198682739769e-01L),
+u2u = LD80C(0xe330c8ead4130733, 0,  1.77492629495841234275e+00L),
+u3u = LD80C(0xd4a213f1a002ec52, 0,  1.66119622514818078064e+00L),
+u4u = LD80C(0xa5a9ca6f5bc62163, -1,  6.47122051417476492989e-01L),
+u5u = LD80C(0xc980e49cd5b019e6, -4,  9.83903751718671509455e-02L),
+u6u = LD80C(0xff636a8bdce7025b, -9,  3.89691687802305743450e-03L),
+v1u = LD80C(0xbd109c533a19fbf5, 1,  2.95413883330948556544e+00L),
+v2u = LD80C(0xd295cbf96f31f099, 1,  3.29039286955665403176e+00L),
+v3u = LD80C(0xdab8bcfee40496cb, 0,  1.70876276441416471410e+00L),
+v4u = LD80C(0xd2f2dc3638567e9f, -2,  4.12009126299534668571e-01L),
+v5u = LD80C(0xa07d9b0851070f41, -5,  3.91822868305682491442e-02L),
+v6u = LD80C(0xe3cd8318f7adb2c4, -11,  8.68998648222144351114e-04L);
+#define        u0      (u0u.e)
+#define        u1      (u1u.e)
+#define        u2      (u2u.e)
+#define        u3      (u3u.e)
+#define        u4      (u4u.e)
+#define        u5      (u5u.e)
+#define        u6      (u6u.e)
+#define        v1      (v1u.e)
+#define        v2      (v2u.e)
+#define        v3      (v3u.e)
+#define        v4      (v4u.e)
+#define        v5      (v5u.e)
+#define        v6      (v6u.e)
+/*
+ * Domain x in (2, 3], range ~[-3.3648e-22, 3.4416e-22]:
+ * |(lgamma(y+2) - 0.5 * y) / y - s(y)/r(y)| < 2**-72.3
+ * with y = x - 2.
+ */
+static const union IEEEl2bits
+s0u = LD80C(0x9e233f1bed863d27, -4, -7.72156649015328606095e-02L),
+s1u = LD80C(0xd3ff0dcc7fa91f94, -3,  2.07027640921219389860e-01L),
+s2u = LD80C(0xb2bb62782478ef31, -2,  3.49085881391362090549e-01L),
+s3u = LD80C(0xb49f7438c4611a74, -3,  1.76389518704213357954e-01L),
+s4u = LD80C(0x9a957008fa27ecf9, -5,  3.77401710862930008071e-02L),
+s5u = LD80C(0xda9b389a6ca7a7ac, -9,  3.33566791452943399399e-03L),
+s6u = LD80C(0xbc7a2263faf59c14, -14,  8.98728786745638844395e-05L),
+r1u = LD80C(0xbf5cff5b11477d4d, 0,  1.49502555796294337722e+00L),
+r2u = LD80C(0xd9aec89de08e3da6, -1,  8.50323236984473285866e-01L),
+r3u = LD80C(0xeab7ae5057c443f9, -3,  2.29216312078225806131e-01L),
+r4u = LD80C(0xf29707d9bd2b1e37, -6,  2.96130326586640089145e-02L),
+r5u = LD80C(0xd376c2f09736c5a3, -10,  1.61334161411590662495e-03L),
+r6u = LD80C(0xc985983d0cd34e3d, -16,  2.40232770710953450636e-05L),
+r7u = LD80C(0xe5c7a4f7fc2ef13d, -25, -5.34997929289167573510e-08L);
+#define        s0      (s0u.e)
+#define        s1      (s1u.e)
+#define        s2      (s2u.e)
+#define        s3      (s3u.e)
+#define        s4      (s4u.e)
+#define        s5      (s5u.e)
+#define        s6      (s6u.e)
+#define        r1      (r1u.e)
+#define        r2      (r2u.e)
+#define        r3      (r3u.e)
+#define        r4      (r4u.e)
+#define        r5      (r5u.e)
+#define        r6      (r6u.e)
+#define        r7      (r7u.e)
+/*
+ * Domain z in [8, 0x1p70], range ~[-3.0235e-22, 3.0563e-22]:
+ * |lgamma(x) - (x - 0.5) * (log(x) - 1) - w(1/x)| < 2**-71.7
+ */
+static const union IEEEl2bits
+w0u = LD80C(0xd67f1c864beb4a69, -2,  4.18938533204672741776e-01L),
+w1u = LD80C(0xaaaaaaaaaaaaaaa1, -4,  8.33333333333333332678e-02L),
+w2u = LD80C(0xb60b60b60b5491c9, -9, -2.77777777777760927870e-03L),
+w3u = LD80C(0xd00d00cf58aede4c, -11,  7.93650793490637233668e-04L),
+w4u = LD80C(0x9c09bf626783d4a5, -11, -5.95238023926039051268e-04L),
+w5u = LD80C(0xdca7cadc5baa517b, -11,  8.41733700408000822962e-04L),
+w6u = LD80C(0xfb060e361e1ffd07, -10, -1.91515849570245136604e-03L),
+w7u = LD80C(0xcbd5101bb58d1f2b, -8,  6.22046743903262649294e-03L),
+w8u = LD80C(0xad27a668d32c821b, -6, -2.11370706734662081843e-02L);
+#define        w0      (w0u.e)
+#define        w1      (w1u.e)
+#define        w2      (w2u.e)
+#define        w3      (w3u.e)
+#define        w4      (w4u.e)
+#define        w5      (w5u.e)
+#define        w6      (w6u.e)
+#define        w7      (w7u.e)
+#define        w8      (w8u.e)
+
+static long double
+sin_pil(long double x)
+{
+       volatile long double vz;
+       long double y,z;
+       uint64_t n;
+       uint16_t hx;
+
+       y = -x;
+
+       vz = y+0x1p63;
+       z = vz-0x1p63;
+       if (z == y)
+           return zero;
+
+       vz = y+0x1p61;
+       EXTRACT_LDBL80_WORDS(hx,n,vz);
+       z = vz-0x1p61;
+       if (z > y) {
+           z -= 0.25;                  /* adjust to round down */
+           n--;
+       }
+       n &= 7;                         /* octant of y mod 2 */
+       y = y - z + n * 0.25;           /* y mod 2 */
+
+       switch (n) {
+           case 0:   y =  __kernel_sinl(pi*y,zero,0); break;
+           case 1:
+           case 2:   y =  __kernel_cosl(pi*(0.5-y),zero); break;
+           case 3:
+           case 4:   y =  __kernel_sinl(pi*(one-y),zero,0); break;
+           case 5:
+           case 6:   y = -__kernel_cosl(pi*(y-1.5),zero); break;
+           default:  y =  __kernel_sinl(pi*(y-2.0),zero,0); break;
+           }
+       return -y;
+}
+
+long double
+lgammal_r(long double x, int *signgamp)
+{
+       long double nadj,p,p1,p2,p3,q,r,t,w,y,z;
+       uint64_t lx;
+       int i;
+       uint16_t hx,ix;
+
+       EXTRACT_LDBL80_WORDS(hx,lx,x);
+
+    /* purge +-Inf and NaNs */
+       *signgamp = 1;
+       ix = hx&0x7fff;
+       if(ix==0x7fff) return x*x;
+
+       ENTERI();
+
+    /* purge +-0 and tiny arguments */
+       *signgamp = 1-2*(hx>>15);
+       if(ix<0x3fff-67) {              /* |x|<2**-(p+3), return -log(|x|) */
+           if((ix|lx)==0)
+               RETURNI(one/vzero);
+           RETURNI(-logl(fabsl(x)));
+       }
+
+    /* purge negative integers and start evaluation for other x < 0 */
+       if(hx&0x8000) {
+           *signgamp = 1;
+           if(ix>=0x3fff+63)           /* |x|>=2**(p-1), must be -integer */
+               RETURNI(one/vzero);
+           t = sin_pil(x);
+           if(t==zero) RETURNI(one/vzero); /* -integer */
+           nadj = logl(pi/fabsl(t*x));
+           if(t<zero) *signgamp = -1;
+           x = -x;
+       }
+
+    /* purge 1 and 2 */
+       if((ix==0x3fff || ix==0x4000) && lx==0x8000000000000000ULL) r = 0;
+    /* for x < 2.0 */
+       else if(ix<0x4000) {
+    /*
+     * XXX Supposedly, one can use the following information to replace the
+     * XXX FP rational expressions.  A similar approach is appropriate
+     * XXX for ld128, but one (may need?) needs to consider llx, too.
+     *
+     * 8.9999961853027344e-01 3ffe e666600000000000
+     * 7.3159980773925781e-01 3ffe bb4a200000000000
+     * 2.3163998126983643e-01 3ffc ed33080000000000
+     * 1.7316312789916992e+00 3fff dda6180000000000
+     * 1.2316322326660156e+00 3fff 9da6200000000000
+     */
+           if(x<8.9999961853027344e-01) {
+               r = -logl(x);
+               if(x>=7.3159980773925781e-01) {y = 1-x; i= 0;}
+               else if(x>=2.3163998126983643e-01) {y= x-(tc-1); i=1;}
+               else {y = x; i=2;}
+           } else {
+               r = 0;
+               if(x>=1.7316312789916992e+00) {y=2-x;i=0;}
+               else if(x>=1.2316322326660156e+00) {y=x-tc;i=1;}
+               else {y=x-1;i=2;}
+           }
+           switch(i) {
+             case 0:
+               z = y*y;
+               p1 = a0+z*(a2+z*(a4+z*(a6+z*(a8+z*(a10+z*a12)))));
+               p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*(a11+z*a13))))));
+               p  = y*p1+p2;
+               r  += p-y/2; break;
+             case 1:
+               p = t0+y*t1+tt+y*y*(t2+y*(t3+y*(t4+y*(t5+y*(t6+y*(t7+y*(t8+
+                   y*(t9+y*(t10+y*(t11+y*(t12+y*(t13+y*(t14+y*(t15+y*(t16+
+                   y*(t17+y*t18))))))))))))))));
+               r += tf + p; break;
+             case 2:
+               p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*(u5+y*u6))))));
+               p2 = 1+y*(v1+y*(v2+y*(v3+y*(v4+y*(v5+y*v6)))));
+               r += p1/p2-y/2;
+           }
+       }
+    /* x < 8.0 */
+       else if(ix<0x4002) {
+           i = x;
+           y = x-i;
+           p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6))))));
+           q = 1+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*(r6+y*r7))))));
+           r = y/2+p/q;
+           z = 1;      /* lgamma(1+s) = log(s) + lgamma(s) */
+           switch(i) {
+           case 7: z *= (y+6);         /* FALLTHRU */
+           case 6: z *= (y+5);         /* FALLTHRU */
+           case 5: z *= (y+4);         /* FALLTHRU */
+           case 4: z *= (y+3);         /* FALLTHRU */
+           case 3: z *= (y+2);         /* FALLTHRU */
+                   r += logl(z); break;
+           }
+    /* 8.0 <= x < 2**(p+3) */
+       } else if (ix<0x3fff+67) {
+           t = logl(x);
+           z = one/x;
+           y = z*z;
+           w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*(w6+y*(w7+y*w8)))))));
+           r = (x-half)*(t-one)+w;
+    /* 2**(p+3) <= x <= inf */
+       } else 
+           r =  x*(logl(x)-1);
+       if(hx&0x8000) r = nadj - r;
+       RETURNI(r);
+}
similarity index 57%
copy from lib/libm/ld80/s_expl.c
copy to lib/libm/ld80/k_expl.h
index 3155ed4..150974c 100644 (file)
@@ -1,3 +1,5 @@
+/* from: FreeBSD: head/lib/msun/ld80/s_expl.c 251343 2013-06-03 19:51:32Z kargl */
+
 /*-
  * Copyright (c) 2009-2013 Steven G. Kargl
  * All rights reserved.
  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
  * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  *
- * Optimized by Bruce D. Evans.
+ * $FreeBSD$
  *
- * $FreeBSD: head/lib/msun/ld80/s_expl.c 251343 2013-06-03 19:51:32Z kargl $
+ * Optimized by Bruce D. Evans.
  */
 
-/**
- * Compute the exponential of x for Intel 80-bit format.  This is based on:
- *
- *   PTP Tang, "Table-driven implementation of the exponential function
- *   in IEEE floating-point arithmetic," ACM Trans. Math. Soft., 15,
- *   144-157 (1989).
+
+/*
+ * See s_expl.c for more comments about __k_expl().
  *
- * where the 32 table entries have been expanded to INTERVALS (see below).
+ * See ../src/e_exp.c and ../src/k_exp.h for precision-independent comments
+ * about the secondary kernels.
  */
 
-#include <float.h>
-
-#ifdef __i386__
-#include <ieeefp.h>
-#endif
-
-#include "fpmath.h"
-#include "math.h"
-#include "math_private.h"
-
 #define        INTERVALS       128
 #define        LOG2_INTERVALS  7
 #define        BIAS    (LDBL_MAX_EXP - 1)
 
-static const long double
-huge = 0x1p10000L,
-twom10000 = 0x1p-10000L;
-/* XXX Prevent gcc from erroneously constant folding this: */
-static volatile const long double tiny = 0x1p-10000L;
-
-static const union IEEEl2bits
-/* log(2**16384 - 0.5) rounded towards zero: */
-/* log(2**16384 - 0.5 + 1) rounded towards zero for expm1l() is the same: */
-o_thresholdu = LD80C(0xb17217f7d1cf79ab, 13,  11356.5234062941439488L),
-#define o_threshold     (o_thresholdu.e)
-/* log(2**(-16381-64-1)) rounded towards zero: */
-u_thresholdu = LD80C(0xb21dfe7f09e2baa9, 13, -11399.4985314888605581L);
-#define u_threshold     (u_thresholdu.e)
-
 static const double
 /*
  * ln2/INTERVALS = L1+L2 (hi+lo decomposition for multiplication).  L1 must
@@ -101,6 +76,12 @@ static const struct {
        double  lo;
 } tbl[INTERVALS] = {
        0x1p+0, 0x0p+0,
+       /*
+        * XXX hi is rounded down, and the formatting is not quite normal.
+        * But I rather like both.  The 0x1.*p format is good for 4N+1
+        * mantissa bits.  Rounding down makes the lo terms positive,
+        * so that the columnar formatting can be simpler.
+        */
        0x1.0163da9fb3335p+0, 0x1.b61299ab8cdb7p-54,
        0x1.02c9a3e778060p+0, 0x1.dcdef95949ef4p-53,
        0x1.04315e86e7f84p+0, 0x1.7ae71f3441b49p-53,
@@ -230,239 +211,95 @@ static const struct {
        0x1.fd3c22b8f71f1p+0, 0x1.2eb74966579e7p-57
 };
 
-long double
-expl(long double x)
+/*
+ * Kernel for expl(x).  x must be finite and not tiny or huge.
+ * "tiny" is anything that would make us underflow (|A6*x^6| < ~LDBL_MIN).
+ * "huge" is anything that would make fn*L1 inexact (|x| > ~2**17*ln2).
+ */
+static inline void
+__k_expl(long double x, long double *hip, long double *lop, int *kp)
 {
-       union IEEEl2bits u, v;
-       long double fn, q, r, r1, r2, t, twopk, twopkp10000;
-       long double z;
-       int k, n, n2;
-       uint16_t hx, ix;
-
-       /* Filter out exceptional cases. */
-       u.e = x;
-       hx = u.xbits.expsign;
-       ix = hx & 0x7fff;
-       if (ix >= BIAS + 13) {          /* |x| >= 8192 or x is NaN */
-               if (ix == BIAS + LDBL_MAX_EXP) {
-                       if (hx & 0x8000)  /* x is -Inf, -NaN or unsupported */
-                               return (-1 / x);
-                       return (x + x); /* x is +Inf, +NaN or unsupported */
-               }
-               if (x > o_threshold)
-                       return (huge * huge);
-               if (x < u_threshold)
-                       return (tiny * tiny);
-       } else if (ix < BIAS - 65) {    /* |x| < 0x1p-65 (includes pseudos) */
-               return (1 + x);         /* 1 with inexact iff x != 0 */
-       }
-
-       ENTERI();
+       long double fn, q, r, r1, r2, t, z;
+       int n, n2;
 
        /* Reduce x to (k*ln2 + endpoint[n2] + r1 + r2). */
        /* Use a specialized rint() to get fn.  Assume round-to-nearest. */
        fn = x * INV_L + 0x1.8p63 - 0x1.8p63;
        r = x - fn * L1 - fn * L2;      /* r = r1 + r2 done independently. */
 #if defined(HAVE_EFFICIENT_IRINTL)
-       n  = irintl(fn);
+       n = irintl(fn);
 #elif defined(HAVE_EFFICIENT_IRINT)
-       n  = irint(fn);
+       n = irint(fn);
 #else
-       n  = (int)fn;
+       n = (int)fn;
 #endif
        n2 = (unsigned)n % INTERVALS;
        /* Depend on the sign bit being propagated: */
-       k = n >> LOG2_INTERVALS;
+       *kp = n >> LOG2_INTERVALS;
        r1 = x - fn * L1;
        r2 = fn * -L2;
 
-       /* Prepare scale factors. */
-       v.e = 1;
-       if (k >= LDBL_MIN_EXP) {
-               v.xbits.expsign = BIAS + k;
-               twopk = v.e;
-       } else {
-               v.xbits.expsign = BIAS + k + 10000;
-               twopkp10000 = v.e;
-       }
-
        /* Evaluate expl(endpoint[n2] + r1 + r2) = tbl[n2] * expl(r1 + r2). */
        z = r * r;
+#if 0
        q = r2 + z * (A2 + r * A3) + z * z * (A4 + r * A5) + z * z * z * A6;
+#else
+       q = r2 + z * A2 + z * r * (A3 + r * A4 + z * (A5 + r * A6));
+#endif
        t = (long double)tbl[n2].lo + tbl[n2].hi;
-       t = tbl[n2].lo + t * (q + r1) + tbl[n2].hi;
+       *hip = tbl[n2].hi;
+       *lop = tbl[n2].lo + t * (q + r1);
+}
+
+static inline void
+k_hexpl(long double x, long double *hip, long double *lop)
+{
+       float twopkm1;
+       int k;
 
-       /* Scale by 2**k. */
-       if (k >= LDBL_MIN_EXP) {
-               if (k == LDBL_MAX_EXP)
-                       RETURNI(t * 2 * 0x1p16383L);
-               RETURNI(t * twopk);
-       } else {
-               RETURNI(t * twopkp10000 * twom10000);
-       }
+       __k_expl(x, hip, lop, &k);
+       SET_FLOAT_WORD(twopkm1, 0x3f800000 + ((k - 1) << 23));
+       *hip *= twopkm1;
+       *lop *= twopkm1;
 }
 
-/**
- * Compute expm1l(x) for Intel 80-bit format.  This is based on:
- *
- *   PTP Tang, "Table-driven implementation of the Expm1 function
- *   in IEEE floating-point arithmetic," ACM Trans. Math. Soft., 18,
- *   211-222 (1992).
- */
+static inline long double
+hexpl(long double x)
+{
+       long double hi, lo, twopkm2;
+       int k;
 
-/*
- * Our T1 and T2 are chosen to be approximately the points where method
- * A and method B have the same accuracy.  Tang's T1 and T2 are the
- * points where method A's accuracy changes by a full bit.  For Tang,
- * this drop in accuracy makes method A immediately less accurate than
- * method B, but our larger INTERVALS makes method A 2 bits more
- * accurate so it remains the most accurate method significantly
- * closer to the origin despite losing the full bit in our extended
- * range for it.
- */
-static const double
-T1 = -0.1659,                          /* ~-30.625/128 * log(2) */
-T2 =  0.1659;                          /* ~30.625/128 * log(2) */
+       twopkm2 = 1;
+       __k_expl(x, &hi, &lo, &k);
+       SET_LDBL_EXPSIGN(twopkm2, BIAS + k - 2);
+       return (lo + hi) * 2 * twopkm2;
+}
 
+#ifdef _COMPLEX_H
 /*
- * Domain [-0.1659, 0.1659], range ~[-1.2027e-22, 3.4417e-22]:
- * |(exp(x)-1-x-x**2/2)/x - p(x)| < 2**-71.2
+ * See ../src/k_exp.c for details.
  */
-static const union IEEEl2bits
-B3 = LD80C(0xaaaaaaaaaaaaaaab, -3,  1.66666666666666666671e-1L),
-B4 = LD80C(0xaaaaaaaaaaaaaaac, -5,  4.16666666666666666712e-2L);
-
-static const double
-B5  =  8.3333333333333245e-3,          /*  0x1.111111111110cp-7 */
-B6  =  1.3888888888888861e-3,          /*  0x1.6c16c16c16c0ap-10 */
-B7  =  1.9841269841532042e-4,          /*  0x1.a01a01a0319f9p-13 */
-B8  =  2.4801587302069236e-5,          /*  0x1.a01a01a03cbbcp-16 */
-B9  =  2.7557316558468562e-6,          /*  0x1.71de37fd33d67p-19 */
-B10 =  2.7557315829785151e-7,          /*  0x1.27e4f91418144p-22 */
-B11 =  2.5063168199779829e-8,          /*  0x1.ae94fabdc6b27p-26 */
-B12 =  2.0887164654459567e-9;          /*  0x1.1f122d6413fe1p-29 */
-
-long double
-expm1l(long double x)
+static inline long double complex
+__ldexp_cexpl(long double complex z, int expt)
 {
-       union IEEEl2bits u, v;
-       long double fn, hx2_hi, hx2_lo, q, r, r1, r2, t, twomk, twopk, x_hi;
-       long double x_lo, x2, z;
-       long double x4;
-       int k, n, n2;
-       uint16_t hx, ix;
-
-       /* Filter out exceptional cases. */
-       u.e = x;
-       hx = u.xbits.expsign;
-       ix = hx & 0x7fff;
-       if (ix >= BIAS + 6) {           /* |x| >= 64 or x is NaN */
-               if (ix == BIAS + LDBL_MAX_EXP) {
-                       if (hx & 0x8000)  /* x is -Inf, -NaN or unsupported */
-                               return (-1 / x - 1);
-                       return (x + x); /* x is +Inf, +NaN or unsupported */
-               }
-               if (x > o_threshold)
-                       return (huge * huge);
-               /*
-                * expm1l() never underflows, but it must avoid
-                * unrepresentable large negative exponents.  We used a
-                * much smaller threshold for large |x| above than in
-                * expl() so as to handle not so large negative exponents
-                * in the same way as large ones here.
-                */
-               if (hx & 0x8000)        /* x <= -64 */
-                       return (tiny - 1);      /* good for x < -65ln2 - eps */
-       }
-
-       ENTERI();
-
-       if (T1 < x && x < T2) {
-               if (ix < BIAS - 64) {   /* |x| < 0x1p-64 (includes pseudos) */
-                       /* x (rounded) with inexact if x != 0: */
-                       RETURNI(x == 0 ? x :
-                           (0x1p100 * x + fabsl(x)) * 0x1p-100);
-               }
+       long double exp_x, hi, lo;
+       long double x, y, scale1, scale2;
+       int half_expt, k;
 
-               x2 = x * x;
-               x4 = x2 * x2;
-               q = x4 * (x2 * (x4 *
-                   /*
-                    * XXX the number of terms is no longer good for
-                    * pairwise grouping of all except B3, and the
-                    * grouping is no longer from highest down.
-                    */
-                   (x2 *            B12  + (x * B11 + B10)) +
-                   (x2 * (x * B9 +  B8) +  (x * B7 +  B6))) +
-                         (x * B5 +  B4.e)) + x2 * x * B3.e;
-
-               x_hi = (float)x;
-               x_lo = x - x_hi;
-               hx2_hi = x_hi * x_hi / 2;
-               hx2_lo = x_lo * (x + x_hi) / 2;
-               if (ix >= BIAS - 7)
-                       RETURNI(hx2_lo + x_lo + q + (hx2_hi + x_hi));
-               else
-                       RETURNI(hx2_lo + q + hx2_hi + x);
-       }
-
-       /* Reduce x to (k*ln2 + endpoint[n2] + r1 + r2). */
-       /* Use a specialized rint() to get fn.  Assume round-to-nearest. */
-       fn = x * INV_L + 0x1.8p63 - 0x1.8p63;
-#if defined(HAVE_EFFICIENT_IRINTL)
-       n = irintl(fn);
-#elif defined(HAVE_EFFICIENT_IRINT)
-       n = irint(fn);
-#else
-       n = (int)fn;
-#endif
-       n2 = (unsigned)n % INTERVALS;
-       k = n >> LOG2_INTERVALS;
-       r1 = x - fn * L1;
-       r2 = fn * -L2;
-       r = r1 + r2;
-
-       /* Prepare scale factor. */
-       v.e = 1;
-       v.xbits.expsign = BIAS + k;
-       twopk = v.e;
-
-       /*
-        * Evaluate lower terms of
-        * expl(endpoint[n2] + r1 + r2) = tbl[n2] * expl(r1 + r2).
-        */
-       z = r * r;
-       q = r2 + z * (A2 + r * A3) + z * z * (A4 + r * A5) + z * z * z * A6;
-
-       t = (long double)tbl[n2].lo + tbl[n2].hi;
+       x = creall(z);
+       y = cimagl(z);
+       __k_expl(x, &hi, &lo, &k);
 
-       if (k == 0) {
-               t = tbl[n2].lo * (r1 + 1) + t * q + tbl[n2].hi * r1 +
-                   (tbl[n2].hi - 1);
-               RETURNI(t);
-       }
-       if (k == -1) {
-               t = tbl[n2].lo * (r1 + 1) + t * q + tbl[n2].hi * r1 + 
-                   (tbl[n2].hi - 2);
-               RETURNI(t / 2);
-       }
-       if (k < -7) {
-               t = tbl[n2].lo + t * (q + r1) + tbl[n2].hi;
-               RETURNI(t * twopk - 1);
-       }
-       if (k > 2 * LDBL_MANT_DIG - 1) {
-               t = tbl[n2].lo + t * (q + r1) + tbl[n2].hi;
-               if (k == LDBL_MAX_EXP)
-                       RETURNI(t * 2 * 0x1p16383L - 1);
-               RETURNI(t * twopk - 1);
-       }
+       exp_x = (lo + hi) * 0x1p16382;
+       expt += k - 16382;
 
-       v.xbits.expsign = BIAS - k;
-       twomk = v.e;
+       scale1 = 1;
+       half_expt = expt / 2;
+       SET_LDBL_EXPSIGN(scale1, BIAS + half_expt);
+       scale2 = 1;
+       SET_LDBL_EXPSIGN(scale1, BIAS + expt - half_expt);
 
-       if (k > LDBL_MANT_DIG - 1)
-               t = tbl[n2].lo - twomk + t * (q + r1) + tbl[n2].hi;
-       else
-               t = tbl[n2].lo + t * (q + r1) + (tbl[n2].hi - twomk);
-       RETURNI(t * twopk);
+       return (cpackl(cos(y) * exp_x * scale1 * scale2,
+           sinl(y) * exp_x * scale1 * scale2));
 }
+#endif /* _COMPLEX_H */
diff --git a/lib/libm/ld80/s_erfl.c b/lib/libm/ld80/s_erfl.c
new file mode 100644 (file)
index 0000000..ececc7e
--- /dev/null
@@ -0,0 +1,335 @@
+/* @(#)s_erf.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+
+/*
+ * See s_erf.c for complete comments.
+ *
+ * Converted to long double by Steven G. Kargl.
+ */
+#include <float.h>
+#ifdef __i386__
+#include <ieeefp.h>
+#endif
+
+#include "fpmath.h"
+#include "math.h"
+#include "math_private.h"
+
+/* XXX Prevent compilers from erroneously constant folding: */
+static const volatile long double tiny = 0x1p-10000L;
+
+static const double
+half= 0.5,
+one = 1,
+two = 2;
+/*
+ * In the domain [0, 2**-34], only the first term in the power series
+ * expansion of erf(x) is used.  The magnitude of the first neglected
+ * terms is less than 2**-102.
+ */
+static const union IEEEl2bits
+efxu  = LD80C(0x8375d410a6db446c, -3,  1.28379167095512573902e-1L),
+efx8u = LD80C(0x8375d410a6db446c,  0,  1.02703333676410059122e+0L),
+/*
+ * Domain [0, 0.84375], range ~[-1.423e-22, 1.423e-22]:
+ * |(erf(x) - x)/x - pp(x)/qq(x)| < 2**-72.573
+ */
+pp0u  = LD80C(0x8375d410a6db446c, -3,   1.28379167095512573902e-1L),
+pp1u  = LD80C(0xa46c7d09ec3d0cec, -2,  -3.21140201054840180596e-1L),
+pp2u  = LD80C(0x9b31e66325576f86, -5,  -3.78893851760347812082e-2L),
+pp3u  = LD80C(0x804ac72c9a0b97dd, -7,  -7.83032847030604679616e-3L),
+pp4u  = LD80C(0x9f42bcbc3d5a601d, -12, -3.03765663857082048459e-4L),
+pp5u  = LD80C(0x9ec4ad6193470693, -16, -1.89266527398167917502e-5L),
+qq1u  = LD80C(0xdb4b8eb713188d6b, -2,   4.28310832832310510579e-1L),
+qq2u  = LD80C(0xa5750835b2459bd1, -4,   8.07896272074540216658e-2L),
+qq3u  = LD80C(0x8b85d6bd6a90b51c, -7,   8.51579638189385354266e-3L),
+qq4u  = LD80C(0x87332f82cff4ff96, -11,  5.15746855583604912827e-4L),
+qq5u  = LD80C(0x83466cb6bf9dca00, -16,  1.56492109706256700009e-5L),
+qq6u  = LD80C(0xf5bf98c2f996bf63, -24,  1.14435527803073879724e-7L);
+#define        efx     (efxu.e)
+#define        efx8    (efx8u.e)
+#define        pp0     (pp0u.e)
+#define        pp1     (pp1u.e)
+#define        pp2     (pp2u.e)
+#define        pp3     (pp3u.e)
+#define        pp4     (pp4u.e)
+#define        pp5     (pp5u.e)
+#define        qq1     (qq1u.e)
+#define        qq2     (qq2u.e)
+#define        qq3     (qq3u.e)
+#define        qq4     (qq4u.e)
+#define        qq5     (qq5u.e)
+#define        qq6     (qq6u.e)
+static const union IEEEl2bits
+erxu  = LD80C(0xd7bb3d0000000000, -1,  8.42700779438018798828e-1L),
+/*
+ * Domain [0.84375, 1.25], range ~[-8.132e-22, 8.113e-22]:
+ * |(erf(x) - erx) - pa(x)/qa(x)| < 2**-71.762
+ */
+pa0u  = LD80C(0xe8211158da02c692, -27,  1.35116960705131296711e-8L),
+pa1u  = LD80C(0xd488f89f36988618, -2,   4.15107507167065612570e-1L),
+pa2u  = LD80C(0xece74f8c63fa3942, -4,  -1.15675565215949226989e-1L),
+pa3u  = LD80C(0xc8d31e020727c006, -4,   9.80589241379624665791e-2L),
+pa4u  = LD80C(0x985d5d5fafb0551f, -5,   3.71984145558422368847e-2L),
+pa5u  = LD80C(0xa5b6c4854d2f5452, -8,  -5.05718799340957673661e-3L),
+pa6u  = LD80C(0x85c8d58fe3993a47, -8,   4.08277919612202243721e-3L),
+pa7u  = LD80C(0xddbfbc23677b35cf, -13,  2.11476292145347530794e-4L),
+qa1u  = LD80C(0xb8a977896f5eff3f, -1,   7.21335860303380361298e-1L),
+qa2u  = LD80C(0x9fcd662c3d4eac86, -1,   6.24227891731886593333e-1L),
+qa3u  = LD80C(0x9d0b618eac67ba07, -2,   3.06727455774491855801e-1L),
+qa4u  = LD80C(0x881a4293f6d6c92d, -3,   1.32912674218195890535e-1L),
+qa5u  = LD80C(0xbab144f07dea45bf, -5,   4.55792134233613027584e-2L),
+qa6u  = LD80C(0xa6c34ba438bdc900, -7,   1.01783980070527682680e-2L),
+qa7u  = LD80C(0x8fa866dc20717a91, -9,   2.19204436518951438183e-3L);
+#define erx    (erxu.e)
+#define pa0    (pa0u.e)
+#define pa1    (pa1u.e)
+#define pa2    (pa2u.e)
+#define pa3    (pa3u.e)
+#define pa4    (pa4u.e)
+#define pa5    (pa5u.e)
+#define pa6    (pa6u.e)
+#define pa7    (pa7u.e)
+#define qa1    (qa1u.e)
+#define qa2    (qa2u.e)
+#define qa3    (qa3u.e)
+#define qa4    (qa4u.e)
+#define qa5    (qa5u.e)
+#define qa6    (qa6u.e)
+#define qa7    (qa7u.e)
+static const union IEEEl2bits
+/*
+ * Domain [1.25,2.85715], range ~[-2.334e-22,2.334e-22]:
+ * |log(x*erfc(x)) + x**2 + 0.5625 - ra(x)/sa(x)| < 2**-71.860
+ */
+ra0u  = LD80C(0xa1a091e0fb4f335a, -7, -9.86494298915814308249e-3L),
+ra1u  = LD80C(0xc2b0d045ae37df6b, -1, -7.60510460864878271275e-1L),
+ra2u  = LD80C(0xf2cec3ee7da636c5, 3,  -1.51754798236892278250e+1L),
+ra3u  = LD80C(0x813cc205395adc7d, 7,  -1.29237335516455333420e+2L),
+ra4u  = LD80C(0x8737c8b7b4062c2f, 9,  -5.40871625829510494776e+2L),
+ra5u  = LD80C(0x8ffe5383c08d4943, 10, -1.15194769466026108551e+3L),
+ra6u  = LD80C(0x983573e64d5015a9, 10, -1.21767039790249025544e+3L),
+ra7u  = LD80C(0x92a794e763a6d4db, 9,  -5.86618463370624636688e+2L),
+ra8u  = LD80C(0xd5ad1fae77c3d9a3, 6,  -1.06838132335777049840e+2L),
+ra9u  = LD80C(0x934c1a247807bb9c, 2,  -4.60303980944467334806e+0L),
+sa1u  = LD80C(0xd342f90012bb1189, 4,   2.64077014928547064865e+1L),
+sa2u  = LD80C(0x839be13d9d5da883, 8,   2.63217811300123973067e+2L),
+sa3u  = LD80C(0x9f8cba6d1ae1b24b, 10,  1.27639775710344617587e+3L),
+sa4u  = LD80C(0xcaa83f403713e33e, 11,  3.24251544209971162003e+3L),
+sa5u  = LD80C(0x8796aff2f3c47968, 12,  4.33883591261332837874e+3L),
+sa6u  = LD80C(0xb6ef97f9c753157b, 11,  2.92697460344182158454e+3L),
+sa7u  = LD80C(0xe02aee5f83773d1c, 9,   8.96670799139389559818e+2L),
+sa8u  = LD80C(0xc82b83855b88e07e, 6,   1.00084987800048510018e+2L),
+sa9u  = LD80C(0x92f030aefadf28ad, 1,   2.29591004455459083843e+0L);
+#define ra0    (ra0u.e)
+#define ra1    (ra1u.e)
+#define ra2    (ra2u.e)
+#define ra3    (ra3u.e)
+#define ra4    (ra4u.e)
+#define ra5    (ra5u.e)
+#define ra6    (ra6u.e)
+#define ra7    (ra7u.e)
+#define ra8    (ra8u.e)
+#define ra9    (ra9u.e)
+#define sa1    (sa1u.e)
+#define sa2    (sa2u.e)
+#define sa3    (sa3u.e)
+#define sa4    (sa4u.e)
+#define sa5    (sa5u.e)
+#define sa6    (sa6u.e)
+#define sa7    (sa7u.e)
+#define sa8    (sa8u.e)
+#define sa9    (sa9u.e)
+/*
+ * Domain [2.85715,7], range ~[-8.323e-22,8.390e-22]:
+ * |log(x*erfc(x)) + x**2 + 0.5625 - rb(x)/sb(x)| < 2**-70.326
+ */
+static const union IEEEl2bits
+rb0u = LD80C(0xa1a091cf43abcd26, -7, -9.86494292470284646962e-3L),
+rb1u = LD80C(0xd19d2df1cbb8da0a, -1, -8.18804618389296662837e-1L),
+rb2u = LD80C(0x9a4dd1383e5daf5b, 4,  -1.92879967111618594779e+1L),
+rb3u = LD80C(0xbff0ae9fc0751de6, 7,  -1.91940164551245394969e+2L),
+rb4u = LD80C(0xdde08465310b472b, 9,  -8.87508080766577324539e+2L),
+rb5u = LD80C(0xe796e1d38c8c70a9, 10, -1.85271506669474503781e+3L),
+rb6u = LD80C(0xbaf655a76e0ab3b5, 10, -1.49569795581333675349e+3L),
+rb7u = LD80C(0x95d21e3e75503c21, 8,  -2.99641547972948019157e+2L),
+sb1u = LD80C(0x814487ed823c8cbd, 5,   3.23169247732868256569e+1L),
+sb2u = LD80C(0xbe4bfbb1301304be, 8,   3.80593618534539961773e+2L),
+sb3u = LD80C(0x809c4ade46b927c7, 11,  2.05776827838541292848e+3L),
+sb4u = LD80C(0xa55284359f3395a8, 12,  5.29031455540062116327e+3L),
+sb5u = LD80C(0xbcfa72da9b820874, 12,  6.04730608102312640462e+3L),
+sb6u = LD80C(0x9d09a35988934631, 11,  2.51260238030767176221e+3L),
+sb7u = LD80C(0xd675bbe542c159fa, 7,   2.14459898308561015684e+2L);
+#define rb0    (rb0u.e)
+#define rb1    (rb1u.e)
+#define rb2    (rb2u.e)
+#define rb3    (rb3u.e)
+#define rb4    (rb4u.e)
+#define rb5    (rb5u.e)
+#define rb6    (rb6u.e)
+#define rb7    (rb7u.e)
+#define sb1    (sb1u.e)
+#define sb2    (sb2u.e)
+#define sb3    (sb3u.e)
+#define sb4    (sb4u.e)
+#define sb5    (sb5u.e)
+#define sb6    (sb6u.e)
+#define sb7    (sb7u.e)
+/*
+ * Domain [7,108], range ~[-4.422e-22,4.422e-22]:
+ * |log(x*erfc(x)) + x**2 + 0.5625 - rc(x)/sc(x)| < 2**-70.938
+ */
+static const union IEEEl2bits
+/* err = -4.422092275318925082e-22 -70.937689 */
+rc0u = LD80C(0xa1a091cf437a17ad, -7, -9.86494292470008707260e-3L),
+rc1u = LD80C(0xbe79c5a978122b00, -1, -7.44045595049165939261e-1L),
+rc2u = LD80C(0xdb26f9bbe31a2794, 3,  -1.36970155085888424425e+1L),
+rc3u = LD80C(0xb5f69a38f5747ac8, 6,  -9.09816453742625888546e+1L),
+rc4u = LD80C(0xd79676d970d0a21a, 7,  -2.15587750997584074147e+2L),
+rc5u = LD80C(0xfe528153c45ec97c, 6,  -1.27161142938347796666e+2L),
+sc1u = LD80C(0xc5e8cd46d5604a96, 4,   2.47386727842204312937e+1L),
+sc2u = LD80C(0xc5f0f5a5484520eb, 7,   1.97941248254913378865e+2L),
+sc3u = LD80C(0x964e3c7b34db9170, 9,   6.01222441484087787522e+2L),
+sc4u = LD80C(0x99be1b89faa0596a, 9,   6.14970430845978077827e+2L),
+sc5u = LD80C(0xf80dfcbf37ffc5ea, 6,   1.24027318931184605891e+2L);
+#define rc0    (rc0u.e)
+#define rc1    (rc1u.e)
+#define rc2    (rc2u.e)
+#define rc3    (rc3u.e)
+#define rc4    (rc4u.e)
+#define rc5    (rc5u.e)
+#define sc1    (sc1u.e)
+#define sc2    (sc2u.e)
+#define sc3    (sc3u.e)
+#define sc4    (sc4u.e)
+#define sc5    (sc5u.e)
+
+long double
+erfl(long double x)
+{
+       long double ax,R,S,P,Q,s,y,z,r;
+       uint64_t lx;
+       int32_t i;
+       uint16_t hx;
+
+       EXTRACT_LDBL80_WORDS(hx, lx, x);
+
+       if((hx & 0x7fff) == 0x7fff) {   /* erfl(nan)=nan */
+               i = (hx>>15)<<1;
+               return (1-i)+one/x;     /* erfl(+-inf)=+-1 */
+       }
+
+       ENTERI();
+
+       ax = fabsl(x);
+       if(ax < 0.84375) {
+           if(ax < 0x1p-34L) {
+               if(ax < 0x1p-16373L)    
+                   RETURNI((8*x+efx8*x)/8);    /* avoid spurious underflow */
+               RETURNI(x + efx*x);
+           }
+           z = x*x;
+           r = pp0+z*(pp1+z*(pp2+z*(pp3+z*(pp4+z*pp5))));
+           s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*(qq5+z*qq6)))));
+           y = r/s;
+           RETURNI(x + x*y);
+       }
+       if(ax < 1.25) {
+           s = ax-one;
+           P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*(pa6+s*pa7))))));
+           Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*(qa6+s*qa7))))));
+           if(x>=0) RETURNI(erx + P/Q); else RETURNI(-erx - P/Q);
+       }
+       if(ax >= 7) {                   /* inf>|x|>= 7 */
+           if(x>=0) RETURNI(one-tiny); else RETURNI(tiny-one);
+       }
+       s = one/(ax*ax);
+       if(ax < 2.85715) {      /* |x| < 2.85715 */
+           R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(ra5+s*(ra6+s*(ra7+
+               s*(ra8+s*ra9))))))));
+           S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(sa5+s*(sa6+s*(sa7+
+               s*(sa8+s*sa9))))))));
+       } else {        /* |x| >= 2.85715 */
+           R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(rb5+s*(rb6+s*rb7))))));
+           S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(sb5+s*(sb6+s*sb7))))));
+       }
+       z=(float)ax;
+       r=expl(-z*z-0.5625)*expl((z-ax)*(z+ax)+R/S);
+       if(x>=0) RETURNI(one-r/ax); else RETURNI(r/ax-one);
+}
+
+long double
+erfcl(long double x)
+{
+       long double ax,R,S,P,Q,s,y,z,r;
+       uint64_t lx;
+       uint16_t hx;
+
+       EXTRACT_LDBL80_WORDS(hx, lx, x);
+
+       if((hx & 0x7fff) == 0x7fff) {   /* erfcl(nan)=nan */
+                                       /* erfcl(+-inf)=0,2 */
+           return ((hx>>15)<<1)+one/x;
+       }
+
+       ENTERI();
+
+       ax = fabsl(x);
+       if(ax < 0.84375L) {
+           if(ax < 0x1p-34L)
+               RETURNI(one-x);
+           z = x*x;
+           r = pp0+z*(pp1+z*(pp2+z*(pp3+z*(pp4+z*pp5))));
+           s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*(qq5+z*qq6)))));
+           y = r/s;
+           if(ax < 0.25L) {    /* x<1/4 */
+               RETURNI(one-(x+x*y));
+           } else {
+               r = x*y;
+               r += (x-half);
+              RETURNI(half - r);
+           }
+       }
+       if(ax < 1.25L) {
+           s = ax-one;
+           P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*(pa6+s*pa7))))));
+           Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*(qa6+s*qa7))))));
+           if(x>=0) {
+               z  = one-erx; RETURNI(z - P/Q);
+           } else {
+               z = (erx+P/Q); RETURNI(one+z);
+           }
+       }
+
+       if(ax < 108) {                  /* |x| < 108 */
+           s = one/(ax*ax);
+           if(ax < 2.85715) {          /* |x| < 2.85715 */
+               R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(ra5+s*(ra6+s*(ra7+
+                   s*(ra8+s*ra9))))))));
+               S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(sa5+s*(sa6+s*(sa7+
+                   s*(sa8+s*sa9))))))));
+           } else if(ax < 7) {         /* | |x| < 7 */
+               R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(rb5+s*(rb6+s*rb7))))));
+               S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(sb5+s*(sb6+s*sb7))))));
+           } else {
+               if(x < -7) RETURNI(two-tiny);/* x < -7 */
+               R=rc0+s*(rc1+s*(rc2+s*(rc3+s*(rc4+s*rc5))));
+               S=one+s*(sc1+s*(sc2+s*(sc3+s*(sc4+s*sc5))));
+           }
+           z = (float)ax;
+           r = expl(-z*z-0.5625)*expl((z-ax)*(z+ax)+R/S);
+           if(x>0) RETURNI(r/ax); else RETURNI(two-r/ax);
+       } else {
+           if(x>0) RETURNI(tiny*tiny); else RETURNI(two-tiny);
+       }
+}
index 3155ed4..824de42 100644 (file)
  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
  * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  *
- * Optimized by Bruce D. Evans.
+ * $FreeBSD: head/lib/msun/ld80/s_expl.c 260066 2013-12-30 00:51:25Z kargl $
  *
- * $FreeBSD: head/lib/msun/ld80/s_expl.c 251343 2013-06-03 19:51:32Z kargl $
+ * Optimized by Bruce D. Evans.
  */
 
+
 /**
  * Compute the exponential of x for Intel 80-bit format.  This is based on:
  *
 #include "fpmath.h"
 #include "math.h"
 #include "math_private.h"
+#include "k_expl.h"
 
-#define        INTERVALS       128
-#define        LOG2_INTERVALS  7
-#define        BIAS    (LDBL_MAX_EXP - 1)
+/* XXX Prevent compilers from erroneously constant folding these: */
+static const volatile long double
+huge = 0x1p10000L,
+tiny = 0x1p-10000L;
 
 static const long double
-huge = 0x1p10000L,
 twom10000 = 0x1p-10000L;
-/* XXX Prevent gcc from erroneously constant folding this: */
-static volatile const long double tiny = 0x1p-10000L;
 
 static const union IEEEl2bits
 /* log(2**16384 - 0.5) rounded towards zero: */
@@ -67,178 +67,16 @@ o_thresholdu = LD80C(0xb17217f7d1cf79ab, 13,  11356.5234062941439488L),
 u_thresholdu = LD80C(0xb21dfe7f09e2baa9, 13, -11399.4985314888605581L);
 #define u_threshold     (u_thresholdu.e)
 
-static const double
-/*
- * ln2/INTERVALS = L1+L2 (hi+lo decomposition for multiplication).  L1 must
- * have at least 22 (= log2(|LDBL_MIN_EXP-extras|) + log2(INTERVALS)) lowest
- * bits zero so that multiplication of it by n is exact.
- */
-INV_L = 1.8466496523378731e+2,         /*  0x171547652b82fe.0p-45 */
-L1 =  5.4152123484527692e-3,           /*  0x162e42ff000000.0p-60 */
-L2 = -3.2819649005320973e-13,          /* -0x1718432a1b0e26.0p-94 */
-/*
- * Domain [-0.002708, 0.002708], range ~[-5.7136e-24, 5.7110e-24]:
- * |exp(x) - p(x)| < 2**-77.2
- * (0.002708 is ln2/(2*INTERVALS) rounded up a little).
- */
-A2 =  0.5,
-A3 =  1.6666666666666119e-1,           /*  0x15555555555490.0p-55 */
-A4 =  4.1666666666665887e-2,           /*  0x155555555554e5.0p-57 */
-A5 =  8.3333354987869413e-3,           /*  0x1111115b789919.0p-59 */
-A6 =  1.3888891738560272e-3;           /*  0x16c16c651633ae.0p-62 */
-
-/*
- * 2^(i/INTERVALS) for i in [0,INTERVALS] is represented by two values where
- * the first 53 bits of the significand are stored in hi and the next 53
- * bits are in lo.  Tang's paper states that the trailing 6 bits of hi must
- * be zero for his algorithm in both single and double precision, because
- * the table is re-used in the implementation of expm1() where a floating
- * point addition involving hi must be exact.  Here hi is double, so
- * converting it to long double gives 11 trailing zero bits.
- */
-static const struct {
-       double  hi;
-       double  lo;
-} tbl[INTERVALS] = {
-       0x1p+0, 0x0p+0,
-       0x1.0163da9fb3335p+0, 0x1.b61299ab8cdb7p-54,
-       0x1.02c9a3e778060p+0, 0x1.dcdef95949ef4p-53,
-       0x1.04315e86e7f84p+0, 0x1.7ae71f3441b49p-53,
-       0x1.059b0d3158574p+0, 0x1.d73e2a475b465p-55,
-       0x1.0706b29ddf6ddp+0, 0x1.8db880753b0f6p-53,
-       0x1.0874518759bc8p+0, 0x1.186be4bb284ffp-57,
-       0x1.09e3ecac6f383p+0, 0x1.1487818316136p-54,
-       0x1.0b5586cf9890fp+0, 0x1.8a62e4adc610bp-54,
-       0x1.0cc922b7247f7p+0, 0x1.01edc16e24f71p-54,
-       0x1.0e3ec32d3d1a2p+0, 0x1.03a1727c57b53p-59,
-       0x1.0fb66affed31ap+0, 0x1.e464123bb1428p-53,
-       0x1.11301d0125b50p+0, 0x1.49d77e35db263p-53,
-       0x1.12abdc06c31cbp+0, 0x1.f72575a649ad2p-53,
-       0x1.1429aaea92ddfp+0, 0x1.66820328764b1p-53,
-       0x1.15a98c8a58e51p+0, 0x1.2406ab9eeab0ap-55,
-       0x1.172b83c7d517ap+0, 0x1.b9bef918a1d63p-53,
-       0x1.18af9388c8de9p+0, 0x1.777ee1734784ap-53,
-       0x1.1a35beb6fcb75p+0, 0x1.e5b4c7b4968e4p-55,
-       0x1.1bbe084045cd3p+0, 0x1.3563ce56884fcp-53,
-       0x1.1d4873168b9aap+0, 0x1.e016e00a2643cp-54,
-       0x1.1ed5022fcd91cp+0, 0x1.71033fec2243ap-53,
-       0x1.2063b88628cd6p+0, 0x1.dc775814a8495p-55,
-       0x1.21f49917ddc96p+0, 0x1.2a97e9494a5eep-55,
-       0x1.2387a6e756238p+0, 0x1.9b07eb6c70573p-54,
-       0x1.251ce4fb2a63fp+0, 0x1.ac155bef4f4a4p-55,
-       0x1.26b4565e27cddp+0, 0x1.2bd339940e9d9p-55,
-       0x1.284dfe1f56380p+0, 0x1.2d9e2b9e07941p-53,
-       0x1.29e9df51fdee1p+0, 0x1.612e8afad1255p-55,
-       0x1.2b87fd0dad98fp+0, 0x1.fbbd48ca71f95p-53,
-       0x1.2d285a6e4030bp+0, 0x1.0024754db41d5p-54,
-       0x1.2ecafa93e2f56p+0, 0x1.1ca0f45d52383p-56,
-       0x1.306fe0a31b715p+0, 0x1.6f46ad23182e4p-55,
-       0x1.32170fc4cd831p+0, 0x1.a9ce78e18047cp-55,
-       0x1.33c08b26416ffp+0, 0x1.32721843659a6p-54,
-       0x1.356c55f929ff0p+0, 0x1.928c468ec6e76p-53,
-       0x1.371a7373aa9cap+0, 0x1.4e28aa05e8a8fp-53,
-       0x1.38cae6d05d865p+0, 0x1.0b53961b37da2p-53,
-       0x1.3a7db34e59ff6p+0, 0x1.d43792533c144p-53,
-       0x1.3c32dc313a8e4p+0, 0x1.08003e4516b1ep-53,
-       0x1.3dea64c123422p+0, 0x1.ada0911f09ebcp-55,
-       0x1.3fa4504ac801bp+0, 0x1.417ee03548306p-53,
-       0x1.4160a21f72e29p+0, 0x1.f0864b71e7b6cp-53,
-       0x1.431f5d950a896p+0, 0x1.b8e088728219ap-53,
-       0x1.44e086061892dp+0, 0x1.89b7a04ef80d0p-59,
-       0x1.46a41ed1d0057p+0, 0x1.c944bd1648a76p-54,
-       0x1.486a2b5c13cd0p+0, 0x1.3c1a3b69062f0p-56,
-       0x1.4a32af0d7d3dep+0, 0x1.9cb62f3d1be56p-54,
-       0x1.4bfdad5362a27p+0, 0x1.d4397afec42e2p-56,
-       0x1.4dcb299fddd0dp+0, 0x1.8ecdbbc6a7833p-54,
-       0x1.4f9b2769d2ca6p+0, 0x1.5a67b16d3540ep-53,
-       0x1.516daa2cf6641p+0, 0x1.8225ea5909b04p-53,
-       0x1.5342b569d4f81p+0, 0x1.be1507893b0d5p-53,
-       0x1.551a4ca5d920ep+0, 0x1.8a5d8c4048699p-53,
-       0x1.56f4736b527dap+0, 0x1.9bb2c011d93adp-54,
-       0x1.58d12d497c7fdp+0, 0x1.295e15b9a1de8p-55,
-       0x1.5ab07dd485429p+0, 0x1.6324c054647adp-54,
-       0x1.5c9268a5946b7p+0, 0x1.c4b1b816986a2p-60,
-       0x1.5e76f15ad2148p+0, 0x1.ba6f93080e65ep-54,
-       0x1.605e1b976dc08p+0, 0x1.60edeb25490dcp-53,
-       0x1.6247eb03a5584p+0, 0x1.63e1f40dfa5b5p-53,
-       0x1.6434634ccc31fp+0, 0x1.8edf0e2989db3p-53,
-       0x1.6623882552224p+0, 0x1.224fb3c5371e6p-53,
-       0x1.68155d44ca973p+0, 0x1.038ae44f73e65p-57,
-       0x1.6a09e667f3bccp+0, 0x1.21165f626cdd5p-53,
-       0x1.6c012750bdabep+0, 0x1.daed533001e9ep-53,
-       0x1.6dfb23c651a2ep+0, 0x1.e441c597c3775p-53,
-       0x1.6ff7df9519483p+0, 0x1.9f0fc369e7c42p-53,
-       0x1.71f75e8ec5f73p+0, 0x1.ba46e1e5de15ap-53,
-       0x1.73f9a48a58173p+0, 0x1.7ab9349cd1562p-53,
-       0x1.75feb564267c8p+0, 0x1.7edd354674916p-53,
-       0x1.780694fde5d3fp+0, 0x1.866b80a02162dp-54,
-       0x1.7a11473eb0186p+0, 0x1.afaa2047ed9b4p-53,
-       0x1.7c1ed0130c132p+0, 0x1.f124cd1164dd6p-54,
-       0x1.7e2f336cf4e62p+0, 0x1.05d02ba15797ep-56,
-       0x1.80427543e1a11p+0, 0x1.6c1bccec9346bp-53,
-       0x1.82589994cce12p+0, 0x1.159f115f56694p-53,
-       0x1.8471a4623c7acp+0, 0x1.9ca5ed72f8c81p-53,
-       0x1.868d99b4492ecp+0, 0x1.01c83b21584a3p-53,
-       0x1.88ac7d98a6699p+0, 0x1.994c2f37cb53ap-54,
-       0x1.8ace5422aa0dbp+0, 0x1.6e9f156864b27p-54,
-       0x1.8cf3216b5448bp+0, 0x1.de55439a2c38bp-53,
-       0x1.8f1ae99157736p+0, 0x1.5cc13a2e3976cp-55,
-       0x1.9145b0b91ffc5p+0, 0x1.114c368d3ed6ep-53,
-       0x1.93737b0cdc5e4p+0, 0x1.e8a0387e4a814p-53,
-       0x1.95a44cbc8520ep+0, 0x1.d36906d2b41f9p-53,
-       0x1.97d829fde4e4fp+0, 0x1.173d241f23d18p-53,
-       0x1.9a0f170ca07b9p+0, 0x1.7462137188ce7p-53,
-       0x1.9c49182a3f090p+0, 0x1.c7c46b071f2bep-56,
-       0x1.9e86319e32323p+0, 0x1.824ca78e64c6ep-56,
-       0x1.a0c667b5de564p+0, 0x1.6535b51719567p-53,
-       0x1.a309bec4a2d33p+0, 0x1.6305c7ddc36abp-54,
-       0x1.a5503b23e255cp+0, 0x1.1684892395f0fp-53,
-       0x1.a799e1330b358p+0, 0x1.bcb7ecac563c7p-54,
-       0x1.a9e6b5579fdbfp+0, 0x1.0fac90ef7fd31p-54,
-       0x1.ac36bbfd3f379p+0, 0x1.81b72cd4624ccp-53,
-       0x1.ae89f995ad3adp+0, 0x1.7a1cd345dcc81p-54,
-       0x1.b0e07298db665p+0, 0x1.2108559bf8deep-53,
-       0x1.b33a2b84f15fap+0, 0x1.ed7fa1cf7b290p-53,
-       0x1.b59728de55939p+0, 0x1.1c7102222c90ep-53,
-       0x1.b7f76f2fb5e46p+0, 0x1.d54f610356a79p-53,
-       0x1.ba5b030a10649p+0, 0x1.0819678d5eb69p-53,
-       0x1.bcc1e904bc1d2p+0, 0x1.23dd07a2d9e84p-55,
-       0x1.bf2c25bd71e08p+0, 0x1.0811ae04a31c7p-53,
-       0x1.c199bdd85529cp+0, 0x1.11065895048ddp-55,
-       0x1.c40ab5fffd07ap+0, 0x1.b4537e083c60ap-54,
-       0x1.c67f12e57d14bp+0, 0x1.2884dff483cadp-54,
-       0x1.c8f6d9406e7b5p+0, 0x1.1acbc48805c44p-56,
-       0x1.cb720dcef9069p+0, 0x1.503cbd1e949dbp-56,
-       0x1.cdf0b555dc3f9p+0, 0x1.889f12b1f58a3p-53,
-       0x1.d072d4a07897bp+0, 0x1.1a1e45e4342b2p-53,
-       0x1.d2f87080d89f1p+0, 0x1.15bc247313d44p-53,
-       0x1.d5818dcfba487p+0, 0x1.2ed02d75b3707p-55,
-       0x1.d80e316c98397p+0, 0x1.7709f3a09100cp-53,
-       0x1.da9e603db3285p+0, 0x1.c2300696db532p-54,
-       0x1.dd321f301b460p+0, 0x1.2da5778f018c3p-54,
-       0x1.dfc97337b9b5ep+0, 0x1.72d195873da52p-53,
-       0x1.e264614f5a128p+0, 0x1.424ec3f42f5b5p-53,
-       0x1.e502ee78b3ff6p+0, 0x1.39e8980a9cc8fp-55,
-       0x1.e7a51fbc74c83p+0, 0x1.2d522ca0c8de2p-54,
-       0x1.ea4afa2a490d9p+0, 0x1.0b1ee7431ebb6p-53,
-       0x1.ecf482d8e67f0p+0, 0x1.1b60625f7293ap-53,
-       0x1.efa1bee615a27p+0, 0x1.dc7f486a4b6b0p-54,
-       0x1.f252b376bba97p+0, 0x1.3a1a5bf0d8e43p-54,
-       0x1.f50765b6e4540p+0, 0x1.9d3e12dd8a18bp-54,
-       0x1.f7bfdad9cbe13p+0, 0x1.1227697fce57bp-53,
-       0x1.fa7c1819e90d8p+0, 0x1.74853f3a5931ep-55,
-       0x1.fd3c22b8f71f1p+0, 0x1.2eb74966579e7p-57
-};
-
 long double
 expl(long double x)
 {
-       union IEEEl2bits u, v;
-       long double fn, q, r, r1, r2, t, twopk, twopkp10000;
-       long double z;
-       int k, n, n2;
+       union IEEEl2bits u;
+       long double hi, lo, t, twopk;
+       int k;
        uint16_t hx, ix;
 
+       DOPRINT_START(&x);
+
        /* Filter out exceptional cases. */
        u.e = x;
        hx = u.xbits.expsign;
@@ -246,59 +84,32 @@ expl(long double x)
        if (ix >= BIAS + 13) {          /* |x| >= 8192 or x is NaN */
                if (ix == BIAS + LDBL_MAX_EXP) {
                        if (hx & 0x8000)  /* x is -Inf, -NaN or unsupported */
-                               return (-1 / x);
-                       return (x + x); /* x is +Inf, +NaN or unsupported */
+                               RETURNP(-1 / x);
+                       RETURNP(x + x); /* x is +Inf, +NaN or unsupported */
                }
                if (x > o_threshold)
-                       return (huge * huge);
+                       RETURNP(huge * huge);
                if (x < u_threshold)
-                       return (tiny * tiny);
-       } else if (ix < BIAS - 65) {    /* |x| < 0x1p-65 (includes pseudos) */
-               return (1 + x);         /* 1 with inexact iff x != 0 */
+                       RETURNP(tiny * tiny);
+       } else if (ix < BIAS - 75) {    /* |x| < 0x1p-75 (includes pseudos) */
+               RETURN2P(1, x);         /* 1 with inexact iff x != 0 */
        }
 
        ENTERI();
 
-       /* Reduce x to (k*ln2 + endpoint[n2] + r1 + r2). */
-       /* Use a specialized rint() to get fn.  Assume round-to-nearest. */
-       fn = x * INV_L + 0x1.8p63 - 0x1.8p63;
-       r = x - fn * L1 - fn * L2;      /* r = r1 + r2 done independently. */
-#if defined(HAVE_EFFICIENT_IRINTL)
-       n  = irintl(fn);
-#elif defined(HAVE_EFFICIENT_IRINT)
-       n  = irint(fn);
-#else
-       n  = (int)fn;
-#endif
-       n2 = (unsigned)n % INTERVALS;
-       /* Depend on the sign bit being propagated: */
-       k = n >> LOG2_INTERVALS;
-       r1 = x - fn * L1;
-       r2 = fn * -L2;
-
-       /* Prepare scale factors. */
-       v.e = 1;
-       if (k >= LDBL_MIN_EXP) {
-               v.xbits.expsign = BIAS + k;
-               twopk = v.e;
-       } else {
-               v.xbits.expsign = BIAS + k + 10000;
-               twopkp10000 = v.e;
-       }
-
-       /* Evaluate expl(endpoint[n2] + r1 + r2) = tbl[n2] * expl(r1 + r2). */
-       z = r * r;
-       q = r2 + z * (A2 + r * A3) + z * z * (A4 + r * A5) + z * z * z * A6;
-       t = (long double)tbl[n2].lo + tbl[n2].hi;
-       t = tbl[n2].lo + t * (q + r1) + tbl[n2].hi;
+       twopk = 1;
+       __k_expl(x, &hi, &lo, &k);
+       t = SUM2P(hi, lo);
 
        /* Scale by 2**k. */
        if (k >= LDBL_MIN_EXP) {
                if (k == LDBL_MAX_EXP)
                        RETURNI(t * 2 * 0x1p16383L);
+               SET_LDBL_EXPSIGN(twopk, BIAS + k);
                RETURNI(t * twopk);
        } else {
-               RETURNI(t * twopkp10000 * twom10000);
+               SET_LDBL_EXPSIGN(twopk, BIAS + k + 10000);
+               RETURNI(t * twopk * twom10000);
        }
 }
 
@@ -325,8 +136,11 @@ T1 = -0.1659,                              /* ~-30.625/128 * log(2) */
 T2 =  0.1659;                          /* ~30.625/128 * log(2) */
 
 /*
- * Domain [-0.1659, 0.1659], range ~[-1.2027e-22, 3.4417e-22]:
- * |(exp(x)-1-x-x**2/2)/x - p(x)| < 2**-71.2
+ * Domain [-0.1659, 0.1659], range ~[-2.6155e-22, 2.5507e-23]:
+ * |(exp(x)-1-x-x**2/2)/x - p(x)| < 2**-71.6
+ *
+ * XXX the coeffs aren't very carefully rounded, and I get 2.8 more bits,
+ * but unlike for ld128 we can't drop any terms.
  */
 static const union IEEEl2bits
 B3 = LD80C(0xaaaaaaaaaaaaaaab, -3,  1.66666666666666666671e-1L),
@@ -352,6 +166,8 @@ expm1l(long double x)
        int k, n, n2;
        uint16_t hx, ix;
 
+       DOPRINT_START(&x);
+
        /* Filter out exceptional cases. */
        u.e = x;
        hx = u.xbits.expsign;
@@ -359,11 +175,11 @@ expm1l(long double x)
        if (ix >= BIAS + 6) {           /* |x| >= 64 or x is NaN */
                if (ix == BIAS + LDBL_MAX_EXP) {
                        if (hx & 0x8000)  /* x is -Inf, -NaN or unsupported */
-                               return (-1 / x - 1);
-                       return (x + x); /* x is +Inf, +NaN or unsupported */
+                               RETURNP(-1 / x - 1);
+                       RETURNP(x + x); /* x is +Inf, +NaN or unsupported */
                }
                if (x > o_threshold)
-                       return (huge * huge);
+                       RETURNP(huge * huge);
                /*
                 * expm1l() never underflows, but it must avoid
                 * unrepresentable large negative exponents.  We used a
@@ -372,15 +188,15 @@ expm1l(long double x)
                 * in the same way as large ones here.
                 */
                if (hx & 0x8000)        /* x <= -64 */
-                       return (tiny - 1);      /* good for x < -65ln2 - eps */
+                       RETURN2P(tiny, -1);     /* good for x < -65ln2 - eps */
        }
 
        ENTERI();
 
        if (T1 < x && x < T2) {
-               if (ix < BIAS - 64) {   /* |x| < 0x1p-64 (includes pseudos) */
+               if (ix < BIAS - 74) {   /* |x| < 0x1p-74 (includes pseudos) */
                        /* x (rounded) with inexact if x != 0: */
-                       RETURNI(x == 0 ? x :
+                       RETURNPI(x == 0 ? x :
                            (0x1p100 * x + fabsl(x)) * 0x1p-100);
                }
 
@@ -401,9 +217,9 @@ expm1l(long double x)
                hx2_hi = x_hi * x_hi / 2;
                hx2_lo = x_lo * (x + x_hi) / 2;
                if (ix >= BIAS - 7)
-                       RETURNI(hx2_lo + x_lo + q + (hx2_hi + x_hi));
+                       RETURN2PI(hx2_hi + x_hi, hx2_lo + x_lo + q);
                else
-                       RETURNI(hx2_lo + q + hx2_hi + x);
+                       RETURN2PI(x, hx2_lo + q + hx2_hi);
        }
 
        /* Reduce x to (k*ln2 + endpoint[n2] + r1 + r2). */
@@ -437,21 +253,21 @@ expm1l(long double x)
        t = (long double)tbl[n2].lo + tbl[n2].hi;
 
        if (k == 0) {
-               t = tbl[n2].lo * (r1 + 1) + t * q + tbl[n2].hi * r1 +
-                   (tbl[n2].hi - 1);
+               t = SUM2P(tbl[n2].hi - 1, tbl[n2].lo * (r1 + 1) + t * q +
+                   tbl[n2].hi * r1);
                RETURNI(t);
        }
        if (k == -1) {
-               t = tbl[n2].lo * (r1 + 1) + t * q + tbl[n2].hi * r1 + 
-                   (tbl[n2].hi - 2);
+               t = SUM2P(tbl[n2].hi - 2, tbl[n2].lo * (r1 + 1) + t * q +
+                   tbl[n2].hi * r1);
                RETURNI(t / 2);
        }
        if (k < -7) {
-               t = tbl[n2].lo + t * (q + r1) + tbl[n2].hi;
+               t = SUM2P(tbl[n2].hi, tbl[n2].lo + t * (q + r1));
                RETURNI(t * twopk - 1);
        }
        if (k > 2 * LDBL_MANT_DIG - 1) {
-               t = tbl[n2].lo + t * (q + r1) + tbl[n2].hi;
+               t = SUM2P(tbl[n2].hi, tbl[n2].lo + t * (q + r1));
                if (k == LDBL_MAX_EXP)
                        RETURNI(t * 2 * 0x1p16383L - 1);
                RETURNI(t * twopk - 1);
@@ -461,8 +277,8 @@ expm1l(long double x)
        twomk = v.e;
 
        if (k > LDBL_MANT_DIG - 1)
-               t = tbl[n2].lo - twomk + t * (q + r1) + tbl[n2].hi;
+               t = SUM2P(tbl[n2].hi, tbl[n2].lo - twomk + t * (q + r1));
        else
-               t = tbl[n2].lo + t * (q + r1) + (tbl[n2].hi - twomk);
+               t = SUM2P(tbl[n2].hi - twomk, tbl[n2].lo + t * (q + r1));
        RETURNI(t * twopk);
 }
index bc95df4..fd9d85e 100644 (file)
 .\"     from: @(#)cosh.3       5.1 (Berkeley) 5/2/91
 .\" $FreeBSD: head/lib/msun/man/cosh.3 226458 2011-10-17 05:41:03Z das $
 .\"
-.Dd January 14, 2005
+.Dd August 17, 2013
 .Dt COSH 3
 .Os
 .Sh NAME
 .Nm cosh ,
-.Nm coshf
+.Nm coshf ,
+.Nm coshl
 .Nd hyperbolic cosine functions
 .Sh LIBRARY
 .Lb libm
 .Fn cosh "double x"
 .Ft float
 .Fn coshf "float x"
+.Ft long double
+.Fn coshl "long double x"
 .Sh DESCRIPTION
 The
-.Fn cosh
-and the
-.Fn coshf
+.Fn cosh ,
+.Fn coshf ,
+and
+.Fn coshl
 functions compute the hyperbolic cosine of
 .Fa x .
 .Sh SEE ALSO
index c34234d..a139895 100644 (file)
 .\"     from: @(#)erf.3        6.4 (Berkeley) 4/20/91
 .\" $FreeBSD: head/lib/msun/man/erf.3 165906 2007-01-09 01:02:06Z imp $
 .\"
-.Dd April 20, 1991
+.Dd July 13, 2014
 .Dt ERF 3
 .Os
 .Sh NAME
 .Nm erf ,
 .Nm erff ,
+.Nm erfl ,
 .Nm erfc ,
-.Nm erfcf
+.Nm erfcf ,
+.Nm erfcl
 .Nd error function operators
 .Sh LIBRARY
 .Lb libm
 .Fn erf "double x"
 .Ft float
 .Fn erff "float x"
+.Ft "long double"
+.Fn erfl "long double x"
 .Ft double
 .Fn erfc "double x"
 .Ft float
 .Fn erfcf "float x"
+.Ft "long double"
+.Fn erfcl "long double x"
 .Sh DESCRIPTION
 These functions calculate the error function of
 .Fa x .
 .Pp
 The
-.Fn erf
-and the
-.Fn erff
+.Fn erf ,
+.Fn erff ,
+and
+.Fn erfl
 functions calculate the error function of x; where
 .Bd -ragged -offset indent
 .if n \{\
@@ -69,9 +76,10 @@ erf\|(x) :=
 .Ed
 .Pp
 The
-.Fn erfc
-and the
-.Fn erfcf
+.Fn erfc ,
+.Fn erfcf ,
+and
+.Fn erfcl
 functions calculate the complementary error function of
 .Fa x ;
 that is
@@ -79,9 +87,6 @@ that is
 subtracts the result of the error function
 .Fn erf x
 from 1.0.
-This is useful, since for large
-.Fa x
-places disappear.
 .Sh SEE ALSO
 .Xr math 3
 .Sh HISTORY
index c477e56..bc269d9 100644 (file)
@@ -28,7 +28,7 @@
 .\"     from: @(#)lgamma.3     6.6 (Berkeley) 12/3/92
 .\" $FreeBSD: head/lib/msun/man/lgamma.3 176388 2008-02-18 17:27:11Z das $
 .\"
-.Dd December 21, 2011
+.Dd September 12, 2014
 .Dt LGAMMA 3
 .Os
 .Sh NAME
@@ -36,6 +36,8 @@
 .Nm lgamma_r ,
 .Nm lgammaf ,
 .Nm lgammaf_r ,
+.Nm lgammal ,
+.Nm lgammal_r ,
 .Nm gamma ,
 .Nm gamma_r ,
 .Nm gammaf ,
 .Fn lgammaf "float x"
 .Ft float
 .Fn lgammaf_r "float x" "int *signgamp"
+.Ft "long double"
+.Fn lgammal "long double x"
+.Ft "long double"
+.Fn lgammal_r "long double x" "int *signgamp"
 .Ft double
 .Fn gamma "double x"
 .Ft double
 .Fn gammaf "float x"
 .Ft float
 .Fn gammaf_r "float x" "int *signgamp"
-.Ft double
+.Ft "long double"
 .Fn tgamma "double x"
 .Ft float
 .Fn tgammaf "float x"
 .Sh DESCRIPTION
-.Fn lgamma x
+.Fn lgamma x ,
+.Fn lgammaf x ,
 and
-.Fn lgammaf x
+.Fn lgammal x
 .if t \{\
 return ln\||\(*G(x)| where
 .Bd -unfilled -offset indent
@@ -87,13 +94,15 @@ The external integer
 .Fa signgam
 returns the sign of \(*G(x).
 .Pp
-.Fn lgamma_r x signgamp
+.Fn lgamma_r x signgamp ,
+.Fn lgammaf_r x signgamp ,
 and
-.Fn lgammaf_r x signgamp
+.Fn lgammal_r x signgamp
 provide the same functionality as
-.Fn lgamma x
+.Fn lgamma x , 
+.Fn lgammaf x ,
 and
-.Fn lgammaf x
+.Fn lgammal x ,
 but the caller must provide an integer to store the sign of \(*G(x).
 .Pp
 The
@@ -115,6 +124,7 @@ are deprecated aliases for
 and
 .Fn lgammaf_r ,
 respectively.
+
 .Sh IDIOSYNCRASIES
 Do not use the expression
 .Dq Li signgam\(**exp(lgamma(x))
@@ -139,14 +149,18 @@ Exponentiation of
 will lose up to 10 significant bits.
 .Sh RETURN VALUES
 .Fn gamma ,
-.Fn gamma_r ,
 .Fn gammaf ,
+.Fn gammal ,
+.Fn gamma_r ,
 .Fn gammaf_r ,
+.Fn gammal_r ,
 .Fn lgamma ,
-.Fn lgamma_r ,
 .Fn lgammaf ,
+.Fn lgammal ,
+.Fn lgamma_r ,
+.Fn lgammaf_r ,
 and
-.Fn lgammaf_r
+.Fn lgammal_r
 return appropriate values unless an argument is out of range.
 Overflow will occur for sufficiently large positive values, and
 non-positive integers.
@@ -159,6 +173,7 @@ will underflow.
 The
 .Fn lgamma ,
 .Fn lgammaf ,
+.Fn lgammal ,
 .Fn tgamma ,
 and
 .Fn tgammaf
index 66ade3f..1f85f8e 100644 (file)
 .\"
 .\"    from: @(#)sinh.3        6.6 (Berkeley) 4/19/91
 .\" $FreeBSD: head/lib/msun/man/sinh.3 226458 2011-10-17 05:41:03Z das $
-.Dd December 21, 2011
+.\"
+.Dd August 17, 2013
 .Dt SINH 3
 .Os
 .Sh NAME
 .Nm sinh ,
-.Nm sinhf
+.Nm sinhf ,
+.Nm sinhl
 .Nd hyperbolic sine function
 .Sh LIBRARY
 .Lb libm
 .Fn sinh "double x"
 .Ft float
 .Fn sinhf "float x"
+.Ft long double
+.Fn sinhl "long double x"
 .Sh DESCRIPTION
 The
-.Fn sinh
-and the
-.Fn sinhf
+.Fn sinh ,
+.Fn sinhf ,
+and
+.Fn sinhl
 functions compute the hyperbolic sine of
 .Fa x .
 .Sh SEE ALSO
index c15c833..bd3c66a 100644 (file)
 .\"     from: @(#)tanh.3       5.1 (Berkeley) 5/2/91
 .\" $FreeBSD: head/lib/msun/man/tanh.3 226458 2011-10-17 05:41:03Z das $
 .\"
-.Dd September 18, 2011
+.Dd August 17, 2013
 .Dt TANH 3
 .Os
 .Sh NAME
 .Nm tanh ,
-.Nm tanhf
+.Nm tanhf ,
+.Nm tanhl
 .Nd hyperbolic tangent functions
 .Sh LIBRARY
 .Lb libm
 .Fn tanh "double x"
 .Ft float
 .Fn tanhf "float x"
+.Ft long double
+.Fn tanhl "long double x"
 .Sh DESCRIPTION
 The
-.Fn tanh
-and the
-.Fn tanhf
+.Fn tanh ,
+.Fn tanhf ,
+and
+.Fn tanhl
 functions compute the hyperbolic tangent of
 .Fa x .
 For a discussion of error due to roundoff, see
 .Xr math 3 .
 .Sh RETURN VALUES
 The
-.Fn tanh
+.Fn tanh ,
+.Fn tanhf ,
 and the
-.Fn tanhf
+.Fn tanhl
 functions return the hyperbolic tangent value.
 .Sh SEE ALSO
 .Xr acos 3 ,
index f407f10..7d45f8f 100644 (file)
@@ -150,13 +150,13 @@ f(double a, double b, double hypot_a_b)
  */
 static inline void
 do_hard_work(double x, double y, double *rx, int *B_is_usable, double *B,
-            double *sqrt_A2my2, double *new_y)
+    double *sqrt_A2my2, double *new_y)
 {
        double R, S, A; /* A, B, R, and S are as in Hull et al. */
        double Am1, Amy; /* A-1, A-y. */
 
-       R = hypot(x, y + 1); /* |z+I| */
-       S = hypot(x, y - 1); /* |z-I| */
+       R = hypot(x, y + 1);            /* |z+I| */
+       S = hypot(x, y - 1);            /* |z-I| */
 
        /* A = (|z+I| + |z-I|) / 2 */
        A = (R + S) / 2;
@@ -250,7 +250,7 @@ do_hard_work(double x, double y, double *rx, int *B_is_usable, double *B,
                         * scaling should avoid any underflow problems.
                         */
                        *sqrt_A2my2 = x * (4 / DBL_EPSILON / DBL_EPSILON) * y /
-                               sqrt((y + 1) * (y - 1));
+                           sqrt((y + 1) * (y - 1));
                        *new_y = y * (4 / DBL_EPSILON / DBL_EPSILON);
                } else {                /* if (y < 1) */
                        /*
@@ -587,7 +587,7 @@ catanh(double complex z)
                /* catanh(NaN + I*+-Inf) = sign(NaN)0 + I*+-PI/2 */
                if (isinf(y))
                        return (cpack(copysign(0, x),
-                                     copysign(pio2_hi + pio2_lo, y)));
+                           copysign(pio2_hi + pio2_lo, y)));
                /*
                 * All other cases involving NaN return NaN + I*NaN.
                 * C99 leaves it optional whether to raise invalid if one of
@@ -598,7 +598,7 @@ catanh(double complex z)
 
        if (ax > RECIP_EPSILON || ay > RECIP_EPSILON)
                return (cpack(real_part_reciprocal(x, y),
-                             copysign(pio2_hi + pio2_lo, y)));
+                   copysign(pio2_hi + pio2_lo, y)));
 
        if (ax < SQRT_3_EPSILON / 2 && ay < SQRT_3_EPSILON / 2) {
                /*
index 5c787c9..cdca3c6 100644 (file)
@@ -85,7 +85,7 @@ f(float a, float b, float hypot_a_b)
 
 static inline void
 do_hard_work(float x, float y, float *rx, int *B_is_usable, float *B,
-            float *sqrt_A2my2, float *new_y)
+    float *sqrt_A2my2, float *new_y)
 {
        float R, S, A;
        float Am1, Amy;
@@ -133,7 +133,7 @@ do_hard_work(float x, float y, float *rx, int *B_is_usable, float *B,
                        *sqrt_A2my2 = sqrtf(Amy * (A + y));
                } else if (y > 1) {
                        *sqrt_A2my2 = x * (4 / FLT_EPSILON / FLT_EPSILON) * y /
-                               sqrtf((y + 1) * (y - 1));
+                           sqrtf((y + 1) * (y - 1));
                        *new_y = y * (4 / FLT_EPSILON / FLT_EPSILON);
                } else {
                        *sqrt_A2my2 = sqrtf((1 - y) * (1 + y));
@@ -169,7 +169,7 @@ casinhf(float complex z)
                else
                        w = clog_for_large_values(-z) + m_ln2;
                return (cpackf(copysignf(crealf(w), x),
-                              copysignf(cimagf(w), y)));
+                   copysignf(cimagf(w), y)));
        }
 
        if (x == 0 && y == 0)
@@ -291,7 +291,7 @@ clog_for_large_values(float complex z)
 
        if (ax > FLT_MAX / 2)
                return (cpackf(logf(hypotf(x / m_e, y / m_e)) + 1,
-                              atan2f(y, x)));
+                   atan2f(y, x)));
 
        if (ax > QUARTER_SQRT_MAX || ay < SQRT_MIN)
                return (cpackf(logf(hypotf(x, y)), atan2f(y, x)));
@@ -355,13 +355,13 @@ catanhf(float complex z)
                        return (cpackf(copysignf(0, x), y + y));
                if (isinf(y))
                        return (cpackf(copysignf(0, x),
-                                      copysignf(pio2_hi + pio2_lo, y)));
+                           copysignf(pio2_hi + pio2_lo, y)));
                return (cpackf(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
        }
 
        if (ax > RECIP_EPSILON || ay > RECIP_EPSILON)
                return (cpackf(real_part_reciprocal(x, y),
-                              copysignf(pio2_hi + pio2_lo, y)));
+                   copysignf(pio2_hi + pio2_lo, y)));
 
        if (ax < SQRT_3_EPSILON / 2 && ay < SQRT_3_EPSILON / 2) {
                raise_inexact();
index b7d8cfc..a82d7c6 100644 (file)
@@ -1,6 +1,5 @@
 
 /* @(#)e_cosh.c 1.3 95/01/18 */
-/* $FreeBSD: head/lib/msun/src/e_cosh.c 226598 2011-10-21 06:28:47Z das $ */
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
@@ -12,6 +11,7 @@
  * ====================================================
  */
 
+
 /* __ieee754_cosh(x)
  * Method : 
  * mathematically cosh(x) if defined to be (exp(x)+exp(-x))/2
@@ -33,6 +33,8 @@
  *     only cosh(0)=1 is exact for finite x.
  */
 
+#include <float.h>
+
 #include "math.h"
 #include "math_private.h"
 
@@ -75,3 +77,7 @@ __ieee754_cosh(double x)
     /* |x| > overflowthresold, cosh(x) overflow */
        return huge*huge;
 }
+
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(cosh, coshl);
+#endif
diff --git a/lib/libm/src/e_coshl.c b/lib/libm/src/e_coshl.c
new file mode 100644 (file)
index 0000000..57beba8
--- /dev/null
@@ -0,0 +1,128 @@
+/* from: FreeBSD: head/lib/msun/src/e_coshl.c XXX */
+
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+
+/*
+ * See e_cosh.c for complete comments.
+ *
+ * Converted to long double by Bruce D. Evans.
+ */
+
+#include <float.h>
+#ifdef __i386__
+#include <ieeefp.h>
+#endif
+
+#include "fpmath.h"
+#include "math.h"
+#include "math_private.h"
+#include "k_expl.h"
+
+#if LDBL_MAX_EXP != 0x4000
+/* We also require the usual expsign encoding. */
+#error "Unsupported long double format"
+#endif
+
+#define        BIAS    (LDBL_MAX_EXP - 1)
+
+static const volatile long double huge = 0x1p10000L, tiny = 0x1p-10000L;
+#if LDBL_MANT_DIG == 64
+/*
+ * Domain [-1, 1], range ~[-1.8211e-21, 1.8211e-21]:
+ * |cosh(x) - c(x)| < 2**-68.8
+ */
+static const union IEEEl2bits
+C4u = LD80C(0xaaaaaaaaaaaaac78, -5,  4.16666666666666682297e-2L);
+#define        C4      C4u.e
+static const double
+C2  =  0.5,
+C6  =  1.3888888888888616e-3,          /*  0x16c16c16c16b99.0p-62 */
+C8  =  2.4801587301767953e-5,          /*  0x1a01a01a027061.0p-68 */
+C10 =  2.7557319163300398e-7,          /*  0x127e4fb6c9b55f.0p-74 */
+C12 =  2.0876768371393075e-9,          /*  0x11eed99406a3f4.0p-81 */
+C14 =  1.1469537039374480e-11,         /*  0x1938c67cd18c48.0p-89 */
+C16 =  4.8473490896852041e-14;         /*  0x1b49c429701e45.0p-97 */
+#elif LDBL_MANT_DIG == 113
+/*
+ * Domain [-1, 1], range ~[-2.3194e-37, 2.3194e-37]:
+ * |cosh(x) - c(x)| < 2**-121.69
+ */
+static const long double
+C4  =  4.16666666666666666666666666666666225e-2L,      /*  0x1555555555555555555555555554e.0p-117L */
+C6  =  1.38888888888888888888888888889434831e-3L,      /*  0x16c16c16c16c16c16c16c16c1dd7a.0p-122L */
+C8  =  2.48015873015873015873015871870962089e-5L,      /*  0x1a01a01a01a01a01a01a017af2756.0p-128L */
+C10 =  2.75573192239858906525574318600800201e-7L,      /*  0x127e4fb7789f5c72ef01c8a040640.0p-134L */
+C12 =  2.08767569878680989791444691755468269e-9L,      /*  0x11eed8eff8d897b543d0679607399.0p-141L */
+C14=  1.14707455977297247387801189650495351e-11L,      /*  0x193974a8c07c9d24ae169a7fa9b54.0p-149L */
+C16 =  4.77947733238737883626416876486279985e-14L;     /*  0x1ae7f3e733b814d4e1b90f5727fe4.0p-157L */
+static const double
+C2  =  0.5,
+C18 =  1.5619206968597871e-16,         /*  0x16827863b9900b.0p-105 */
+C20 =  4.1103176218528049e-19,         /*  0x1e542ba3d3c269.0p-114 */
+C22 =  8.8967926401641701e-22,         /*  0x10ce399542a014.0p-122 */
+C24 =  1.6116681626523904e-24,         /*  0x1f2c981d1f0cb7.0p-132 */
+C26 =  2.5022374732804632e-27;         /*  0x18c7ecf8b2c4a0.0p-141 */
+#else
+#error "Unsupported long double format"
+#endif /* LDBL_MANT_DIG == 64 */
+
+/* log(2**16385 - 0.5) rounded up: */
+static const float
+o_threshold =  1.13572168e4;           /*  0xb174de.0p-10 */
+
+long double
+coshl(long double x)
+{
+       long double hi,lo,x2,x4;
+       double dx2;
+       uint16_t ix;
+
+       GET_LDBL_EXPSIGN(ix,x);
+       ix &= 0x7fff;
+
+    /* x is INF or NaN */
+       if(ix>=0x7fff) return x*x;
+
+       ENTERI();
+
+    /* |x| < 1, return 1 or c(x) */
+       if(ix<0x3fff) {
+           if (ix<BIAS-(LDBL_MANT_DIG+1)/2)    /* |x| < TINY */
+               RETURNI(1+tiny);        /* cosh(tiny) = 1(+) with inexact */
+           x2 = x*x;
+#if LDBL_MANT_DIG == 64
+           x4 = x2*x2;
+           RETURNI(((C16*x2 + C14)*x4 + (C12*x2 + C10))*(x4*x4*x2) +
+               ((C8*x2 + C6)*x2 + C4)*x4 + C2*x2 + 1);
+#elif LDBL_MANT_DIG == 113
+           dx2 = x2;
+           RETURNI((((((((((((C26*dx2 + C24)*dx2 + C22)*dx2 +
+               C20)*x2 + C18)*x2 +
+               C16)*x2 + C14)*x2 + C12)*x2 + C10)*x2 + C8)*x2 + C6)*x2 +
+               C4)*(x2*x2) + C2*x2 + 1);
+#endif
+       }
+
+    /* |x| in [1, 64), return accurate exp(|x|)/2+1/exp(|x|)/2 */
+       if (ix < 0x4005) {
+           k_hexpl(fabsl(x), &hi, &lo);
+           RETURNI(lo + 0.25/(hi + lo) + hi);
+       }
+
+    /* |x| in [64, o_threshold], return correctly-overflowing exp(|x|)/2 */
+       if (fabsl(x) <= o_threshold)
+           RETURNI(hexpl(fabsl(x)));
+
+    /* |x| > o_threshold, cosh(x) overflow */
+       RETURNI(huge*huge);
+}
index e9cfde4..4451362 100644 (file)
@@ -1,6 +1,5 @@
 
 /* @(#)e_gamma.c 1.3 95/01/18 */
-/* $FreeBSD: head/lib/msun/src/e_gamma.c 176451 2008-02-22 02:30:36Z das $ */
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
@@ -13,6 +12,7 @@
  *
  */
 
+
 /* __ieee754_gamma(x)
  * Return the logarithm of the Gamma function of x.
  *
index 0259a38..d689320 100644 (file)
@@ -1,6 +1,5 @@
 
 /* @(#)e_lgamma.c 1.3 95/01/18 */
-/* $FreeBSD: head/lib/msun/src/e_lgamma.c 176451 2008-02-22 02:30:36Z das $ */
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  *
  */
 
+
 /* __ieee754_lgamma(x)
  * Return the logarithm of the Gamma function of x.
  *
  * Method: call __ieee754_lgamma_r
  */
 
+#include <float.h>
+
 #include "math.h"
 #include "math_private.h"
 
@@ -29,3 +31,7 @@ __ieee754_lgamma(double x)
 {
        return __ieee754_lgamma_r(x,&signgam);
 }
+
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(lgamma, lgammal);
+#endif
index 79f3074..75996e3 100644 (file)
@@ -1,25 +1,23 @@
-
 /* @(#)e_lgamma_r.c 1.3 95/01/18 */
-/* $FreeBSD: head/lib/msun/src/e_lgamma_r.c 226380 2011-10-15 07:00:28Z 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 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
- *
  */
 
+
 /* __ieee754_lgamma_r(x, signgamp)
- * Reentrant version of the logarithm of the Gamma function 
- * with user provide pointer for the sign of Gamma(x). 
+ * Reentrant version of the logarithm of the Gamma function
+ * with user provide pointer for the sign of Gamma(x).
  *
  * Method:
  *   1. Argument Reduction for 0 < x <= 8
- *     Since gamma(1+s)=s*gamma(s), for x in [0,8], we may 
+ *     Since gamma(1+s)=s*gamma(s), for x in [0,8], we may
  *     reduce x to a number in [1.5,2.5] by
  *             lgamma(1+s) = log(s) + lgamma(s)
  *     for example,
  *     by
  *                                 3       5             11
  *             w = w0 + w1*z + w2*z  + w3*z  + ... + w6*z
- *     where 
+ *     where
  *             |w - f(z)| < 2**-58.74
- *             
+ *
  *   4. For negative x, since (G is gamma function)
  *             -x*G(-x)*G(x) = pi/sin(pi*x),
  *     we have
  *             G(x) = pi/(sin(pi*x)*(-x)*G(-x))
  *     since G(-x) is positive, sign(G(x)) = sign(sin(pi*x)) for x<0
- *     Hence, for x<0, signgam = sign(sin(pi*x)) and 
+ *     Hence, for x<0, signgam = sign(sin(pi*x)) and
  *             lgamma(x) = log(|Gamma(x)|)
  *                       = log(pi/(|x*sin(pi*x)|)) - lgamma(-x);
- *     Note: one should avoid compute pi*(-x) directly in the 
+ *     Note: one should avoid compute pi*(-x) directly in the
  *           computation of sin(pi*(-x)).
- *             
+ *
  *   5. Special Cases
  *             lgamma(2+s) ~ s*(1-Euler) for tiny s
  *             lgamma(1) = lgamma(2) = 0
  *             lgamma(0) = lgamma(neg.integer) = inf and raise divide-by-zero
  *             lgamma(inf) = inf
  *             lgamma(-inf) = inf (bug for bug compatible with C99!?)
- *     
  */
 
+#include <float.h>
+
 #include "math.h"
 #include "math_private.h"
 
+static const volatile double vzero = 0;
+
 static const double
-two52=  4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
+zero=  0.00000000000000000000e+00,
 half=  5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
 one =  1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
 pi  =  3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */
@@ -152,44 +153,40 @@ w4  = -5.95187557450339963135e-04, /* 0xBF4380CB, 0x8C0FE741 */
 w5  =  8.36339918996282139126e-04, /* 0x3F4B67BA, 0x4CDAD5D1 */
 w6  = -1.63092934096575273989e-03; /* 0xBF5AB89D, 0x0B9E43E4 */
 
-static const double zero=  0.00000000000000000000e+00;
-
-       static double sin_pi(double x)
+/*
+ * Compute sin(pi*x) without actually doing the pi*x multiplication.
+ * sin_pi(x) is only called for x < 0 and |x| < 2**(p-1) where p is
+ * the precision of x.
+ */
+static double
+sin_pi(double x)
 {
+       volatile double vz;
        double y,z;
-       int n,ix;
+       int n;
+
+       y = -x;
 
-       GET_HIGH_WORD(ix,x);
-       ix &= 0x7fffffff;
+       vz = y+0x1p52;                  /* depend on 0 <= y < 0x1p52 */
+       z = vz-0x1p52;                  /* rint(y) for the above range */
+       if (z == y)
+           return zero;
 
-       if(ix<0x3fd00000) return __kernel_sin(pi*x,zero,0);
-       y = -x;         /* x is assume negative */
+       vz = y+0x1p50;
+       GET_LOW_WORD(n,vz);             /* bits for rounded y (units 0.25) */
+       z = vz-0x1p50;                  /* y rounded to a multiple of 0.25 */
+       if (z > y) {
+           z -= 0.25;                  /* adjust to round down */
+           n--;
+       }
+       n &= 7;                         /* octant of y mod 2 */
+       y = y - z + n * 0.25;           /* y mod 2 */
 
-    /*
-     * argument reduction, make sure inexact flag not raised if input
-     * is an integer
-     */
-       z = floor(y);
-       if(z!=y) {                              /* inexact anyway */
-           y  *= 0.5;
-           y   = 2.0*(y - floor(y));           /* y = |x| mod 2.0 */
-           n   = (int) (y*4.0);
-       } else {
-            if(ix>=0x43400000) {
-                y = zero; n = 0;                 /* y must be even */
-            } else {
-                if(ix<0x43300000) z = y+two52; /* exact */
-               GET_LOW_WORD(n,z);
-               n &= 1;
-                y  = n;
-                n<<= 2;
-            }
-        }
        switch (n) {
            case 0:   y =  __kernel_sin(pi*y,zero,0); break;
-           case 1:   
+           case 1:
            case 2:   y =  __kernel_cos(pi*(0.5-y),zero); break;
-           case 3:  
+           case 3:
            case 4:   y =  __kernel_sin(pi*(one-y),zero,0); break;
            case 5:
            case 6:   y = -__kernel_cos(pi*(y-1.5),zero); break;
@@ -202,34 +199,38 @@ static const double zero=  0.00000000000000000000e+00;
 double
 __ieee754_lgamma_r(double x, int *signgamp)
 {
-       double t,y,z,nadj,p,p1,p2,p3,q,r,w;
+       double nadj,p,p1,p2,p3,q,r,t,w,y,z;
        int32_t hx;
-       int i,lx,ix;
+       int i,ix,lx;
 
        EXTRACT_WORDS(hx,lx,x);
 
-    /* purge off +-inf, NaN, +-0, tiny and negative arguments */
+    /* purge +-Inf and NaNs */
        *signgamp = 1;
        ix = hx&0x7fffffff;
        if(ix>=0x7ff00000) return x*x;
-       if((ix|lx)==0) return one/zero;
-       if(ix<0x3b900000) {     /* |x|<2**-70, return -log(|x|) */
-           if(hx<0) {
-               *signgamp = -1;
-               return -__ieee754_log(-x);
-           } else return -__ieee754_log(x);
+
+    /* purge +-0 and tiny arguments */
+       *signgamp = 1-2*((uint32_t)hx>>31);
+       if(ix<0x3c700000) {     /* |x|<2**-56, return -log(|x|) */
+           if((ix|lx)==0)
+               return one/vzero;
+           return -__ieee754_log(fabs(x));
        }
+
+    /* purge negative integers and start evaluation for other x < 0 */
        if(hx<0) {
+           *signgamp = 1;
            if(ix>=0x43300000)  /* |x|>=2**52, must be -integer */
-               return one/zero;
+               return one/vzero;
            t = sin_pi(x);
-           if(t==zero) return one/zero; /* -integer */
+           if(t==zero) return one/vzero; /* -integer */
            nadj = __ieee754_log(pi/fabs(t*x));
            if(t<zero) *signgamp = -1;
            x = -x;
        }
 
-    /* purge off 1 and 2 */
+    /* purge 1 and 2 */
        if((((ix-0x3ff00000)|lx)==0)||(((ix-0x40000000)|lx)==0)) r = 0;
     /* for x < 2.0 */
        else if(ix<0x40000000) {
@@ -250,7 +251,7 @@ __ieee754_lgamma_r(double x, int *signgamp)
                p1 = a0+z*(a2+z*(a4+z*(a6+z*(a8+z*a10))));
                p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*a11)))));
                p  = y*p1+p2;
-               r  += (p-0.5*y); break;
+               r  += p-y/2; break;
              case 1:
                z = y*y;
                w = z*y;
@@ -258,38 +259,43 @@ __ieee754_lgamma_r(double x, int *signgamp)
                p2 = t1+w*(t4+w*(t7+w*(t10+w*t13)));
                p3 = t2+w*(t5+w*(t8+w*(t11+w*t14)));
                p  = z*p1-(tt-w*(p2+y*p3));
-               r += (tf + p); break;
-             case 2:   
+               r += tf + p; break;
+             case 2:
                p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*u5)))));
                p2 = one+y*(v1+y*(v2+y*(v3+y*(v4+y*v5))));
-               r += (-0.5*y + p1/p2);
+               r += p1/p2-y/2;
            }
        }
-       else if(ix<0x40200000) {                        /* x < 8.0 */
-           i = (int)x;
-           y = x-(double)i;
+    /* x < 8.0 */
+       else if(ix<0x40200000) {
+           i = x;
+           y = x-i;
            p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6))))));
            q = one+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))));
-           r = half*y+p/q;
+           r = y/2+p/q;
            z = one;    /* lgamma(1+s) = log(s) + lgamma(s) */
            switch(i) {
-           case 7: z *= (y+6.0);       /* FALLTHRU */
-           case 6: z *= (y+5.0);       /* FALLTHRU */
-           case 5: z *= (y+4.0);       /* FALLTHRU */
-           case 4: z *= (y+3.0);       /* FALLTHRU */
-           case 3: z *= (y+2.0);       /* FALLTHRU */
+           case 7: z *= (y+6);         /* FALLTHRU */
+           case 6: z *= (y+5);         /* FALLTHRU */
+           case 5: z *= (y+4);         /* FALLTHRU */
+           case 4: z *= (y+3);         /* FALLTHRU */
+           case 3: z *= (y+2);         /* FALLTHRU */
                    r += __ieee754_log(z); break;
            }
-    /* 8.0 <= x < 2**58 */
-       } else if (ix < 0x43900000) {
+    /* 8.0 <= x < 2**56 */
+       } else if (ix < 0x43700000) {
            t = __ieee754_log(x);
            z = one/x;
            y = z*z;
            w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6)))));
            r = (x-half)*(t-one)+w;
-       } else 
-    /* 2**58 <= x <= inf */
+       } else
+    /* 2**56 <= x <= inf */
            r =  x*(__ieee754_log(x)-one);
        if(hx<0) r = nadj - r;
        return r;
 }
+
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(lgamma_r, lgammal_r);
+#endif
index 2f12b40..a722265 100644 (file)
@@ -1,5 +1,6 @@
 /* e_lgammaf_r.c -- float version of e_lgamma_r.c.
  * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ * Conversion to float fixed By Steven G. Kargl.
  */
 
 /*
  * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
- *
- * $FreeBSD: head/lib/msun/src/e_lgammaf_r.c 226380 2011-10-15 07:00:28Z das $
  */
 
+
 #include "math.h"
 #include "math_private.h"
 
+static const volatile float vzero = 0;
+
 static const float
-two23=  8.3886080000e+06, /* 0x4b000000 */
-half=  5.0000000000e-01, /* 0x3f000000 */
-one =  1.0000000000e+00, /* 0x3f800000 */
+zero=  0,
+half=  0.5,
+one =  1,
 pi  =  3.1415927410e+00, /* 0x40490fdb */
-a0  =  7.7215664089e-02, /* 0x3d9e233f */
-a1  =  3.2246702909e-01, /* 0x3ea51a66 */
-a2  =  6.7352302372e-02, /* 0x3d89f001 */
-a3  =  2.0580807701e-02, /* 0x3ca89915 */
-a4  =  7.3855509982e-03, /* 0x3bf2027e */
-a5  =  2.8905137442e-03, /* 0x3b3d6ec6 */
-a6  =  1.1927076848e-03, /* 0x3a9c54a1 */
-a7  =  5.1006977446e-04, /* 0x3a05b634 */
-a8  =  2.2086278477e-04, /* 0x39679767 */
-a9  =  1.0801156895e-04, /* 0x38e28445 */
-a10 =  2.5214456400e-05, /* 0x37d383a2 */
-a11 =  4.4864096708e-05, /* 0x383c2c75 */
-tc  =  1.4616321325e+00, /* 0x3fbb16c3 */
-tf  = -1.2148628384e-01, /* 0xbdf8cdcd */
-/* tt = -(tail of tf) */
-tt  =  6.6971006518e-09, /* 0x31e61c52 */
-t0  =  4.8383611441e-01, /* 0x3ef7b95e */
-t1  = -1.4758771658e-01, /* 0xbe17213c */
-t2  =  6.4624942839e-02, /* 0x3d845a15 */
-t3  = -3.2788541168e-02, /* 0xbd064d47 */
-t4  =  1.7970675603e-02, /* 0x3c93373d */
-t5  = -1.0314224288e-02, /* 0xbc28fcfe */
-t6  =  6.1005386524e-03, /* 0x3bc7e707 */
-t7  = -3.6845202558e-03, /* 0xbb7177fe */
-t8  =  2.2596477065e-03, /* 0x3b141699 */
-t9  = -1.4034647029e-03, /* 0xbab7f476 */
-t10 =  8.8108185446e-04, /* 0x3a66f867 */
-t11 = -5.3859531181e-04, /* 0xba0d3085 */
-t12 =  3.1563205994e-04, /* 0x39a57b6b */
-t13 = -3.1275415677e-04, /* 0xb9a3f927 */
-t14 =  3.3552918467e-04, /* 0x39afe9f7 */
-u0  = -7.7215664089e-02, /* 0xbd9e233f */
-u1  =  6.3282704353e-01, /* 0x3f2200f4 */
-u2  =  1.4549225569e+00, /* 0x3fba3ae7 */
-u3  =  9.7771751881e-01, /* 0x3f7a4bb2 */
-u4  =  2.2896373272e-01, /* 0x3e6a7578 */
-u5  =  1.3381091878e-02, /* 0x3c5b3c5e */
-v1  =  2.4559779167e+00, /* 0x401d2ebe */
-v2  =  2.1284897327e+00, /* 0x4008392d */
-v3  =  7.6928514242e-01, /* 0x3f44efdf */
-v4  =  1.0422264785e-01, /* 0x3dd572af */
-v5  =  3.2170924824e-03, /* 0x3b52d5db */
-s0  = -7.7215664089e-02, /* 0xbd9e233f */
-s1  =  2.1498242021e-01, /* 0x3e5c245a */
-s2  =  3.2577878237e-01, /* 0x3ea6cc7a */
-s3  =  1.4635047317e-01, /* 0x3e15dce6 */
-s4  =  2.6642270386e-02, /* 0x3cda40e4 */
-s5  =  1.8402845599e-03, /* 0x3af135b4 */
-s6  =  3.1947532989e-05, /* 0x3805ff67 */
-r1  =  1.3920053244e+00, /* 0x3fb22d3b */
-r2  =  7.2193557024e-01, /* 0x3f38d0c5 */
-r3  =  1.7193385959e-01, /* 0x3e300f6e */
-r4  =  1.8645919859e-02, /* 0x3c98bf54 */
-r5  =  7.7794247773e-04, /* 0x3a4beed6 */
-r6  =  7.3266842264e-06, /* 0x36f5d7bd */
-w0  =  4.1893854737e-01, /* 0x3ed67f1d */
-w1  =  8.3333335817e-02, /* 0x3daaaaab */
-w2  = -2.7777778450e-03, /* 0xbb360b61 */
-w3  =  7.9365057172e-04, /* 0x3a500cfd */
-w4  = -5.9518753551e-04, /* 0xba1c065c */
-w5  =  8.3633989561e-04, /* 0x3a5b3dd2 */
-w6  = -1.6309292987e-03; /* 0xbad5c4e8 */
-
-static const float zero=  0.0000000000e+00;
-
-       static float sin_pif(float x)
+/*
+ * Domain y in [0x1p-27, 0.27], range ~[-3.4599e-10, 3.4590e-10]:
+ * |(lgamma(2 - y) + 0.5 * y) / y - a(y)| < 2**-31.4
+ */
+a0  =  7.72156641e-02, /* 0x3d9e233f */
+a1  =  3.22467119e-01, /* 0x3ea51a69 */
+a2  =  6.73484802e-02, /* 0x3d89ee00 */
+a3  =  2.06395667e-02, /* 0x3ca9144f */
+a4  =  6.98275631e-03, /* 0x3be4cf9b */
+a5  =  4.11768444e-03, /* 0x3b86eda4 */
+/*
+ * Domain x in [tc-0.24, tc+0.28], range ~[-5.6577e-10, 5.5677e-10]:
+ * |(lgamma(x) - tf) - t(x - tc)| < 2**-30.8.
+ */
+tc  =  1.46163213e+00, /* 0x3fbb16c3 */
+tf  = -1.21486291e-01, /* 0xbdf8cdce */
+t0  = -2.94064460e-11, /* 0xae0154b7 */
+t1  = -2.35939837e-08, /* 0xb2caabb8 */
+t2  =  4.83836412e-01, /* 0x3ef7b968 */
+t3  = -1.47586212e-01, /* 0xbe1720d7 */
+t4  =  6.46013096e-02, /* 0x3d844db1 */
+t5  = -3.28450352e-02, /* 0xbd068884 */
+t6  =  1.86483748e-02, /* 0x3c98c47a */
+t7  = -9.89206228e-03, /* 0xbc221251 */
+/*
+ * Domain y in [-0.1, 0.232], range ~[-8.4931e-10, 8.7794e-10]:
+ * |(lgamma(1 + y) + 0.5 * y) / y - u(y) / v(y)| < 2**-31.2
+ */
+u0  = -7.72156641e-02, /* 0xbd9e233f */
+u1  =  7.36789703e-01, /* 0x3f3c9e40 */
+u2  =  4.95649040e-01, /* 0x3efdc5b6 */
+v1  =  1.10958421e+00, /* 0x3f8e06db */
+v2  =  2.10598111e-01, /* 0x3e57a708 */
+v3  = -1.02995494e-02, /* 0xbc28bf71 */
+/*
+ * Domain x in (2, 3], range ~[-5.5189e-11, 5.2317e-11]:
+ * |(lgamma(y+2) - 0.5 * y) / y - s(y)/r(y)| < 2**-35.0
+ * with y = x - 2.
+ */
+s0 = -7.72156641e-02, /* 0xbd9e233f */
+s1 =  2.69987404e-01, /* 0x3e8a3bca */
+s2 =  1.42851010e-01, /* 0x3e124789 */
+s3 =  1.19389519e-02, /* 0x3c439b98 */
+r1 =  6.79650068e-01, /* 0x3f2dfd8c */
+r2 =  1.16058730e-01, /* 0x3dedb033 */
+r3 =  3.75673687e-03, /* 0x3b763396 */
+/*
+ * Domain z in [8, 0x1p24], range ~[-1.2640e-09, 1.2640e-09]:
+ * |lgamma(x) - (x - 0.5) * (log(x) - 1) - w(1/x)| < 2**-29.6.
+ */
+w0 =  4.18938547e-01, /* 0x3ed67f1d */
+w1 =  8.33332464e-02, /* 0x3daaaa9f */
+w2 = -2.76129087e-03; /* 0xbb34f6c6 */
+
+static float
+sin_pif(float x)
 {
+       volatile float vz;
        float y,z;
-       int n,ix;
-
-       GET_FLOAT_WORD(ix,x);
-       ix &= 0x7fffffff;
-
-       if(ix<0x3e800000) return __kernel_sindf(pi*x);
-       y = -x;         /* x is assume negative */
-
-    /*
-     * argument reduction, make sure inexact flag not raised if input
-     * is an integer
-     */
-       z = floorf(y);
-       if(z!=y) {                              /* inexact anyway */
-           y  *= (float)0.5;
-           y   = (float)2.0*(y - floorf(y));   /* y = |x| mod 2.0 */
-           n   = (int) (y*(float)4.0);
-       } else {
-            if(ix>=0x4b800000) {
-                y = zero; n = 0;                 /* y must be even */
-            } else {
-                if(ix<0x4b000000) z = y+two23; /* exact */
-               GET_FLOAT_WORD(n,z);
-               n &= 1;
-                y  = n;
-                n<<= 2;
-            }
-        }
+       int n;
+
+       y = -x;
+
+       vz = y+0x1p23F;                 /* depend on 0 <= y < 0x1p23 */
+       z = vz-0x1p23F;                 /* rintf(y) for the above range */
+       if (z == y)
+           return zero;
+
+       vz = y+0x1p21F;
+       GET_FLOAT_WORD(n,vz);           /* bits for rounded y (units 0.25) */
+       z = vz-0x1p21F;                 /* y rounded to a multiple of 0.25 */
+       if (z > y) {
+           z -= 0.25F;                 /* adjust to round down */
+           n--;
+       }
+       n &= 7;                         /* octant of y mod 2 */
+       y = y - z + n * 0.25F;          /* y mod 2 */
+
        switch (n) {
            case 0:   y =  __kernel_sindf(pi*y); break;
            case 1:
@@ -136,34 +120,38 @@ static const float zero=  0.0000000000e+00;
 float
 __ieee754_lgammaf_r(float x, int *signgamp)
 {
-       float t,y,z,nadj,p,p1,p2,p3,q,r,w;
+       float nadj,p,p1,p2,p3,q,r,t,w,y,z;
        int32_t hx;
        int i,ix;
 
        GET_FLOAT_WORD(hx,x);
 
-    /* purge off +-inf, NaN, +-0, tiny and negative arguments */
+    /* purge +-Inf and NaNs */
        *signgamp = 1;
        ix = hx&0x7fffffff;
        if(ix>=0x7f800000) return x*x;
-       if(ix==0) return one/zero;
-       if(ix<0x35000000) {     /* |x|<2**-21, return -log(|x|) */
-           if(hx<0) {
-               *signgamp = -1;
-               return -__ieee754_logf(-x);
-           } else return -__ieee754_logf(x);
+
+    /* purge +-0 and tiny arguments */
+       *signgamp = 1-2*((uint32_t)hx>>31);
+       if(ix<0x32000000) {             /* |x|<2**-27, return -log(|x|) */
+           if(ix==0)
+               return one/vzero;
+           return -__ieee754_logf(fabsf(x));
        }
+
+    /* purge negative integers and start evaluation for other x < 0 */
        if(hx<0) {
-           if(ix>=0x4b000000)  /* |x|>=2**23, must be -integer */
-               return one/zero;
+           *signgamp = 1;
+           if(ix>=0x4b000000)          /* |x|>=2**23, must be -integer */
+               return one/vzero;
            t = sin_pif(x);
-           if(t==zero) return one/zero; /* -integer */
+           if(t==zero) return one/vzero; /* -integer */
            nadj = __ieee754_logf(pi/fabsf(t*x));
            if(t<zero) *signgamp = -1;
            x = -x;
        }
 
-    /* purge off 1 and 2 */
+    /* purge 1 and 2 */
        if (ix==0x3f800000||ix==0x40000000) r = 0;
     /* for x < 2.0 */
        else if(ix<0x40000000) {
@@ -174,55 +162,51 @@ __ieee754_lgammaf_r(float x, int *signgamp)
                else {y = x; i=2;}
            } else {
                r = zero;
-               if(ix>=0x3fdda618) {y=(float)2.0-x;i=0;} /* [1.7316,2] */
+               if(ix>=0x3fdda618) {y=2-x;i=0;} /* [1.7316,2] */
                else if(ix>=0x3F9da620) {y=x-tc;i=1;} /* [1.23,1.73] */
                else {y=x-one;i=2;}
            }
            switch(i) {
              case 0:
                z = y*y;
-               p1 = a0+z*(a2+z*(a4+z*(a6+z*(a8+z*a10))));
-               p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*a11)))));
+               p1 = a0+z*(a2+z*a4);
+               p2 = z*(a1+z*(a3+z*a5));
                p  = y*p1+p2;
-               r  += (p-(float)0.5*y); break;
+               r  += p-y/2; break;
              case 1:
-               z = y*y;
-               w = z*y;
-               p1 = t0+w*(t3+w*(t6+w*(t9 +w*t12)));    /* parallel comp */
-               p2 = t1+w*(t4+w*(t7+w*(t10+w*t13)));
-               p3 = t2+w*(t5+w*(t8+w*(t11+w*t14)));
-               p  = z*p1-(tt-w*(p2+y*p3));
-               r += (tf + p); break;
+               p = t0+y*t1+y*y*(t2+y*(t3+y*(t4+y*(t5+y*(t6+y*t7)))));
+               r += tf + p; break;
              case 2:
-               p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*u5)))));
-               p2 = one+y*(v1+y*(v2+y*(v3+y*(v4+y*v5))));
-               r += (-(float)0.5*y + p1/p2);
+               p1 = y*(u0+y*(u1+y*u2));
+               p2 = one+y*(v1+y*(v2+y*v3));
+               r += p1/p2-y/2;
            }
        }
-       else if(ix<0x41000000) {                        /* x < 8.0 */
-           i = (int)x;
-           y = x-(float)i;
-           p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6))))));
-           q = one+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))));
-           r = half*y+p/q;
+    /* x < 8.0 */
+       else if(ix<0x41000000) {
+           i = x;
+           y = x-i;
+           p = y*(s0+y*(s1+y*(s2+y*s3)));
+           q = one+y*(r1+y*(r2+y*r3));
+           r = y/2+p/q;
            z = one;    /* lgamma(1+s) = log(s) + lgamma(s) */
            switch(i) {
-           case 7: z *= (y+(float)6.0);        /* FALLTHRU */
-           case 6: z *= (y+(float)5.0);        /* FALLTHRU */
-           case 5: z *= (y+(float)4.0);        /* FALLTHRU */
-           case 4: z *= (y+(float)3.0);        /* FALLTHRU */
-           case 3: z *= (y+(float)2.0);        /* FALLTHRU */
+           case 7: z *= (y+6);         /* FALLTHRU */
+           case 6: z *= (y+5);         /* FALLTHRU */
+           case 5: z *= (y+4);         /* FALLTHRU */
+           case 4: z *= (y+3);         /* FALLTHRU */
+           case 3: z *= (y+2);         /* FALLTHRU */
                    r += __ieee754_logf(z); break;
            }
-    /* 8.0 <= x < 2**58 */
-       } else if (ix < 0x5c800000) {
+    /* 8.0 <= x < 2**27 */
+       } else if (ix < 0x4d000000) {
            t = __ieee754_logf(x);
            z = one/x;
            y = z*z;
-           w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6)))));
+           w = w0+z*(w1+y*w2);
            r = (x-half)*(t-one)+w;
        } else
-    /* 2**58 <= x <= inf */
+    /* 2**27 <= x <= inf */
            r =  x*(__ieee754_logf(x)-one);
        if(hx<0) r = nadj - r;
        return r;
similarity index 56%
copy from lib/libm/src/e_lgamma.c
copy to lib/libm/src/e_lgammal.c
index 0259a38..4632724 100644 (file)
@@ -1,31 +1,23 @@
-
 /* @(#)e_lgamma.c 1.3 95/01/18 */
-/* $FreeBSD: head/lib/msun/src/e_lgamma.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 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
- *
  */
 
-/* __ieee754_lgamma(x)
- * Return the logarithm of the Gamma function of x.
- *
- * Method: call __ieee754_lgamma_r
- */
 
 #include "math.h"
 #include "math_private.h"
 
 extern int signgam;
 
-double
-__ieee754_lgamma(double x)
+long double
+lgammal(long double x)
 {
-       return __ieee754_lgamma_r(x,&signgam);
+       return lgammal_r(x,&signgam);
 }
index 2e6189c..43f96e2 100644 (file)
  *     1. Compute and return log2(x) in two pieces:
  *             log2(x) = w1 + w2,
  *        where w1 has 53-24 = 29 bit trailing zeros.
- *     2. Perform y*log2(x) = n+y' by simulating muti-precision 
+ *     2. Perform y*log2(x) = n+y' by simulating multi-precision 
  *        arithmetic, where |y'|<=0.5.
  *     3. Return x**y = 2**n*exp(y'*log2)
  *
  * Special cases:
  *     1.  (anything) ** 0  is 1
  *     2.  (anything) ** 1  is itself
- *     3.  (anything) ** NAN is NAN
+ *     3.  (anything) ** NAN is NAN except 1 ** NAN = 1
  *     4.  NAN ** (anything except 0) is NAN
  *     5.  +-(|x| > 1) **  +INF is +INF
  *     6.  +-(|x| > 1) **  -INF is +0
  *     7.  +-(|x| < 1) **  +INF is +0
  *     8.  +-(|x| < 1) **  -INF is +INF
- *     9.  +-1         ** +-INF is NAN
+ *     9.  +-1         ** +-INF is 1
  *     10. +0 ** (+anything except 0, NAN)               is +0
  *     11. -0 ** (+anything except 0, NAN, odd integer)  is +0
  *     12. +0 ** (-anything except 0, NAN)               is +INF
@@ -139,7 +139,7 @@ __ieee754_pow(double x, double y)
        if(ly==0) {     
            if (iy==0x7ff00000) {       /* y is +-inf */
                if(((ix-0x3ff00000)|lx)==0)
-                   return  one;        /* (-1)**+-inf is NaN */
+                   return  one;        /* (-1)**+-inf is 1 */
                else if (ix >= 0x3ff00000)/* (|x|>1)**+-inf = inf,0 */
                    return (hy>=0)? y: zero;
                else                    /* (|x|<1)**-,+inf = inf,0 */
diff --git a/lib/libm/src/e_remainder.c b/lib/libm/src/e_remainder.c
new file mode 100644 (file)
index 0000000..27c55a0
--- /dev/null
@@ -0,0 +1,77 @@
+
+/* @(#)e_remainder.c 1.3 95/01/18 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+
+/* __ieee754_remainder(x,p)
+ * Return :                  
+ *     returns  x REM p  =  x - [x/p]*p as if in infinite 
+ *     precise arithmetic, where [x/p] is the (infinite bit) 
+ *     integer nearest x/p (in half way case choose the even one).
+ * Method : 
+ *     Based on fmod() return x-[x/p]chopped*p exactlp.
+ */
+
+#include <float.h>
+
+#include "math.h"
+#include "math_private.h"
+
+static const double zero = 0.0;
+
+
+double
+__ieee754_remainder(double x, double p)
+{
+       int32_t hx,hp;
+       u_int32_t sx,lx,lp;
+       double p_half;
+
+       EXTRACT_WORDS(hx,lx,x);
+       EXTRACT_WORDS(hp,lp,p);
+       sx = hx&0x80000000;
+       hp &= 0x7fffffff;
+       hx &= 0x7fffffff;
+
+    /* purge off exception values */
+       if((hp|lp)==0) return (x*p)/(x*p);      /* p = 0 */
+       if((hx>=0x7ff00000)||                   /* x not finite */
+         ((hp>=0x7ff00000)&&                   /* p is NaN */
+         (((hp-0x7ff00000)|lp)!=0)))
+           return ((long double)x*p)/((long double)x*p);
+
+
+       if (hp<=0x7fdfffff) x = __ieee754_fmod(x,p+p);  /* now x < 2p */
+       if (((hx-hp)|(lx-lp))==0) return zero*x;
+       x  = fabs(x);
+       p  = fabs(p);
+       if (hp<0x00200000) {
+           if(x+x>p) {
+               x-=p;
+               if(x+x>=p) x -= p;
+           }
+       } else {
+           p_half = 0.5*p;
+           if(x>p_half) {
+               x-=p;
+               if(x>=p_half) x -= p;
+           }
+       }
+       GET_HIGH_WORD(hx,x);
+       if ((hx&0x7fffffff)==0) hx = 0;
+       SET_HIGH_WORD(x,hx^sx);
+       return x;
+}
+
+#if LDBL_MANT_DIG == 53
+__weak_reference(remainder, remainderl);
+#endif
diff --git a/lib/libm/src/e_remainderf.c b/lib/libm/src/e_remainderf.c
new file mode 100644 (file)
index 0000000..6537d0a
--- /dev/null
@@ -0,0 +1,63 @@
+/* e_remainderf.c -- float version of e_remainder.c.
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ */
+
+/*
+ * ====================================================
+ * 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 "math.h"
+#include "math_private.h"
+
+static const float zero = 0.0;
+
+
+float
+__ieee754_remainderf(float x, float p)
+{
+       int32_t hx,hp;
+       u_int32_t sx;
+       float p_half;
+
+       GET_FLOAT_WORD(hx,x);
+       GET_FLOAT_WORD(hp,p);
+       sx = hx&0x80000000;
+       hp &= 0x7fffffff;
+       hx &= 0x7fffffff;
+
+    /* purge off exception values */
+       if(hp==0) return (x*p)/(x*p);           /* p = 0 */
+       if((hx>=0x7f800000)||                   /* x not finite */
+         ((hp>0x7f800000)))                    /* p is NaN */
+           return ((long double)x*p)/((long double)x*p);
+
+
+       if (hp<=0x7effffff) x = __ieee754_fmodf(x,p+p); /* now x < 2p */
+       if ((hx-hp)==0) return zero*x;
+       x  = fabsf(x);
+       p  = fabsf(p);
+       if (hp<0x01000000) {
+           if(x+x>p) {
+               x-=p;
+               if(x+x>=p) x -= p;
+           }
+       } else {
+           p_half = (float)0.5*p;
+           if(x>p_half) {
+               x-=p;
+               if(x>=p_half) x -= p;
+           }
+       }
+       GET_FLOAT_WORD(hx,x);
+       if ((hx&0x7fffffff)==0) hx = 0;
+       SET_FLOAT_WORD(x,hx^sx);
+       return x;
+}
similarity index 56%
copy from lib/libm/src/imprecise.c
copy to lib/libm/src/e_remainderl.c
index 2498be3..4c25496 100644 (file)
@@ -1,5 +1,5 @@
 /*-
- * Copyright (c) 2013 David Chisnall
+ * Copyright (c) 2008 David Schultz <das@FreeBSD.ORG>
  * All rights reserved.
  *
  * Redistribution and use in source and binary forms, with or without
  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  * SUCH DAMAGE.
  *
- * $FreeBSD: head/lib/msun/src/imprecise.c 255294 2013-09-06 07:58:23Z theraven $
+ * $FreeBSD: head/lib/msun/src/e_remainderl.c 177765 2008-03-30 20:47:42Z das $
  */
 
-#include <float.h>
-#include <math.h>
 
-/*
- * If long double is not the same size as double, then these will lose
- * precision and we should emit a warning whenever something links against
- * them.
- */
-#if (LDBL_MANT_DIG > 53)
-#define WARN_IMPRECISE(x) \
-       __warn_references(x, # x " has lower than advertised precision");
-#else
-#define WARN_IMPRECISE(x)
-#endif
-/*
- * Declare the functions as weak variants so that other libraries providing
- * real versions can override them.
- */
-#define        DECLARE_WEAK(x)\
-       __weak_reference(imprecise_## x, x);\
-       WARN_IMPRECISE(x)
+#include <math.h>
 
 long double
-imprecise_powl(long double x, long double y)
+remainderl(long double x, long double y)
 {
+       int quo;
 
-       return pow(x, y);
+       return (remquol(x, y, &quo));
 }
-DECLARE_WEAK(powl);
-
-#define DECLARE_IMPRECISE(f) \
-       long double imprecise_ ## f ## l(long double v) { return f(v); }\
-       DECLARE_WEAK(f ## l)
-
-DECLARE_IMPRECISE(cosh);
-DECLARE_IMPRECISE(erfc);
-DECLARE_IMPRECISE(erf);
-DECLARE_IMPRECISE(lgamma);
-DECLARE_IMPRECISE(sinh);
-DECLARE_IMPRECISE(tanh);
-DECLARE_IMPRECISE(tgamma);
index 0e46787..34b601a 100644 (file)
@@ -1,6 +1,5 @@
 
 /* @(#)e_sinh.c 1.3 95/01/18 */
-/* $FreeBSD: head/lib/msun/src/e_sinh.c 226598 2011-10-21 06:28:47Z das $ */
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
@@ -12,6 +11,7 @@
  * ====================================================
  */
 
+
 /* __ieee754_sinh(x)
  * Method : 
  * mathematically sinh(x) if defined to be (exp(x)-exp(-x))/2
@@ -30,6 +30,8 @@
  *     only sinh(0)=0 is exact for finite x.
  */
 
+#include <float.h>
+
 #include "math.h"
 #include "math_private.h"
 
@@ -69,3 +71,7 @@ __ieee754_sinh(double x)
     /* |x| > overflowthresold, sinh(x) overflow */
        return x*shuge;
 }
+
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(sinh, sinhl);
+#endif
diff --git a/lib/libm/src/e_sinhl.c b/lib/libm/src/e_sinhl.c
new file mode 100644 (file)
index 0000000..866b3cc
--- /dev/null
@@ -0,0 +1,129 @@
+/* from: FreeBSD: head/lib/msun/src/e_sinhl.c XXX */
+
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+
+/*
+ * See e_sinh.c for complete comments.
+ *
+ * Converted to long double by Bruce D. Evans.
+ */
+
+#include <float.h>
+#ifdef __i386__
+#include <ieeefp.h>
+#endif
+
+#include "fpmath.h"
+#include "math.h"
+#include "math_private.h"
+#include "k_expl.h"
+
+#if LDBL_MAX_EXP != 0x4000
+/* We also require the usual expsign encoding. */
+#error "Unsupported long double format"
+#endif
+
+#define        BIAS    (LDBL_MAX_EXP - 1)
+
+static const long double shuge = 0x1p16383L;
+#if LDBL_MANT_DIG == 64
+/*
+ * Domain [-1, 1], range ~[-6.6749e-22, 6.6749e-22]:
+ * |sinh(x)/x - s(x)| < 2**-70.3
+ */
+static const union IEEEl2bits
+S3u = LD80C(0xaaaaaaaaaaaaaaaa, -3,  1.66666666666666666658e-1L);
+#define        S3      S3u.e
+static const double
+S5  =  8.3333333333333332e-3,          /*  0x11111111111111.0p-59 */
+S7  =  1.9841269841270074e-4,          /*  0x1a01a01a01a070.0p-65 */
+S9  =  2.7557319223873889e-6,          /*  0x171de3a5565fe6.0p-71 */
+S11 =  2.5052108406704084e-8,          /*  0x1ae6456857530f.0p-78 */
+S13 =  1.6059042748655297e-10,         /*  0x161245fa910697.0p-85 */
+S15 =  7.6470006914396920e-13,         /*  0x1ae7ce4eff2792.0p-93 */
+S17 =  2.8346142308424267e-15;         /*  0x19882ce789ffc6.0p-101 */
+#elif LDBL_MANT_DIG == 113
+/*
+ * Domain [-1, 1], range ~[-2.9673e-36, 2.9673e-36]:
+ * |sinh(x)/x - s(x)| < 2**-118.0
+ */
+static const long double
+S3  =  1.66666666666666666666666666666666033e-1L,      /*  0x1555555555555555555555555553b.0p-115L */
+S5  =  8.33333333333333333333333333337643193e-3L,      /*  0x111111111111111111111111180f5.0p-119L */
+S7  =  1.98412698412698412698412697391263199e-4L,      /*  0x1a01a01a01a01a01a01a0176aad11.0p-125L */
+S9  =  2.75573192239858906525574406205464218e-6L,      /*  0x171de3a556c7338faac243aaa9592.0p-131L */
+S11 =  2.50521083854417187749675637460977997e-8L,      /*  0x1ae64567f544e38fe59b3380d7413.0p-138L */
+S13 =  1.60590438368216146368737762431552702e-10L,     /*  0x16124613a86d098059c7620850fc2.0p-145L */
+S15 =  7.64716373181980539786802470969096440e-13L,     /*  0x1ae7f3e733b814193af09ce723043.0p-153L */
+S17 =  2.81145725434775409870584280722701574e-15L;     /*  0x1952c77030c36898c3fd0b6dfc562.0p-161L */
+static const double
+S19=  8.2206352435411005e-18,          /*  0x12f49b4662b86d.0p-109 */
+S21=  1.9572943931418891e-20,          /*  0x171b8f2fab9628.0p-118 */
+S23 =  3.8679983530666939e-23,         /*  0x17617002b73afc.0p-127 */
+S25 =  6.5067867911512749e-26;         /*  0x1423352626048a.0p-136 */
+#else
+#error "Unsupported long double format"
+#endif /* LDBL_MANT_DIG == 64 */
+
+/* log(2**16385 - 0.5) rounded up: */
+static const float
+o_threshold =  1.13572168e4;           /*  0xb174de.0p-10 */
+
+long double
+sinhl(long double x)
+{
+       long double hi,lo,x2,x4;
+       double dx2,s;
+       int16_t ix,jx;
+
+       GET_LDBL_EXPSIGN(jx,x);
+       ix = jx&0x7fff;
+
+    /* x is INF or NaN */
+       if(ix>=0x7fff) return x+x;
+
+       ENTERI();
+
+       s = 1;
+       if (jx<0) s = -1;
+
+    /* |x| < 64, return x, s(x), or accurate s*(exp(|x|)/2-1/exp(|x|)/2) */
+       if (ix<0x4005) {                /* |x|<64 */
+           if (ix<BIAS-(LDBL_MANT_DIG+1)/2)    /* |x|<TINY */
+               if(shuge+x>1) RETURNI(x);  /* sinh(tiny) = tiny with inexact */
+           if (ix<0x3fff) {            /* |x|<1 */
+               x2 = x*x;
+#if LDBL_MANT_DIG == 64
+               x4 = x2*x2;
+               RETURNI(((S17*x2 + S15)*x4 + (S13*x2 + S11))*(x2*x*x4*x4) +
+                   ((S9*x2 + S7)*x2 + S5)*(x2*x*x2) + S3*(x2*x) + x);
+#elif LDBL_MANT_DIG == 113
+               dx2 = x2;
+               RETURNI(((((((((((S25*dx2 + S23)*dx2 +
+                   S21)*x2 + S19)*x2 +
+                   S17)*x2 + S15)*x2 + S13)*x2 + S11)*x2 + S9)*x2 + S7)*x2 +
+                   S5)* (x2*x*x2) +
+                   S3*(x2*x) + x);
+#endif
+           }
+           k_hexpl(fabsl(x), &hi, &lo);
+           RETURNI(s*(lo - 0.25/(hi + lo) + hi));
+       }
+
+    /* |x| in [64, o_threshold], return correctly-overflowing s*exp(|x|)/2 */
+       if (fabsl(x) <= o_threshold)
+           RETURNI(s*hexpl(fabsl(x)));
+
+    /* |x| > o_threshold, sinh(x) overflow */
+       return x*shuge;
+}
diff --git a/lib/libm/src/e_sqrt.c b/lib/libm/src/e_sqrt.c
new file mode 100644 (file)
index 0000000..5e4d540
--- /dev/null
@@ -0,0 +1,449 @@
+
+/* @(#)e_sqrt.c 1.3 95/01/18 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+
+/* __ieee754_sqrt(x)
+ * Return correctly rounded sqrt.
+ *           ------------------------------------------
+ *          |  Use the hardware sqrt if you have one |
+ *           ------------------------------------------
+ * Method: 
+ *   Bit by bit method using integer arithmetic. (Slow, but portable) 
+ *   1. Normalization
+ *     Scale x to y in [1,4) with even powers of 2: 
+ *     find an integer k such that  1 <= (y=x*2^(2k)) < 4, then
+ *             sqrt(x) = 2^k * sqrt(y)
+ *   2. Bit by bit computation
+ *     Let q  = sqrt(y) truncated to i bit after binary point (q = 1),
+ *          i                                                   0
+ *                                     i+1         2
+ *         s  = 2*q , and      y  =  2   * ( y - q  ).         (1)
+ *          i      i            i                 i
+ *                                                        
+ *     To compute q    from q , one checks whether 
+ *                 i+1       i                       
+ *
+ *                           -(i+1) 2
+ *                     (q + 2      ) <= y.                     (2)
+ *                               i
+ *                                                           -(i+1)
+ *     If (2) is false, then q   = q ; otherwise q   = q  + 2      .
+ *                            i+1   i             i+1   i
+ *
+ *     With some algebric manipulation, it is not difficult to see
+ *     that (2) is equivalent to 
+ *                             -(i+1)
+ *                     s  +  2       <= y                      (3)
+ *                      i                i
+ *
+ *     The advantage of (3) is that s  and y  can be computed by 
+ *                                   i      i
+ *     the following recurrence formula:
+ *         if (3) is false
+ *
+ *         s     =  s  ,       y    = y   ;                    (4)
+ *          i+1      i          i+1    i
+ *
+ *         otherwise,
+ *                         -i                     -(i+1)
+ *         s     =  s  + 2  ,  y    = y  -  s  - 2             (5)
+ *           i+1      i          i+1    i     i
+ *                             
+ *     One may easily use induction to prove (4) and (5). 
+ *     Note. Since the left hand side of (3) contain only i+2 bits,
+ *           it does not necessary to do a full (53-bit) comparison 
+ *           in (3).
+ *   3. Final rounding
+ *     After generating the 53 bits result, we compute one more bit.
+ *     Together with the remainder, we can decide whether the
+ *     result is exact, bigger than 1/2ulp, or less than 1/2ulp
+ *     (it will never equal to 1/2ulp).
+ *     The rounding mode can be detected by checking whether
+ *     huge + tiny is equal to huge, and whether huge - tiny is
+ *     equal to huge for some floating point number "huge" and "tiny".
+ *             
+ * Special cases:
+ *     sqrt(+-0) = +-0         ... exact
+ *     sqrt(inf) = inf
+ *     sqrt(-ve) = NaN         ... with invalid signal
+ *     sqrt(NaN) = NaN         ... with invalid signal for signaling NaN
+ *
+ * Other methods : see the appended file at the end of the program below.
+ *---------------
+ */
+
+#include <float.h>
+
+#include "math.h"
+#include "math_private.h"
+
+static const double    one     = 1.0, tiny=1.0e-300;
+
+double
+__ieee754_sqrt(double x)
+{
+       double z;
+       int32_t sign = (int)0x80000000;
+       int32_t ix0,s0,q,m,t,i;
+       u_int32_t r,t1,s1,ix1,q1;
+
+       EXTRACT_WORDS(ix0,ix1,x);
+
+    /* take care of Inf and NaN */
+       if((ix0&0x7ff00000)==0x7ff00000) {                      
+           return x*x+x;               /* sqrt(NaN)=NaN, sqrt(+inf)=+inf
+                                          sqrt(-inf)=sNaN */
+       } 
+    /* take care of zero */
+       if(ix0<=0) {
+           if(((ix0&(~sign))|ix1)==0) return x;/* sqrt(+-0) = +-0 */
+           else if(ix0<0)
+               return (x-x)/(x-x);             /* sqrt(-ve) = sNaN */
+       }
+    /* normalize x */
+       m = (ix0>>20);
+       if(m==0) {                              /* subnormal x */
+           while(ix0==0) {
+               m -= 21;
+               ix0 |= (ix1>>11); ix1 <<= 21;
+           }
+           for(i=0;(ix0&0x00100000)==0;i++) ix0<<=1;
+           m -= i-1;
+           ix0 |= (ix1>>(32-i));
+           ix1 <<= i;
+       }
+       m -= 1023;      /* unbias exponent */
+       ix0 = (ix0&0x000fffff)|0x00100000;
+       if(m&1){        /* odd m, double x to make it even */
+           ix0 += ix0 + ((ix1&sign)>>31);
+           ix1 += ix1;
+       }
+       m >>= 1;        /* m = [m/2] */
+
+    /* generate sqrt(x) bit by bit */
+       ix0 += ix0 + ((ix1&sign)>>31);
+       ix1 += ix1;
+       q = q1 = s0 = s1 = 0;   /* [q,q1] = sqrt(x) */
+       r = 0x00200000;         /* r = moving bit from right to left */
+
+       while(r!=0) {
+           t = s0+r; 
+           if(t<=ix0) { 
+               s0   = t+r; 
+               ix0 -= t; 
+               q   += r; 
+           } 
+           ix0 += ix0 + ((ix1&sign)>>31);
+           ix1 += ix1;
+           r>>=1;
+       }
+
+       r = sign;
+       while(r!=0) {
+           t1 = s1+r; 
+           t  = s0;
+           if((t<ix0)||((t==ix0)&&(t1<=ix1))) { 
+               s1  = t1+r;
+               if(((t1&sign)==sign)&&(s1&sign)==0) s0 += 1;
+               ix0 -= t;
+               if (ix1 < t1) ix0 -= 1;
+               ix1 -= t1;
+               q1  += r;
+           }
+           ix0 += ix0 + ((ix1&sign)>>31);
+           ix1 += ix1;
+           r>>=1;
+       }
+
+    /* use floating add to find out rounding direction */
+       if((ix0|ix1)!=0) {
+           z = one-tiny; /* trigger inexact flag */
+           if (z>=one) {
+               z = one+tiny;
+               if (q1==(u_int32_t)0xffffffff) { q1=0; q += 1;}
+               else if (z>one) {
+                   if (q1==(u_int32_t)0xfffffffe) q+=1;
+                   q1+=2; 
+               } else
+                   q1 += (q1&1);
+           }
+       }
+       ix0 = (q>>1)+0x3fe00000;
+       ix1 =  q1>>1;
+       if ((q&1)==1) ix1 |= sign;
+       ix0 += (m <<20);
+       INSERT_WORDS(z,ix0,ix1);
+       return z;
+}
+
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(sqrt, sqrtl);
+#endif
+
+/*
+Other methods  (use floating-point arithmetic)
+-------------
+(This is a copy of a drafted paper by Prof W. Kahan 
+and K.C. Ng, written in May, 1986)
+
+       Two algorithms are given here to implement sqrt(x) 
+       (IEEE double precision arithmetic) in software.
+       Both supply sqrt(x) correctly rounded. The first algorithm (in
+       Section A) uses newton iterations and involves four divisions.
+       The second one uses reciproot iterations to avoid division, but
+       requires more multiplications. Both algorithms need the ability
+       to chop results of arithmetic operations instead of round them, 
+       and the INEXACT flag to indicate when an arithmetic operation
+       is executed exactly with no roundoff error, all part of the 
+       standard (IEEE 754-1985). The ability to perform shift, add,
+       subtract and logical AND operations upon 32-bit words is needed
+       too, though not part of the standard.
+
+A.  sqrt(x) by Newton Iteration
+
+   (1) Initial approximation
+
+       Let x0 and x1 be the leading and the trailing 32-bit words of
+       a floating point number x (in IEEE double format) respectively 
+
+           1    11                  52                           ...widths
+          ------------------------------------------------------
+       x: |s|    e     |             f                         |
+          ------------------------------------------------------
+             msb    lsb  msb                                 lsb ...order
+
+            ------------------------        ------------------------
+       x0:  |s|   e    |    f1     |    x1: |          f2           |
+            ------------------------        ------------------------
+
+       By performing shifts and subtracts on x0 and x1 (both regarded
+       as integers), we obtain an 8-bit approximation of sqrt(x) as
+       follows.
+
+               k  := (x0>>1) + 0x1ff80000;
+               y0 := k - T1[31&(k>>15)].       ... y ~ sqrt(x) to 8 bits
+       Here k is a 32-bit integer and T1[] is an integer array containing
+       correction terms. Now magically the floating value of y (y's
+       leading 32-bit word is y0, the value of its trailing word is 0)
+       approximates sqrt(x) to almost 8-bit.
+
+       Value of T1:
+       static int T1[32]= {
+       0,      1024,   3062,   5746,   9193,   13348,  18162,  23592,
+       29598,  36145,  43202,  50740,  58733,  67158,  75992,  85215,
+       83599,  71378,  60428,  50647,  41945,  34246,  27478,  21581,
+       16499,  12183,  8588,   5674,   3403,   1742,   661,    130,};
+
+    (2)        Iterative refinement
+
+       Apply Heron's rule three times to y, we have y approximates 
+       sqrt(x) to within 1 ulp (Unit in the Last Place):
+
+               y := (y+x/y)/2          ... almost 17 sig. bits
+               y := (y+x/y)/2          ... almost 35 sig. bits
+               y := y-(y-x/y)/2        ... within 1 ulp
+
+
+       Remark 1.
+           Another way to improve y to within 1 ulp is:
+
+               y := (y+x/y)            ... almost 17 sig. bits to 2*sqrt(x)
+               y := y - 0x00100006     ... almost 18 sig. bits to sqrt(x)
+
+                               2
+                           (x-y )*y
+               y := y + 2* ----------  ...within 1 ulp
+                              2
+                            3y  + x
+
+
+       This formula has one division fewer than the one above; however,
+       it requires more multiplications and additions. Also x must be
+       scaled in advance to avoid spurious overflow in evaluating the
+       expression 3y*y+x. Hence it is not recommended uless division
+       is slow. If division is very slow, then one should use the 
+       reciproot algorithm given in section B.
+
+    (3) Final adjustment
+
+       By twiddling y's last bit it is possible to force y to be 
+       correctly rounded according to the prevailing rounding mode
+       as follows. Let r and i be copies of the rounding mode and
+       inexact flag before entering the square root program. Also we
+       use the expression y+-ulp for the next representable floating
+       numbers (up and down) of y. Note that y+-ulp = either fixed
+       point y+-1, or multiply y by nextafter(1,+-inf) in chopped
+       mode.
+
+               I := FALSE;     ... reset INEXACT flag I
+               R := RZ;        ... set rounding mode to round-toward-zero
+               z := x/y;       ... chopped quotient, possibly inexact
+               If(not I) then {        ... if the quotient is exact
+                   if(z=y) {
+                       I := i;  ... restore inexact flag
+                       R := r;  ... restore rounded mode
+                       return sqrt(x):=y.
+                   } else {
+                       z := z - ulp;   ... special rounding
+                   }
+               }
+               i := TRUE;              ... sqrt(x) is inexact
+               If (r=RN) then z=z+ulp  ... rounded-to-nearest
+               If (r=RP) then {        ... round-toward-+inf
+                   y = y+ulp; z=z+ulp;
+               }
+               y := y+z;               ... chopped sum
+               y0:=y0-0x00100000;      ... y := y/2 is correctly rounded.
+               I := i;                 ... restore inexact flag
+               R := r;                 ... restore rounded mode
+               return sqrt(x):=y.
+                   
+    (4)        Special cases
+
+       Square root of +inf, +-0, or NaN is itself;
+       Square root of a negative number is NaN with invalid signal.
+
+
+B.  sqrt(x) by Reciproot Iteration
+
+   (1) Initial approximation
+
+       Let x0 and x1 be the leading and the trailing 32-bit words of
+       a floating point number x (in IEEE double format) respectively
+       (see section A). By performing shifs and subtracts on x0 and y0,
+       we obtain a 7.8-bit approximation of 1/sqrt(x) as follows.
+
+           k := 0x5fe80000 - (x0>>1);
+           y0:= k - T2[63&(k>>14)].    ... y ~ 1/sqrt(x) to 7.8 bits
+
+       Here k is a 32-bit integer and T2[] is an integer array 
+       containing correction terms. Now magically the floating
+       value of y (y's leading 32-bit word is y0, the value of
+       its trailing word y1 is set to zero) approximates 1/sqrt(x)
+       to almost 7.8-bit.
+
+       Value of T2:
+       static int T2[64]= {
+       0x1500, 0x2ef8, 0x4d67, 0x6b02, 0x87be, 0xa395, 0xbe7a, 0xd866,
+       0xf14a, 0x1091b,0x11fcd,0x13552,0x14999,0x15c98,0x16e34,0x17e5f,
+       0x18d03,0x19a01,0x1a545,0x1ae8a,0x1b5c4,0x1bb01,0x1bfde,0x1c28d,
+       0x1c2de,0x1c0db,0x1ba73,0x1b11c,0x1a4b5,0x1953d,0x18266,0x16be0,
+       0x1683e,0x179d8,0x18a4d,0x19992,0x1a789,0x1b445,0x1bf61,0x1c989,
+       0x1d16d,0x1d77b,0x1dddf,0x1e2ad,0x1e5bf,0x1e6e8,0x1e654,0x1e3cd,
+       0x1df2a,0x1d635,0x1cb16,0x1be2c,0x1ae4e,0x19bde,0x1868e,0x16e2e,
+       0x1527f,0x1334a,0x11051,0xe951, 0xbe01, 0x8e0d, 0x5924, 0x1edd,};
+
+    (2)        Iterative refinement
+
+       Apply Reciproot iteration three times to y and multiply the
+       result by x to get an approximation z that matches sqrt(x)
+       to about 1 ulp. To be exact, we will have 
+               -1ulp < sqrt(x)-z<1.0625ulp.
+       
+       ... set rounding mode to Round-to-nearest
+          y := y*(1.5-0.5*x*y*y)       ... almost 15 sig. bits to 1/sqrt(x)
+          y := y*((1.5-2^-30)+0.5*x*y*y)... about 29 sig. bits to 1/sqrt(x)
+       ... special arrangement for better accuracy
+          z := x*y                     ... 29 bits to sqrt(x), with z*y<1
+          z := z + 0.5*z*(1-z*y)       ... about 1 ulp to sqrt(x)
+
+       Remark 2. The constant 1.5-2^-30 is chosen to bias the error so that
+       (a) the term z*y in the final iteration is always less than 1; 
+       (b) the error in the final result is biased upward so that
+               -1 ulp < sqrt(x) - z < 1.0625 ulp
+           instead of |sqrt(x)-z|<1.03125ulp.
+
+    (3)        Final adjustment
+
+       By twiddling y's last bit it is possible to force y to be 
+       correctly rounded according to the prevailing rounding mode
+       as follows. Let r and i be copies of the rounding mode and
+       inexact flag before entering the square root program. Also we
+       use the expression y+-ulp for the next representable floating
+       numbers (up and down) of y. Note that y+-ulp = either fixed
+       point y+-1, or multiply y by nextafter(1,+-inf) in chopped
+       mode.
+
+       R := RZ;                ... set rounding mode to round-toward-zero
+       switch(r) {
+           case RN:            ... round-to-nearest
+              if(x<= z*(z-ulp)...chopped) z = z - ulp; else
+              if(x<= z*(z+ulp)...chopped) z = z; else z = z+ulp;
+              break;
+           case RZ:case RM:    ... round-to-zero or round-to--inf
+              R:=RP;           ... reset rounding mod to round-to-+inf
+              if(x<z*z ... rounded up) z = z - ulp; else
+              if(x>=(z+ulp)*(z+ulp) ...rounded up) z = z+ulp;
+              break;
+           case RP:            ... round-to-+inf
+              if(x>(z+ulp)*(z+ulp)...chopped) z = z+2*ulp; else
+              if(x>z*z ...chopped) z = z+ulp;
+              break;
+       }
+
+       Remark 3. The above comparisons can be done in fixed point. For
+       example, to compare x and w=z*z chopped, it suffices to compare
+       x1 and w1 (the trailing parts of x and w), regarding them as
+       two's complement integers.
+
+       ...Is z an exact square root?
+       To determine whether z is an exact square root of x, let z1 be the
+       trailing part of z, and also let x0 and x1 be the leading and
+       trailing parts of x.
+
+       If ((z1&0x03ffffff)!=0) ... not exact if trailing 26 bits of z!=0
+           I := 1;             ... Raise Inexact flag: z is not exact
+       else {
+           j := 1 - [(x0>>20)&1]       ... j = logb(x) mod 2
+           k := z1 >> 26;              ... get z's 25-th and 26-th 
+                                           fraction bits
+           I := i or (k&j) or ((k&(j+j+1))!=(x1&3));
+       }
+       R:= r           ... restore rounded mode
+       return sqrt(x):=z.
+
+       If multiplication is cheaper then the foregoing red tape, the 
+       Inexact flag can be evaluated by
+
+           I := i;
+           I := (z*z!=x) or I.
+
+       Note that z*z can overwrite I; this value must be sensed if it is 
+       True.
+
+       Remark 4. If z*z = x exactly, then bit 25 to bit 0 of z1 must be
+       zero.
+
+                   --------------------
+               z1: |        f2        | 
+                   --------------------
+               bit 31             bit 0
+
+       Further more, bit 27 and 26 of z1, bit 0 and 1 of x1, and the odd
+       or even of logb(x) have the following relations:
+
+       -------------------------------------------------
+       bit 27,26 of z1         bit 1,0 of x1   logb(x)
+       -------------------------------------------------
+       00                      00              odd and even
+       01                      01              even
+       10                      10              odd
+       10                      00              even
+       11                      01              even
+       -------------------------------------------------
+
+    (4)        Special cases (see (4) of Section A).   
+ */
diff --git a/lib/libm/src/e_sqrtf.c b/lib/libm/src/e_sqrtf.c
new file mode 100644 (file)
index 0000000..40a2138
--- /dev/null
@@ -0,0 +1,88 @@
+/* e_sqrtf.c -- float version of e_sqrt.c.
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ */
+
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ *
+ * $FreeBSD: head/lib/msun/src/e_sqrtf.c 97413 2002-05-28 18:15:04Z alfred $
+ */
+
+
+#include "math.h"
+#include "math_private.h"
+
+static const float     one     = 1.0, tiny=1.0e-30;
+
+float
+__ieee754_sqrtf(float x)
+{
+       float z;
+       int32_t sign = (int)0x80000000;
+       int32_t ix,s,q,m,t,i;
+       u_int32_t r;
+
+       GET_FLOAT_WORD(ix,x);
+
+    /* take care of Inf and NaN */
+       if((ix&0x7f800000)==0x7f800000) {
+           return x*x+x;               /* sqrt(NaN)=NaN, sqrt(+inf)=+inf
+                                          sqrt(-inf)=sNaN */
+       }
+    /* take care of zero */
+       if(ix<=0) {
+           if((ix&(~sign))==0) return x;/* sqrt(+-0) = +-0 */
+           else if(ix<0)
+               return (x-x)/(x-x);             /* sqrt(-ve) = sNaN */
+       }
+    /* normalize x */
+       m = (ix>>23);
+       if(m==0) {                              /* subnormal x */
+           for(i=0;(ix&0x00800000)==0;i++) ix<<=1;
+           m -= i-1;
+       }
+       m -= 127;       /* unbias exponent */
+       ix = (ix&0x007fffff)|0x00800000;
+       if(m&1) /* odd m, double x to make it even */
+           ix += ix;
+       m >>= 1;        /* m = [m/2] */
+
+    /* generate sqrt(x) bit by bit */
+       ix += ix;
+       q = s = 0;              /* q = sqrt(x) */
+       r = 0x01000000;         /* r = moving bit from right to left */
+
+       while(r!=0) {
+           t = s+r;
+           if(t<=ix) {
+               s    = t+r;
+               ix  -= t;
+               q   += r;
+           }
+           ix += ix;
+           r>>=1;
+       }
+
+    /* use floating add to find out rounding direction */
+       if(ix!=0) {
+           z = one-tiny; /* trigger inexact flag */
+           if (z>=one) {
+               z = one+tiny;
+               if (z>one)
+                   q += 2;
+               else
+                   q += (q&1);
+           }
+       }
+       ix = (q>>1)+0x3f000000;
+       ix += (m <<23);
+       SET_FLOAT_WORD(z,ix);
+       return z;
+}
diff --git a/lib/libm/src/e_sqrtl.c b/lib/libm/src/e_sqrtl.c
new file mode 100644 (file)
index 0000000..047a581
--- /dev/null
@@ -0,0 +1,159 @@
+/*-
+ * Copyright (c) 2007 Steven G. Kargl
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice unmodified, this list of conditions, and the following
+ *    disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
+ * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
+ * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
+ * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
+ * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
+ * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+ * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+ * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
+ * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * $FreeBSD: head/lib/msun/src/e_sqrtl.c 176720 2008-03-02 01:47:58Z das $
+ */
+
+
+#include <fenv.h>
+#include <float.h>
+
+#include "fpmath.h"
+#include "math.h"
+
+/* Return (x + ulp) for normal positive x. Assumes no overflow. */
+static inline long double
+inc(long double x)
+{
+       union IEEEl2bits u;
+
+       u.e = x;
+       if (++u.bits.manl == 0) {
+               if (++u.bits.manh == 0) {
+                       u.bits.exp++;
+                       u.bits.manh |= LDBL_NBIT;
+               }
+       }
+       return (u.e);
+}
+
+/* Return (x - ulp) for normal positive x. Assumes no underflow. */
+static inline long double
+dec(long double x)
+{
+       union IEEEl2bits u;
+
+       u.e = x;
+       if (u.bits.manl-- == 0) {
+               if (u.bits.manh-- == LDBL_NBIT) {
+                       u.bits.exp--;
+                       u.bits.manh |= LDBL_NBIT;
+               }
+       }
+       return (u.e);
+}
+
+#pragma STDC FENV_ACCESS ON
+
+/*
+ * This is slow, but simple and portable. You should use hardware sqrt
+ * if possible.
+ */
+
+long double
+sqrtl(long double x)
+{
+       union IEEEl2bits u;
+       int k, r;
+       long double lo, xn;
+       fenv_t env;
+
+       u.e = x;
+
+       /* If x = NaN, then sqrt(x) = NaN. */
+       /* If x = Inf, then sqrt(x) = Inf. */
+       /* If x = -Inf, then sqrt(x) = NaN. */
+       if (u.bits.exp == LDBL_MAX_EXP * 2 - 1)
+               return (x * x + x);
+
+       /* If x = +-0, then sqrt(x) = +-0. */
+       if ((u.bits.manh | u.bits.manl | u.bits.exp) == 0)
+               return (x);
+
+       /* If x < 0, then raise invalid and return NaN */
+       if (u.bits.sign)
+               return ((x - x) / (x - x));
+
+       feholdexcept(&env);
+
+       if (u.bits.exp == 0) {
+               /* Adjust subnormal numbers. */
+               u.e *= 0x1.0p514;
+               k = -514;
+       } else {
+               k = 0;
+       }
+       /*
+        * u.e is a normal number, so break it into u.e = e*2^n where
+        * u.e = (2*e)*2^2k for odd n and u.e = (4*e)*2^2k for even n.
+        */
+       if ((u.bits.exp - 0x3ffe) & 1) {        /* n is odd.     */
+               k += u.bits.exp - 0x3fff;       /* 2k = n - 1.   */
+               u.bits.exp = 0x3fff;            /* u.e in [1,2). */
+       } else {
+               k += u.bits.exp - 0x4000;       /* 2k = n - 2.   */
+               u.bits.exp = 0x4000;            /* u.e in [2,4). */
+       }
+
+       /*
+        * Newton's iteration.
+        * Split u.e into a high and low part to achieve additional precision.
+        */
+       xn = sqrt(u.e);                 /* 53-bit estimate of sqrtl(x). */
+#if LDBL_MANT_DIG > 100
+       xn = (xn + (u.e / xn)) * 0.5;   /* 106-bit estimate. */
+#endif
+       lo = u.e;
+       u.bits.manl = 0;                /* Zero out lower bits. */
+       lo = (lo - u.e) / xn;           /* Low bits divided by xn. */
+       xn = xn + (u.e / xn);           /* High portion of estimate. */
+       u.e = xn + lo;                  /* Combine everything. */
+       u.bits.exp += (k >> 1) - 1;
+
+       feclearexcept(FE_INEXACT);
+       r = fegetround();
+       fesetround(FE_TOWARDZERO);      /* Set to round-toward-zero. */
+       xn = x / u.e;                   /* Chopped quotient (inexact?). */
+
+       if (!fetestexcept(FE_INEXACT)) { /* Quotient is exact. */
+               if (xn == u.e) {
+                       fesetenv(&env);
+                       return (u.e);
+               }
+               /* Round correctly for inputs like x = y**2 - ulp. */
+               xn = dec(xn);           /* xn = xn - ulp. */
+       }
+
+       if (r == FE_TONEAREST) {
+               xn = inc(xn);           /* xn = xn + ulp. */
+       } else if (r == FE_UPWARD) {
+               u.e = inc(u.e);         /* u.e = u.e + ulp. */
+               xn = inc(xn);           /* xn  = xn + ulp. */
+       }
+       u.e = u.e + xn;                         /* Chopped sum. */
+       feupdateenv(&env);      /* Restore env and raise inexact */
+       u.bits.exp--;
+       return (u.e);
+}
diff --git a/lib/libm/src/fenv-softfloat.h b/lib/libm/src/fenv-softfloat.h
new file mode 100644 (file)
index 0000000..48c1277
--- /dev/null
@@ -0,0 +1,184 @@
+/*-
+ * Copyright (c) 2004-2011 David Schultz <das@FreeBSD.ORG>
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ *
+ * $FreeBSD$
+ */
+
+#ifndef        _FENV_H_
+#error "This file is meant to be included only by <fenv.h>."
+#endif
+
+/*
+ * This file implements the functionality of <fenv.h> on platforms that
+ * lack an FPU and use softfloat in libc for floating point.  To use it,
+ * you must write an <fenv.h> that provides the following:
+ *
+ *   - a typedef for fenv_t, which may be an integer or struct type
+ *   - a typedef for fexcept_t (XXX This file assumes fexcept_t is a
+ *     simple integer type containing the exception mask.)
+ *   - definitions of FE_* constants for the five exceptions and four
+ *     rounding modes in IEEE 754, as described in fenv(3)
+ *   - a definition, and the corresponding external symbol, for FE_DFL_ENV
+ *   - a macro __set_env(env, flags, mask, rnd), which sets the given fenv_t
+ *     from the exception flags, mask, and rounding mode
+ *   - macros __env_flags(env), __env_mask(env), and __env_round(env), which
+ *     extract fields from an fenv_t
+ *   - a definition of __fenv_static
+ *
+ * If the architecture supports an optional FPU, it's recommended that you
+ * define fenv_t and fexcept_t to match the hardware ABI.  Otherwise, it
+ * doesn't matter how you define them.
+ */
+
+extern int __softfloat_float_exception_flags;
+extern int __softfloat_float_exception_mask;
+extern int __softfloat_float_rounding_mode;
+void __softfloat_float_raise(int);
+
+__fenv_static inline int
+feclearexcept(int __excepts)
+{
+
+       __softfloat_float_exception_flags &= ~__excepts;
+       return (0);
+}
+
+__fenv_static inline int
+fegetexceptflag(fexcept_t *__flagp, int __excepts)
+{
+
+       *__flagp = __softfloat_float_exception_flags & __excepts;
+       return (0);
+}
+
+__fenv_static inline int
+fesetexceptflag(const fexcept_t *__flagp, int __excepts)
+{
+
+       __softfloat_float_exception_flags &= ~__excepts;
+       __softfloat_float_exception_flags |= *__flagp & __excepts;
+       return (0);
+}
+
+__fenv_static inline int
+feraiseexcept(int __excepts)
+{
+
+       __softfloat_float_raise(__excepts);
+       return (0);
+}
+
+__fenv_static inline int
+fetestexcept(int __excepts)
+{
+
+       return (__softfloat_float_exception_flags & __excepts);
+}
+
+__fenv_static inline int
+fegetround(void)
+{
+
+       return (__softfloat_float_rounding_mode);
+}
+
+__fenv_static inline int
+fesetround(int __round)
+{
+
+       __softfloat_float_rounding_mode = __round;
+       return (0);
+}
+
+__fenv_static inline int
+fegetenv(fenv_t *__envp)
+{
+
+       __set_env(*__envp, __softfloat_float_exception_flags,
+           __softfloat_float_exception_mask, __softfloat_float_rounding_mode);
+       return (0);
+}
+
+__fenv_static inline int
+feholdexcept(fenv_t *__envp)
+{
+       fenv_t __env;
+
+       fegetenv(__envp);
+       __softfloat_float_exception_flags = 0;
+       __softfloat_float_exception_mask = 0;
+       return (0);
+}
+
+__fenv_static inline int
+fesetenv(const fenv_t *__envp)
+{
+
+       __softfloat_float_exception_flags = __env_flags(*__envp);
+       __softfloat_float_exception_mask = __env_mask(*__envp);
+       __softfloat_float_rounding_mode = __env_round(*__envp);
+       return (0);
+}
+
+__fenv_static inline int
+feupdateenv(const fenv_t *__envp)
+{
+       int __oflags = __softfloat_float_exception_flags;
+
+       fesetenv(__envp);
+       feraiseexcept(__oflags);
+       return (0);
+}
+
+#if __BSD_VISIBLE
+
+/* We currently provide no external definitions of the functions below. */
+
+__fenv_static inline int
+feenableexcept(int __mask)
+{
+       int __omask = __softfloat_float_exception_mask;
+
+       __softfloat_float_exception_mask |= __mask;
+       return (__omask);
+}
+
+__fenv_static inline int
+fedisableexcept(int __mask)
+{
+       int __omask = __softfloat_float_exception_mask;
+
+       __softfloat_float_exception_mask &= ~__mask;
+       return (__omask);
+}
+
+__fenv_static inline int
+fegetexcept(void)
+{
+
+       return (__softfloat_float_exception_mask);
+}
+
+#endif /* __BSD_VISIBLE */
index 2498be3..c8e7b47 100644 (file)
@@ -56,14 +56,18 @@ imprecise_powl(long double x, long double y)
 }
 DECLARE_WEAK(powl);
 
+#define DECLARE_FORMER_IMPRECISE(f) \
+       __asm__(".symver imprecise_" #f "l," #f "l@DF306.1"); \
+       long double imprecise_ ## f ## l(long double v) { return f(v); }
+
 #define DECLARE_IMPRECISE(f) \
        long double imprecise_ ## f ## l(long double v) { return f(v); }\
        DECLARE_WEAK(f ## l)
 
-DECLARE_IMPRECISE(cosh);
-DECLARE_IMPRECISE(erfc);
-DECLARE_IMPRECISE(erf);
 DECLARE_IMPRECISE(lgamma);
-DECLARE_IMPRECISE(sinh);
-DECLARE_IMPRECISE(tanh);
-DECLARE_IMPRECISE(tgamma);
+DECLARE_FORMER_IMPRECISE(cosh);
+DECLARE_FORMER_IMPRECISE(erfc);
+DECLARE_FORMER_IMPRECISE(erf);
+DECLARE_FORMER_IMPRECISE(sinh);
+DECLARE_FORMER_IMPRECISE(tanh);
+DECLARE_FORMER_IMPRECISE(tgamma);
index 22c2553..5cd239b 100644 (file)
@@ -11,7 +11,7 @@
 
 /*
  * from: @(#)fdlibm.h 5.1 93/09/24
- * $FreeBSD: head/lib/msun/src/math.h 253766 2013-07-29 12:33:03Z theraven $
+ * $FreeBSD: head/lib/msun/src/math.h 271651 2014-09-15 23:21:57Z kargl $
  */
 
 #ifndef _MATH_H_
@@ -452,7 +452,10 @@ long double        atanl(long double);
 long double    cbrtl(long double);
 long double    ceill(long double);
 long double    copysignl(long double, long double) __pure2;
+long double    coshl(long double);
 long double    cosl(long double);
+long double    erfcl(long double);
+long double    erfl(long double);
 long double    exp2l(long double);
 long double    expl(long double);
 long double    expm1l(long double);
@@ -467,6 +470,7 @@ long double frexpl(long double value, int *); /* fundamentally !__pure2 */
 long double    hypotl(long double, long double);
 int            ilogbl(long double) __pure2;
 long double    ldexpl(long double, int);
+long double    lgammal(long double);
 long long      llrintl(long double);
 long long      llroundl(long double);
 long double    log10l(long double);
@@ -483,45 +487,26 @@ long double       nextafterl(long double, long double);
 double         nexttoward(double, long double);
 float          nexttowardf(float, long double);
 long double    nexttowardl(long double, long double);
+long double    powl(long double, long double);
 long double    remainderl(long double, long double);
 long double    remquol(long double, long double, int *);
 long double    rintl(long double);
 long double    roundl(long double);
 long double    scalblnl(long double, long);
 long double    scalbnl(long double, int);
+long double    sinhl(long double);
 long double    sinl(long double);
 long double    sqrtl(long double);
+long double    tanhl(long double);
 long double    tanl(long double);
+long double    tgammal(long double);
 long double    truncl(long double);
-
 #endif /* __ISO_C_VISIBLE >= 1999 */
-__END_DECLS
-
-#endif /* !_MATH_H_ */
 
-/* separate header for cmath */
-#ifndef _MATH_EXTRA_H_
-#if __ISO_C_VISIBLE >= 1999
-#if _DECLARE_C99_LDBL_MATH
-
-#define _MATH_EXTRA_H_
-
-/*
- * extra long double versions of math functions for C99 and cmath
- */
-__BEGIN_DECLS
-
-long double    coshl(long double);
-long double    erfcl(long double);
-long double    erfl(long double);
-long double    lgammal(long double);
-long double    powl(long double, long double);
-long double    sinhl(long double);
-long double    tanhl(long double);
-long double    tgammal(long double);
+#if __BSD_VISIBLE
+long double    lgammal_r(long double, int *);
+#endif
 
 __END_DECLS
 
-#endif /* !_DECLARE_C99_LDBL_MATH */
-#endif /* __ISO_C_VISIBLE >= 1999 */
-#endif /* !_MATH_EXTRA_H_ */
+#endif /* !_MATH_H_ */
index d018c3c..f4207b1 100644 (file)
@@ -7,8 +7,6 @@
  * Permission to use, copy, modify, and distribute this
  * software is freely granted, provided that this notice
  * is preserved.
- *
- * $FreeBSD: head/lib/msun/src/s_erf.c 254994 2013-08-28 16:59:55Z kargl $
  * ====================================================
  */
 
 #include "math.h"
 #include "math_private.h"
 
+/* XXX Prevent compilers from erroneously constant folding: */
+static const volatile double tiny= 1e-300;
+
 static const double
-tiny       = 1e-300,
-half=  5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
-one =  1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
-two =  2.00000000000000000000e+00, /* 0x40000000, 0x00000000 */
-       /* c = (float)0.84506291151 */
+half= 0.5,
+one = 1,
+two = 2,
+/* c = (float)0.84506291151 */
 erx =  8.45062911510467529297e-01, /* 0x3FEB0AC1, 0x60000000 */
 /*
- * Coefficients for approximation to  erf on [0,0.84375]
+ * In the domain [0, 2**-28], only the first term in the power series
+ * expansion of erf(x) is used.  The magnitude of the first neglected
+ * terms is less than 2**-84.
  */
 efx =  1.28379167095512586316e-01, /* 0x3FC06EBA, 0x8214DB69 */
 efx8=  1.02703333676410069053e+00, /* 0x3FF06EBA, 0x8214DB69 */
+/*
+ * Coefficients for approximation to erf on [0,0.84375]
+ */
 pp0  =  1.28379167095512558561e-01, /* 0x3FC06EBA, 0x8214DB68 */
 pp1  = -3.25042107247001499370e-01, /* 0xBFD4CD7D, 0x691CB913 */
 pp2  = -2.84817495755985104766e-02, /* 0xBF9D2A51, 0xDBD7194F */
@@ -134,7 +139,7 @@ qq3  =  5.08130628187576562776e-03, /* 0x3F74D022, 0xC4D36B0F */
 qq4  =  1.32494738004321644526e-04, /* 0x3F215DC9, 0x221C1A10 */
 qq5  = -3.96022827877536812320e-06, /* 0xBED09C43, 0x42A26120 */
 /*
- * Coefficients for approximation to  erf  in [0.84375,1.25]
+ * Coefficients for approximation to erf in [0.84375,1.25]
  */
 pa0  = -2.36211856075265944077e-03, /* 0xBF6359B8, 0xBEF77538 */
 pa1  =  4.14856118683748331666e-01, /* 0x3FDA8D00, 0xAD92B34D */
@@ -150,7 +155,7 @@ qa4  =  1.26171219808761642112e-01, /* 0x3FC02660, 0xE763351F */
 qa5  =  1.36370839120290507362e-02, /* 0x3F8BEDC2, 0x6B51DD1C */
 qa6  =  1.19844998467991074170e-02, /* 0x3F888B54, 0x5735151D */
 /*
- * Coefficients for approximation to  erfc in [1.25,1/0.35]
+ * Coefficients for approximation to erfc in [1.25,1/0.35]
  */
 ra0  = -9.86494403484714822705e-03, /* 0xBF843412, 0x600D6435 */
 ra1  = -6.93858572707181764372e-01, /* 0xBFE63416, 0xE4BA7360 */
@@ -169,7 +174,7 @@ sa6  =  1.08635005541779435134e+02, /* 0x405B28A3, 0xEE48AE2C */
 sa7  =  6.57024977031928170135e+00, /* 0x401A47EF, 0x8E484A93 */
 sa8  = -6.04244152148580987438e-02, /* 0xBFAEEFF2, 0xEE749A62 */
 /*
- * Coefficients for approximation to  erfc in [1/.35,28]
+ * Coefficients for approximation to erfc in [1/.35,28]
  */
 rb0  = -9.86494292470009928597e-03, /* 0xBF843412, 0x39E86F4A */
 rb1  = -7.99283237680523006574e-01, /* 0xBFE993BA, 0x70C285DE */
@@ -222,15 +227,12 @@ erf(double x)
        x = fabs(x);
        s = one/(x*x);
        if(ix< 0x4006DB6E) {    /* |x| < 1/0.35 */
-           R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
-                               ra5+s*(ra6+s*ra7))))));
-           S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
-                               sa5+s*(sa6+s*(sa7+s*sa8)))))));
+           R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(ra5+s*(ra6+s*ra7))))));
+           S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(sa5+s*(sa6+s*(sa7+
+               s*sa8)))))));
        } else {        /* |x| >= 1/0.35 */
-           R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
-                               rb5+s*rb6)))));
-           S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
-                               sb5+s*(sb6+s*sb7))))));
+           R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(rb5+s*rb6)))));
+           S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(sb5+s*(sb6+s*sb7))))));
        }
        z  = x;
        SET_LOW_WORD(z,0);
@@ -238,6 +240,10 @@ erf(double x)
        if(hx>=0) return one-r/x; else return  r/x-one;
 }
 
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(erf, erfl);
+#endif
+
 double
 erfc(double x)
 {
@@ -279,23 +285,23 @@ erfc(double x)
            x = fabs(x);
            s = one/(x*x);
            if(ix< 0x4006DB6D) {        /* |x| < 1/.35 ~ 2.857143*/
-               R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
-                               ra5+s*(ra6+s*ra7))))));
-               S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
-                               sa5+s*(sa6+s*(sa7+s*sa8)))))));
+               R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(ra5+s*(ra6+s*ra7))))));
+               S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(sa5+s*(sa6+s*(sa7+
+                   s*sa8)))))));
            } else {                    /* |x| >= 1/.35 ~ 2.857143 */
                if(hx<0&&ix>=0x40180000) return two-tiny;/* x < -6 */
-               R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
-                               rb5+s*rb6)))));
-               S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
-                               sb5+s*(sb6+s*sb7))))));
+               R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(rb5+s*rb6)))));
+               S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(sb5+s*(sb6+s*sb7))))));
            }
            z  = x;
            SET_LOW_WORD(z,0);
-           r  =  __ieee754_exp(-z*z-0.5625)*
-                       __ieee754_exp((z-x)*(z+x)+R/S);
+           r  =  __ieee754_exp(-z*z-0.5625)*__ieee754_exp((z-x)*(z+x)+R/S);
            if(hx>0) return r/x; else return two-r/x;
        } else {
            if(hx>0) return tiny*tiny; else return two-tiny;
        }
 }
+
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(erfc, erfcl);
+#endif
index eece494..d71f8c4 100644 (file)
  * is preserved.
  * ====================================================
  *
- * $FreeBSD: head/lib/msun/src/s_erff.c 254969 2013-08-27 19:46:56Z kargl $
+ * $FreeBSD: head/lib/msun/src/s_erff.c 268590 2014-07-13 16:24:16Z kargl $
  */
 
 
 #include "math.h"
 #include "math_private.h"
 
+/* XXX Prevent compilers from erroneously constant folding: */
+static const volatile float tiny = 1e-30;
+
 static const float
-tiny       = 1e-30,
-half=  5.0000000000e-01, /* 0x3F000000 */
-one =  1.0000000000e+00, /* 0x3F800000 */
-two =  2.0000000000e+00, /* 0x40000000 */
+half= 0.5,
+one = 1,
+two = 2,
+erx = 8.42697144e-01,                  /* 0x3f57bb00 */
 /*
- * Coefficients for approximation to erf on [0,0.84375]
+ * In the domain [0, 2**-14], only the first term in the power series
+ * expansion of erf(x) is used.  The magnitude of the first neglected
+ * terms is less than 2**-42.
  */
-efx =  1.2837916613e-01, /* 0x3e0375d4 */
-efx8=  1.0270333290e+00, /* 0x3f8375d4 */
+efx = 1.28379166e-01, /* 0x3e0375d4 */
+efx8= 1.02703333e+00, /* 0x3f8375d4 */
 /*
- *  Domain [0, 0.84375], range ~[-5.4446e-10,5.5197e-10]:
- *  |(erf(x) - x)/x - p(x)/q(x)| < 2**-31.
+ * Domain [0, 0.84375], range ~[-5.4419e-10, 5.5179e-10]:
+ * |(erf(x) - x)/x - pp(x)/qq(x)| < 2**-31
  */
-pp0  =  1.28379166e-01F, /*  0x1.06eba8p-3 */
-pp1  = -3.36030394e-01F, /* -0x1.58185ap-2 */
-pp2  = -1.86260219e-03F, /* -0x1.e8451ep-10 */
-qq1  =  3.12324286e-01F, /*  0x1.3fd1f0p-2 */
-qq2  =  2.16070302e-02F, /*  0x1.620274p-6 */
-qq3  = -1.98859419e-03F, /* -0x1.04a626p-9 */
+pp0  =  1.28379166e-01, /* 0x3e0375d4 */
+pp1  = -3.36030394e-01, /* 0xbeac0c2d */
+pp2  = -1.86261395e-03, /* 0xbaf422f4 */
+qq1  =  3.12324315e-01, /* 0x3e9fe8f9 */
+qq2  =  2.16070414e-02, /* 0x3cb10140 */
+qq3  = -1.98859372e-03, /* 0xbb025311 */
 /*
- * Domain [0.84375, 1.25], range ~[-1.953e-11,1.940e-11]:
- * |(erf(x) - erx) - p(x)/q(x)| < 2**-36.
+ * Domain [0.84375, 1.25], range ~[-1.023e-9, 1.023e-9]:
+ * |(erf(x) - erx) - pa(x)/qa(x)| < 2**-31
  */
-erx  =  8.42697144e-01F, /*  0x1.af7600p-1.  erf(1) rounded to 16 bits. */
-pa0  =  3.64939137e-06F, /*  0x1.e9d022p-19 */
-pa1  =  4.15109694e-01F, /*  0x1.a91284p-2 */
-pa2  = -1.65179938e-01F, /* -0x1.5249dcp-3 */
-pa3  =  1.10914491e-01F, /*  0x1.c64e46p-4 */
-qa1  =  6.02074385e-01F, /*  0x1.344318p-1 */
-qa2  =  5.35934687e-01F, /*  0x1.126608p-1 */
-qa3  =  1.68576106e-01F, /*  0x1.593e6ep-3 */
-qa4  =  5.62181212e-02F, /*  0x1.cc89f2p-5 */
+pa0  =  3.65041046e-06, /* 0x3674f993 */
+pa1  =  4.15109307e-01, /* 0x3ed48935 */
+pa2  = -2.09395722e-01, /* 0xbe566bd5 */
+pa3  =  8.67677554e-02, /* 0x3db1b34b */
+qa1  =  4.95560974e-01, /* 0x3efdba2b */
+qa2  =  3.71248513e-01, /* 0x3ebe1449 */
+qa3  =  3.92478965e-02, /* 0x3d20c267 */
 /*
- * Domain [1.25,1/0.35], range ~[-7.043e-10,7.457e-10]:
- * |log(x*erfc(x)) + x**2 + 0.5625 - r(x)/s(x)| < 2**-30
+ * Domain [1.25,1/0.35], range ~[-4.821e-9, 4.927e-9]:
+ * |log(x*erfc(x)) + x**2 + 0.5625 - ra(x)/sa(x)| < 2**-28
  */
-ra0  = -9.87132732e-03F, /* -0x1.4376b2p-7 */
-ra1  = -5.53605914e-01F, /* -0x1.1b723cp-1 */
-ra2  = -2.17589188e+00F, /* -0x1.1683a0p+1 */
-ra3  = -1.43268085e+00F, /* -0x1.6ec42cp+0 */
-sa1  =  5.45995426e+00F, /*  0x1.5d6fe4p+2 */
-sa2  =  6.69798088e+00F, /*  0x1.acabb8p+2 */
-sa3  =  1.43113089e+00F, /*  0x1.6e5e98p+0 */
-sa4  = -5.77397496e-02F, /* -0x1.d90108p-5 */
+ra0  = -9.88156721e-03, /* 0xbc21e64c */
+ra1  = -5.43658376e-01, /* 0xbf0b2d32 */
+ra2  = -1.66828310e+00, /* 0xbfd58a4d */
+ra3  = -6.91554189e-01, /* 0xbf3109b2 */
+sa1  =  4.48581553e+00, /* 0x408f8bcd */
+sa2  =  4.10799170e+00, /* 0x408374ab */
+sa3  =  5.53855181e-01, /* 0x3f0dc974 */
 /*
- * Domain [1/0.35, 11], range ~[-2.264e-13,2.336e-13]:
- * |log(x*erfc(x)) + x**2 + 0.5625 - r(x)/s(x)| < 2**-42
+ * Domain [2.85715, 11], range ~[-1.484e-9, 1.505e-9]:
+ * |log(x*erfc(x)) + x**2 + 0.5625 - rb(x)/sb(x)| < 2**-30
  */
-rb0  = -9.86494310e-03F, /* -0x1.434124p-7 */
-rb1  = -6.25171244e-01F, /* -0x1.401672p-1 */
-rb2  = -6.16498327e+00F, /* -0x1.8a8f16p+2 */
-rb3  = -1.66696873e+01F, /* -0x1.0ab70ap+4 */
-rb4  = -9.53764343e+00F, /* -0x1.313460p+3 */
-sb1  =  1.26884899e+01F, /*  0x1.96081cp+3 */
-sb2  =  4.51839523e+01F, /*  0x1.6978bcp+5 */
-sb3  =  4.72810211e+01F, /*  0x1.7a3f88p+5 */
-sb4  =  8.93033314e+00F; /*  0x1.1dc54ap+3 */
+rb0  = -9.86496918e-03, /* 0xbc21a0ae */
+rb1  = -5.48049808e-01, /* 0xbf0c4cfe */
+rb2  = -1.84115684e+00, /* 0xbfebab07 */
+sb1  =  4.87132740e+00, /* 0x409be1ea */
+sb2  =  3.04982710e+00, /* 0x4043305e */
+sb3  = -7.61900663e-01; /* 0xbf430bec */
 
 float
 erff(float x)
@@ -85,9 +84,9 @@ erff(float x)
        float R,S,P,Q,s,y,z,r;
        GET_FLOAT_WORD(hx,x);
        ix = hx&0x7fffffff;
-       if(ix>=0x7f800000) {            /* erf(nan)=nan */
+       if(ix>=0x7f800000) {            /* erff(nan)=nan */
            i = ((u_int32_t)hx>>31)<<1;
-           return (float)(1-i)+one/x;  /* erf(+-inf)=+-1 */
+           return (float)(1-i)+one/x;  /* erff(+-inf)=+-1 */
        }
 
        if(ix < 0x3f580000) {           /* |x|<0.84375 */
@@ -105,7 +104,7 @@ erff(float x)
        if(ix < 0x3fa00000) {           /* 0.84375 <= |x| < 1.25 */
            s = fabsf(x)-one;
            P = pa0+s*(pa1+s*(pa2+s*pa3));
-           Q = one+s*(qa1+s*(qa2+s*(qa3+s*qa4)));
+           Q = one+s*(qa1+s*(qa2+s*qa3));
            if(hx>=0) return erx + P/Q; else return -erx - P/Q;
        }
        if (ix >= 0x40800000) {         /* inf>|x|>=4 */
@@ -113,12 +112,12 @@ erff(float x)
        }
        x = fabsf(x);
        s = one/(x*x);
-       if(ix< 0x4036DB6E) {    /* |x| < 1/0.35 */
+       if(ix< 0x4036db8c) {    /* |x| < 2.85715 ~ 1/0.35 */
            R=ra0+s*(ra1+s*(ra2+s*ra3));
-           S=one+s*(sa1+s*(sa2+s*(sa3+s*sa4)));
-       } else {        /* |x| >= 1/0.35 */
-           R=rb0+s*(rb1+s*(rb2+s*(rb3+s*rb4)));
-           S=one+s*(sb1+s*(sb2+s*(sb3+s*sb4)));
+           S=one+s*(sa1+s*(sa2+s*sa3));
+       } else {        /* |x| >= 2.85715 ~ 1/0.35 */
+           R=rb0+s*(rb1+s*rb2);
+           S=one+s*(sb1+s*(sb2+s*sb3));
        }
        SET_FLOAT_WORD(z,hx&0xffffe000);
        r  = expf(-z*z-0.5625F)*expf((z-x)*(z+x)+R/S);
@@ -132,8 +131,8 @@ erfcf(float x)
        float R,S,P,Q,s,y,z,r;
        GET_FLOAT_WORD(hx,x);
        ix = hx&0x7fffffff;
-       if(ix>=0x7f800000) {                    /* erfc(nan)=nan */
-                                               /* erfc(+-inf)=0,2 */
+       if(ix>=0x7f800000) {                    /* erfcf(nan)=nan */
+                                               /* erfcf(+-inf)=0,2 */
            return (float)(((u_int32_t)hx>>31)<<1)+one/x;
        }
 
@@ -155,7 +154,7 @@ erfcf(float x)
        if(ix < 0x3fa00000) {           /* 0.84375 <= |x| < 1.25 */
            s = fabsf(x)-one;
            P = pa0+s*(pa1+s*(pa2+s*pa3));
-           Q = one+s*(qa1+s*(qa2+s*(qa3+s*qa4)));
+           Q = one+s*(qa1+s*(qa2+s*qa3));
            if(hx>=0) {
                z  = one-erx; return z - P/Q;
            } else {
@@ -165,13 +164,13 @@ erfcf(float x)
        if (ix < 0x41300000) {          /* |x|<11 */
            x = fabsf(x);
            s = one/(x*x);
-           if(ix< 0x4036DB6D) {        /* |x| < 1/.35 ~ 2.857143*/
-               R=ra0+s*(ra1+s*(ra2+s*ra3));
-               S=one+s*(sa1+s*(sa2+s*(sa3+s*sa4)));
-           } else {                    /* |x| >= 1/.35 ~ 2.857143 */
+           if(ix< 0x4036db8c) {        /* |x| < 2.85715 ~ 1/.35 */
+               R=ra0+s*(ra1+s*(ra2+s*ra3));
+               S=one+s*(sa1+s*(sa2+s*sa3));
+           } else {                    /* |x| >= 2.85715 ~ 1/.35 */
                if(hx<0&&ix>=0x40a00000) return two-tiny;/* x < -5 */
-               R=rb0+s*(rb1+s*(rb2+s*(rb3+s*rb4)));
-               S=one+s*(sb1+s*(sb2+s*(sb3+s*sb4)));
+               R=rb0+s*(rb1+s*rb2);
+               S=one+s*(sb1+s*(sb2+s*sb3));
            }
            SET_FLOAT_WORD(z,hx&0xffffe000);
            r  = expf(-z*z-0.5625F)*expf((z-x)*(z+x)+R/S);
diff --git a/lib/libm/src/s_llrint.c b/lib/libm/src/s_llrint.c
new file mode 100644 (file)
index 0000000..b1a9447
--- /dev/null
@@ -0,0 +1,7 @@
+
+#define type           double
+#define        roundit         rint
+#define dtype          long long
+#define        fn              llrint
+
+#include "s_lrint.c"
diff --git a/lib/libm/src/s_llrintf.c b/lib/libm/src/s_llrintf.c
new file mode 100644 (file)
index 0000000..1ffd3ce
--- /dev/null
@@ -0,0 +1,7 @@
+
+#define type           float
+#define        roundit         rintf
+#define dtype          long long
+#define        fn              llrintf
+
+#include "s_lrint.c"
diff --git a/lib/libm/src/s_llrintl.c b/lib/libm/src/s_llrintl.c
new file mode 100644 (file)
index 0000000..6b1c4be
--- /dev/null
@@ -0,0 +1,7 @@
+
+#define type           long double
+#define        roundit         rintl
+#define dtype          long long
+#define        fn              llrintl
+
+#include "s_lrint.c"
diff --git a/lib/libm/src/s_logbl.c b/lib/libm/src/s_logbl.c
new file mode 100644 (file)
index 0000000..e06f328
--- /dev/null
@@ -0,0 +1,53 @@
+/*
+ * From: @(#)s_ilogb.c 5.1 93/09/24
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ * $FreeBSD: head/lib/msun/src/s_logbl.c 174698 2007-12-17 03:53:38Z das $
+ */
+
+
+#include <float.h>
+#include <limits.h>
+#include <math.h>
+
+#include "fpmath.h"
+
+long double
+logbl(long double x)
+{
+       union IEEEl2bits u;
+       unsigned long m;
+       int b;
+
+       u.e = x;
+       if (u.bits.exp == 0) {
+               if ((u.bits.manl | u.bits.manh) == 0) { /* x == 0 */
+                       u.bits.sign = 1;
+                       return (1.0L / u.e);
+               }
+               /* denormalized */
+               if (u.bits.manh == 0) {
+                       m = 1lu << (LDBL_MANL_SIZE - 1);
+                       for (b = LDBL_MANH_SIZE; !(u.bits.manl & m); m >>= 1)
+                               b++;
+               } else {
+                       m = 1lu << (LDBL_MANH_SIZE - 1);
+                       for (b = 0; !(u.bits.manh & m); m >>= 1)
+                               b++;
+               }
+#ifdef LDBL_IMPLICIT_NBIT
+               b++;
+#endif
+               return ((long double)(LDBL_MIN_EXP - b - 1));
+       }
+       if (u.bits.exp < (LDBL_MAX_EXP << 1) - 1)       /* normal */
+               return ((long double)(u.bits.exp - LDBL_MAX_EXP + 1));
+       else                                            /* +/- inf or nan */
+               return (x * x);
+}
similarity index 56%
copy from lib/libm/src/imprecise.c
copy to lib/libm/src/s_lrint.c
index 2498be3..f8cd420 100644 (file)
@@ -1,5 +1,5 @@
 /*-
- * Copyright (c) 2013 David Chisnall
+ * Copyright (c) 2005 David Schultz <das@FreeBSD.ORG>
  * All rights reserved.
  *
  * Redistribution and use in source and binary forms, with or without
  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  * SUCH DAMAGE.
  *
- * $FreeBSD: head/lib/msun/src/imprecise.c 255294 2013-09-06 07:58:23Z theraven $
+ * $FreeBSD: head/lib/msun/src/s_lrint.c 140088 2005-01-11 23:12:55Z das $
  */
 
-#include <float.h>
+#include <fenv.h>
 #include <math.h>
 
-/*
- * If long double is not the same size as double, then these will lose
- * precision and we should emit a warning whenever something links against
- * them.
- */
-#if (LDBL_MANT_DIG > 53)
-#define WARN_IMPRECISE(x) \
-       __warn_references(x, # x " has lower than advertised precision");
-#else
-#define WARN_IMPRECISE(x)
+#ifndef type
+#define type           double
+#define        roundit         rint
+#define dtype          long
+#define        fn              lrint
 #endif
+
 /*
- * Declare the functions as weak variants so that other libraries providing
- * real versions can override them.
+ * C99 says we should not raise a spurious inexact exception when an
+ * invalid exception is raised.  Unfortunately, the set of inputs
+ * that overflows depends on the rounding mode when 'dtype' has more
+ * significant bits than 'type'.  Hence, we bend over backwards for the
+ * sake of correctness; an MD implementation could be more efficient.
  */
-#define        DECLARE_WEAK(x)\
-       __weak_reference(imprecise_## x, x);\
-       WARN_IMPRECISE(x)
-
-long double
-imprecise_powl(long double x, long double y)
+dtype
+fn(type x)
 {
+       fenv_t env;
+       dtype d;
 
-       return pow(x, y);
+       feholdexcept(&env);
+       d = (dtype)roundit(x);
+       if (fetestexcept(FE_INVALID))
+               feclearexcept(FE_INEXACT);
+       feupdateenv(&env);
+       return (d);
 }
-DECLARE_WEAK(powl);
-
-#define DECLARE_IMPRECISE(f) \
-       long double imprecise_ ## f ## l(long double v) { return f(v); }\
-       DECLARE_WEAK(f ## l)
-
-DECLARE_IMPRECISE(cosh);
-DECLARE_IMPRECISE(erfc);
-DECLARE_IMPRECISE(erf);
-DECLARE_IMPRECISE(lgamma);
-DECLARE_IMPRECISE(sinh);
-DECLARE_IMPRECISE(tanh);
-DECLARE_IMPRECISE(tgamma);
diff --git a/lib/libm/src/s_lrintf.c b/lib/libm/src/s_lrintf.c
new file mode 100644 (file)
index 0000000..5dc7873
--- /dev/null
@@ -0,0 +1,7 @@
+
+#define type           float
+#define        roundit         rintf
+#define dtype          long
+#define        fn              lrintf
+
+#include "s_lrint.c"
diff --git a/lib/libm/src/s_lrintl.c b/lib/libm/src/s_lrintl.c
new file mode 100644 (file)
index 0000000..c7edaca
--- /dev/null
@@ -0,0 +1,7 @@
+
+#define type           long double
+#define        roundit         rintl
+#define dtype          long
+#define        fn              lrintl
+
+#include "s_lrint.c"
diff --git a/lib/libm/src/s_remquo.c b/lib/libm/src/s_remquo.c
new file mode 100644 (file)
index 0000000..d9d9032
--- /dev/null
@@ -0,0 +1,157 @@
+/* @(#)e_fmod.c 1.3 95/01/18 */
+/*-
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+
+#include <float.h>
+
+#include "math.h"
+#include "math_private.h"
+
+static const double Zero[] = {0.0, -0.0,};
+
+/*
+ * Return the IEEE remainder and set *quo to the last n bits of the
+ * quotient, rounded to the nearest integer.  We choose n=31 because
+ * we wind up computing all the integer bits of the quotient anyway as
+ * a side-effect of computing the remainder by the shift and subtract
+ * method.  In practice, this is far more bits than are needed to use
+ * remquo in reduction algorithms.
+ */
+double
+remquo(double x, double y, int *quo)
+{
+       int32_t n,hx,hy,hz,ix,iy,sx,i;
+       u_int32_t lx,ly,lz,q,sxy;
+
+       EXTRACT_WORDS(hx,lx,x);
+       EXTRACT_WORDS(hy,ly,y);
+       sxy = (hx ^ hy) & 0x80000000;
+       sx = hx&0x80000000;             /* sign of x */
+       hx ^=sx;                /* |x| */
+       hy &= 0x7fffffff;       /* |y| */
+
+    /* purge off exception values */
+       if((hy|ly)==0||(hx>=0x7ff00000)||       /* y=0,or x not finite */
+         ((hy|((ly|-ly)>>31))>0x7ff00000))     /* or y is NaN */
+           return (x*y)/(x*y);
+       if(hx<=hy) {
+           if((hx<hy)||(lx<ly)) {
+               q = 0;
+               goto fixup;     /* |x|<|y| return x or x-y */
+           }
+           if(lx==ly) {
+               *quo = (sxy ? -1 : 1);
+               return Zero[(u_int32_t)sx>>31]; /* |x|=|y| return x*0*/
+           }
+       }
+
+    /* determine ix = ilogb(x) */
+       if(hx<0x00100000) {     /* subnormal x */
+           if(hx==0) {
+               for (ix = -1043, i=lx; i>0; i<<=1) ix -=1;
+           } else {
+               for (ix = -1022,i=(hx<<11); i>0; i<<=1) ix -=1;
+           }
+       } else ix = (hx>>20)-1023;
+
+    /* determine iy = ilogb(y) */
+       if(hy<0x00100000) {     /* subnormal y */
+           if(hy==0) {
+               for (iy = -1043, i=ly; i>0; i<<=1) iy -=1;
+           } else {
+               for (iy = -1022,i=(hy<<11); i>0; i<<=1) iy -=1;
+           }
+       } else iy = (hy>>20)-1023;
+
+    /* set up {hx,lx}, {hy,ly} and align y to x */
+       if(ix >= -1022) 
+           hx = 0x00100000|(0x000fffff&hx);
+       else {          /* subnormal x, shift x to normal */
+           n = -1022-ix;
+           if(n<=31) {
+               hx = (hx<<n)|(lx>>(32-n));
+               lx <<= n;
+           } else {
+               hx = lx<<(n-32);
+               lx = 0;
+           }
+       }
+       if(iy >= -1022) 
+           hy = 0x00100000|(0x000fffff&hy);
+       else {          /* subnormal y, shift y to normal */
+           n = -1022-iy;
+           if(n<=31) {
+               hy = (hy<<n)|(ly>>(32-n));
+               ly <<= n;
+           } else {
+               hy = ly<<(n-32);
+               ly = 0;
+           }
+       }
+
+    /* fix point fmod */
+       n = ix - iy;
+       q = 0;
+       while(n--) {
+           hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
+           if(hz<0){hx = hx+hx+(lx>>31); lx = lx+lx;}
+           else {hx = hz+hz+(lz>>31); lx = lz+lz; q++;}
+           q <<= 1;
+       }
+       hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
+       if(hz>=0) {hx=hz;lx=lz;q++;}
+
+    /* convert back to floating value and restore the sign */
+       if((hx|lx)==0) {                        /* return sign(x)*0 */
+           q &= 0x7fffffff;
+           *quo = (sxy ? -q : q);
+           return Zero[(u_int32_t)sx>>31];
+       }
+       while(hx<0x00100000) {          /* normalize x */
+           hx = hx+hx+(lx>>31); lx = lx+lx;
+           iy -= 1;
+       }
+       if(iy>= -1022) {        /* normalize output */
+           hx = ((hx-0x00100000)|((iy+1023)<<20));
+       } else {                /* subnormal output */
+           n = -1022 - iy;
+           if(n<=20) {
+               lx = (lx>>n)|((u_int32_t)hx<<(32-n));
+               hx >>= n;
+           } else if (n<=31) {
+               lx = (hx<<(32-n))|(lx>>n); hx = 0;
+           } else {
+               lx = hx>>(n-32); hx = 0;
+           }
+       }
+fixup:
+       INSERT_WORDS(x,hx,lx);
+       y = fabs(y);
+       if (y < 0x1p-1021) {
+           if (x+x>y || (x+x==y && (q & 1))) {
+               q++;
+               x-=y;
+           }
+       } else if (x>0.5*y || (x==0.5*y && (q & 1))) {
+           q++;
+           x-=y;
+       }
+       GET_HIGH_WORD(hx,x);
+       SET_HIGH_WORD(x,hx^sx);
+       q &= 0x7fffffff;
+       *quo = (sxy ? -q : q);
+       return x;
+}
+
+#if LDBL_MANT_DIG == 53
+__weak_reference(remquo, remquol);
+#endif
diff --git a/lib/libm/src/s_remquof.c b/lib/libm/src/s_remquof.c
new file mode 100644 (file)
index 0000000..7bc66e1
--- /dev/null
@@ -0,0 +1,120 @@
+/* @(#)e_fmod.c 1.3 95/01/18 */
+/*-
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+
+#include "math.h"
+#include "math_private.h"
+
+static const float Zero[] = {0.0, -0.0,};
+
+/*
+ * Return the IEEE remainder and set *quo to the last n bits of the
+ * quotient, rounded to the nearest integer.  We choose n=31 because
+ * we wind up computing all the integer bits of the quotient anyway as
+ * a side-effect of computing the remainder by the shift and subtract
+ * method.  In practice, this is far more bits than are needed to use
+ * remquo in reduction algorithms.
+ */
+float
+remquof(float x, float y, int *quo)
+{
+       int32_t n,hx,hy,hz,ix,iy,sx,i;
+       u_int32_t q,sxy;
+
+       GET_FLOAT_WORD(hx,x);
+       GET_FLOAT_WORD(hy,y);
+       sxy = (hx ^ hy) & 0x80000000;
+       sx = hx&0x80000000;             /* sign of x */
+       hx ^=sx;                /* |x| */
+       hy &= 0x7fffffff;       /* |y| */
+
+    /* purge off exception values */
+       if(hy==0||hx>=0x7f800000||hy>0x7f800000) /* y=0,NaN;or x not finite */
+           return (x*y)/(x*y);
+       if(hx<hy) {
+           q = 0;
+           goto fixup; /* |x|<|y| return x or x-y */
+       } else if(hx==hy) {
+           *quo = (sxy ? -1 : 1);
+           return Zero[(u_int32_t)sx>>31];     /* |x|=|y| return x*0*/
+       }
+
+    /* determine ix = ilogb(x) */
+       if(hx<0x00800000) {     /* subnormal x */
+           for (ix = -126,i=(hx<<8); i>0; i<<=1) ix -=1;
+       } else ix = (hx>>23)-127;
+
+    /* determine iy = ilogb(y) */
+       if(hy<0x00800000) {     /* subnormal y */
+           for (iy = -126,i=(hy<<8); i>0; i<<=1) iy -=1;
+       } else iy = (hy>>23)-127;
+
+    /* set up {hx,lx}, {hy,ly} and align y to x */
+       if(ix >= -126)
+           hx = 0x00800000|(0x007fffff&hx);
+       else {          /* subnormal x, shift x to normal */
+           n = -126-ix;
+           hx <<= n;
+       }
+       if(iy >= -126)
+           hy = 0x00800000|(0x007fffff&hy);
+       else {          /* subnormal y, shift y to normal */
+           n = -126-iy;
+           hy <<= n;
+       }
+
+    /* fix point fmod */
+       n = ix - iy;
+       q = 0;
+       while(n--) {
+           hz=hx-hy;
+           if(hz<0) hx = hx << 1;
+           else {hx = hz << 1; q++;}
+           q <<= 1;
+       }
+       hz=hx-hy;
+       if(hz>=0) {hx=hz;q++;}
+
+    /* convert back to floating value and restore the sign */
+       if(hx==0) {                             /* return sign(x)*0 */
+           q &= 0x7fffffff;
+           *quo = (sxy ? -q : q);
+           return Zero[(u_int32_t)sx>>31];
+       }
+       while(hx<0x00800000) {          /* normalize x */
+           hx <<= 1;
+           iy -= 1;
+       }
+       if(iy>= -126) {         /* normalize output */
+           hx = ((hx-0x00800000)|((iy+127)<<23));
+       } else {                /* subnormal output */
+           n = -126 - iy;
+           hx >>= n;
+       }
+fixup:
+       SET_FLOAT_WORD(x,hx);
+       y = fabsf(y);
+       if (y < 0x1p-125f) {
+           if (x+x>y || (x+x==y && (q & 1))) {
+               q++;
+               x-=y;
+           }
+       } else if (x>0.5f*y || (x==0.5f*y && (q & 1))) {
+           q++;
+           x-=y;
+       }
+       GET_FLOAT_WORD(hx,x);
+       SET_FLOAT_WORD(x,hx^sx);
+       q &= 0x7fffffff;
+       *quo = (sxy ? -q : q);
+       return x;
+}
diff --git a/lib/libm/src/s_remquol.c b/lib/libm/src/s_remquol.c
new file mode 100644 (file)
index 0000000..330583a
--- /dev/null
@@ -0,0 +1,176 @@
+/* @(#)e_fmod.c 1.3 95/01/18 */
+/*-
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+
+#include <float.h>
+#include <stdint.h>
+
+#include "fpmath.h"
+#include "math.h"
+#include "math_private.h"
+
+#define        BIAS (LDBL_MAX_EXP - 1)
+
+#if LDBL_MANL_SIZE > 32
+typedef        uint64_t manl_t;
+#else
+typedef        uint32_t manl_t;
+#endif
+
+#if LDBL_MANH_SIZE > 32
+typedef        uint64_t manh_t;
+#else
+typedef        uint32_t manh_t;
+#endif
+
+/*
+ * These macros add and remove an explicit integer bit in front of the
+ * fractional mantissa, if the architecture doesn't have such a bit by
+ * default already.
+ */
+#ifdef LDBL_IMPLICIT_NBIT
+#define        SET_NBIT(hx)    ((hx) | (1ULL << LDBL_MANH_SIZE))
+#define        HFRAC_BITS      LDBL_MANH_SIZE
+#else
+#define        SET_NBIT(hx)    (hx)
+#define        HFRAC_BITS      (LDBL_MANH_SIZE - 1)
+#endif
+
+#define        MANL_SHIFT      (LDBL_MANL_SIZE - 1)
+
+static const long double Zero[] = {0.0L, -0.0L};
+
+/*
+ * Return the IEEE remainder and set *quo to the last n bits of the
+ * quotient, rounded to the nearest integer.  We choose n=31 because
+ * we wind up computing all the integer bits of the quotient anyway as
+ * a side-effect of computing the remainder by the shift and subtract
+ * method.  In practice, this is far more bits than are needed to use
+ * remquo in reduction algorithms.
+ *
+ * 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
+remquol(long double x, long double y, int *quo)
+{
+       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,q,sx,sxy;
+
+       ux.e = x;
+       uy.e = y;
+       sx = ux.bits.sign;
+       sxy = sx ^ uy.bits.sign;
+       ux.bits.sign = 0;       /* |x| */
+       uy.bits.sign = 0;       /* |y| */
+       x = ux.e;
+
+    /* purge off exception values */
+       if((uy.bits.exp|uy.bits.manh|uy.bits.manl)==0 || /* y=0 */
+          (ux.bits.exp == BIAS + LDBL_MAX_EXP) ||       /* or x not finite */
+          (uy.bits.exp == BIAS + LDBL_MAX_EXP &&
+           ((uy.bits.manh&~LDBL_NBIT)|uy.bits.manl)!=0)) /* or y is NaN */
+           return (x*y)/(x*y);
+       if(ux.bits.exp<=uy.bits.exp) {
+           if((ux.bits.exp<uy.bits.exp) ||
+              (ux.bits.manh<=uy.bits.manh &&
+               (ux.bits.manh<uy.bits.manh ||
+                ux.bits.manl<uy.bits.manl))) {
+               q = 0;
+               goto fixup;     /* |x|<|y| return x or x-y */
+           }
+           if(ux.bits.manh==uy.bits.manh && ux.bits.manl==uy.bits.manl) {
+               *quo = (sxy ? -1 : 1);
+               return Zero[sx];        /* |x|=|y| return x*0*/
+           }
+       }
+
+    /* determine ix = ilogb(x) */
+       if(ux.bits.exp == 0) {  /* subnormal x */
+           ux.e *= 0x1.0p512;
+           ix = ux.bits.exp - (BIAS + 512);
+       } else {
+           ix = ux.bits.exp - BIAS;
+       }
+
+    /* determine iy = ilogb(y) */
+       if(uy.bits.exp == 0) {  /* subnormal y */
+           uy.e *= 0x1.0p512;
+           iy = uy.bits.exp - (BIAS + 512);
+       } else {
+           iy = uy.bits.exp - BIAS;
+       }
+
+    /* set up {hx,lx}, {hy,ly} and align y to x */
+       hx = SET_NBIT(ux.bits.manh);
+       hy = SET_NBIT(uy.bits.manh);
+       lx = ux.bits.manl;
+       ly = uy.bits.manl;
+
+    /* fix point fmod */
+       n = ix - iy;
+       q = 0;
+
+       while(n--) {
+           hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
+           if(hz<0){hx = hx+hx+(lx>>MANL_SHIFT); lx = lx+lx;}
+           else {hx = hz+hz+(lz>>MANL_SHIFT); lx = lz+lz; q++;}
+           q <<= 1;
+       }
+       hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
+       if(hz>=0) {hx=hz;lx=lz;q++;}
+
+    /* convert back to floating value and restore the sign */
+       if((hx|lx)==0) {                        /* return sign(x)*0 */
+           q &= 0x7fffffff;
+           *quo = (sxy ? -q : q);
+           return Zero[sx];
+       }
+       while(hx<(1ULL<<HFRAC_BITS)) {  /* normalize x */
+           hx = hx+hx+(lx>>MANL_SHIFT); lx = lx+lx;
+           iy -= 1;
+       }
+       ux.bits.manh = hx; /* The integer bit is truncated here if needed. */
+       ux.bits.manl = lx;
+       if (iy < LDBL_MIN_EXP) {
+           ux.bits.exp = iy + (BIAS + 512);
+           ux.e *= 0x1p-512;
+       } else {
+           ux.bits.exp = iy + BIAS;
+       }
+       ux.bits.sign = 0;
+       x = ux.e;
+fixup:
+       y = fabsl(y);
+       if (y < LDBL_MIN * 2) {
+           if (x+x>y || (x+x==y && (q & 1))) {
+               q++;
+               x-=y;
+           }
+       } else if (x>0.5*y || (x==0.5*y && (q & 1))) {
+           q++;
+           x-=y;
+       }
+
+       ux.e = x;
+       ux.bits.sign ^= sx;
+       x = ux.e;
+
+       q &= 0x7fffffff;
+       *quo = (sxy ? -q : q);
+       return x;
+}
diff --git a/lib/libm/src/s_rintl.c b/lib/libm/src/s_rintl.c
new file mode 100644 (file)
index 0000000..324f17a
--- /dev/null
@@ -0,0 +1,90 @@
+/*-
+ * Copyright (c) 2008 David Schultz <das@FreeBSD.ORG>
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ *
+ * $FreeBSD: head/lib/msun/src/s_rintl.c 176461 2008-02-22 11:59:05Z bde $
+ */
+
+
+#include <float.h>
+#include <math.h>
+
+#include "fpmath.h"
+
+#if LDBL_MAX_EXP != 0x4000
+/* We also require the usual bias, min exp and expsign packing. */
+#error "Unsupported long double format"
+#endif
+
+#define        BIAS    (LDBL_MAX_EXP - 1)
+
+static const float
+shift[2] = {
+#if LDBL_MANT_DIG == 64
+       0x1.0p63, -0x1.0p63
+#elif LDBL_MANT_DIG == 113
+       0x1.0p112, -0x1.0p112
+#else
+#error "Unsupported long double format"
+#endif
+};
+static const float zero[2] = { 0.0, -0.0 };
+
+long double
+rintl(long double x)
+{
+       union IEEEl2bits u;
+       uint32_t expsign;
+       int ex, sign;
+
+       u.e = x;
+       expsign = u.xbits.expsign;
+       ex = expsign & 0x7fff;
+
+       if (ex >= BIAS + LDBL_MANT_DIG - 1) {
+               if (ex == BIAS + LDBL_MAX_EXP)
+                       return (x + x); /* Inf, NaN, or unsupported format */
+               return (x);             /* finite and already an integer */
+       }
+       sign = expsign >> 15;
+
+       /*
+        * The following code assumes that intermediate results are
+        * evaluated in long double precision. If they are evaluated in
+        * greater precision, double rounding may occur, and if they are
+        * evaluated in less precision (as on i386), results will be
+        * wildly incorrect.
+        */
+       x += shift[sign];
+       x -= shift[sign];
+
+       /*
+        * If the result is +-0, then it must have the same sign as x, but
+        * the above calculation doesn't always give this.  Fix up the sign.
+        */
+       if (ex < BIAS && x == 0.0L)
+               return (zero[sign]);
+
+       return (x);
+}
index 3bdc7c7..8798a6f 100644 (file)
  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
  * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  *
- * $FreeBSD: head/lib/msun/src/s_round.c 153017 2005-12-02 13:45:06Z bde $
+ * $FreeBSD: head/lib/msun/src/s_round.c 257823 2013-11-07 22:46:13Z kargl $
  */
 
-#include <math.h>
+
+#include <float.h>
+
+#include "math.h"
+#include "math_private.h"
 
 double
 round(double x)
 {
        double t;
+       uint32_t hx;
 
-       if (!isfinite(x))
-               return (x);
+       GET_HIGH_WORD(hx, x);
+       if ((hx & 0x7fffffff) == 0x7ff00000)
+               return (x + x);
 
-       if (x >= 0.0) {
+       if (!(hx & 0x80000000)) {
                t = floor(x);
                if (t - x <= -0.5)
-                       t += 1.0;
+                       t += 1;
                return (t);
        } else {
                t = floor(-x);
                if (t + x <= -0.5)
-                       t += 1.0;
+                       t += 1;
                return (-t);
        }
 }
+
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(round, roundl);
+#endif
index c602f78..5196b70 100644 (file)
  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
  * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  *
- * $FreeBSD: head/lib/msun/src/s_roundf.c 153017 2005-12-02 13:45:06Z bde $
+ * $FreeBSD: head/lib/msun/src/s_roundf.c 257770 2013-11-06 23:44:52Z kargl $
  */
 
-#include <math.h>
+
+#include "math.h"
+#include "math_private.h"
 
 float
 roundf(float x)
 {
        float t;
+       uint32_t hx;
 
-       if (!isfinite(x))
-               return (x);
+       GET_FLOAT_WORD(hx, x);
+       if ((hx & 0x7fffffff) == 0x7f800000)
+               return (x + x);
 
-       if (x >= 0.0) {
+       if (!(hx & 0x80000000)) {
                t = floorf(x);
-               if (t - x <= -0.5)
-                       t += 1.0;
+               if (t - x <= -0.5F)
+                       t += 1;
                return (t);
        } else {
                t = floorf(-x);
-               if (t + x <= -0.5)
-                       t += 1.0;
+               if (t + x <= -0.5F)
+                       t += 1;
                return (-t);
        }
 }
index 07bba13..c22b64c 100644 (file)
  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
  * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  *
- * $FreeBSD: head/lib/msun/src/s_roundl.c 153017 2005-12-02 13:45:06Z bde $
+ * $FreeBSD: head/lib/msun/src/s_roundl.c 257770 2013-11-06 23:44:52Z kargl $
  */
 
-#include <math.h>
+
+#include <float.h>
+#ifdef __i386__
+#include <ieeefp.h>
+#endif
+
+#include "fpmath.h"
+#include "math.h"
+#include "math_private.h"
 
 long double
 roundl(long double x)
 {
        long double t;
+       uint16_t hx;
+
+       GET_LDBL_EXPSIGN(hx, x);
+       if ((hx & 0x7fff) == 0x7fff)
+               return (x + x);
 
-       if (!isfinite(x))
-               return (x);
+       ENTERI();
 
-       if (x >= 0.0) {
+       if (!(hx & 0x8000)) {
                t = floorl(x);
-               if (t - x <= -0.5)
-                       t += 1.0;
-               return (t);
+               if (t - x <= -0.5L)
+                       t += 1;
+               RETURNI(t);
        } else {
                t = floorl(-x);
-               if (t + x <= -0.5)
-                       t += 1.0;
-               return (-t);
+               if (t + x <= -0.5L)
+                       t += 1;
+               RETURNI(-t);
        }
 }
diff --git a/lib/libm/src/s_scalbn.c b/lib/libm/src/s_scalbn.c
new file mode 100644 (file)
index 0000000..8df777b
--- /dev/null
@@ -0,0 +1,64 @@
+/* @(#)s_scalbn.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ *
+ * $FreeBSD: head/lib/msun/src/s_scalbn.c 143264 2005-03-07 21:27:37Z das $
+ */
+
+
+/*
+ * scalbn (double x, int n)
+ * scalbn(x,n) returns x* 2**n  computed by  exponent
+ * manipulation rather than by actually performing an
+ * exponentiation or a multiplication.
+ */
+
+#include <float.h>
+
+#include "math.h"
+#include "math_private.h"
+
+static const double
+two54   =  1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
+twom54  =  5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */
+huge   = 1.0e+300,
+tiny   = 1.0e-300;
+
+double
+scalbn (double x, int n)
+{
+       int32_t k,hx,lx;
+       EXTRACT_WORDS(hx,lx,x);
+        k = (hx&0x7ff00000)>>20;               /* extract exponent */
+        if (k==0) {                            /* 0 or subnormal x */
+            if ((lx|(hx&0x7fffffff))==0) return x; /* +-0 */
+           x *= two54;
+           GET_HIGH_WORD(hx,x);
+           k = ((hx&0x7ff00000)>>20) - 54;
+            if (n< -50000) return tiny*x;      /*underflow*/
+           }
+        if (k==0x7ff) return x+x;              /* NaN or Inf */
+        k = k+n;
+        if (k >  0x7fe) return huge*copysign(huge,x); /* overflow  */
+        if (k > 0)                             /* normal result */
+           {SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20)); return x;}
+        if (k <= -54)
+            if (n > 50000)     /* in case integer overflow in n+k */
+               return huge*copysign(huge,x);   /*overflow*/
+           else return tiny*copysign(tiny,x);  /*underflow*/
+        k += 54;                               /* subnormal result */
+       SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20));
+        return x*twom54;
+}
+
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(scalbn, ldexpl);
+__weak_reference(scalbn, scalbnl);
+#endif
diff --git a/lib/libm/src/s_scalbnf.c b/lib/libm/src/s_scalbnf.c
new file mode 100644 (file)
index 0000000..794d6de
--- /dev/null
@@ -0,0 +1,55 @@
+/* s_scalbnf.c -- float version of s_scalbn.c.
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ */
+
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * 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.
+ * ====================================================
+ *
+ * $FreeBSD: head/lib/msun/src/s_scalbnf.c 143205 2005-03-07 04:52:43Z das $
+ */
+
+
+#include "math.h"
+#include "math_private.h"
+
+static const float
+two25   =  3.355443200e+07,    /* 0x4c000000 */
+twom25  =  2.9802322388e-08,   /* 0x33000000 */
+huge   = 1.0e+30,
+tiny   = 1.0e-30;
+
+float
+scalbnf (float x, int n)
+{
+       int32_t k,ix;
+       GET_FLOAT_WORD(ix,x);
+        k = (ix&0x7f800000)>>23;               /* extract exponent */
+        if (k==0) {                            /* 0 or subnormal x */
+            if ((ix&0x7fffffff)==0) return x; /* +-0 */
+           x *= two25;
+           GET_FLOAT_WORD(ix,x);
+           k = ((ix&0x7f800000)>>23) - 25;
+            if (n< -50000) return tiny*x;      /*underflow*/
+           }
+        if (k==0xff) return x+x;               /* NaN or Inf */
+        k = k+n;
+        if (k >  0xfe) return huge*copysignf(huge,x); /* overflow  */
+        if (k > 0)                             /* normal result */
+           {SET_FLOAT_WORD(x,(ix&0x807fffff)|(k<<23)); return x;}
+        if (k <= -25)
+            if (n > 50000)     /* in case integer overflow in n+k */
+               return huge*copysignf(huge,x);  /*overflow*/
+           else return tiny*copysignf(tiny,x); /*underflow*/
+        k += 25;                               /* subnormal result */
+       SET_FLOAT_WORD(x,(ix&0x807fffff)|(k<<23));
+        return x*twom25;
+}
+
+__strong_reference(scalbnf, ldexpf);
diff --git a/lib/libm/src/s_scalbnl.c b/lib/libm/src/s_scalbnl.c
new file mode 100644 (file)
index 0000000..57f0d81
--- /dev/null
@@ -0,0 +1,69 @@
+/* @(#)s_scalbn.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ *
+ * $FreeBSD: head/lib/msun/src/s_scalbnl.c 143206 2005-03-07 04:52:58Z das
+ */
+
+
+/*
+ * scalbnl (long double x, int n)
+ * scalbnl(x,n) returns x* 2**n  computed by  exponent
+ * manipulation rather than by actually performing an
+ * exponentiation or a multiplication.
+ */
+
+/*
+ * We assume that a long double has a 15-bit exponent.  On systems
+ * where long double is the same as double, scalbnl() is an alias
+ * for scalbn(), so we don't use this routine.
+ */
+
+#include <float.h>
+#include <math.h>
+
+#include "fpmath.h"
+
+#if LDBL_MAX_EXP != 0x4000
+#error "Unsupported long double format"
+#endif
+
+static const long double
+huge = 0x1p16000L,
+tiny = 0x1p-16000L;
+
+long double
+scalbnl (long double x, int n)
+{
+       union IEEEl2bits u;
+       int k;
+       u.e = x;
+        k = u.bits.exp;                                /* extract exponent */
+        if (k==0) {                            /* 0 or subnormal x */
+            if ((u.bits.manh|u.bits.manl)==0) return x;        /* +-0 */
+           u.e *= 0x1p+128;
+           k = u.bits.exp - 128;
+            if (n< -50000) return tiny*x;      /*underflow*/
+           }
+        if (k==0x7fff) return x+x;             /* NaN or Inf */
+        k = k+n;
+        if (k >= 0x7fff) return huge*copysignl(huge,x); /* overflow  */
+        if (k > 0)                             /* normal result */
+           {u.bits.exp = k; return u.e;}
+        if (k <= -128)
+            if (n > 50000)     /* in case integer overflow in n+k */
+               return huge*copysign(huge,x);   /*overflow*/
+           else return tiny*copysign(tiny,x);  /*underflow*/
+        k += 128;                              /* subnormal result */
+       u.bits.exp = k;
+        return u.e*0x1p-128;
+}
+
+__strong_reference(scalbnl, ldexpl);
index 1514b35..62f268c 100644 (file)
@@ -1,5 +1,4 @@
 /* @(#)s_tanh.c 5.1 93/09/24 */
-/* $FreeBSD: head/lib/msun/src/s_tanh.c 176451 2008-02-22 02:30:36Z das $ */
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
@@ -11,6 +10,7 @@
  * ====================================================
  */
 
+
 /* Tanh(x)
  * Return the Hyperbolic Tangent of x
  *
  *     only tanh(0)=0 is exact for finite argument.
  */
 
+#include <float.h>
+
 #include "math.h"
 #include "math_private.h"
 
-static const double one = 1.0, two = 2.0, tiny = 1.0e-300, huge = 1.0e300;
+static const volatile double tiny = 1.0e-300;
+static const double one = 1.0, two = 2.0, huge = 1.0e300;
 
 double
 tanh(double x)
@@ -73,3 +76,7 @@ tanh(double x)
        }
        return (jx>=0)? z: -z;
 }
+
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(tanh, tanhl);
+#endif
index acbdf3f..24a8dc2 100644 (file)
  * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
- *
- * $FreeBSD: head/lib/msun/src/s_tanhf.c 176451 2008-02-22 02:30:36Z das $
  */
 
+
 #include "math.h"
 #include "math_private.h"
 
-static const float one=1.0, two=2.0, tiny = 1.0e-30, huge = 1.0e30;
+static const volatile float tiny = 1.0e-30;
+static const float one=1.0, two=2.0, huge = 1.0e30;
+
 float
 tanhf(float x)
 {
diff --git a/lib/libm/src/s_tanhl.c b/lib/libm/src/s_tanhl.c
new file mode 100644 (file)
index 0000000..c32b925
--- /dev/null
@@ -0,0 +1,170 @@
+/* from: FreeBSD: head/lib/msun/src/s_tanhl.c XXX */
+
+/* @(#)s_tanh.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+
+/*
+ * See s_tanh.c for complete comments.
+ *
+ * Converted to long double by Bruce D. Evans.
+ */
+
+#include <float.h>
+#ifdef __i386__
+#include <ieeefp.h>
+#endif
+
+#include "math.h"
+#include "math_private.h"
+#include "fpmath.h"
+#include "k_expl.h"
+
+#if LDBL_MAX_EXP != 0x4000
+/* We also require the usual expsign encoding. */
+#error "Unsupported long double format"
+#endif
+
+#define        BIAS    (LDBL_MAX_EXP - 1)
+
+static const volatile double tiny = 1.0e-300;
+static const double one = 1.0;
+#if LDBL_MANT_DIG == 64
+/*
+ * Domain [-0.25, 0.25], range ~[-1.6304e-22, 1.6304e-22]:
+ * |tanh(x)/x - t(x)| < 2**-72.3
+ */
+static const union IEEEl2bits
+T3u = LD80C(0xaaaaaaaaaaaaaa9f, -2, -3.33333333333333333017e-1L);
+#define        T3      T3u.e
+static const double
+T5  =  1.3333333333333314e-1,          /*  0x1111111111110a.0p-55 */
+T7  = -5.3968253968210485e-2,          /* -0x1ba1ba1ba1a1a1.0p-57 */
+T9  =  2.1869488531393817e-2,          /*  0x1664f488172022.0p-58 */
+T11 = -8.8632352345964591e-3,          /* -0x1226e34bc138d5.0p-59 */
+T13 =  3.5921169709993771e-3,          /*  0x1d6d371d3e400f.0p-61 */
+T15 = -1.4555786415756001e-3,          /* -0x17d923aa63814d.0p-62 */
+T17 =  5.8645267876296793e-4,          /*  0x13378589b85aa7.0p-63 */
+T19 = -2.1121033571392224e-4;          /* -0x1baf0af80c4090.0p-65 */
+#elif LDBL_MANT_DIG == 113
+/*
+ * Domain [-0.25, 0.25], range ~[-2.4211e-37, 2.4211e-37]:
+ * |tanh(x)/x - t(x)| < 2**121.6
+ */
+static const long double
+T3 = -3.33333333333333333333333333333332980e-1L,       /* -0x1555555555555555555555555554e.0p-114L */
+T5  =  1.33333333333333333333333333332707260e-1L,      /*  0x1111111111111111111111110ab7b.0p-115L */
+T7  = -5.39682539682539682539682535723482314e-2L,      /* -0x1ba1ba1ba1ba1ba1ba1ba17b5fc98.0p-117L */
+T9  =  2.18694885361552028218693591149061717e-2L,      /*  0x1664f4882c10f9f32d6b1a12a25e5.0p-118L */
+T11 = -8.86323552990219656883762347736381851e-3L,      /* -0x1226e355e6c23c8f5a5a0f386cb4d.0p-119L */
+T13 =  3.59212803657248101358314398220822722e-3L,      /*  0x1d6d3d0e157ddfb403ad3637442c6.0p-121L */
+T15 = -1.45583438705131796512568010348874662e-3L;      /* -0x17da36452b75e150c44cc34253b34.0p-122L */
+static const double
+T17 =  5.9002744094556621e-4,          /*  0x1355824803668e.0p-63 */
+T19 = -2.3912911424260516e-4,          /* -0x1f57d7734c8dde.0p-65 */
+T21 =  9.6915379535512898e-5,          /*  0x1967e18ad6a6ca.0p-66 */
+T23 = -3.9278322983156353e-5,          /* -0x1497d8e6b75729.0p-67 */
+T25 =  1.5918887220143869e-5,          /*  0x10b1319998cafa.0p-68 */
+T27 = -6.4514295231630956e-6,          /* -0x1b0f2b71b218eb.0p-70 */
+T29 =  2.6120754043964365e-6,          /*  0x15e963a3cf3a39.0p-71 */
+T31 = -1.0407567231003314e-6,          /* -0x1176041e656869.0p-72 */
+T33 =  3.4744117554063574e-7;          /*  0x1750fe732cab9c.0p-74 */
+#endif /* LDBL_MANT_DIG == 64 */
+
+static inline long double
+divl(long double a, long double b, long double c, long double d,
+    long double e, long double f)
+{
+       long double inv, r;
+       float fr, fw;
+
+       _2sumF(a, c);
+       b = b + c;
+       _2sumF(d, f);
+       e = e + f;
+
+       inv = 1 / (d + e);
+
+       r = (a + b) * inv;
+       fr = r;
+       r = fr;
+
+       fw = d + e;
+       e = d - fw + e;
+       d = fw;
+
+       r = r + (a - d * r + b - e * r) * inv;
+
+       return r;
+}
+
+long double
+tanhl(long double x)
+{
+       long double hi,lo,s,x2,x4,z;
+       double dx2;
+       int16_t jx,ix;
+
+       GET_LDBL_EXPSIGN(jx,x);
+       ix = jx&0x7fff;
+
+    /* x is INF or NaN */
+       if(ix>=0x7fff) {
+           if (jx>=0) return one/x+one;    /* tanh(+-inf)=+-1 */
+           else       return one/x-one;    /* tanh(NaN) = NaN */
+       }
+
+       ENTERI();
+
+    /* |x| < 40 */
+       if (ix < 0x4004 || fabsl(x) < 40) {     /* |x|<40 */
+           if (__predict_false(ix<BIAS-(LDBL_MANT_DIG+1)/2)) { /* |x|<TINY */
+               /* tanh(+-0) = +0; tanh(tiny) = tiny(-+) with inexact: */
+               return (x == 0 ? x : (0x1p200 * x - x) * 0x1p-200);
+           }
+           if (ix<0x3ffd) {            /* |x|<0.25 */
+               x2 = x*x;
+#if LDBL_MANT_DIG == 64
+               x4 = x2*x2;
+               RETURNI(((T19*x2 + T17)*x4 + (T15*x2 + T13))*(x2*x*x2*x4*x4) +
+                   ((T11*x2 + T9)*x4 + (T7*x2 + T5))*(x2*x*x2) +
+                   T3*(x2*x) + x);
+#elif LDBL_MANT_DIG == 113
+               dx2 = x2;
+#if 0
+               RETURNI(((((((((((((((T33*dx2 + T31)*dx2 + T29)*dx2 + T27)*dx2 +
+                   T25)*x2 + T23)*x2 + T21)*x2 + T19)*x2 + T17)*x2 +
+                   T15)*x2 + T13)*x2 + T11)*x2 + T9)*x2 + T7)*x2 + T5)*
+                   (x2*x*x2) +
+                   T3*(x2*x) + x);
+#else
+               long double q = ((((((((((((((T33*dx2 + T31)*dx2 + T29)*dx2 + T27)*dx2 +
+                   T25)*x2 + T23)*x2 + T21)*x2 + T19)*x2 + T17)*x2 +
+                   T15)*x2 + T13)*x2 + T11)*x2 + T9)*x2 + T7)*x2 + T5)*
+                   (x2*x*x2);
+               RETURNI(q + T3*(x2*x) + x);
+#endif
+#endif
+           }
+           k_hexpl(2*fabsl(x), &hi, &lo);
+           if (ix<0x4001 && fabsl(x) < 1.5)    /* |x|<1.5 */
+               z = divl(hi, lo, -0.5, hi, lo, 0.5);
+           else
+               z = one - one/(lo+0.5+hi);
+    /* |x| >= 40, return +-1 */
+       } else {
+           z = one - tiny;             /* raise inexact flag */
+       }
+       s = 1;
+       if (jx<0) s = -1;
+       RETURNI(s*z);
+}