libm: Add several new functions and symbol versioning
authorJohn Marino <draco@marino.st>
Tue, 11 Jun 2013 22:28:04 +0000 (00:28 +0200)
committerJohn Marino <draco@marino.st>
Wed, 12 Jun 2013 06:15:53 +0000 (08:15 +0200)
The following long double functions were added to the math library:
  logl
  log2l
  log10l
  log1pl
  expm1l
  acoshl
  asinhl
  atanhl

In addition, the FreeBSD functionality that creates symbol versioning
for libraries was adapted for FreeBSD.  The first version is called
"DFLY36.0".  If it is necessary to create a new version of the 3.5 or
3.6 branch, the number after the decimal will be incremented.  The 3.7
branch will start with "DFLY38.0" if it needs its own version.

libm was baselined with all symbols being the same version: DFLY36.0.
With symbol versioning, it will not be necessary to increment the major
version anymore, so this library shall always be known as libm.so.4 from
this point on.

33 files changed:
lib/libc/Versions.def [new file with mode: 0644]
lib/libm/Makefile
lib/libm/Symbol.map [new file with mode: 0644]
lib/libm/i386/Makefile.inc
lib/libm/i386/Symbol.map [new file with mode: 0644]
lib/libm/ld80/s_expl.c
lib/libm/ld80/s_logl.c [new file with mode: 0644]
lib/libm/man/acosh.3
lib/libm/man/asinh.3
lib/libm/man/atanh.3
lib/libm/man/exp.3
lib/libm/man/log.3
lib/libm/man/sin.3
lib/libm/src/catrig.c
lib/libm/src/catrigf.c
lib/libm/src/e_acosh.c
lib/libm/src/e_acoshl.c [new file with mode: 0644]
lib/libm/src/e_atanh.c
lib/libm/src/e_atanhl.c [new file with mode: 0644]
lib/libm/src/e_log.c
lib/libm/src/e_log10.c
lib/libm/src/e_log2.c
lib/libm/src/math.h
lib/libm/src/math_private.h
lib/libm/src/s_asinh.c
lib/libm/src/s_asinhl.c [new file with mode: 0644]
lib/libm/src/s_expm1.c
lib/libm/src/s_log1p.c
lib/libm/x86_64/Makefile.inc
lib/libm/x86_64/Symbol.map [new file with mode: 0644]
share/mk/bsd.lib.mk
share/mk/bsd.symver.mk [new file with mode: 0644]
share/mk/version_gen.awk [new file with mode: 0644]

diff --git a/lib/libc/Versions.def b/lib/libc/Versions.def
new file mode 100644 (file)
index 0000000..d504116
--- /dev/null
@@ -0,0 +1,19 @@
+#
+# Note: Whenever bumping the DFLY version, always make
+#       DFLYprivate_1.0 depend on the new DFLY version.
+#       This will keep it at the end of the dependency chain.
+#
+
+# This is our first version; it depends on no other.
+# This version was first added to DragonFly 3.5 
+DFLY36.0 {
+};
+
+# 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
+# simple ABI-checking tools.
+#
+# Please do NOT increment the version of this namespace.
+DFLYprivate_1.0 {
+} DFLY36.0;
index 240fef6..a6700da 100644 (file)
@@ -62,16 +62,21 @@ COMMON_SRCS= b_exp.c b_log.c b_tgamma.c \
 # Location of fpmath.h and _fpmath.h
 LIBCDIR=${.CURDIR}/../libc
 CFLAGS+=-I${.CURDIR}/src -I${LIBCDIR}/include -I${LIBCDIR}/${MACHINE_ARCH}
+SYM_MAPS+=     ${.CURDIR}/Symbol.map
+
+VERSION_DEF=   ${LIBCDIR}/Versions.def
+SYMBOL_MAPS=   ${SYM_MAPS}
 
 # C99 long double functions
 COMMON_SRCS+=  s_copysignl.c s_fabsl.c s_modfl.c
 # If long double != double use these; otherwise, we alias the double versions.
-COMMON_SRCS+=  e_acosl.c e_asinl.c e_atan2l.c e_fmodl.c \
-       e_hypotl.c \
+COMMON_SRCS+=  e_acoshl.c e_acosl.c e_asinl.c e_atan2l.c e_atanhl.c \
+       e_fmodl.c e_hypotl.c \
        invtrig.c k_cosl.c k_sinl.c k_tanl.c \
-       s_atanl.c s_cbrtl.c s_ceill.c s_cosl.c s_cprojl.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_nanl.c s_nextafterl.c s_nexttoward.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
 
 # C99 complex functions
@@ -107,10 +112,12 @@ MAN=      acos.3 acosh.3 asin.3 asinh.3 atan.3 atan2.3 atanh.3 \
 
 MLINKS+=acos.3 acosf.3 \
        acos.3 acosl.3
-MLINKS+=acosh.3 acoshf.3
+MLINKS+=acosh.3 acoshf.3 \
+       acosh.3 acoshl.3
 MLINKS+=asin.3 asinf.3 \
        asin.3 asinl.3
-MLINKS+=asinh.3 asinhf.3
+MLINKS+=asinh.3 asinhf.3 \
+       asinh.3 asinhl.3
 MLINKS+=atan.3 atanf.3 \
        atan.3 atanl.3
 MLINKS+=atan2.3 atan2f.3 \
@@ -118,7 +125,8 @@ MLINKS+=atan2.3 atan2f.3 \
        atan2.3 carg.3 \
        atan2.3 cargf.3 \
        atan2.3 cargl.3
-MLINKS+=atanh.3 atanhf.3
+MLINKS+=atanh.3 atanhf.3 \
+       atanh.3 atanhl.3
 MLINKS+=cacos.3 cacosf.3 \
        cacos.3 cacosh.3 \
        cacos.3 cacoshf.3 \
@@ -171,6 +179,7 @@ MLINKS+=exp.3 expf.3 \
        exp.3 exp2l.3 \
        exp.3 expm1.3 \
        exp.3 expm1f.3 \
+       exp.3 expm1l.3 \
        exp.3 pow.3 \
        exp.3 powf.3
 MLINKS+=fabs.3 fabsf.3 \
@@ -237,10 +246,13 @@ MLINKS+=log.3 logf.3 \
        log.3 logl.3 \
        log.3 log10.3 \
        log.3 log10f.3 \
+       log.3 log10l.3 \
        log.3 log2.3 \
        log.3 log2f.3 \
+       log.3 log2l.e \
        log.3 log1p.3 \
-       log.3 log1pf.3
+       log.3 log1pf.3 \
+       log.3 log1pl.3
 MLINKS+=lrint.3 llrint.3 \
        lrint.3 llrintf.3 \
        lrint.3 llrintl.3 \
diff --git a/lib/libm/Symbol.map b/lib/libm/Symbol.map
new file mode 100644 (file)
index 0000000..6129c8e
--- /dev/null
@@ -0,0 +1,260 @@
+/*
+ * $FreeBSD: head/lib/msun/Symbol.map 251599 2013-06-10 06:04:58Z das $
+ */
+
+DFLY36.0 {
+       __fe_dfl_env;
+       tgamma;
+       acos;
+       acosf;
+       acosh;
+       acoshf;
+       asin;
+       asinf;
+       atan2;
+       atan2f;
+       atanh;
+       atanhf;
+       cosh;
+       coshf;
+       exp;
+       expf;
+       fmod;
+       fmodf;
+       gamma;
+       gamma_r;
+       gammaf;
+       gammaf_r;
+       hypot;
+       hypotf;
+       j0;
+       y0;
+       j0f;
+       y0f;
+       j1;
+       y1;
+       j1f;
+       y1f;
+       jn;
+       yn;
+       jnf;
+       ynf;
+       lgamma;
+       lgamma_r;
+       lgammaf;
+       lgammaf_r;
+       log;
+       log10;
+       log10f;
+       logf;
+       pow;
+       powf;
+       remainder;
+       remainderf;
+       scalb;
+       scalbf;
+       sinh;
+       sinhf;
+       sqrt;
+       sqrtf;
+       asinh;
+       asinhf;
+       atan;
+       atanf;
+       cbrt;
+       cbrtf;
+       ceil;
+       ceilf;
+       ceill;
+       cimag;
+       cimagf;
+       cimagl;
+       conj;
+       conjf;
+       conjl;
+       copysign;
+       copysignf;
+       copysignl;
+       cos;
+       cosf;
+       creal;
+       crealf;
+       creall;
+       erf;
+       erfc;
+       erff;
+       erfcf;
+       exp2;
+       exp2f;
+       expm1;
+       expm1f;
+       fabs;
+       fabsf;
+       fabsl;
+       fdim;
+       fdimf;
+       fdiml;
+       finite;
+       finitef;
+       floor;
+       floorf;
+       floorl;
+       fma;
+       fmaf;
+       fmal;
+       fmax;
+       fmaxf;
+       fmaxl;
+       fmin;
+       fminf;
+       fminl;
+       frexp;
+       frexpf;
+       frexpl;
+       ilogb;
+       ilogbf;
+       ilogbl;
+       __isfinite;
+       __isfinitef;
+       __isfinitel;
+       isnanf;
+       __isnanl;
+       __isnormal;
+       __isnormalf;
+       __isnormall;
+       llrint;
+       llrintf;
+       llround;
+       llroundf;
+       llroundl;
+       log1p;
+       log1pf;
+       logb;
+       logbf;
+       lrint;
+       lrintf;
+       lround;
+       lroundf;
+       lroundl;
+       modff;
+       modfl;
+       nearbyint;
+       nearbyintf;
+       nextafter;
+       nexttoward;
+       nexttowardl;
+       nextafterl;
+       nextafterf;
+       nexttowardf;
+       remquo;
+       remquof;
+       rint;
+       rintf;
+       round;
+       roundf;
+       roundl;
+       scalbln;
+       scalblnf;
+       scalblnl;
+       scalbn;
+       scalbnl;
+       scalbnf;
+       ldexpf;
+       ldexpl;
+       __signbit;
+       __signbitf;
+       __signbitl;
+       signgam;
+       significand;
+       significandf;
+       sin;
+       sinf;
+       tan;
+       tanf;
+       tanh;
+       tanhf;
+       trunc;
+       truncf;
+       truncl;
+       cabs;
+       cabsf;
+       drem;
+       dremf;
+       carg;
+       cargf;
+       csqrt;
+       csqrtf;
+       logbl;
+       nan;
+       nanf;
+       nanl;
+       llrintl;
+       lrintl;
+       nearbyintl;
+       rintl;
+       exp2l;
+       sinl;
+       cosl;
+       tanl;
+       tgammaf;
+       sqrtl;
+       hypotl;
+       cabsl;
+       csqrtl;
+       remquol;
+       remainderl;
+       fmodl;
+       acosl;
+       asinl;
+       atan2l;
+       atanl;
+       cargl;
+       cproj;
+       cprojf;
+       cprojl;
+       __isnanf;
+       cbrtl;
+       cexp;
+       cexpf;
+       log2;
+       log2f;
+       feclearexcept;
+       fegetexceptflag;
+       fetestexcept;
+       fegetround;
+       fesetround;
+       fesetenv;
+       acoshl;
+       asinhl;
+       atanhl;
+       cacos;
+       cacosf;
+       cacosh;
+       cacoshf;
+       casin;
+       casinf;
+       casinh;
+       casinhf;
+       catan;
+       catanf;
+       catanh;
+       catanhf;
+       csin;
+       csinf;
+       csinh;
+       csinhf;
+       ccos;
+       ccosf;
+       ccosh;
+       ccoshf;
+       ctan;
+       ctanf;
+       ctanh;
+       ctanhf;
+       expl;
+       expm1l;
+       log10l;
+       log1pl;
+       log2l;
+       logl;
+};
index 6769cee..9397a29 100644 (file)
@@ -16,3 +16,5 @@ ARCH_SRCS+= e_log10f.S e_logf.S e_remainderf.S \
 ARCH_SRCS+= e_remainderl.S e_sqrtl.S s_ceill.S s_copysignl.S \
            s_floorl.S s_llrintl.S \
            s_logbl.S s_lrintl.S s_remquol.S s_rintl.S s_scalbnl.S s_truncl.S
+
+SYM_MAPS += ${.CURDIR}/i386/Symbol.map
diff --git a/lib/libm/i386/Symbol.map b/lib/libm/i386/Symbol.map
new file mode 100644 (file)
index 0000000..aa0ebc5
--- /dev/null
@@ -0,0 +1,14 @@
+/*
+ * $FreeBSD: head/lib/msun/i387/Symbol.map 226218 2011-10-10 15:43:09Z das $
+ */
+DFLY36.0 {
+       __has_sse;
+       __test_sse;
+       fesetexceptflag;
+       feraiseexcept;
+       fegetenv;
+       feholdexcept;
+       feupdateenv;
+       feenableexcept;
+       fedisableexcept;
+};
index 48a92ab..3155ed4 100644 (file)
@@ -1,5 +1,5 @@
 /*-
- * Copyright (c) 2009-2012 Steven G. Kargl
+ * Copyright (c) 2009-2013 Steven G. Kargl
  * All rights reserved.
  *
  * Redistribution and use in source and binary forms, with or without
  *
  * Optimized by Bruce D. Evans.
  *
- * $FreeBSD: head/lib/msun/ld80/s_expl.c 241516 2012-10-13 19:53:11Z kargl $
+ * $FreeBSD: head/lib/msun/ld80/s_expl.c 251343 2013-06-03 19:51:32Z kargl $
  */
 
-/*-
+/**
  * Compute the exponential of x for Intel 80-bit format.  This is based on:
  *
  *   PTP Tang, "Table-driven implementation of the exponential function
@@ -49,6 +49,7 @@
 #include "math_private.h"
 
 #define        INTERVALS       128
+#define        LOG2_INTERVALS  7
 #define        BIAS    (LDBL_MAX_EXP - 1)
 
 static const long double
@@ -59,9 +60,12 @@ static volatile const long double tiny = 0x1p-10000L;
 
 static const union IEEEl2bits
 /* log(2**16384 - 0.5) rounded towards zero: */
-o_threshold = LD80C(0xb17217f7d1cf79ab, 13,  11356.5234062941439488L),
+/* 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_threshold = LD80C(0xb21dfe7f09e2baa9, 13, -11399.4985314888605581L);
+u_thresholdu = LD80C(0xb21dfe7f09e2baa9, 13, -11399.4985314888605581L);
+#define u_threshold     (u_thresholdu.e)
 
 static const double
 /*
@@ -77,11 +81,11 @@ L2 = -3.2819649005320973e-13,               /* -0x1718432a1b0e26.0p-94 */
  * |exp(x) - p(x)| < 2**-77.2
  * (0.002708 is ln2/(2*INTERVALS) rounded up a little).
  */
-P2 =  0.5,
-P3 =  1.6666666666666119e-1,           /*  0x15555555555490.0p-55 */
-P4 =  4.1666666666665887e-2,           /*  0x155555555554e5.0p-57 */
-P5 =  8.3333354987869413e-3,           /*  0x1111115b789919.0p-59 */
-P6 =  1.3888891738560272e-3;           /*  0x16c16c651633ae.0p-62 */
+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
@@ -95,8 +99,7 @@ P6 =  1.3888891738560272e-3;          /*  0x16c16c651633ae.0p-62 */
 static const struct {
        double  hi;
        double  lo;
-/* XXX should rename 's'. */
-} s[INTERVALS] = {
+} tbl[INTERVALS] = {
        0x1p+0, 0x0p+0,
        0x1.0163da9fb3335p+0, 0x1.b61299ab8cdb7p-54,
        0x1.02c9a3e778060p+0, 0x1.dcdef95949ef4p-53,
@@ -231,7 +234,8 @@ long double
 expl(long double x)
 {
        union IEEEl2bits u, v;
-       long double fn, q, r, r1, r2, t, t23, t45, twopk, twopkp10000, z;
+       long double fn, q, r, r1, r2, t, twopk, twopkp10000;
+       long double z;
        int k, n, n2;
        uint16_t hx, ix;
 
@@ -241,23 +245,21 @@ expl(long double x)
        ix = hx & 0x7fff;
        if (ix >= BIAS + 13) {          /* |x| >= 8192 or x is NaN */
                if (ix == BIAS + LDBL_MAX_EXP) {
-                       if (hx & 0x8000 && u.xbits.man == 1ULL << 63)
-                               return (0.0L);  /* x is -Inf */
-                       return (x + x); /* x is +Inf, NaN or unsupported */
+                       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.e)
+               if (x > o_threshold)
                        return (huge * huge);
-               if (x < u_threshold.e)
+               if (x < u_threshold)
                        return (tiny * tiny);
-       } else if (ix < BIAS - 66) {    /* |x| < 0x1p-66 */
-                                       /* includes pseudo-denormals */
-               if (huge + x > 1.0L)    /* trigger inexact iff x != 0 */
-                       return (1.0L + x);
+       } else if (ix < BIAS - 65) {    /* |x| < 0x1p-65 (includes pseudos) */
+               return (1 + x);         /* 1 with inexact iff x != 0 */
        }
 
        ENTERI();
 
-       /* Reduce x to (k*ln2 + midpoint[n2] + r1 + r2). */
+       /* 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. */
@@ -269,12 +271,13 @@ expl(long double x)
        n  = (int)fn;
 #endif
        n2 = (unsigned)n % INTERVALS;
-       k = (n - n2) / INTERVALS;
+       /* Depend on the sign bit being propagated: */
+       k = n >> LOG2_INTERVALS;
        r1 = x - fn * L1;
-       r2 = -fn * L2;
+       r2 = fn * -L2;
 
        /* Prepare scale factors. */
-       v.xbits.man = 1ULL << 63;
+       v.e = 1;
        if (k >= LDBL_MIN_EXP) {
                v.xbits.expsign = BIAS + k;
                twopk = v.e;
@@ -283,21 +286,183 @@ expl(long double x)
                twopkp10000 = v.e;
        }
 
-       /* Evaluate expl(midpoint[n2] + r1 + r2) = s[n2] * expl(r1 + r2). */
-       /* Here q = q(r), not q(r1), since r1 is lopped like L1. */
-       t45 = r * P5 + P4;
+       /* Evaluate expl(endpoint[n2] + r1 + r2) = tbl[n2] * expl(r1 + r2). */
        z = r * r;
-       t23 = r * P3 + P2;
-       q = r2 + z * t23 + z * z * t45 + z * z * z * P6;
-       t = (long double)s[n2].lo + s[n2].hi;
-       t = s[n2].lo + t * (q + r1) + s[n2].hi;
+       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;
 
        /* Scale by 2**k. */
        if (k >= LDBL_MIN_EXP) {
                if (k == LDBL_MAX_EXP)
-                       RETURNI(t * 2.0L * 0x1p16383L);
+                       RETURNI(t * 2 * 0x1p16383L);
                RETURNI(t * twopk);
        } else {
                RETURNI(t * twopkp10000 * twom10000);
        }
 }
+
+/**
+ * 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).
+ */
+
+/*
+ * 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) */
+
+/*
+ * 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
+ */
+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)
+{
+       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);
+               }
+
+               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;
+
+       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);
+       }
+
+       v.xbits.expsign = BIAS - k;
+       twomk = v.e;
+
+       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);
+}
diff --git a/lib/libm/ld80/s_logl.c b/lib/libm/ld80/s_logl.c
new file mode 100644 (file)
index 0000000..e8b41ca
--- /dev/null
@@ -0,0 +1,716 @@
+/*-
+ * Copyright (c) 2007-2013 Bruce D. Evans
+ * 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/ld80/s_logl.c 251292 2013-06-03 09:14:31Z das $
+ */
+
+/**
+ * Implementation of the natural logarithm of x for Intel 80-bit format.
+ *
+ * First decompose x into its base 2 representation:
+ *
+ *    log(x) = log(X * 2**k), where X is in [1, 2)
+ *           = log(X) + k * log(2).
+ *
+ * Let X = X_i + e, where X_i is the center of one of the intervals
+ * [-1.0/256, 1.0/256), [1.0/256, 3.0/256), .... [2.0-1.0/256, 2.0+1.0/256)
+ * and X is in this interval.  Then
+ *
+ *    log(X) = log(X_i + e)
+ *           = log(X_i * (1 + e / X_i))
+ *           = log(X_i) + log(1 + e / X_i).
+ *
+ * The values log(X_i) are tabulated below.  Let d = e / X_i and use
+ *
+ *    log(1 + d) = p(d)
+ *
+ * where p(d) = d - 0.5*d*d + ... is a special minimax polynomial of
+ * suitably high degree.
+ *
+ * To get sufficiently small roundoff errors, k * log(2), log(X_i), and
+ * sometimes (if |k| is not large) the first term in p(d) must be evaluated
+ * and added up in extra precision.  Extra precision is not needed for the
+ * rest of p(d).  In the worst case when k = 0 and log(X_i) is 0, the final
+ * error is controlled mainly by the error in the second term in p(d).  The
+ * error in this term itself is at most 0.5 ulps from the d*d operation in
+ * it.  The error in this term relative to the first term is thus at most
+ * 0.5 * |-0.5| * |d| < 1.0/1024 ulps.  We aim for an accumulated error of
+ * at most twice this at the point of the final rounding step.  Thus the
+ * final error should be at most 0.5 + 1.0/512 = 0.5020 ulps.  Exhaustive
+ * testing of a float variant of this function showed a maximum final error
+ * of 0.5008 ulps.  Non-exhaustive testing of a double variant of this
+ * function showed a maximum final error of 0.5078 ulps (near 1+1.0/256).
+ *
+ * We made the maximum of |d| (and thus the total relative error and the
+ * degree of p(d)) small by using a large number of intervals.  Using
+ * centers of intervals instead of endpoints reduces this maximum by a
+ * factor of 2 for a given number of intervals.  p(d) is special only
+ * in beginning with the Taylor coefficients 0 + 1*d, which tends to happen
+ * naturally.  The most accurate minimax polynomial of a given degree might
+ * be different, but then we wouldn't want it since we would have to do
+ * extra work to avoid roundoff error (especially for P0*d instead of d).
+ */
+
+#ifdef DEBUG
+#include <assert.h>
+#include <fenv.h>
+#endif
+
+#ifdef __i386__
+#include <ieeefp.h>
+#endif
+
+#include "fpmath.h"
+#include "math.h"
+#define        i386_SSE_GOOD
+#ifndef NO_STRUCT_RETURN
+#define        STRUCT_RETURN
+#endif
+#include "math_private.h"
+
+#if !defined(NO_UTAB) && !defined(NO_UTABL)
+#define        USE_UTAB
+#endif
+
+/*
+ * Domain [-0.005280, 0.004838], range ~[-5.1736e-22, 5.1738e-22]:
+ * |log(1 + d)/d - p(d)| < 2**-70.7
+ */
+static const double
+P2 = -0.5,
+P3 =  3.3333333333333359e-1,           /*  0x1555555555555a.0p-54 */
+P4 = -2.5000000000004424e-1,           /* -0x1000000000031d.0p-54 */
+P5 =  1.9999999992970016e-1,           /*  0x1999999972f3c7.0p-55 */
+P6 = -1.6666666072191585e-1,           /* -0x15555548912c09.0p-55 */
+P7 =  1.4286227413310518e-1,           /*  0x12494f9d9def91.0p-55 */
+P8 = -1.2518388626763144e-1;           /* -0x1006068cc0b97c.0p-55 */
+
+static volatile const double zero = 0;
+
+#define        INTERVALS       128
+#define        LOG2_INTERVALS  7
+#define        TSIZE           (INTERVALS + 1)
+#define        G(i)            (T[(i)].G)
+#define        F_hi(i)         (T[(i)].F_hi)
+#define        F_lo(i)         (T[(i)].F_lo)
+#define        ln2_hi          F_hi(TSIZE - 1)
+#define        ln2_lo          F_lo(TSIZE - 1)
+#define        E(i)            (U[(i)].E)
+#define        H(i)            (U[(i)].H)
+
+static const struct {
+       float   G;                      /* 1/(1 + i/128) rounded to 8/9 bits */
+       float   F_hi;                   /* log(1 / G_i) rounded (see below) */
+       double  F_lo;                   /* next 53 bits for log(1 / G_i) */
+} T[TSIZE] = {
+       /*
+        * ln2_hi and each F_hi(i) are rounded to a number of bits that
+        * makes F_hi(i) + dk*ln2_hi exact for all i and all dk.
+        *
+        * The last entry (for X just below 2) is used to define ln2_hi
+        * and ln2_lo, to ensure that F_hi(i) and F_lo(i) cancel exactly
+        * with dk*ln2_hi and dk*ln2_lo, respectively, when dk = -1.
+        * This is needed for accuracy when x is just below 1.  (To avoid
+        * special cases, such x are "reduced" strangely to X just below
+        * 2 and dk = -1, and then the exact cancellation is needed
+        * because any the error from any non-exactness would be too
+        * large).
+        *
+        * We want to share this table between double precision and ld80,
+        * so the relevant range of dk is the larger one of ld80
+        * ([-16445, 16383]) and the relevant exactness requirement is
+        * the stricter one of double precision.  The maximum number of
+        * bits in F_hi(i) that works is very dependent on i but has
+        * a minimum of 33.  We only need about 12 bits in F_hi(i) for
+        * it to provide enough extra precision in double precision (11
+        * more than that are required for ld80).
+        *
+        * We round F_hi(i) to 24 bits so that it can have type float,
+        * mainly to minimize the size of the table.  Using all 24 bits
+        * in a float for it automatically satisfies the above constraints.
+        */
+        0x800000.0p-23,  0,               0,
+        0xfe0000.0p-24,  0x8080ac.0p-30, -0x14ee431dae6675.0p-84,
+        0xfc0000.0p-24,  0x8102b3.0p-29, -0x1db29ee2d83718.0p-84,
+        0xfa0000.0p-24,  0xc24929.0p-29,  0x1191957d173698.0p-83,
+        0xf80000.0p-24,  0x820aec.0p-28,  0x13ce8888e02e79.0p-82,
+        0xf60000.0p-24,  0xa33577.0p-28, -0x17a4382ce6eb7c.0p-82,
+        0xf48000.0p-24,  0xbc42cb.0p-28, -0x172a21161a1076.0p-83,
+        0xf30000.0p-24,  0xd57797.0p-28, -0x1e09de07cb9589.0p-82,
+        0xf10000.0p-24,  0xf7518e.0p-28,  0x1ae1eec1b036c5.0p-91,
+        0xef0000.0p-24,  0x8cb9df.0p-27, -0x1d7355325d560e.0p-81,
+        0xed8000.0p-24,  0x999ec0.0p-27, -0x1f9f02d256d503.0p-82,
+        0xec0000.0p-24,  0xa6988b.0p-27, -0x16fc0a9d12c17a.0p-83,
+        0xea0000.0p-24,  0xb80698.0p-27,  0x15d581c1e8da9a.0p-81,
+        0xe80000.0p-24,  0xc99af3.0p-27, -0x1535b3ba8f150b.0p-83,
+        0xe70000.0p-24,  0xd273b2.0p-27,  0x163786f5251af0.0p-85,
+        0xe50000.0p-24,  0xe442c0.0p-27,  0x1bc4b2368e32d5.0p-84,
+        0xe38000.0p-24,  0xf1b83f.0p-27,  0x1c6090f684e676.0p-81,
+        0xe20000.0p-24,  0xff448a.0p-27, -0x1890aa69ac9f42.0p-82,
+        0xe08000.0p-24,  0x8673f6.0p-26,  0x1b9985194b6b00.0p-80,
+        0xdf0000.0p-24,  0x8d515c.0p-26, -0x1dc08d61c6ef1e.0p-83,
+        0xdd8000.0p-24,  0x943a9e.0p-26, -0x1f72a2dac729b4.0p-82,
+        0xdc0000.0p-24,  0x9b2fe6.0p-26, -0x1fd4dfd3a0afb9.0p-80,
+        0xda8000.0p-24,  0xa2315d.0p-26, -0x11b26121629c47.0p-82,
+        0xd90000.0p-24,  0xa93f2f.0p-26,  0x1286d633e8e569.0p-81,
+        0xd78000.0p-24,  0xb05988.0p-26,  0x16128eba936770.0p-84,
+        0xd60000.0p-24,  0xb78094.0p-26,  0x16ead577390d32.0p-80,
+        0xd50000.0p-24,  0xbc4c6c.0p-26,  0x151131ccf7c7b7.0p-81,
+        0xd38000.0p-24,  0xc3890a.0p-26, -0x115e2cd714bd06.0p-80,
+        0xd20000.0p-24,  0xcad2d7.0p-26, -0x1847f406ebd3b0.0p-82,
+        0xd10000.0p-24,  0xcfb620.0p-26,  0x1c2259904d6866.0p-81,
+        0xcf8000.0p-24,  0xd71653.0p-26,  0x1ece57a8d5ae55.0p-80,
+        0xce0000.0p-24,  0xde843a.0p-26, -0x1f109d4bc45954.0p-81,
+        0xcd0000.0p-24,  0xe37fde.0p-26,  0x1bc03dc271a74d.0p-81,
+        0xcb8000.0p-24,  0xeb050c.0p-26, -0x1bf2badc0df842.0p-85,
+        0xca0000.0p-24,  0xf29878.0p-26, -0x18efededd89fbe.0p-87,
+        0xc90000.0p-24,  0xf7ad6f.0p-26,  0x1373ff977baa69.0p-81,
+        0xc80000.0p-24,  0xfcc8e3.0p-26,  0x196766f2fb3283.0p-80,
+        0xc68000.0p-24,  0x823f30.0p-25,  0x19bd076f7c434e.0p-79,
+        0xc58000.0p-24,  0x84d52c.0p-25, -0x1a327257af0f46.0p-79,
+        0xc40000.0p-24,  0x88bc74.0p-25,  0x113f23def19c5a.0p-81,
+        0xc30000.0p-24,  0x8b5ae6.0p-25,  0x1759f6e6b37de9.0p-79,
+        0xc20000.0p-24,  0x8dfccb.0p-25,  0x1ad35ca6ed5148.0p-81,
+        0xc10000.0p-24,  0x90a22b.0p-25,  0x1a1d71a87deba4.0p-79,
+        0xbf8000.0p-24,  0x94a0d8.0p-25, -0x139e5210c2b731.0p-80,
+        0xbe8000.0p-24,  0x974f16.0p-25, -0x18f6ebcff3ed73.0p-81,
+        0xbd8000.0p-24,  0x9a00f1.0p-25, -0x1aa268be39aab7.0p-79,
+        0xbc8000.0p-24,  0x9cb672.0p-25, -0x14c8815839c566.0p-79,
+        0xbb0000.0p-24,  0xa0cda1.0p-25,  0x1eaf46390dbb24.0p-81,
+        0xba0000.0p-24,  0xa38c6e.0p-25,  0x138e20d831f698.0p-81,
+        0xb90000.0p-24,  0xa64f05.0p-25, -0x1e8d3c41123616.0p-82,
+        0xb80000.0p-24,  0xa91570.0p-25,  0x1ce28f5f3840b2.0p-80,
+        0xb70000.0p-24,  0xabdfbb.0p-25, -0x186e5c0a424234.0p-79,
+        0xb60000.0p-24,  0xaeadef.0p-25, -0x14d41a0b2a08a4.0p-83,
+        0xb50000.0p-24,  0xb18018.0p-25,  0x16755892770634.0p-79,
+        0xb40000.0p-24,  0xb45642.0p-25, -0x16395ebe59b152.0p-82,
+        0xb30000.0p-24,  0xb73077.0p-25,  0x1abc65c8595f09.0p-80,
+        0xb20000.0p-24,  0xba0ec4.0p-25, -0x1273089d3dad89.0p-79,
+        0xb10000.0p-24,  0xbcf133.0p-25,  0x10f9f67b1f4bbf.0p-79,
+        0xb00000.0p-24,  0xbfd7d2.0p-25, -0x109fab90486409.0p-80,
+        0xaf0000.0p-24,  0xc2c2ac.0p-25, -0x1124680aa43333.0p-79,
+        0xae8000.0p-24,  0xc439b3.0p-25, -0x1f360cc4710fc0.0p-80,
+        0xad8000.0p-24,  0xc72afd.0p-25, -0x132d91f21d89c9.0p-80,
+        0xac8000.0p-24,  0xca20a2.0p-25, -0x16bf9b4d1f8da8.0p-79,
+        0xab8000.0p-24,  0xcd1aae.0p-25,  0x19deb5ce6a6a87.0p-81,
+        0xaa8000.0p-24,  0xd0192f.0p-25,  0x1a29fb48f7d3cb.0p-79,
+        0xaa0000.0p-24,  0xd19a20.0p-25,  0x1127d3c6457f9d.0p-81,
+        0xa90000.0p-24,  0xd49f6a.0p-25, -0x1ba930e486a0ac.0p-81,
+        0xa80000.0p-24,  0xd7a94b.0p-25, -0x1b6e645f31549e.0p-79,
+        0xa70000.0p-24,  0xdab7d0.0p-25,  0x1118a425494b61.0p-80,
+        0xa68000.0p-24,  0xdc40d5.0p-25,  0x1966f24d29d3a3.0p-80,
+        0xa58000.0p-24,  0xdf566d.0p-25, -0x1d8e52eb2248f1.0p-82,
+        0xa48000.0p-24,  0xe270ce.0p-25, -0x1ee370f96e6b68.0p-80,
+        0xa40000.0p-24,  0xe3ffce.0p-25,  0x1d155324911f57.0p-80,
+        0xa30000.0p-24,  0xe72179.0p-25, -0x1fe6e2f2f867d9.0p-80,
+        0xa20000.0p-24,  0xea4812.0p-25,  0x1b7be9add7f4d4.0p-80,
+        0xa18000.0p-24,  0xebdd3d.0p-25,  0x1b3cfb3f7511dd.0p-79,
+        0xa08000.0p-24,  0xef0b5b.0p-25, -0x1220de1f730190.0p-79,
+        0xa00000.0p-24,  0xf0a451.0p-25, -0x176364c9ac81cd.0p-80,
+        0x9f0000.0p-24,  0xf3da16.0p-25,  0x1eed6b9aafac8d.0p-81,
+        0x9e8000.0p-24,  0xf576e9.0p-25,  0x1d593218675af2.0p-79,
+        0x9d8000.0p-24,  0xf8b47c.0p-25, -0x13e8eb7da053e0.0p-84,
+        0x9d0000.0p-24,  0xfa553f.0p-25,  0x1c063259bcade0.0p-79,
+        0x9c0000.0p-24,  0xfd9ac5.0p-25,  0x1ef491085fa3c1.0p-79,
+        0x9b8000.0p-24,  0xff3f8c.0p-25,  0x1d607a7c2b8c53.0p-79,
+        0x9a8000.0p-24,  0x814697.0p-24, -0x12ad3817004f3f.0p-78,
+        0x9a0000.0p-24,  0x821b06.0p-24, -0x189fc53117f9e5.0p-81,
+        0x990000.0p-24,  0x83c5f8.0p-24,  0x14cf15a048907b.0p-79,
+        0x988000.0p-24,  0x849c7d.0p-24,  0x1cbb1d35fb8287.0p-78,
+        0x978000.0p-24,  0x864ba6.0p-24,  0x1128639b814f9c.0p-78,
+        0x970000.0p-24,  0x87244c.0p-24,  0x184733853300f0.0p-79,
+        0x968000.0p-24,  0x87fdaa.0p-24,  0x109d23aef77dd6.0p-80,
+        0x958000.0p-24,  0x89b293.0p-24, -0x1a81ef367a59de.0p-78,
+        0x950000.0p-24,  0x8a8e20.0p-24, -0x121ad3dbb2f452.0p-78,
+        0x948000.0p-24,  0x8b6a6a.0p-24, -0x1cfb981628af72.0p-79,
+        0x938000.0p-24,  0x8d253a.0p-24, -0x1d21730ea76cfe.0p-79,
+        0x930000.0p-24,  0x8e03c2.0p-24,  0x135cc00e566f77.0p-78,
+        0x928000.0p-24,  0x8ee30d.0p-24, -0x10fcb5df257a26.0p-80,
+        0x918000.0p-24,  0x90a3ee.0p-24, -0x16e171b15433d7.0p-79,
+        0x910000.0p-24,  0x918587.0p-24, -0x1d050da07f3237.0p-79,
+        0x908000.0p-24,  0x9267e7.0p-24,  0x1be03669a5268d.0p-79,
+        0x8f8000.0p-24,  0x942f04.0p-24,  0x10b28e0e26c337.0p-79,
+        0x8f0000.0p-24,  0x9513c3.0p-24,  0x1a1d820da57cf3.0p-78,
+        0x8e8000.0p-24,  0x95f950.0p-24, -0x19ef8f13ae3cf1.0p-79,
+        0x8e0000.0p-24,  0x96dfab.0p-24, -0x109e417a6e507c.0p-78,
+        0x8d0000.0p-24,  0x98aed2.0p-24,  0x10d01a2c5b0e98.0p-79,
+        0x8c8000.0p-24,  0x9997a2.0p-24, -0x1d6a50d4b61ea7.0p-78,
+        0x8c0000.0p-24,  0x9a8145.0p-24,  0x1b3b190b83f952.0p-78,
+        0x8b8000.0p-24,  0x9b6bbf.0p-24,  0x13a69fad7e7abe.0p-78,
+        0x8b0000.0p-24,  0x9c5711.0p-24, -0x11cd12316f576b.0p-78,
+        0x8a8000.0p-24,  0x9d433b.0p-24,  0x1c95c444b807a2.0p-79,
+        0x898000.0p-24,  0x9f1e22.0p-24, -0x1b9c224ea698c3.0p-79,
+        0x890000.0p-24,  0xa00ce1.0p-24,  0x125ca93186cf0f.0p-81,
+        0x888000.0p-24,  0xa0fc80.0p-24, -0x1ee38a7bc228b3.0p-79,
+        0x880000.0p-24,  0xa1ed00.0p-24, -0x1a0db876613d20.0p-78,
+        0x878000.0p-24,  0xa2de62.0p-24,  0x193224e8516c01.0p-79,
+        0x870000.0p-24,  0xa3d0a9.0p-24,  0x1fa28b4d2541ad.0p-79,
+        0x868000.0p-24,  0xa4c3d6.0p-24,  0x1c1b5760fb4572.0p-78,
+        0x858000.0p-24,  0xa6acea.0p-24,  0x1fed5d0f65949c.0p-80,
+        0x850000.0p-24,  0xa7a2d4.0p-24,  0x1ad270c9d74936.0p-80,
+        0x848000.0p-24,  0xa899ab.0p-24,  0x199ff15ce53266.0p-79,
+        0x840000.0p-24,  0xa99171.0p-24,  0x1a19e15ccc45d2.0p-79,
+        0x838000.0p-24,  0xaa8a28.0p-24, -0x121a14ec532b36.0p-80,
+        0x830000.0p-24,  0xab83d1.0p-24,  0x1aee319980bff3.0p-79,
+        0x828000.0p-24,  0xac7e6f.0p-24, -0x18ffd9e3900346.0p-80,
+        0x820000.0p-24,  0xad7a03.0p-24, -0x1e4db102ce29f8.0p-80,
+        0x818000.0p-24,  0xae768f.0p-24,  0x17c35c55a04a83.0p-81,
+        0x810000.0p-24,  0xaf7415.0p-24,  0x1448324047019b.0p-78,
+        0x808000.0p-24,  0xb07298.0p-24, -0x1750ee3915a198.0p-78,
+        0x800000.0p-24,  0xb17218.0p-24, -0x105c610ca86c39.0p-81,
+};
+
+#ifdef USE_UTAB
+static const struct {
+       float   H;                      /* 1 + i/INTERVALS (exact) */
+       float   E;                      /* H(i) * G(i) - 1 (exact) */
+} U[TSIZE] = {
+        0x800000.0p-23,  0,
+        0x810000.0p-23, -0x800000.0p-37,
+        0x820000.0p-23, -0x800000.0p-35,
+        0x830000.0p-23, -0x900000.0p-34,
+        0x840000.0p-23, -0x800000.0p-33,
+        0x850000.0p-23, -0xc80000.0p-33,
+        0x860000.0p-23, -0xa00000.0p-36,
+        0x870000.0p-23,  0x940000.0p-33,
+        0x880000.0p-23,  0x800000.0p-35,
+        0x890000.0p-23, -0xc80000.0p-34,
+        0x8a0000.0p-23,  0xe00000.0p-36,
+        0x8b0000.0p-23,  0x900000.0p-33,
+        0x8c0000.0p-23, -0x800000.0p-35,
+        0x8d0000.0p-23, -0xe00000.0p-33,
+        0x8e0000.0p-23,  0x880000.0p-33,
+        0x8f0000.0p-23, -0xa80000.0p-34,
+        0x900000.0p-23, -0x800000.0p-35,
+        0x910000.0p-23,  0x800000.0p-37,
+        0x920000.0p-23,  0x900000.0p-35,
+        0x930000.0p-23,  0xd00000.0p-35,
+        0x940000.0p-23,  0xe00000.0p-35,
+        0x950000.0p-23,  0xc00000.0p-35,
+        0x960000.0p-23,  0xe00000.0p-36,
+        0x970000.0p-23, -0x800000.0p-38,
+        0x980000.0p-23, -0xc00000.0p-35,
+        0x990000.0p-23, -0xd00000.0p-34,
+        0x9a0000.0p-23,  0x880000.0p-33,
+        0x9b0000.0p-23,  0xe80000.0p-35,
+        0x9c0000.0p-23, -0x800000.0p-35,
+        0x9d0000.0p-23,  0xb40000.0p-33,
+        0x9e0000.0p-23,  0x880000.0p-34,
+        0x9f0000.0p-23, -0xe00000.0p-35,
+        0xa00000.0p-23,  0x800000.0p-33,
+        0xa10000.0p-23, -0x900000.0p-36,
+        0xa20000.0p-23, -0xb00000.0p-33,
+        0xa30000.0p-23, -0xa00000.0p-36,
+        0xa40000.0p-23,  0x800000.0p-33,
+        0xa50000.0p-23, -0xf80000.0p-35,
+        0xa60000.0p-23,  0x880000.0p-34,
+        0xa70000.0p-23, -0x900000.0p-33,
+        0xa80000.0p-23, -0x800000.0p-35,
+        0xa90000.0p-23,  0x900000.0p-34,
+        0xaa0000.0p-23,  0xa80000.0p-33,
+        0xab0000.0p-23, -0xac0000.0p-34,
+        0xac0000.0p-23, -0x800000.0p-37,
+        0xad0000.0p-23,  0xf80000.0p-35,
+        0xae0000.0p-23,  0xf80000.0p-34,
+        0xaf0000.0p-23, -0xac0000.0p-33,
+        0xb00000.0p-23, -0x800000.0p-33,
+        0xb10000.0p-23, -0xb80000.0p-34,
+        0xb20000.0p-23, -0x800000.0p-34,
+        0xb30000.0p-23, -0xb00000.0p-35,
+        0xb40000.0p-23, -0x800000.0p-35,
+        0xb50000.0p-23, -0xe00000.0p-36,
+        0xb60000.0p-23, -0x800000.0p-35,
+        0xb70000.0p-23, -0xb00000.0p-35,
+        0xb80000.0p-23, -0x800000.0p-34,
+        0xb90000.0p-23, -0xb80000.0p-34,
+        0xba0000.0p-23, -0x800000.0p-33,
+        0xbb0000.0p-23, -0xac0000.0p-33,
+        0xbc0000.0p-23,  0x980000.0p-33,
+        0xbd0000.0p-23,  0xbc0000.0p-34,
+        0xbe0000.0p-23,  0xe00000.0p-36,
+        0xbf0000.0p-23, -0xb80000.0p-35,
+        0xc00000.0p-23, -0x800000.0p-33,
+        0xc10000.0p-23,  0xa80000.0p-33,
+        0xc20000.0p-23,  0x900000.0p-34,
+        0xc30000.0p-23, -0x800000.0p-35,
+        0xc40000.0p-23, -0x900000.0p-33,
+        0xc50000.0p-23,  0x820000.0p-33,
+        0xc60000.0p-23,  0x800000.0p-38,
+        0xc70000.0p-23, -0x820000.0p-33,
+        0xc80000.0p-23,  0x800000.0p-33,
+        0xc90000.0p-23, -0xa00000.0p-36,
+        0xca0000.0p-23, -0xb00000.0p-33,
+        0xcb0000.0p-23,  0x840000.0p-34,
+        0xcc0000.0p-23, -0xd00000.0p-34,
+        0xcd0000.0p-23,  0x800000.0p-33,
+        0xce0000.0p-23, -0xe00000.0p-35,
+        0xcf0000.0p-23,  0xa60000.0p-33,
+        0xd00000.0p-23, -0x800000.0p-35,
+        0xd10000.0p-23,  0xb40000.0p-33,
+        0xd20000.0p-23, -0x800000.0p-35,
+        0xd30000.0p-23,  0xaa0000.0p-33,
+        0xd40000.0p-23, -0xe00000.0p-35,
+        0xd50000.0p-23,  0x880000.0p-33,
+        0xd60000.0p-23, -0xd00000.0p-34,
+        0xd70000.0p-23,  0x9c0000.0p-34,
+        0xd80000.0p-23, -0xb00000.0p-33,
+        0xd90000.0p-23, -0x800000.0p-38,
+        0xda0000.0p-23,  0xa40000.0p-33,
+        0xdb0000.0p-23, -0xdc0000.0p-34,
+        0xdc0000.0p-23,  0xc00000.0p-35,
+        0xdd0000.0p-23,  0xca0000.0p-33,
+        0xde0000.0p-23, -0xb80000.0p-34,
+        0xdf0000.0p-23,  0xd00000.0p-35,
+        0xe00000.0p-23,  0xc00000.0p-33,
+        0xe10000.0p-23, -0xf40000.0p-34,
+        0xe20000.0p-23,  0x800000.0p-37,
+        0xe30000.0p-23,  0x860000.0p-33,
+        0xe40000.0p-23, -0xc80000.0p-33,
+        0xe50000.0p-23, -0xa80000.0p-34,
+        0xe60000.0p-23,  0xe00000.0p-36,
+        0xe70000.0p-23,  0x880000.0p-33,
+        0xe80000.0p-23, -0xe00000.0p-33,
+        0xe90000.0p-23, -0xfc0000.0p-34,
+        0xea0000.0p-23, -0x800000.0p-35,
+        0xeb0000.0p-23,  0xe80000.0p-35,
+        0xec0000.0p-23,  0x900000.0p-33,
+        0xed0000.0p-23,  0xe20000.0p-33,
+        0xee0000.0p-23, -0xac0000.0p-33,
+        0xef0000.0p-23, -0xc80000.0p-34,
+        0xf00000.0p-23, -0x800000.0p-35,
+        0xf10000.0p-23,  0x800000.0p-35,
+        0xf20000.0p-23,  0xb80000.0p-34,
+        0xf30000.0p-23,  0x940000.0p-33,
+        0xf40000.0p-23,  0xc80000.0p-33,
+        0xf50000.0p-23, -0xf20000.0p-33,
+        0xf60000.0p-23, -0xc80000.0p-33,
+        0xf70000.0p-23, -0xa20000.0p-33,
+        0xf80000.0p-23, -0x800000.0p-33,
+        0xf90000.0p-23, -0xc40000.0p-34,
+        0xfa0000.0p-23, -0x900000.0p-34,
+        0xfb0000.0p-23, -0xc80000.0p-35,
+        0xfc0000.0p-23, -0x800000.0p-35,
+        0xfd0000.0p-23, -0x900000.0p-36,
+        0xfe0000.0p-23, -0x800000.0p-37,
+        0xff0000.0p-23, -0x800000.0p-39,
+        0x800000.0p-22,  0,
+};
+#endif /* USE_UTAB */
+
+#ifdef STRUCT_RETURN
+#define        RETURN1(rp, v) do {     \
+       (rp)->hi = (v);         \
+       (rp)->lo_set = 0;       \
+       return;                 \
+} while (0)
+
+#define        RETURN2(rp, h, l) do {  \
+       (rp)->hi = (h);         \
+       (rp)->lo = (l);         \
+       (rp)->lo_set = 1;       \
+       return;                 \
+} while (0)
+
+struct ld {
+       long double hi;
+       long double lo;
+       int     lo_set;
+};
+#else
+#define        RETURN1(rp, v)  RETURNF(v)
+#define        RETURN2(rp, h, l)       RETURNI((h) + (l))
+#endif
+
+#ifdef STRUCT_RETURN
+static inline __always_inline void
+k_logl(long double x, struct ld *rp)
+#else
+long double
+logl(long double x)
+#endif
+{
+       long double d, dk, val_hi, val_lo, z;
+       uint64_t ix, lx;
+       int i, k;
+       uint16_t hx;
+
+       EXTRACT_LDBL80_WORDS(hx, lx, x);
+       k = -16383;
+#if 0 /* Hard to do efficiently.  Don't do it until we support all modes. */
+       if (x == 1)
+               RETURN1(rp, 0);         /* log(1) = +0 in all rounding modes */
+#endif
+       if (hx == 0 || hx >= 0x8000) {  /* zero, negative or subnormal? */
+               if (((hx & 0x7fff) | lx) == 0)
+                       RETURN1(rp, -1 / zero); /* log(+-0) = -Inf */
+               if (hx != 0)
+                       /* log(neg or [pseudo-]NaN) = qNaN: */
+                       RETURN1(rp, (x - x) / zero);
+               x *= 0x1.0p65;          /* subnormal; scale up x */
+                                       /* including pseudo-subnormals */
+               EXTRACT_LDBL80_WORDS(hx, lx, x);
+               k = -16383 - 65;
+       } else if (hx >= 0x7fff || (lx & 0x8000000000000000ULL) == 0)
+               RETURN1(rp, x + x);     /* log(Inf or NaN) = Inf or qNaN */
+                                       /* log(pseudo-Inf) = qNaN */
+                                       /* log(pseudo-NaN) = qNaN */
+                                       /* log(unnormal) = qNaN */
+#ifndef STRUCT_RETURN
+       ENTERI();
+#endif
+       k += hx;
+       ix = lx & 0x7fffffffffffffffULL;
+       dk = k;
+
+       /* Scale x to be in [1, 2). */
+       SET_LDBL_EXPSIGN(x, 0x3fff);
+
+       /* 0 <= i <= INTERVALS: */
+#define        L2I     (64 - LOG2_INTERVALS)
+       i = (ix + (1LL << (L2I - 2))) >> (L2I - 1);
+
+       /*
+        * -0.005280 < d < 0.004838.  In particular, the infinite-
+        * precision |d| is <= 2**-7.  Rounding of G(i) to 8 bits
+        * ensures that d is representable without extra precision for
+        * this bound on |d| (since when this calculation is expressed
+        * as x*G(i)-1, the multiplication needs as many extra bits as
+        * G(i) has and the subtraction cancels 8 bits).  But for
+        * most i (107 cases out of 129), the infinite-precision |d|
+        * is <= 2**-8.  G(i) is rounded to 9 bits for such i to give
+        * better accuracy (this works by improving the bound on |d|,
+        * which in turn allows rounding to 9 bits in more cases).
+        * This is only important when the original x is near 1 -- it
+        * lets us avoid using a special method to give the desired
+        * accuracy for such x.
+        */
+       if (0)
+               d = x * G(i) - 1;
+       else {
+#ifdef USE_UTAB
+               d = (x - H(i)) * G(i) + E(i);
+#else
+               long double x_hi, x_lo;
+               float fx_hi;
+
+               /*
+                * Split x into x_hi + x_lo to calculate x*G(i)-1 exactly.
+                * G(i) has at most 9 bits, so the splitting point is not
+                * critical.
+                */
+               SET_FLOAT_WORD(fx_hi, (lx >> 40) | 0x3f800000);
+               x_hi = fx_hi;
+               x_lo = x - x_hi;
+               d = x_hi * G(i) - 1 + x_lo * G(i);
+#endif
+       }
+
+       /*
+        * Our algorithm depends on exact cancellation of F_lo(i) and
+        * F_hi(i) with dk*ln_2_lo and dk*ln2_hi when k is -1 and i is
+        * at the end of the table.  This and other technical complications
+        * make it difficult to avoid the double scaling in (dk*ln2) *
+        * log(base) for base != e without losing more accuracy and/or
+        * efficiency than is gained.
+        */
+       z = d * d;
+       val_lo = z * d * z * (z * (d * P8 + P7) + (d * P6 + P5)) +
+           (F_lo(i) + dk * ln2_lo + z * d * (d * P4 + P3)) + z * P2;
+       val_hi = d;
+#ifdef DEBUG
+       if (fetestexcept(FE_UNDERFLOW))
+               breakpoint();
+#endif
+
+       _3sumF(val_hi, val_lo, F_hi(i) + dk * ln2_hi);
+       RETURN2(rp, val_hi, val_lo);
+}
+
+long double
+log1pl(long double x)
+{
+       long double d, d_hi, d_lo, dk, f_lo, val_hi, val_lo, z;
+       long double f_hi, twopminusk;
+       uint64_t ix, lx;
+       int i, k;
+       int16_t ax, hx;
+
+       DOPRINT_START(&x);
+       EXTRACT_LDBL80_WORDS(hx, lx, x);
+       if (hx < 0x3fff) {              /* x < 1, or x neg NaN */
+               ax = hx & 0x7fff;
+               if (ax >= 0x3fff) {     /* x <= -1, or x neg NaN */
+                       if (ax == 0x3fff && lx == 0x8000000000000000ULL)
+                               RETURNP(-1 / zero);     /* log1p(-1) = -Inf */
+                       /* log1p(x < 1, or x [pseudo-]NaN) = qNaN: */
+                       RETURNP((x - x) / (x - x));
+               }
+               if (ax <= 0x3fbe) {     /* |x| < 2**-64 */
+                       if ((int)x == 0)
+                               RETURNP(x);     /* x with inexact if x != 0 */
+               }
+               f_hi = 1;
+               f_lo = x;
+       } else if (hx >= 0x7fff) {      /* x +Inf or non-neg NaN */
+               RETURNP(x + x);         /* log1p(Inf or NaN) = Inf or qNaN */
+                                       /* log1p(pseudo-Inf) = qNaN */
+                                       /* log1p(pseudo-NaN) = qNaN */
+                                       /* log1p(unnormal) = qNaN */
+       } else if (hx < 0x407f) {       /* 1 <= x < 2**128 */
+               f_hi = x;
+               f_lo = 1;
+       } else {                        /* 2**128 <= x < +Inf */
+               f_hi = x;
+               f_lo = 0;               /* avoid underflow of the P5 term */
+       }
+       ENTERI();
+       x = f_hi + f_lo;
+       f_lo = (f_hi - x) + f_lo;
+
+       EXTRACT_LDBL80_WORDS(hx, lx, x);
+       k = -16383;
+
+       k += hx;
+       ix = lx & 0x7fffffffffffffffULL;
+       dk = k;
+
+       SET_LDBL_EXPSIGN(x, 0x3fff);
+       twopminusk = 1;
+       SET_LDBL_EXPSIGN(twopminusk, 0x7ffe - (hx & 0x7fff));
+       f_lo *= twopminusk;
+
+       i = (ix + (1LL << (L2I - 2))) >> (L2I - 1);
+
+       /*
+        * x*G(i)-1 (with a reduced x) can be represented exactly, as
+        * above, but now we need to evaluate the polynomial on d =
+        * (x+f_lo)*G(i)-1 and extra precision is needed for that.
+        * Since x+x_lo is a hi+lo decomposition and subtracting 1
+        * doesn't lose too many bits, an inexact calculation for
+        * f_lo*G(i) is good enough.
+        */
+       if (0)
+               d_hi = x * G(i) - 1;
+       else {
+#ifdef USE_UTAB
+               d_hi = (x - H(i)) * G(i) + E(i);
+#else
+               long double x_hi, x_lo;
+               float fx_hi;
+
+               SET_FLOAT_WORD(fx_hi, (lx >> 40) | 0x3f800000);
+               x_hi = fx_hi;
+               x_lo = x - x_hi;
+               d_hi = x_hi * G(i) - 1 + x_lo * G(i);
+#endif
+       }
+       d_lo = f_lo * G(i);
+
+       /*
+        * This is _2sumF(d_hi, d_lo) inlined.  The condition
+        * (d_hi == 0 || |d_hi| >= |d_lo|) for using _2sumF() is not
+        * always satisifed, so it is not clear that this works, but
+        * it works in practice.  It works even if it gives a wrong
+        * normalized d_lo, since |d_lo| > |d_hi| implies that i is
+        * nonzero and d is tiny, so the F(i) term dominates d_lo.
+        * In float precision:
+        * (By exhaustive testing, the worst case is d_hi = 0x1.bp-25.
+        * And if d is only a little tinier than that, we would have
+        * another underflow problem for the P3 term; this is also ruled
+        * out by exhaustive testing.)
+        */
+       d = d_hi + d_lo;
+       d_lo = d_hi - d + d_lo;
+       d_hi = d;
+
+       z = d * d;
+       val_lo = z * d * z * (z * (d * P8 + P7) + (d * P6 + P5)) +
+           (F_lo(i) + dk * ln2_lo + d_lo + z * d * (d * P4 + P3)) + z * P2;
+       val_hi = d_hi;
+#ifdef DEBUG
+       if (fetestexcept(FE_UNDERFLOW))
+               breakpoint();
+#endif
+
+       _3sumF(val_hi, val_lo, F_hi(i) + dk * ln2_hi);
+       RETURN2PI(val_hi, val_lo);
+}
+
+#ifdef STRUCT_RETURN
+
+long double
+logl(long double x)
+{
+       struct ld r;
+
+       ENTERI();
+       DOPRINT_START(&x);
+       k_logl(x, &r);
+       RETURNSPI(&r);
+}
+
+static const double
+invln10_hi =  4.3429448190317999e-1,           /*  0x1bcb7b1526e000.0p-54 */
+invln10_lo =  7.1842412889749798e-14,          /*  0x1438ca9aadd558.0p-96 */
+invln2_hi =  1.4426950408887933e0,             /*  0x171547652b8000.0p-52 */
+invln2_lo =  1.7010652264631490e-13;           /*  0x17f0bbbe87fed0.0p-95 */
+
+long double
+log10l(long double x)
+{
+       struct ld r;
+       long double hi, lo;
+
+       ENTERI();
+       DOPRINT_START(&x);
+       k_logl(x, &r);
+       if (!r.lo_set)
+               RETURNPI(r.hi);
+       _2sumF(r.hi, r.lo);
+       hi = (float)r.hi;
+       lo = r.lo + (r.hi - hi);
+       RETURN2PI(invln10_hi * hi,
+           (invln10_lo + invln10_hi) * lo + invln10_lo * hi);
+}
+
+long double
+log2l(long double x)
+{
+       struct ld r;
+       long double hi, lo;
+
+       ENTERI();
+       DOPRINT_START(&x);
+       k_logl(x, &r);
+       if (!r.lo_set)
+               RETURNPI(r.hi);
+       _2sumF(r.hi, r.lo);
+       hi = (float)r.hi;
+       lo = r.lo + (r.hi - hi);
+       RETURN2PI(invln2_hi * hi,
+           (invln2_lo + invln2_hi) * lo + invln2_lo * hi);
+}
+
+#endif /* STRUCT_RETURN */
index 8db27ea..2bd787d 100644 (file)
 .\" SUCH DAMAGE.
 .\"
 .\"     from: @(#)acosh.3      5.2 (Berkeley) 5/6/91
-.\" $FreeBSD: head/lib/msun/man/acosh.3 165906 2007-01-09 01:02:06Z imp $
+.\" $FreeBSD: head/lib/msun/man/acosh.3 251599 2013-06-10 06:04:58Z das $
 .\"
-.Dd January 14, 2005
+.Dd June 11, 2013
 .Dt ACOSH 3
 .Os
 .Sh NAME
 .Nm acosh ,
-.Nm acoshf
+.Nm acoshf ,
+.Nm acoshl
 .Nd inverse hyperbolic cosine functions
 .Sh LIBRARY
 .Lb libm
 .Fn acosh "double x"
 .Ft float
 .Fn acoshf "float x"
+.Ft long double
+.Fn acoshl "long double x"
 .Sh DESCRIPTION
 The
-.Fn acosh
-and the
-.Fn acoshf
+.Fn acosh ,
+.Fn acoshf ,
+and
+.Fn acoshl
 functions compute the inverse hyperbolic cosine
 of the real
 argument
@@ -55,11 +59,7 @@ argument
 For a discussion of error due to roundoff, see
 .Xr math 3 .
 .Sh RETURN VALUES
-The
-.Fn acosh
-and the
-.Fn acoshf
-functions
+These functions
 return the inverse hyperbolic cosine of
 .Ar x .
 If the argument is less than 1,
@@ -73,6 +73,13 @@ raises an invalid exception and returns an \*(Na.
 .Xr math 3
 .Sh HISTORY
 The
-.Fn acosh
-function appeared in
-.Bx 4.3 .
+.Fn acosh ,
+.Fn acoshf ,
+and
+.Fn acoshl
+functions appeared in
+.Bx 4.3 ,
+.Fx 2.0 ,
+and
+.Dx 3.5 ,
+respectively.
index 77ac19e..c08ba30 100644 (file)
 .\" SUCH DAMAGE.
 .\"
 .\"     from: @(#)asinh.3      6.4 (Berkeley) 5/6/91
-.\" $FreeBSD: head/lib/msun/man/asinh.3 165906 2007-01-09 01:02:06Z imp $
+.\" $FreeBSD: head/lib/msun/man/asinh.3 251599 2013-06-10 06:04:58Z das $
 .\"
-.Dd May 6, 1991
+.Dd June 11, 2013
 .Dt ASINH 3
 .Os
 .Sh NAME
 .Nm asinh ,
-.Nm asinhf
+.Nm asinhf ,
+.Nm asinhl
 .Nd inverse hyperbolic sine functions
 .Sh LIBRARY
 .Lb libm
 .Fn asinh "double x"
 .Ft float
 .Fn asinhf "float x"
+.Ft long double
+.Fn asinhl "long double x"
 .Sh DESCRIPTION
 The
-.Fn asinh
-and the
-.Fn asinhf
+.Fn asinh ,
+.Fn asinhf ,
+and
+.Fn asinhl
 functions compute the inverse hyperbolic sine
 of the real
 argument
@@ -55,11 +59,7 @@ argument
 For a discussion of error due to roundoff, see
 .Xr math 3 .
 .Sh RETURN VALUES
-The
-.Fn asinh
-and the
-.Fn asinhf
-functions
+These functions
 return the inverse hyperbolic sine of
 .Ar x .
 .Sh SEE ALSO
@@ -69,6 +69,13 @@ return the inverse hyperbolic sine of
 .Xr math 3
 .Sh HISTORY
 The
-.Fn asinh
-function appeared in
-.Bx 4.3 .
+.Fn asinh ,
+.Fn asinhf ,
+and
+.Fn asinhl
+functions appeared in
+.Bx 4.3 ,
+.Fx 2.0 ,
+and
+.Dx 3.5 ,
+respectively.
index 9c1480a..bb6d967 100644 (file)
 .\" SUCH DAMAGE.
 .\"
 .\"     from: @(#)atanh.3      5.2 (Berkeley) 5/6/91
-.\" $FreeBSD: head/lib/msun/man/atanh.3 165906 2007-01-09 01:02:06Z imp $
+.\" $FreeBSD: head/lib/msun/man/atanh.3 251599 2013-06-10 06:04:58Z das $
 .\"
-.Dd January 14, 2005
+.Dd June 11, 2013
 .Dt ATANH 3
 .Os
 .Sh NAME
 .Nm atanh ,
-.Nm atanhf
+.Nm atanhf ,
+.Nm atanhl
 .Nd inverse hyperbolic tangent functions
 .Sh LIBRARY
 .Lb libm
 .Fn atanh "double x"
 .Ft float
 .Fn atanhf "float x"
+.Ft long double
+.Fn atanhl "long double x"
 .Sh DESCRIPTION
 The
-.Fn atanh
-and the
-.Fn atanhf
+.Fn atanh ,
+.Fn atanhf ,
+and
+.Fn atanhl
 functions compute the inverse hyperbolic tangent
 of the real
 argument
@@ -55,11 +59,7 @@ argument
 For a discussion of error due to roundoff, see
 .Xr math 3 .
 .Sh RETURN VALUES
-The
-.Fn atanh
-and the
-.Fn atanhf
-functions
+These functions
 return the inverse hyperbolic tangent of
 .Ar x
 if successful.
@@ -76,6 +76,13 @@ If
 .Xr math 3
 .Sh HISTORY
 The
-.Fn atanh
-function appeared in
-.Bx 4.3 .
+.Fn atanh ,
+.Fn atanhf ,
+and
+.Fn atanhl
+functions appeared in
+.Bx 4.3 ,
+.Fx 2.0 ,
+and
+.Dx 3.5 ,
+respectively.
index 143aae2..29af562 100644 (file)
@@ -26,9 +26,9 @@
 .\" SUCH DAMAGE.
 .\"
 .\"     from: @(#)exp.3        6.12 (Berkeley) 7/31/91
-.\" $FreeBSD: head/lib/msun/man/exp.3 238722 2012-07-23 19:13:55Z kargl $
+.\" $FreeBSD: head/lib/msun/man/exp.3 251343 2013-06-03 19:51:32Z kargl $
 .\"
-.Dd July 10, 2012
+.Dd June 11, 2013
 .Dt EXP 3
 .Os
 .Sh NAME
@@ -41,6 +41,7 @@
 .Nm exp2l ,
 .Nm expm1 ,
 .Nm expm1f ,
+.Nm expm1l ,
 .Nm pow ,
 .Nm powf
 .Nd exponential and power functions
@@ -64,6 +65,8 @@
 .Fn expm1 "double x"
 .Ft float
 .Fn expm1f "float x"
+.Ft long double
+.Fn expm1l "long double x"
 .Ft double
 .Fn pow "double x" "double y"
 .Ft float
@@ -88,9 +91,10 @@ functions compute the base 2 exponential of the given argument
 .Fa x .
 .Pp
 The
-.Fn expm1
+.Fn expm1 ,
+.Fn expm1f ,
 and the
-.Fn expm1f
+.Fn expm1l
 functions compute the value exp(x)\-1 accurately even for tiny argument
 .Fa x .
 .Pp
index c068b0a..a18e732 100644 (file)
@@ -22,9 +22,9 @@
 .\" OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
 .\" SUCH DAMAGE.
 .\"
-.\" $FreeBSD: head/lib/msun/man/log.3 216211 2010-12-05 22:11:22Z das $
+.\" $FreeBSD: head/lib/msun/man/log.3 251292 2013-06-03 09:14:31Z das $
 .\"
-.Dd December 21, 2011
+.Dd June 11, 2013
 .Dt LOG 3
 .Os
 .Sh NAME
 .Nm logl ,
 .Nm log10 ,
 .Nm log10f ,
+.Nm log10l ,
 .Nm log2 ,
 .Nm log2f ,
+.Nm log2l ,
 .Nm log1p ,
-.Nm log1pf
+.Nm log1pf ,
+.Nm log1pl
 .Nd logarithm functions
 .Sh LIBRARY
 .Lb libm
 .Fn log "double x"
 .Ft float
 .Fn logf "float x"
+.Ft long double
+.Fn logl "long double x"
 .Ft double
 .Fn log10 "double x"
 .Ft float
 .Fn log10f "float x"
+.Ft long double
+.Fn log10l "long double x"
 .Ft double
 .Fn log2 "double x"
 .Ft float
 .Fn log2f "float x"
+.Ft long double
+.Fn log2l "long double x"
 .Ft double
 .Fn log1p "double x"
 .Ft float
 .Fn log1pf "float x"
+.Ft long double
+.Fn log1pl "long double x"
 .Sh DESCRIPTION
 The
-.Fn log
+.Fn log ,
+.Fn logf ,
 and
-.Fn logf
+.Fn logl
 functions compute the natural logarithm of
 .Fa x .
 .Pp
 The
-.Fn log10
+.Fn log10 ,
+.Fn log10f ,
 and
-.Fn log10f
+.Fn log10l
 functions compute the logarithm base 10 of
 .Fa x ,
 while
-.Fn log2
+.Fn log2 ,
+.Fn log2f ,
 and
-.Fn log2f
+.Fn log2l
 compute the logarithm base 2 of
 .Fa x .
 .Pp
 The
-.Fn log1p
+.Fn log1p ,
+.Fn log1pf ,
 and
-.Fn log1pf
+.Fn log1pl
 functions compute the natural logarithm of
 .No "1 + x" .
 Computing the natural logarithm as
@@ -107,12 +122,16 @@ results in an invalid exception and a return value of \*(Na.
 The
 .Fn log ,
 .Fn logf ,
+.Fn logl ,
 .Fn log10 ,
 .Fn log10f ,
+.Fn log10l ,
 .Fn log2 ,
 .Fn log2f ,
+.Fn log2l ,
 .Fn log1p ,
+.Fn log1pf ,
 and
-.Fn log1pf
+.Fn log1pl
 functions conform to
 .St -isoC-99 .
index d5497ac..049cca5 100644 (file)
@@ -70,9 +70,9 @@ functions return the sine value.
 .Xr asin 3 ,
 .Xr atan 3 ,
 .Xr atan2 3 ,
+.Xr csin 3 ,
 .Xr cos 3 ,
 .Xr cosh 3 ,
-.Xr csin 3 ,
 .Xr math 3 ,
 .Xr sinh 3 ,
 .Xr tan 3 ,
index fd1312d..f407f10 100644 (file)
@@ -23,7 +23,7 @@
  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  * SUCH DAMAGE.
  *
- * $FreeBSD: head/lib/msun/src/catrig.c 251121 2013-05-30 04:49:26Z das $
+ * $FreeBSD: head/lib/msun/src/catrig.c 251404 2013-06-05 05:33:01Z das $
  */
 
 #include <complex.h>
@@ -173,7 +173,7 @@ do_hard_work(double x, double y, double *rx, int *B_is_usable, double *B,
                 * Am1 = fp + fm, where fp = f(x, 1+y), and fm = f(x, 1-y).
                 * rx = log1p(Am1 + sqrt(Am1*(A+1)))
                 */
-               if (y == 1 && x < DBL_EPSILON*DBL_EPSILON / 128) {
+               if (y == 1 && x < DBL_EPSILON * DBL_EPSILON / 128) {
                        /*
                         * fp is of order x^2, and fm = x/2.
                         * A = 1 (inexactly).
@@ -192,7 +192,7 @@ do_hard_work(double x, double y, double *rx, int *B_is_usable, double *B,
                         * A = 1 (inexactly).
                         */
                        *rx = x / sqrt((1 - y) * (1 + y));
-               } else /* if (y > 1) */ {
+               } else {                /* if (y > 1) */
                        /*
                         * A-1 = y-1 (inexactly).
                         */
@@ -252,7 +252,7 @@ do_hard_work(double x, double y, double *rx, int *B_is_usable, double *B,
                        *sqrt_A2my2 = x * (4 / DBL_EPSILON / DBL_EPSILON) * y /
                                sqrt((y + 1) * (y - 1));
                        *new_y = y * (4 / DBL_EPSILON / DBL_EPSILON);
-               } else /* if (y < 1) */ {
+               } else {                /* if (y < 1) */
                        /*
                         * fm = 1-y >= DBL_EPSILON, fp is of order x^2, and
                         * A = 1 (inexactly).
@@ -297,7 +297,6 @@ casinh(double complex z)
                 * C99 leaves it optional whether to raise invalid if one of
                 * the arguments is not NaN, so we opt not to raise it.
                 */
-               /* Bruce Evans tells me this is the way to do this: */
                return (cpack(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
        }
 
@@ -336,6 +335,7 @@ double complex
 casin(double complex z)
 {
        double complex w = casinh(cpack(cimag(z), creal(z)));
+
        return (cpack(cimag(w), creal(w)));
 }
 
@@ -401,17 +401,17 @@ cacos(double complex z)
        /* All remaining cases are inexact. */
        raise_inexact();
 
-       if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON/4)
+       if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4)
                return (cpack(pio2_hi - (x - pio2_lo), -y));
 
        do_hard_work(ay, ax, &ry, &B_is_usable, &B, &sqrt_A2mx2, &new_x);
        if (B_is_usable) {
-               if (sx==0)
+               if (sx == 0)
                        rx = acos(B);
                else
                        rx = acos(-B);
        } else {
-               if (sx==0)
+               if (sx == 0)
                        rx = atan2(sqrt_A2mx2, new_x);
                else
                        rx = atan2(sqrt_A2mx2, -new_x);
@@ -486,10 +486,6 @@ clog_for_large_values(double complex z)
        return (cpack(log(ax * ax + ay * ay) / 2, atan2(y, x)));
 }
 
-/*
- *=============================================================================
- */
-
 /*
  *                             =================
  *                             | catanh, catan |
@@ -510,6 +506,7 @@ sum_squares(double x, double y)
        /* Avoid underflow when y is small. */
        if (y < SQRT_MIN)
                return (x * x);
+
        return (x * x + y * y);
 }
 
@@ -588,10 +585,9 @@ catanh(double complex z)
                if (isinf(x))
                        return (cpack(copysign(0, x), y + y));
                /* catanh(NaN + I*+-Inf) = sign(NaN)0 + I*+-PI/2 */
-               if (isinf(y)) {
+               if (isinf(y))
                        return (cpack(copysign(0, x),
                                      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
@@ -600,10 +596,9 @@ catanh(double complex z)
                return (cpack(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
        }
 
-       if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
+       if (ax > RECIP_EPSILON || ay > RECIP_EPSILON)
                return (cpack(real_part_reciprocal(x, y),
                              copysign(pio2_hi + pio2_lo, y)));
-       }
 
        if (ax < SQRT_3_EPSILON / 2 && ay < SQRT_3_EPSILON / 2) {
                /*
@@ -616,7 +611,7 @@ catanh(double complex z)
        }
 
        if (ax == 1 && ay < DBL_EPSILON)
-               rx = (log(ay) - m_ln2) / -2;
+               rx = (m_ln2 - log(ay)) / 2;
        else
                rx = log1p(4 * ax / sum_squares(ax - 1, ay)) / 4;
 
@@ -638,5 +633,6 @@ double complex
 catan(double complex z)
 {
        double complex w = catanh(cpack(cimag(z), creal(z)));
+
        return (cpack(cimag(w), creal(w)));
 }
index f44a319..5c787c9 100644 (file)
@@ -23,7 +23,7 @@
  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  * SUCH DAMAGE.
  *
- * $FreeBSD: head/lib/msun/src/catrigf.c 251121 2013-05-30 04:49:26Z das $
+ * $FreeBSD: head/lib/msun/src/catrigf.c 251404 2013-06-05 05:33:01Z das $
  */
 
 /*
  * Mathematical Software, Volume 23 Issue 3, 1997, Pages 299-335,
  * http://dl.acm.org/citation.cfm?id=275324.
  *
- * The code for catrig.c contains complete comments.
+ * See catrig.c for complete comments.
+ *
+ * XXX comments were removed automatically, and even short ones on the right
+ * of statements were removed (all of them), contrary to normal style.  Only
+ * a few comments on the right of declarations remain.
  */
 
 #include <complex.h>
@@ -100,7 +104,7 @@ do_hard_work(float x, float y, float *rx, int *B_is_usable, float *B,
                        Am1 = f(x, 1 + y, R) + f(x, 1 - y, S);
                        *rx = log1pf(Am1 + sqrtf(Am1 * (A + 1)));
                } else if (y < 1) {
-                       *rx = x / sqrtf((1 - y)*(1 + y));
+                       *rx = x / sqrtf((1 - y) * (1 + y));
                } else {
                        *rx = log1pf((y - 1) + sqrtf((y - 1) * (y + 1)));
                }
@@ -188,6 +192,7 @@ float complex
 casinf(float complex z)
 {
        float complex w = casinhf(cpackf(cimagf(z), crealf(z)));
+
        return (cpackf(cimagf(w), crealf(w)));
 }
 
@@ -211,7 +216,8 @@ cacosf(float complex z)
                        return (cpackf(y + y, -INFINITY));
                if (isinf(y))
                        return (cpackf(x + x, -y));
-               if (x == 0) return (cpackf(pio2_hi + pio2_lo, y + y));
+               if (x == 0)
+                       return (cpackf(pio2_hi + pio2_lo, y + y));
                return (cpackf(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
        }
 
@@ -234,17 +240,17 @@ cacosf(float complex z)
 
        do_hard_work(ay, ax, &ry, &B_is_usable, &B, &sqrt_A2mx2, &new_x);
        if (B_is_usable) {
-               if (sx==0)
+               if (sx == 0)
                        rx = acosf(B);
                else
                        rx = acosf(-B);
        } else {
-               if (sx==0)
+               if (sx == 0)
                        rx = atan2f(sqrt_A2mx2, new_x);
                else
                        rx = atan2f(sqrt_A2mx2, -new_x);
        }
-       if (sy==0)
+       if (sy == 0)
                ry = -ry;
        return (cpackf(rx, ry));
 }
@@ -283,10 +289,9 @@ clog_for_large_values(float complex z)
                ay = t;
        }
 
-       if (ax > FLT_MAX / 2) {
+       if (ax > FLT_MAX / 2)
                return (cpackf(logf(hypotf(x / m_e, y / m_e)) + 1,
                               atan2f(y, x)));
-       }
 
        if (ax > QUARTER_SQRT_MAX || ay < SQRT_MIN)
                return (cpackf(logf(hypotf(x, y)), atan2f(y, x)));
@@ -299,8 +304,9 @@ sum_squares(float x, float y)
 {
 
        if (y < SQRT_MIN)
-               return (x*x);
-       return (x*x + y*y);
+               return (x * x);
+
+       return (x * x + y * y);
 }
 
 static inline float
@@ -317,9 +323,9 @@ real_part_reciprocal(float x, float y)
 #define        BIAS    (FLT_MAX_EXP - 1)
 #define        CUTOFF  (FLT_MANT_DIG / 2 + 1)
        if (ix - iy >= CUTOFF << 23 || isinf(x))
-               return (1/x);
+               return (1 / x);
        if (iy - ix >= CUTOFF << 23)
-               return (x/y/y);
+               return (x / y / y);
        if (ix <= (BIAS + FLT_MAX_EXP / 2 - CUTOFF) << 23)
                return (x / (x * x + y * y));
        SET_FLOAT_WORD(scale, 0x7f800000 - ix);
@@ -346,18 +352,16 @@ catanhf(float complex z)
 
        if (isnan(x) || isnan(y)) {
                if (isinf(x))
-                       return (cpackf(copysignf(0, x), y+y));
-               if (isinf(y)) {
+                       return (cpackf(copysignf(0, x), y + y));
+               if (isinf(y))
                        return (cpackf(copysignf(0, x),
                                       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) {
+       if (ax > RECIP_EPSILON || ay > RECIP_EPSILON)
                return (cpackf(real_part_reciprocal(x, y),
                               copysignf(pio2_hi + pio2_lo, y)));
-       }
 
        if (ax < SQRT_3_EPSILON / 2 && ay < SQRT_3_EPSILON / 2) {
                raise_inexact();
@@ -365,7 +369,7 @@ catanhf(float complex z)
        }
 
        if (ax == 1 && ay < FLT_EPSILON)
-               rx = (logf(ay) - m_ln2) / -2;
+               rx = (m_ln2 - logf(ay)) / 2;
        else
                rx = log1pf(4 * ax / sum_squares(ax - 1, ay)) / 4;
 
@@ -383,5 +387,6 @@ float complex
 catanf(float complex z)
 {
        float complex w = catanhf(cpackf(cimagf(z), crealf(z)));
+
        return (cpackf(cimagf(w), crealf(w)));
 }
index a14e09c..8de9689 100644 (file)
@@ -1,6 +1,6 @@
 
 /* @(#)e_acosh.c 1.3 95/01/18 */
-/* $FreeBSD: head/lib/msun/src/e_acosh.c 176451 2008-02-22 02:30:36Z das $ */
+/* $FreeBSD: head/lib/msun/src/e_acosh.c 251599 2013-06-10 06:04:58Z das $
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
@@ -10,6 +10,7 @@
  * software is freely granted, provided that this notice 
  * is preserved.
  * ====================================================
+ *
  */
 
 /* __ieee754_acosh(x)
@@ -26,6 +27,8 @@
  *     acosh(NaN) is NaN without signal.
  */
 
+#include <float.h>
+
 #include "math.h"
 #include "math_private.h"
 
@@ -57,3 +60,7 @@ __ieee754_acosh(double x)
            return log1p(t+sqrt(2.0*t+t*t));
        }
 }
+
+#if LDBL_MANT_DIG == 53
+__weak_reference(acosh, acoshl);
+#endif
diff --git a/lib/libm/src/e_acoshl.c b/lib/libm/src/e_acoshl.c
new file mode 100644 (file)
index 0000000..e25cd9c
--- /dev/null
@@ -0,0 +1,87 @@
+/* from: FreeBSD: head/lib/msun/src/e_acosh.c 176451 2008-02-22 02:30:36Z das */
+
+/* @(#)e_acosh.c 1.3 95/01/18 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ *
+ * $FreeBSD: head/lib/msun/src/e_acoshl.c 251599 2013-06-10 06:04:58Z das $
+ */
+
+/*
+ * See e_acosh.c for complete comments.
+ *
+ * Converted to long double by David Schultz <das@FreeBSD.ORG> and
+ * Bruce D. Evans.
+ */
+
+#include <float.h>
+#ifdef __i386__
+#include <ieeefp.h>
+#endif
+
+#include "fpmath.h"
+#include "math.h"
+#include "math_private.h"
+
+/* EXP_LARGE is the threshold above which we use acosh(x) ~= log(2x). */
+#if LDBL_MANT_DIG == 64
+#define        EXP_LARGE       34
+#elif LDBL_MANT_DIG == 113
+#define        EXP_LARGE       58
+#else
+#error "Unsupported long double format"
+#endif
+
+#if LDBL_MAX_EXP != 0x4000
+/* We also require the usual expsign encoding. */
+#error "Unsupported long double format"
+#endif
+
+#define        BIAS    (LDBL_MAX_EXP - 1)
+
+static const double
+one    = 1.0;
+
+#if LDBL_MANT_DIG == 64
+static const union IEEEl2bits
+u_ln2 =  LD80C(0xb17217f7d1cf79ac, -1, 6.93147180559945309417e-1L);
+#define        ln2     u_ln2.e
+#elif LDBL_MANT_DIG == 113
+static const long double
+ln2 =  6.93147180559945309417232121458176568e-1L;      /* 0x162e42fefa39ef35793c7673007e6.0p-113 */
+#else
+#error "Unsupported long double format"
+#endif
+
+long double
+acoshl(long double x)
+{
+       long double t;
+       int16_t hx;
+
+       ENTERI();
+       GET_LDBL_EXPSIGN(hx, x);
+       if (hx < 0x3fff) {              /* x < 1, or misnormal */
+           RETURNI((x-x)/(x-x));
+       } else if (hx >= BIAS + EXP_LARGE) { /* x >= LARGE */
+           if (hx >= 0x7fff) {         /* x is inf, NaN or misnormal */
+               RETURNI(x+x);
+           } else 
+               RETURNI(logl(x)+ln2);   /* acosh(huge)=log(2x), or misnormal */
+       } else if (hx == 0x3fff && x == 1) {
+           RETURNI(0.0);               /* acosh(1) = 0 */
+       } else if (hx >= 0x4000) {      /* LARGE > x >= 2, or misnormal */
+           t=x*x;
+           RETURNI(logl(2.0*x-one/(x+sqrtl(t-one))));
+       } else {                        /* 1<x<2 */
+           t = x-one;
+           RETURNI(log1pl(t+sqrtl(2.0*t+t*t)));
+       }
+}
index c91f1ad..12a6f87 100644 (file)
@@ -1,6 +1,6 @@
 
 /* @(#)e_atanh.c 1.3 95/01/18 */
-/* $FreeBSD: head/lib/msun/src/e_atanh.c 176451 2008-02-22 02:30:36Z das $ */
+/* $FreeBSD: head/lib/msun/src/e_atanh.c 251599 2013-06-10 06:04:58Z das $ */
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
@@ -31,6 +31,8 @@
  *
  */
 
+#include <float.h>
+
 #include "math.h"
 #include "math_private.h"
 
@@ -58,3 +60,7 @@ __ieee754_atanh(double x)
            t = 0.5*log1p((x+x)/(one-x));
        if(hx>=0) return t; else return -t;
 }
+
+#if LDBL_MANT_DIG == 53
+__weak_reference(atanh, atanhl);
+#endif
diff --git a/lib/libm/src/e_atanhl.c b/lib/libm/src/e_atanhl.c
new file mode 100644 (file)
index 0000000..be26e0c
--- /dev/null
@@ -0,0 +1,72 @@
+/* from: FreeBSD: head/lib/msun/src/e_atanh.c 176451 2008-02-22 02:30:36Z das */
+
+/* @(#)e_atanh.c 1.3 95/01/18 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ *
+ * $FreeBSD: head/lib/msun/src/e_atanhl.c 251599 2013-06-10 06:04:58Z das $
+ */
+
+/*
+ * See e_atanh.c for complete comments.
+ *
+ * Converted to long double by David Schultz <das@FreeBSD.ORG> and
+ * Bruce D. Evans.
+ */
+
+#include <float.h>
+#ifdef __i386__
+#include <ieeefp.h>
+#endif
+
+#include "fpmath.h"
+#include "math.h"
+#include "math_private.h"
+
+/* EXP_TINY is the threshold below which we use atanh(x) ~= x. */
+#if LDBL_MANT_DIG == 64
+#define        EXP_TINY        -34
+#elif LDBL_MANT_DIG == 113
+#define        EXP_TINY        -58
+#else
+#error "Unsupported long double format"
+#endif
+
+#if LDBL_MAX_EXP != 0x4000
+/* We also require the usual expsign encoding. */
+#error "Unsupported long double format"
+#endif
+
+#define        BIAS    (LDBL_MAX_EXP - 1)
+
+static const double one = 1.0, huge = 1e300;
+static const double zero = 0.0;
+
+long double
+atanhl(long double x)
+{
+       long double t;
+       uint16_t hx, ix;
+
+       ENTERI();
+       GET_LDBL_EXPSIGN(hx, x);
+       ix = hx & 0x7fff;
+       if (ix >= 0x3fff)               /* |x| >= 1, or NaN or misnormal */
+           RETURNI(fabsl(x) == 1 ? x / zero : (x - x) / (x - x));
+       if (ix < BIAS + EXP_TINY && (huge + x) > zero)
+           RETURNI(x);                 /* x is tiny */
+       SET_LDBL_EXPSIGN(x, ix);
+       if (ix < 0x3ffe) {              /* |x| < 0.5, or misnormal */
+           t = x+x;
+           t = 0.5*log1pl(t+t*x/(one-x));
+       } else 
+           t = 0.5*log1pl((x+x)/(one-x));
+       RETURNI((hx & 0x8000) == 0 ? t : -t);
+}
index dd875c1..2fc3f06 100644 (file)
@@ -1,6 +1,6 @@
 
 /* @(#)e_log.c 1.3 95/01/18 */
-/* $FreeBSD: head/lib/msun/src/e_log.c 251024 2013-05-27 08:50:10Z das $ */
+/* $FreeBSD: head/lib/msun/src/e_log.c 251292 2013-06-03 09:14:31Z das $
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
@@ -63,6 +63,8 @@
  * to produce the hexadecimal values shown.
  */
 
+#include <float.h>
+
 #include "math.h"
 #include "math_private.h"
 
@@ -137,3 +139,7 @@ __ieee754_log(double x)
                     return dk*ln2_hi-((s*(f-R)-dk*ln2_lo)-f);
        }
 }
+
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(log, logl);
+#endif
index d744fb5..ccb6373 100644 (file)
@@ -1,6 +1,6 @@
 
 /* @(#)e_log10.c 1.3 95/01/18 */
-/* $FreeBSD: head/lib/msun/src/e_log10.c 251024 2013-05-27 08:50:10Z das $ */
+/* $FreeBSD: head/lib/msun/src/e_log10.c 251292 2013-06-03 09:14:31Z das $ */
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
@@ -20,6 +20,8 @@
  * in not-quite-routine extra precision.
  */
 
+#include <float.h>
+
 #include "math.h"
 #include "math_private.h"
 #include "k_log.h"
@@ -84,3 +86,7 @@ __ieee754_log10(double x)
 
        return val_lo + val_hi;
 }
+
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(log10, log10l);
+#endif
index a7388e5..080e46b 100644 (file)
@@ -1,6 +1,6 @@
 
 /* @(#)e_log10.c 1.3 95/01/18 */
-/* $FreeBSD: head/lib/msun/src/e_log2.c 251024 2013-05-27 08:50:10Z das $ */
+/* $FreeBSD: head/lib/msun/src/e_log2.c 251404 2013-06-05 05:33:01Z das $ */
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
@@ -22,6 +22,8 @@
  * in not-quite-routine extra precision.
  */
 
+#include <float.h>
+
 #include "math.h"
 #include "math_private.h"
 #include "k_log.h"
@@ -107,3 +109,7 @@ __ieee754_log2(double x)
 
        return val_lo + val_hi;
 }
+
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(log2, log2l);
+#endif
index 3768ac3..f9befc7 100644 (file)
@@ -11,7 +11,7 @@
 
 /*
  * from: @(#)fdlibm.h 5.1 93/09/24
- * $FreeBSD: head/lib/msun/src/math.h 238722 2012-07-23 19:13:55Z kargl $
+ * $FreeBSD: head/lib/msun/src/math.h 251599 2013-06-10 06:04:58Z das $
  */
 
 #ifndef _MATH_H_
@@ -396,9 +396,12 @@ float      significandf(float);
  * long double versions of ISO/POSIX math functions
  */
 #if __ISO_C_VISIBLE >= 1999
+long double    acoshl(long double);
 long double    acosl(long double);
+long double    asinhl(long double);
 long double    asinl(long double);
 long double    atan2l(long double, long double);
+long double    atanhl(long double);
 long double    atanl(long double);
 long double    cbrtl(long double);
 long double    ceill(long double);
@@ -406,6 +409,7 @@ long double copysignl(long double, long double) __pure2;
 long double    cosl(long double);
 long double    exp2l(long double);
 long double    expl(long double);
+long double    expm1l(long double);
 long double    fabsl(long double) __pure2;
 long double    fdiml(long double, long double);
 long double    floorl(long double);
@@ -419,7 +423,11 @@ int                ilogbl(long double) __pure2;
 long double    ldexpl(long double, int);
 long long      llrintl(long double);
 long long      llroundl(long double);
+long double    log10l(long double);
+long double    log1pl(long double);
+long double    log2l(long double);
 long double    logbl(long double);
+long double    logl(long double);
 long           lrintl(long double);
 long           lroundl(long double);
 long double    modfl(long double, long double *); /* fundamentally !__pure2 */
@@ -457,18 +465,10 @@ __END_DECLS
  */
 __BEGIN_DECLS
 
-long double    acoshl(long double);
-long double    asinhl(long double);
-long double    atanhl(long double);
 long double    coshl(long double);
 long double    erfcl(long double);
 long double    erfl(long double);
-long double    expm1l(long double);
 long double    lgammal(long double);
-long double    log10l(long double);
-long double    log1pl(long double);
-long double    log2l(long double);
-long double    logl(long double);
 long double    powl(long double, long double);
 long double    sinhl(long double);
 long double    tanhl(long double);
index 34e8f59..d4623af 100644 (file)
@@ -11,7 +11,7 @@
 
 /*
  * from: @(#)fdlibm.h 5.1 93/09/24
- * $FreeBSD: head/lib/msun/src/math_private.h 241051 2012-09-29 16:40:12Z kargl $
+ * $FreeBSD: head/lib/msun/src/math_private.h 251292 2013-06-03 09:14:31Z das $
  */
 
 #ifndef _MATH_PRIVATE_H_
@@ -180,6 +180,33 @@ do {                                                               \
   (d) = sf_u.value;                                            \
 } while (0)
 
+/*
+ * Get expsign and mantissa as 16 bit and 64 bit ints from an 80 bit long
+ * double.
+ */
+
+#define        EXTRACT_LDBL80_WORDS(ix0,ix1,d)                         \
+do {                                                           \
+  union IEEEl2bits ew_u;                                       \
+  ew_u.e = (d);                                                        \
+  (ix0) = ew_u.xbits.expsign;                                  \
+  (ix1) = ew_u.xbits.man;                                      \
+} while (0)
+
+/*
+ * Get expsign and mantissa as one 16 bit and two 64 bit ints from a 128 bit
+ * long double.
+ */
+
+#define        EXTRACT_LDBL128_WORDS(ix0,ix1,ix2,d)                    \
+do {                                                           \
+  union IEEEl2bits ew_u;                                       \
+  ew_u.e = (d);                                                        \
+  (ix0) = ew_u.xbits.expsign;                                  \
+  (ix1) = ew_u.xbits.manh;                                     \
+  (ix2) = ew_u.xbits.manl;                                     \
+} while (0)
+
 /* Get expsign as a 16 bit int from a long double.  */
 
 #define        GET_LDBL_EXPSIGN(i,d)                                   \
@@ -189,6 +216,33 @@ do {                                                               \
   (i) = ge_u.xbits.expsign;                                    \
 } while (0)
 
+/*
+ * Set an 80 bit long double from a 16 bit int expsign and a 64 bit int
+ * mantissa.
+ */
+
+#define        INSERT_LDBL80_WORDS(d,ix0,ix1)                          \
+do {                                                           \
+  union IEEEl2bits iw_u;                                       \
+  iw_u.xbits.expsign = (ix0);                                  \
+  iw_u.xbits.man = (ix1);                                      \
+  (d) = iw_u.e;                                                        \
+} while (0)
+
+/*
+ * Set a 128 bit long double from a 16 bit int expsign and two 64 bit ints
+ * comprising the mantissa.
+ */
+
+#define        INSERT_LDBL128_WORDS(d,ix0,ix1,ix2)                     \
+do {                                                           \
+  union IEEEl2bits iw_u;                                       \
+  iw_u.xbits.expsign = (ix0);                                  \
+  iw_u.xbits.manh = (ix1);                                     \
+  iw_u.xbits.manl = (ix2);                                     \
+  (d) = iw_u.e;                                                        \
+} while (0)
+
 /* Set expsign of a long double from a 16 bit int.  */
 
 #define        SET_LDBL_EXPSIGN(d,v)                                   \
@@ -252,6 +306,110 @@ do {                                                              \
 /* Default return statement if hack*_t() is not used. */
 #define      RETURNF(v)      return (v)
 
+/*
+ * 2sum gives the same result as 2sumF without requiring |a| >= |b| or
+ * a == 0, but is slower.
+ */
+#define        _2sum(a, b) do {        \
+       __typeof(a) __s, __w;   \
+                               \
+       __w = (a) + (b);        \
+       __s = __w - (a);        \
+       (b) = ((a) - (__w - __s)) + ((b) - __s); \
+       (a) = __w;              \
+} while (0)
+
+/*
+ * 2sumF algorithm.
+ *
+ * "Normalize" the terms in the infinite-precision expression a + b for
+ * the sum of 2 floating point values so that b is as small as possible
+ * relative to 'a'.  (The resulting 'a' is the value of the expression in
+ * the same precision as 'a' and the resulting b is the rounding error.)
+ * |a| must be >= |b| or 0, b's type must be no larger than 'a's type, and
+ * exponent overflow or underflow must not occur.  This uses a Theorem of
+ * Dekker (1971).  See Knuth (1981) 4.2.2 Theorem C.  The name "TwoSum"
+ * is apparently due to Skewchuk (1997).
+ *
+ * For this to always work, assignment of a + b to 'a' must not retain any
+ * extra precision in a + b.  This is required by C standards but broken
+ * in many compilers.  The brokenness cannot be worked around using
+ * STRICT_ASSIGN() like we do elsewhere, since the efficiency of this
+ * algorithm would be destroyed by non-null strict assignments.  (The
+ * compilers are correct to be broken -- the efficiency of all floating
+ * point code calculations would be destroyed similarly if they forced the
+ * conversions.)
+ *
+ * Fortunately, a case that works well can usually be arranged by building
+ * any extra precision into the type of 'a' -- 'a' should have type float_t,
+ * double_t or long double.  b's type should be no larger than 'a's type.
+ * Callers should use these types with scopes as large as possible, to
+ * reduce their own extra-precision and efficiciency problems.  In
+ * particular, they shouldn't convert back and forth just to call here.
+ */
+#ifdef DEBUG
+#define        _2sumF(a, b) do {                               \
+       __typeof(a) __w;                                \
+       volatile __typeof(a) __ia, __ib, __r, __vw;     \
+                                                       \
+       __ia = (a);                                     \
+       __ib = (b);                                     \
+       assert(__ia == 0 || fabsl(__ia) >= fabsl(__ib));        \
+                                                       \
+       __w = (a) + (b);                                \
+       (b) = ((a) - __w) + (b);                        \
+       (a) = __w;                                      \
+                                                       \
+       /* The next 2 assertions are weak if (a) is already long double. */ \
+       assert((long double)__ia + __ib == (long double)(a) + (b));     \
+       __vw = __ia + __ib;                             \
+       __r = __ia - __vw;                              \
+       __r += __ib;                                    \
+       assert(__vw == (a) && __r == (b));              \
+} while (0)
+#else /* !DEBUG */
+#define        _2sumF(a, b) do {       \
+       __typeof(a) __w;        \
+                               \
+       __w = (a) + (b);        \
+       (b) = ((a) - __w) + (b); \
+       (a) = __w;              \
+} while (0)
+#endif /* DEBUG */
+
+/*
+ * Set x += c, where x is represented in extra precision as a + b.
+ * x must be sufficiently normalized and sufficiently larger than c,
+ * and the result is then sufficiently normalized.
+ *
+ * The details of ordering are that |a| must be >= |c| (so that (a, c)
+ * can be normalized without extra work to swap 'a' with c).  The details of
+ * the normalization are that b must be small relative to the normalized 'a'.
+ * Normalization of (a, c) makes the normalized c tiny relative to the
+ * normalized a, so b remains small relative to 'a' in the result.  However,
+ * b need not ever be tiny relative to 'a'.  For example, b might be about
+ * 2**20 times smaller than 'a' to give about 20 extra bits of precision.
+ * That is usually enough, and adding c (which by normalization is about
+ * 2**53 times smaller than a) cannot change b significantly.  However,
+ * cancellation of 'a' with c in normalization of (a, c) may reduce 'a'
+ * significantly relative to b.  The caller must ensure that significant
+ * cancellation doesn't occur, either by having c of the same sign as 'a',
+ * or by having |c| a few percent smaller than |a|.  Pre-normalization of
+ * (a, b) may help.
+ *
+ * This is is a variant of an algorithm of Kahan (see Knuth (1981) 4.2.2
+ * exercise 19).  We gain considerable efficiency by requiring the terms to
+ * be sufficiently normalized and sufficiently increasing.
+ */
+#define        _3sumF(a, b, c) do {    \
+       __typeof(a) __tmp;      \
+                               \
+       __tmp = (c);            \
+       _2sumF(__tmp, (a));     \
+       (b) += (a);             \
+       (a) = __tmp;            \
+} while (0)
+
 /*
  * Common routine to process the arguments to nan(), nanf(), and nanl().
  */
@@ -348,6 +506,7 @@ irint(double x)
 #define        HAVE_EFFICIENT_IRINT
 #endif
 
+#if defined(__x86_64__) || defined(__i386__)
 static __inline int
 irintl(long double x)
 {
@@ -357,9 +516,144 @@ irintl(long double x)
        return (n);
 }
 #define        HAVE_EFFICIENT_IRINTL
+#endif
 
 #endif /* __GNUC__ */
 
+#ifdef DEBUG
+#if defined(__amd64__) || defined(__i386__)
+#define        breakpoint()    asm("int $3")
+#else
+#include <signal.h>
+
+#define        breakpoint()    raise(SIGTRAP)
+#endif
+#endif
+
+/* Write a pari script to test things externally. */
+#ifdef DOPRINT
+#include <stdio.h>
+
+#ifndef DOPRINT_SWIZZLE
+#define        DOPRINT_SWIZZLE         0
+#endif
+
+#ifdef DOPRINT_LD80
+
+#define        DOPRINT_START(xp) do {                                          \
+       uint64_t __lx;                                                  \
+       uint16_t __hx;                                                  \
+                                                                       \
+       /* Hack to give more-problematic args. */                       \
+       EXTRACT_LDBL80_WORDS(__hx, __lx, *xp);                          \
+       __lx ^= DOPRINT_SWIZZLE;                                        \
+       INSERT_LDBL80_WORDS(*xp, __hx, __lx);                           \
+       printf("x = %.21Lg; ", (long double)*xp);                       \
+} while (0)
+#define        DOPRINT_END1(v)                                                 \
+       printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v))
+#define        DOPRINT_END2(hi, lo)                                            \
+       printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n",              \
+           (long double)(hi), (long double)(lo))
+
+#elif defined(DOPRINT_D64)
+
+#define        DOPRINT_START(xp) do {                                          \
+       uint32_t __hx, __lx;                                            \
+                                                                       \
+       EXTRACT_WORDS(__hx, __lx, *xp);                                 \
+       __lx ^= DOPRINT_SWIZZLE;                                        \
+       INSERT_WORDS(*xp, __hx, __lx);                                  \
+       printf("x = %.21Lg; ", (long double)*xp);                       \
+} while (0)
+#define        DOPRINT_END1(v)                                                 \
+       printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v))
+#define        DOPRINT_END2(hi, lo)                                            \
+       printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n",              \
+           (long double)(hi), (long double)(lo))
+
+#elif defined(DOPRINT_F32)
+
+#define        DOPRINT_START(xp) do {                                          \
+       uint32_t __hx;                                                  \
+                                                                       \
+       GET_FLOAT_WORD(__hx, *xp);                                      \
+       __hx ^= DOPRINT_SWIZZLE;                                        \
+       SET_FLOAT_WORD(*xp, __hx);                                      \
+       printf("x = %.21Lg; ", (long double)*xp);                       \
+} while (0)
+#define        DOPRINT_END1(v)                                                 \
+       printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v))
+#define        DOPRINT_END2(hi, lo)                                            \
+       printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n",              \
+           (long double)(hi), (long double)(lo))
+
+#else /* !DOPRINT_LD80 && !DOPRINT_D64 (LD128 only) */
+
+#ifndef DOPRINT_SWIZZLE_HIGH
+#define        DOPRINT_SWIZZLE_HIGH    0
+#endif
+
+#define        DOPRINT_START(xp) do {                                          \
+       uint64_t __lx, __llx;                                           \
+       uint16_t __hx;                                                  \
+                                                                       \
+       EXTRACT_LDBL128_WORDS(__hx, __lx, __llx, *xp);                  \
+       __llx ^= DOPRINT_SWIZZLE;                                       \
+       __lx ^= DOPRINT_SWIZZLE_HIGH;                                   \
+       INSERT_LDBL128_WORDS(*xp, __hx, __lx, __llx);                   \
+       printf("x = %.36Lg; ", (long double)*xp);                                       \
+} while (0)
+#define        DOPRINT_END1(v)                                                 \
+       printf("y = %.36Lg; z = 0; show(x, y, z);\n", (long double)(v))
+#define        DOPRINT_END2(hi, lo)                                            \
+       printf("y = %.36Lg; z = %.36Lg; show(x, y, z);\n",              \
+           (long double)(hi), (long double)(lo))
+
+#endif /* DOPRINT_LD80 */
+
+#else /* !DOPRINT */
+#define        DOPRINT_START(xp)
+#define        DOPRINT_END1(v)
+#define        DOPRINT_END2(hi, lo)
+#endif /* DOPRINT */
+
+#define        RETURNP(x) do {                 \
+       DOPRINT_END1(x);                \
+       RETURNF(x);                     \
+} while (0)
+#define        RETURNPI(x) do {                \
+       DOPRINT_END1(x);                \
+       RETURNI(x);                     \
+} while (0)
+#define        RETURN2P(x, y) do {             \
+       DOPRINT_END2((x), (y));         \
+       RETURNF((x) + (y));             \
+} while (0)
+#define        RETURN2PI(x, y) do {            \
+       DOPRINT_END2((x), (y));         \
+       RETURNI((x) + (y));             \
+} while (0)
+#ifdef STRUCT_RETURN
+#define        RETURNSP(rp) do {               \
+       if (!(rp)->lo_set)              \
+               RETURNP((rp)->hi);      \
+       RETURN2P((rp)->hi, (rp)->lo);   \
+} while (0)
+#define        RETURNSPI(rp) do {              \
+       if (!(rp)->lo_set)              \
+               RETURNPI((rp)->hi);     \
+       RETURN2PI((rp)->hi, (rp)->lo);  \
+} while (0)
+#endif
+#define        SUM2P(x, y) ({                  \
+       const __typeof (x) __x = (x);   \
+       const __typeof (y) __y = (y);   \
+                                       \
+       DOPRINT_END2(__x, __y);         \
+       __x + __y;                      \
+})
+
 /*
  * ieee style elementary functions
  *
index 50720bc..38ac00e 100644 (file)
@@ -1,5 +1,5 @@
 /* @(#)s_asinh.c 5.1 93/09/24 */
-/* $FreeBSD: head/lib/msun/src/s_asinh.c 176451 2008-02-22 02:30:36Z das $ */
+/* $FreeBSD: head/lib/msun/src/s_asinh.c 251599 2013-06-10 06:04:58Z das $
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
@@ -22,6 +22,8 @@
  *              := sign(x)*log1p(|x| + x^2/(1 + sqrt(1+x^2)))
  */
 
+#include <float.h>
+
 #include "math.h"
 #include "math_private.h"
 
@@ -52,3 +54,7 @@ asinh(double x)
        }
        if(hx>0) return w; else return -w;
 }
+
+#if LDBL_MANT_DIG == 53
+__weak_reference(asinh, asinhl);
+#endif
diff --git a/lib/libm/src/s_asinhl.c b/lib/libm/src/s_asinhl.c
new file mode 100644 (file)
index 0000000..d6577b8
--- /dev/null
@@ -0,0 +1,90 @@
+/* from: FreeBSD: head/lib/msun/src/e_acosh.c 176451 2008-02-22 02:30:36Z das */
+
+/* @(#)s_asinh.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ *
+ * $FreeBSD: head/lib/msun/src/s_asinhl.c 251599 2013-06-10 06:04:58Z das $
+ */
+
+/*
+ * See s_asinh.c for complete comments.
+ *
+ * Converted to long double by David Schultz <das@FreeBSD.ORG> and
+ * Bruce D. Evans.
+ */
+
+#include <float.h>
+#ifdef __i386__
+#include <ieeefp.h>
+#endif
+
+#include "fpmath.h"
+#include "math.h"
+#include "math_private.h"
+
+/* EXP_LARGE is the threshold above which we use asinh(x) ~= log(2x). */
+/* EXP_TINY is the threshold below which we use asinh(x) ~= x. */
+#if LDBL_MANT_DIG == 64
+#define        EXP_LARGE       34
+#define        EXP_TINY        -34
+#elif LDBL_MANT_DIG == 113
+#define        EXP_LARGE       58
+#define        EXP_TINY        -58
+#else
+#error "Unsupported long double format"
+#endif
+
+#if LDBL_MAX_EXP != 0x4000
+/* We also require the usual expsign encoding. */
+#error "Unsupported long double format"
+#endif
+
+#define        BIAS    (LDBL_MAX_EXP - 1)
+
+static const double
+one =  1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
+huge=  1.00000000000000000000e+300;
+
+#if LDBL_MANT_DIG == 64
+static const union IEEEl2bits
+u_ln2 =  LD80C(0xb17217f7d1cf79ac, -1, 6.93147180559945309417e-1L);
+#define        ln2     u_ln2.e
+#elif LDBL_MANT_DIG == 113
+static const long double
+ln2 =  6.93147180559945309417232121458176568e-1L;      /* 0x162e42fefa39ef35793c7673007e6.0p-113 */
+#else
+#error "Unsupported long double format"
+#endif
+
+long double
+asinhl(long double x)
+{
+       long double t, w;
+       uint16_t hx, ix;
+
+       ENTERI();
+       GET_LDBL_EXPSIGN(hx, x);
+       ix = hx & 0x7fff;
+       if (ix >= 0x7fff) RETURNI(x+x); /* x is inf, NaN or misnormal */
+       if (ix < BIAS + EXP_TINY) {     /* |x| < TINY, or misnormal */
+           if (huge + x > one) RETURNI(x);     /* return x inexact except 0 */
+       }
+       if (ix >= BIAS + EXP_LARGE) {   /* |x| >= LARGE, or misnormal */
+           w = logl(fabsl(x))+ln2;
+       } else if (ix >= 0x4000) {      /* LARGE > |x| >= 2.0, or misnormal */
+           t = fabsl(x);
+           w = logl(2.0*t+one/(sqrtl(x*x+one)+t));
+       } else {                /* 2.0 > |x| >= TINY, or misnormal */
+           t = x*x;
+           w =log1pl(fabsl(x)+t/(one+sqrtl(one+t)));
+       }
+       RETURNI((hx & 0x8000) == 0 ? w : -w);
+}
index dfae2af..89ac674 100644 (file)
@@ -1,5 +1,5 @@
 /* @(#)s_expm1.c 5.1 93/09/24 */
-/* $FreeBSD: head/lib/msun/src/s_expm1.c 251024 2013-05-27 08:50:10Z das $ */
+/* $FreeBSD: head/lib/msun/src/s_expm1.c 251343 2013-06-03 19:51:32Z kargl $
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
@@ -214,3 +214,7 @@ expm1(double x)
        }
        return y;
 }
+
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(expm1, expm1l);
+#endif
index 7d68280..7cf2a40 100644 (file)
@@ -1,5 +1,5 @@
 /* @(#)s_log1p.c 5.1 93/09/24 */
-/* $FreeBSD: head/lib/msun/src/s_log1p.c 251024 2013-05-27 08:50:10Z das $ */
+/* $FreeBSD: head/lib/msun/src/s_log1p.c 251292 2013-06-03 09:14:31Z das $ */
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
@@ -172,3 +172,7 @@ log1p(double x)
        if(k==0) return f-(hfsq-s*(hfsq+R)); else
                 return k*ln2_hi-((hfsq-(s*(hfsq+R)+(k*ln2_lo+c)))-f);
 }
+
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(log1p, log1pl);
+#endif
index abd2e85..b69d369 100644 (file)
@@ -6,3 +6,5 @@ ARCH_SRCS = e_remainder.S e_remainderf.S e_remainderl.S \
            s_logbl.S s_lrint.S s_lrintf.S s_lrintl.S \
            s_remquo.S s_remquof.S s_remquol.S \
            s_rintl.S s_scalbn.S s_scalbnf.S s_scalbnl.S
+
+SYM_MAPS += ${.CURDIR}/x86_64/Symbol.map
diff --git a/lib/libm/x86_64/Symbol.map b/lib/libm/x86_64/Symbol.map
new file mode 100644 (file)
index 0000000..e889f36
--- /dev/null
@@ -0,0 +1,12 @@
+/*
+ * $FreeBSD: head/lib/msun/amd64/Symbol.map 226218 2011-10-10 15:43:09Z das $
+ */
+DFLY36.0 {
+       fesetexceptflag;
+       feraiseexcept;
+       fegetenv;
+       feholdexcept;
+       feupdateenv;
+       feenableexcept;
+       fedisableexcept;
+};
index 336a470..3800071 100644 (file)
@@ -109,6 +109,15 @@ PO_CXXFLAGS=${CXXFLAGS:N-ffunction-sections}
 
 all: objwarn
 
+.include <bsd.symver.mk>
+
+# Allow libraries to specify their own version map or have it
+# automatically generated (see bsd.symver.mk above).
+.if !defined(NO_SYMVER) && !empty(VERSION_MAP)
+${SHLIB_NAME}: ${VERSION_MAP}
+LDFLAGS+=      -Wl,--version-script=${VERSION_MAP}
+.endif
+
 .include <bsd.patch.mk>
 
 .if defined(LIB) && !empty(LIB) || defined(SHLIB_NAME)
diff --git a/share/mk/bsd.symver.mk b/share/mk/bsd.symver.mk
new file mode 100644 (file)
index 0000000..29d3471
--- /dev/null
@@ -0,0 +1,29 @@
+# $FreeBSD$
+
+.if !target(__<bsd.symver.mk>__)
+__<bsd.symver.mk>__:
+
+.include <bsd.init.mk>
+
+# Generate the version map given the version definitions
+# and symbol maps.
+.if !defined(NO_SYMVER) && !empty(VERSION_DEF) && !empty(SYMBOL_MAPS)
+# Find the awk script that generates the version map.
+VERSION_GEN?=  version_gen.awk
+VERSION_MAP?=  Version.map
+
+CLEANFILES+=   ${VERSION_MAP}
+
+.if exists(${.PARSEDIR}/${VERSION_GEN})
+_vgen:=        ${.PARSEDIR}/${VERSION_GEN}
+.else
+.error ${VERSION_GEN} not found in ${.PARSEDIR}
+.endif
+
+# Run the symbol maps through the C preprocessor before passing
+# them to the symbol version generator.
+${VERSION_MAP}: ${VERSION_DEF} ${_vgen} ${SYMBOL_MAPS}
+       cat ${SYMBOL_MAPS} | ${CPP} - - \
+           | awk -v vfile=${VERSION_DEF} -f ${_vgen} > ${.TARGET}
+.endif # !empty(VERSION_DEF) && !empty(SYMBOL_MAPS)
+.endif  # !target(__<bsd.symver.mk>__)
diff --git a/share/mk/version_gen.awk b/share/mk/version_gen.awk
new file mode 100644 (file)
index 0000000..93fbc4f
--- /dev/null
@@ -0,0 +1,249 @@
+#
+# Copyright (C) 2006 Daniel M. Eischen.  All rights reserved.
+# 
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions
+# are met:
+# 1. Redistributions of source code must retain the above copyright
+#    notice, this list of conditions and the following disclaimer.
+# 2. Redistributions in binary form must reproduce the above copyright
+#    notice, this list of conditions and the following disclaimer in the
+#    documentation and/or other materials provided with the distribution.
+# 
+# THIS SOFTWARE IS PROVIDED BY AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+# ARE DISCLAIMED.  IN NO EVENT SHALL AUTHOR OR CONTRIBUTORS BE LIABLE
+# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+# OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+# HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+# OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+# SUCH DAMAGE.
+#
+# $FreeBSD$
+#
+
+#
+# Make a list of all the library versions listed in the master file.
+#
+#   versions[] - array indexed by version name, contains number
+#                of symbols (+ 1) found for each version.
+#   successors[] - array index by version name, contains successor
+#                  version name.
+#   symbols[][] - array index by [version name, symbol index], contains
+#                 names of symbols defined for each version.
+#   names[] - array index is symbol name and value is its first version seen,
+#            used to check for duplicate symbols and warn about them.
+#
+BEGIN {
+       brackets = 0;
+       errors = warns = 0;
+       version_count = 0;
+       current_version = "";
+       stderr = "/dev/stderr";
+       while (getline < vfile) {
+               # Strip comments.
+               sub("#.*$", "", $0);
+
+               # Strip leading and trailing whitespace.
+               sub("^[ \t]+", "", $0);
+               sub("[ \t]+$", "", $0);
+
+               if (/^[a-zA-Z0-9._]+[ \t]*{$/) {
+                       # Strip brace.
+                       sub("{", "", $1);
+                       brackets++;
+                       symver = $1;
+                       versions[symver] = 1;
+                       successors[symver] = "";
+                       generated[symver] = 0;
+                       version_count++;
+               }
+               else if (/^}[ \t]*[a-zA-Z0-9._]+[ \t]*;$/) {
+                       v = $1 != "}" ? $1 : $2;
+                       # Strip brace.
+                       sub("}", "", v);
+                       # Strip semicolon.
+                       sub(";", "", v);
+                       if (symver == "") {
+                               printf("File %s: Unmatched bracket.\n",
+                               vfile) > stderr;
+                               errors++;
+                       }
+                       else if (versions[v] != 1) {
+                               printf("File %s: `%s' has unknown " \
+                                   "successor `%s'.\n",
+                                   vfile, symver, v) > stderr;
+                               errors++;
+                       }
+                       else
+                               successors[symver] = v;
+                       brackets--;
+               }
+               else if (/^}[ \t]*;$/) {
+                       if (symver == "") {
+                               printf("File %s: Unmatched bracket.\n",
+                                   vfile) > stderr;
+                               errors++;
+                       }
+                       # No successor
+                       brackets--;
+               }
+               else if (/^}$/) {
+                       printf("File %s: Missing final semicolon.\n",
+                           vfile) > stderr;
+                       errors++;
+               }
+               else if (/^$/)
+                       ;  # Ignore blank lines.
+               else {
+                       printf("File %s: Unknown directive: `%s'.\n",
+                           vfile, $0) > stderr;
+                       errors++;
+               }
+       }
+       brackets = 0;
+}
+
+{
+       # Set meaningful filename for diagnostics.
+       filename = FILENAME != "" ? FILENAME : "<stdin>";
+
+       # Delete comments, preceding and trailing whitespace, then
+       # consume blank lines.
+       sub("#.*$", "", $0);
+       sub("^[ \t]+", "", $0);
+       sub("[ \t]+$", "", $0);
+       if ($0 == "")
+               next;
+}
+
+/^[a-zA-Z0-9._]+[ \t]*{$/ {
+       # Strip bracket from version name.
+       sub("{", "", $1);
+       if (current_version != "") {
+               printf("File %s, line %d: Illegal nesting detected.\n",
+                   filename, FNR) > stderr;
+               errors++;
+       }
+       else if (versions[$1] == 0) {
+               printf("File %s, line %d: Undefined " \
+                   "library version `%s'.\n", filename, FNR, $1) > stderr;
+               errors++;
+               # Remove this entry from the versions.
+               delete versions[$1];
+       }
+       else
+               current_version = $1;
+       brackets++;
+       next;
+}
+
+/^[a-zA-Z0-9._]+[ \t]*;$/ {
+       # Strip semicolon.
+       sub(";", "", $1);
+       if (current_version != "") {
+               count = versions[current_version];
+               versions[current_version]++;
+               symbols[current_version, count] = $1;
+               if ($1 in names && names[$1] != current_version) {
+                       #
+                       # A graver case when a dup symbol appears under
+                       # different versions in the map.  That can result
+                       # in subtle problems with the library later.
+                       #
+                       printf("File %s, line %d: Duplicated symbol `%s' " \
+                           "in version `%s', first seen in `%s'. " \
+                           "Did you forget to move it to ObsoleteVersions?\n",
+                           filename, FNR, $1,
+                           current_version, names[$1]) > stderr;
+                       errors++;
+               }
+               else if (names[$1] == current_version) {
+                       #
+                       # A harmless case: a dup symbol with the same version.
+                       #
+                       printf("File %s, line %d: warning: " \
+                           "Duplicated symbol `%s' in version `%s'.\n",
+                           filename, FNR, $1, current_version) > stderr;
+                       warns++;
+               }
+               else
+                       names[$1] = current_version;
+       }
+       else {
+               printf("File %s, line %d: Symbol `%s' outside version scope.\n",
+                   filename, FNR, $1) > stderr;
+               errors++;
+       }
+       next;
+}
+
+/^}[ \t]*;$/ {
+       brackets--;
+       if (brackets < 0) {
+               printf("File %s, line %d: Unmatched bracket.\n",
+                   filename, FNR, $1) > stderr;
+               errors++;
+               brackets = 0;   # Reset
+       }
+       current_version = "";
+       next;
+}
+
+
+{
+       printf("File %s, line %d: Unknown directive: `%s'.\n",
+           filename, FNR, $0) > stderr;
+       errors++;
+}
+
+function print_version(v)
+{
+       # This function is recursive, so return if this version
+       # has already been printed.  Otherwise, if there is an
+       # ancestral version, recursively print its symbols before
+       # printing the symbols for this version.
+       #
+       if (generated[v] == 1)
+               return;
+       if (successors[v] != "")
+               print_version(successors[v]);
+
+       printf("%s {\n", v);
+
+       # The version count is always one more that actual,
+       # so the loop ranges from 1 to n-1.
+       #
+       for (i = 1; i < versions[v]; i++) {
+               if (i == 1)
+                       printf("global:\n");
+               printf("\t%s;\n", symbols[v, i]);
+       }
+
+       version_count--;
+       if (version_count == 0) {
+               printf("local:\n");
+               printf("\t*;\n");
+       }
+       if (successors[v] == "")
+               printf("};\n");
+       else
+               printf("} %s;\n", successors[v]);
+       printf("\n");
+
+       generated[v] = 1;
+    }
+
+END {
+       if (errors) {
+               printf("%d error(s) total.\n", errors) > stderr;
+               exit(1);
+       }
+       # OK, no errors.
+       for (v in versions) {
+               print_version(v);
+       }
+}