libm: Update erf, add imprecise versions of missing c++11 functions
authorJohn Marino <draco@marino.st>
Sun, 29 Sep 2013 23:29:20 +0000 (01:29 +0200)
committerJohn Marino <draco@marino.st>
Sun, 29 Sep 2013 23:32:21 +0000 (01:32 +0200)
Provide updates to erf and erff functions.

Also add weak versions of the missing c++11 long double functions by
using taking arguments of type double.  Use of these versions will
result in a linker warning to discourage program that really need
extra precision from using them.

Note that since the c/c++ specs only guarantee that long double has
precision equal to double, code that relies on these functions having
greater precision is unportable at best and broken at worst.

Taken-from: FreeBSD

lib/libc/Versions.def
lib/libm/Makefile
lib/libm/Symbol.map
lib/libm/src/imprecise.c [new file with mode: 0644]
lib/libm/src/math.h
lib/libm/src/s_erf.c
lib/libm/src/s_erff.c

index db7cf87..0c9cd87 100644 (file)
@@ -9,6 +9,10 @@
 DF306.0 {
 };
 
+# Second version for DragonFly 3.5 on per library basis
+DF306.1 {
+} DF306.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
@@ -16,4 +20,4 @@ DF306.0 {
 #
 # Please do NOT increment the version of this namespace.
 DFprivate_1.0 {
-} DF306.0;
+} DF306.1;
index 1d1f7c1..1ea5083 100644 (file)
@@ -36,6 +36,7 @@ COMMON_SRCS= b_exp.c b_log.c b_tgamma.c \
        e_pow.c e_powf.c e_rem_pio2.c \
        e_rem_pio2f.c e_scalb.c e_scalbf.c \
        e_sinh.c e_sinhf.c fenv.c \
+       imprecise.c \
        k_cos.c k_cosf.c k_exp.c k_expf.c k_rem_pio2.c k_sin.c k_sinf.c \
        k_tan.c k_tanf.c \
        s_asinh.c s_asinhf.c s_atan.c s_atanf.c s_carg.c s_cargf.c s_cargl.c \
index bea5739..6418f90 100644 (file)
@@ -258,3 +258,18 @@ DF306.0 {
        log2l;
        logl;
 };
+
+/*
+ * Implemented as weak aliases for imprecise versions
+ */
+
+DF306.1 {
+       coshl;
+       erfcl;
+       erfl;
+       lgammal;
+       powl;
+       sinhl;
+       tanhl;
+       tgammal;
+};
diff --git a/lib/libm/src/imprecise.c b/lib/libm/src/imprecise.c
new file mode 100644 (file)
index 0000000..2498be3
--- /dev/null
@@ -0,0 +1,69 @@
+/*-
+ * Copyright (c) 2013 David Chisnall
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ *
+ * $FreeBSD: head/lib/msun/src/imprecise.c 255294 2013-09-06 07:58:23Z theraven $
+ */
+
+#include <float.h>
+#include <math.h>
+
+/*
+ * If long double is not the same size as double, then these will lose
+ * precision and we should emit a warning whenever something links against
+ * them.
+ */
+#if (LDBL_MANT_DIG > 53)
+#define WARN_IMPRECISE(x) \
+       __warn_references(x, # x " has lower than advertised precision");
+#else
+#define WARN_IMPRECISE(x)
+#endif
+/*
+ * Declare the functions as weak variants so that other libraries providing
+ * real versions can override them.
+ */
+#define        DECLARE_WEAK(x)\
+       __weak_reference(imprecise_## x, x);\
+       WARN_IMPRECISE(x)
+
+long double
+imprecise_powl(long double x, long double y)
+{
+
+       return pow(x, y);
+}
+DECLARE_WEAK(powl);
+
+#define DECLARE_IMPRECISE(f) \
+       long double imprecise_ ## f ## l(long double v) { return f(v); }\
+       DECLARE_WEAK(f ## l)
+
+DECLARE_IMPRECISE(cosh);
+DECLARE_IMPRECISE(erfc);
+DECLARE_IMPRECISE(erf);
+DECLARE_IMPRECISE(lgamma);
+DECLARE_IMPRECISE(sinh);
+DECLARE_IMPRECISE(tanh);
+DECLARE_IMPRECISE(tgamma);
index f9befc7..22c2553 100644 (file)
@@ -11,7 +11,7 @@
 
 /*
  * from: @(#)fdlibm.h 5.1 93/09/24
- * $FreeBSD: head/lib/msun/src/math.h 251599 2013-06-10 06:04:58Z das $
+ * $FreeBSD: head/lib/msun/src/math.h 253766 2013-07-29 12:33:03Z theraven $
  */
 
 #ifndef _MATH_H_
@@ -76,27 +76,43 @@ extern const union __nan_un {
 #define        FP_NORMAL       0x04
 #define        FP_SUBNORMAL    0x08
 #define        FP_ZERO         0x10
+
+#if (__STDC_VERSION__ >= 201112L && defined(__clang__)) || \
+    __has_extension(c_generic_selections)
+#define        __fp_type_select(x, f, d, ld) _Generic((x),                     \
+    float: f(x),                                                       \
+    double: d(x),                                                      \
+    long double: ld(x),                                                        \
+    volatile float: f(x),                                              \
+    volatile double: d(x),                                             \
+    volatile long double: ld(x),                                       \
+    volatile const float: f(x),                                                \
+    volatile const double: d(x),                                       \
+    volatile const long double: ld(x),                                 \
+    const float: f(x),                                                 \
+    const double: d(x),                                                        \
+    const long double: ld(x))
+#elif __GNUC_PREREQ__(3, 1) && !defined(__cplusplus)
+#define        __fp_type_select(x, f, d, ld) __builtin_choose_expr(            \
+    __builtin_types_compatible_p(__typeof(x), long double), ld(x),     \
+    __builtin_choose_expr(                                             \
+    __builtin_types_compatible_p(__typeof(x), double), d(x),           \
+    __builtin_choose_expr(                                             \
+    __builtin_types_compatible_p(__typeof(x), float), f(x), (void)0)))
+#else
+#define         __fp_type_select(x, f, d, ld)                                  \
+    ((sizeof(x) == sizeof(float)) ? f(x)                               \
+    : (sizeof(x) == sizeof(double)) ? d(x)                             \
+    : ld(x))
+#endif
+
 #define        fpclassify(x) \
-    ((sizeof (x) == sizeof (float)) ? __fpclassifyf(x) \
-    : (sizeof (x) == sizeof (double)) ? __fpclassifyd(x) \
-    : __fpclassifyl(x))
-
-#define        isfinite(x)                                     \
-    ((sizeof (x) == sizeof (float)) ? __isfinitef(x)   \
-    : (sizeof (x) == sizeof (double)) ? __isfinite(x)  \
-    : __isfinitel(x))
-#define        isinf(x)                                        \
-    ((sizeof (x) == sizeof (float)) ? __isinff(x)      \
-    : (sizeof (x) == sizeof (double)) ? isinf(x)       \
-    : __isinfl(x))
-#define        isnan(x)                                        \
-    ((sizeof (x) == sizeof (float)) ? __isnanf(x)      \
-    : (sizeof (x) == sizeof (double)) ? isnan(x)       \
-    : __isnanl(x))
-#define        isnormal(x)                                     \
-    ((sizeof (x) == sizeof (float)) ? __isnormalf(x)   \
-    : (sizeof (x) == sizeof (double)) ? __isnormal(x)  \
-    : __isnormall(x))
+       __fp_type_select(x, __fpclassifyf, __fpclassifyd, __fpclassifyl)
+#define        isfinite(x) __fp_type_select(x, __isfinitef, __isfinite, __isfinitel)
+#define        isinf(x) __fp_type_select(x, __isinff, __isinf, __isinfl)
+#define        isnan(x) \
+       __fp_type_select(x, __inline_isnanf, __inline_isnan, __inline_isnanl)
+#define        isnormal(x) __fp_type_select(x, __isnormalf, __isnormal, __isnormall)
 
 #ifdef __MATH_BUILTIN_RELOPS
 #define        isgreater(x, y)         __builtin_isgreater((x), (y))
@@ -115,10 +131,7 @@ extern const union __nan_un {
 #define        isunordered(x, y)       (isnan(x) || isnan(y))
 #endif /* __MATH_BUILTIN_RELOPS */
 
-#define        signbit(x)                                      \
-    ((sizeof (x) == sizeof (float)) ? __signbitf(x)    \
-    : (sizeof (x) == sizeof (double)) ? __signbit(x)   \
-    : __signbitl(x))
+#define        signbit(x) __fp_type_select(x, __signbitf, __signbit, __signbitl)
 
 #ifdef __LP64__
 typedef        double          double_t;
@@ -176,9 +189,8 @@ int __isfinitef(float) __pure2;
 int    __isfinite(double) __pure2;
 int    __isfinitel(long double) __pure2;
 int    __isinff(float) __pure2;
+int    __isinf(double) __pure2;
 int    __isinfl(long double) __pure2;
-int    __isnanf(float) __pure2;
-int    __isnanl(long double) __pure2;
 int    __isnormalf(float) __pure2;
 int    __isnormal(double) __pure2;
 int    __isnormall(long double) __pure2;
@@ -186,6 +198,42 @@ int        __signbit(double) __pure2;
 int    __signbitf(float) __pure2;
 int    __signbitl(long double) __pure2;
 
+static __inline int
+__inline_isnan(__const double __x)
+{
+
+       return (__x != __x);
+}
+
+static __inline int
+__inline_isnanf(__const float __x)
+{
+
+       return (__x != __x);
+}
+
+static __inline int
+__inline_isnanl(__const long double __x)
+{
+
+       return (__x != __x);
+}
+
+/*
+ * Version 2 of the Single UNIX Specification (UNIX98) defined isnan() and
+ * isinf() as functions taking double.  C99, and the subsequent POSIX revisions
+ * (SUSv3, POSIX.1-2001, define it as a macro that accepts any real floating
+ * point type.  If we are targeting SUSv2 and C99 or C11 (or C++11) then we
+ * expose the newer definition, assuming that the language spec takes
+ * precedence over the operating system interface spec.
+ */
+#if    __XSI_VISIBLE > 0 && __XSI_VISIBLE < 600 && __ISO_C_VISIBLE < 1999
+#undef isinf
+#undef isnan
+int    isinf(double);
+int    isnan(double);
+#endif
+
 double acos(double);
 double asin(double);
 double atan(double);
@@ -228,8 +276,6 @@ double      expm1(double);
 double fma(double, double, double);
 double hypot(double, double);
 int    ilogb(double) __pure2;
-int    (isinf)(double) __pure2;
-int    (isnan)(double) __pure2;
 double lgamma(double);
 long long llrint(double);
 long long llround(double);
index fb0204d..d018c3c 100644 (file)
@@ -1,5 +1,4 @@
 /* @(#)s_erf.c 5.1 93/09/24 */
-/* $FreeBSD: head/lib/msun/src/s_erf.c 176451 2008-02-22 02:30:36Z das $ */
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
@@ -8,9 +7,12 @@
  * Permission to use, copy, modify, and distribute this
  * software is freely granted, provided that this notice
  * is preserved.
+ *
+ * $FreeBSD: head/lib/msun/src/s_erf.c 254994 2013-08-28 16:59:55Z kargl $
  * ====================================================
  */
 
+
 /* double erf(double x)
  * double erfc(double x)
  *                          x
@@ -199,7 +201,7 @@ erf(double x)
        if(ix < 0x3feb0000) {           /* |x|<0.84375 */
            if(ix < 0x3e300000) {       /* |x|<2**-28 */
                if (ix < 0x00800000)
-                   return 0.125*(8.0*x+efx8*x);  /*avoid underflow */
+                   return (8*x+efx8*x)/8;      /* avoid spurious underflow */
                return x + efx*x;
            }
            z = x*x;
index e0f180d..eece494 100644 (file)
  * is preserved.
  * ====================================================
  *
- * $FreeBSD: head/lib/msun/src/s_erff.c 176451 2008-02-22 02:30:36Z das $
+ * $FreeBSD: head/lib/msun/src/s_erff.c 254969 2013-08-27 19:46:56Z kargl $
  */
 
+
 #include "math.h"
 #include "math_private.h"
 
@@ -23,75 +24,59 @@ tiny            = 1e-30,
 half=  5.0000000000e-01, /* 0x3F000000 */
 one =  1.0000000000e+00, /* 0x3F800000 */
 two =  2.0000000000e+00, /* 0x40000000 */
-       /* c = (subfloat)0.84506291151 */
-erx =  8.4506291151e-01, /* 0x3f58560b */
 /*
- * Coefficients for approximation to  erf on [0,0.84375]
+ * Coefficients for approximation to erf on [0,0.84375]
  */
 efx =  1.2837916613e-01, /* 0x3e0375d4 */
 efx8=  1.0270333290e+00, /* 0x3f8375d4 */
-pp0  =  1.2837916613e-01, /* 0x3e0375d4 */
-pp1  = -3.2504209876e-01, /* 0xbea66beb */
-pp2  = -2.8481749818e-02, /* 0xbce9528f */
-pp3  = -5.7702702470e-03, /* 0xbbbd1489 */
-pp4  = -2.3763017452e-05, /* 0xb7c756b1 */
-qq1  =  3.9791721106e-01, /* 0x3ecbbbce */
-qq2  =  6.5022252500e-02, /* 0x3d852a63 */
-qq3  =  5.0813062117e-03, /* 0x3ba68116 */
-qq4  =  1.3249473704e-04, /* 0x390aee49 */
-qq5  = -3.9602282413e-06, /* 0xb684e21a */
 /*
- * Coefficients for approximation to  erf  in [0.84375,1.25]
+ *  Domain [0, 0.84375], range ~[-5.4446e-10,5.5197e-10]:
+ *  |(erf(x) - x)/x - p(x)/q(x)| < 2**-31.
+ */
+pp0  =  1.28379166e-01F, /*  0x1.06eba8p-3 */
+pp1  = -3.36030394e-01F, /* -0x1.58185ap-2 */
+pp2  = -1.86260219e-03F, /* -0x1.e8451ep-10 */
+qq1  =  3.12324286e-01F, /*  0x1.3fd1f0p-2 */
+qq2  =  2.16070302e-02F, /*  0x1.620274p-6 */
+qq3  = -1.98859419e-03F, /* -0x1.04a626p-9 */
+/*
+ * Domain [0.84375, 1.25], range ~[-1.953e-11,1.940e-11]:
+ * |(erf(x) - erx) - p(x)/q(x)| < 2**-36.
  */
-pa0  = -2.3621185683e-03, /* 0xbb1acdc6 */
-pa1  =  4.1485610604e-01, /* 0x3ed46805 */
-pa2  = -3.7220788002e-01, /* 0xbebe9208 */
-pa3  =  3.1834661961e-01, /* 0x3ea2fe54 */
-pa4  = -1.1089469492e-01, /* 0xbde31cc2 */
-pa5  =  3.5478305072e-02, /* 0x3d1151b3 */
-pa6  = -2.1663755178e-03, /* 0xbb0df9c0 */
-qa1  =  1.0642088205e-01, /* 0x3dd9f331 */
-qa2  =  5.4039794207e-01, /* 0x3f0a5785 */
-qa3  =  7.1828655899e-02, /* 0x3d931ae7 */
-qa4  =  1.2617121637e-01, /* 0x3e013307 */
-qa5  =  1.3637083583e-02, /* 0x3c5f6e13 */
-qa6  =  1.1984500103e-02, /* 0x3c445aa3 */
+erx  =  8.42697144e-01F, /*  0x1.af7600p-1.  erf(1) rounded to 16 bits. */
+pa0  =  3.64939137e-06F, /*  0x1.e9d022p-19 */
+pa1  =  4.15109694e-01F, /*  0x1.a91284p-2 */
+pa2  = -1.65179938e-01F, /* -0x1.5249dcp-3 */
+pa3  =  1.10914491e-01F, /*  0x1.c64e46p-4 */
+qa1  =  6.02074385e-01F, /*  0x1.344318p-1 */
+qa2  =  5.35934687e-01F, /*  0x1.126608p-1 */
+qa3  =  1.68576106e-01F, /*  0x1.593e6ep-3 */
+qa4  =  5.62181212e-02F, /*  0x1.cc89f2p-5 */
 /*
- * Coefficients for approximation to  erfc in [1.25,1/0.35]
+ * Domain [1.25,1/0.35], range ~[-7.043e-10,7.457e-10]:
+ * |log(x*erfc(x)) + x**2 + 0.5625 - r(x)/s(x)| < 2**-30
  */
-ra0  = -9.8649440333e-03, /* 0xbc21a093 */
-ra1  = -6.9385856390e-01, /* 0xbf31a0b7 */
-ra2  = -1.0558626175e+01, /* 0xc128f022 */
-ra3  = -6.2375331879e+01, /* 0xc2798057 */
-ra4  = -1.6239666748e+02, /* 0xc322658c */
-ra5  = -1.8460508728e+02, /* 0xc3389ae7 */
-ra6  = -8.1287437439e+01, /* 0xc2a2932b */
-ra7  = -9.8143291473e+00, /* 0xc11d077e */
-sa1  =  1.9651271820e+01, /* 0x419d35ce */
-sa2  =  1.3765776062e+02, /* 0x4309a863 */
-sa3  =  4.3456588745e+02, /* 0x43d9486f */
-sa4  =  6.4538726807e+02, /* 0x442158c9 */
-sa5  =  4.2900814819e+02, /* 0x43d6810b */
-sa6  =  1.0863500214e+02, /* 0x42d9451f */
-sa7  =  6.5702495575e+00, /* 0x40d23f7c */
-sa8  = -6.0424413532e-02, /* 0xbd777f97 */
+ra0  = -9.87132732e-03F, /* -0x1.4376b2p-7 */
+ra1  = -5.53605914e-01F, /* -0x1.1b723cp-1 */
+ra2  = -2.17589188e+00F, /* -0x1.1683a0p+1 */
+ra3  = -1.43268085e+00F, /* -0x1.6ec42cp+0 */
+sa1  =  5.45995426e+00F, /*  0x1.5d6fe4p+2 */
+sa2  =  6.69798088e+00F, /*  0x1.acabb8p+2 */
+sa3  =  1.43113089e+00F, /*  0x1.6e5e98p+0 */
+sa4  = -5.77397496e-02F, /* -0x1.d90108p-5 */
 /*
- * Coefficients for approximation to  erfc in [1/.35,28]
+ * Domain [1/0.35, 11], range ~[-2.264e-13,2.336e-13]:
+ * |log(x*erfc(x)) + x**2 + 0.5625 - r(x)/s(x)| < 2**-42
  */
-rb0  = -9.8649431020e-03, /* 0xbc21a092 */
-rb1  = -7.9928326607e-01, /* 0xbf4c9dd4 */
-rb2  = -1.7757955551e+01, /* 0xc18e104b */
-rb3  = -1.6063638306e+02, /* 0xc320a2ea */
-rb4  = -6.3756646729e+02, /* 0xc41f6441 */
-rb5  = -1.0250950928e+03, /* 0xc480230b */
-rb6  = -4.8351919556e+02, /* 0xc3f1c275 */
-sb1  =  3.0338060379e+01, /* 0x41f2b459 */
-sb2  =  3.2579251099e+02, /* 0x43a2e571 */
-sb3  =  1.5367296143e+03, /* 0x44c01759 */
-sb4  =  3.1998581543e+03, /* 0x4547fdbb */
-sb5  =  2.5530502930e+03, /* 0x451f90ce */
-sb6  =  4.7452853394e+02, /* 0x43ed43a7 */
-sb7  = -2.2440952301e+01; /* 0xc1b38712 */
+rb0  = -9.86494310e-03F, /* -0x1.434124p-7 */
+rb1  = -6.25171244e-01F, /* -0x1.401672p-1 */
+rb2  = -6.16498327e+00F, /* -0x1.8a8f16p+2 */
+rb3  = -1.66696873e+01F, /* -0x1.0ab70ap+4 */
+rb4  = -9.53764343e+00F, /* -0x1.313460p+3 */
+sb1  =  1.26884899e+01F, /*  0x1.96081cp+3 */
+sb2  =  4.51839523e+01F, /*  0x1.6978bcp+5 */
+sb3  =  4.72810211e+01F, /*  0x1.7a3f88p+5 */
+sb4  =  8.93033314e+00F; /*  0x1.1dc54ap+3 */
 
 float
 erff(float x)
@@ -106,43 +91,37 @@ erff(float x)
        }
 
        if(ix < 0x3f580000) {           /* |x|<0.84375 */
-           if(ix < 0x31800000) {       /* |x|<2**-28 */
-               if (ix < 0x04000000)
-                   /*avoid underflow */
-                   return (float)0.125*((float)8.0*x+efx8*x);
+           if(ix < 0x38800000) {       /* |x|<2**-14 */
+               if (ix < 0x04000000)    /* |x|<0x1p-119 */
+                   return (8*x+efx8*x)/8;      /* avoid spurious underflow */
                return x + efx*x;
            }
            z = x*x;
-           r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
-           s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
+           r = pp0+z*(pp1+z*pp2);
+           s = one+z*(qq1+z*(qq2+z*qq3));
            y = r/s;
            return x + x*y;
        }
        if(ix < 0x3fa00000) {           /* 0.84375 <= |x| < 1.25 */
            s = fabsf(x)-one;
-           P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
-           Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
+           P = pa0+s*(pa1+s*(pa2+s*pa3));
+           Q = one+s*(qa1+s*(qa2+s*(qa3+s*qa4)));
            if(hx>=0) return erx + P/Q; else return -erx - P/Q;
        }
-       if (ix >= 0x40c00000) {         /* inf>|x|>=6 */
+       if (ix >= 0x40800000) {         /* inf>|x|>=4 */
            if(hx>=0) return one-tiny; else return tiny-one;
        }
        x = fabsf(x);
        s = one/(x*x);
        if(ix< 0x4036DB6E) {    /* |x| < 1/0.35 */
-           R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
-                               ra5+s*(ra6+s*ra7))))));
-           S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
-                               sa5+s*(sa6+s*(sa7+s*sa8)))))));
+           R=ra0+s*(ra1+s*(ra2+s*ra3));
+           S=one+s*(sa1+s*(sa2+s*(sa3+s*sa4)));
        } else {        /* |x| >= 1/0.35 */
-           R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
-                               rb5+s*rb6)))));
-           S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
-                               sb5+s*(sb6+s*sb7))))));
+           R=rb0+s*(rb1+s*(rb2+s*(rb3+s*rb4)));
+           S=one+s*(sb1+s*(sb2+s*(sb3+s*sb4)));
        }
-       GET_FLOAT_WORD(ix,x);
-       SET_FLOAT_WORD(z,ix&0xfffff000);
-       r  =  __ieee754_expf(-z*z-(float)0.5625)*__ieee754_expf((z-x)*(z+x)+R/S);
+       SET_FLOAT_WORD(z,hx&0xffffe000);
+       r  = expf(-z*z-0.5625F)*expf((z-x)*(z+x)+R/S);
        if(hx>=0) return one-r/x; else return  r/x-one;
 }
 
@@ -159,11 +138,11 @@ erfcf(float x)
        }
 
        if(ix < 0x3f580000) {           /* |x|<0.84375 */
-           if(ix < 0x23800000)         /* |x|<2**-56 */
+           if(ix < 0x33800000)         /* |x|<2**-24 */
                return one-x;
            z = x*x;
-           r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
-           s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
+           r = pp0+z*(pp1+z*pp2);
+           s = one+z*(qq1+z*(qq2+z*qq3));
            y = r/s;
            if(hx < 0x3e800000) {       /* x<1/4 */
                return one-(x+x*y);
@@ -175,33 +154,27 @@ erfcf(float x)
        }
        if(ix < 0x3fa00000) {           /* 0.84375 <= |x| < 1.25 */
            s = fabsf(x)-one;
-           P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
-           Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
+           P = pa0+s*(pa1+s*(pa2+s*pa3));
+           Q = one+s*(qa1+s*(qa2+s*(qa3+s*qa4)));
            if(hx>=0) {
                z  = one-erx; return z - P/Q;
            } else {
                z = erx+P/Q; return one+z;
            }
        }
-       if (ix < 0x41e00000) {          /* |x|<28 */
+       if (ix < 0x41300000) {          /* |x|<11 */
            x = fabsf(x);
            s = one/(x*x);
            if(ix< 0x4036DB6D) {        /* |x| < 1/.35 ~ 2.857143*/
-               R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
-                               ra5+s*(ra6+s*ra7))))));
-               S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
-                               sa5+s*(sa6+s*(sa7+s*sa8)))))));
+               R=ra0+s*(ra1+s*(ra2+s*ra3));
+               S=one+s*(sa1+s*(sa2+s*(sa3+s*sa4)));
            } else {                    /* |x| >= 1/.35 ~ 2.857143 */
-               if(hx<0&&ix>=0x40c00000) return two-tiny;/* x < -6 */
-               R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
-                               rb5+s*rb6)))));
-               S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
-                               sb5+s*(sb6+s*sb7))))));
+               if(hx<0&&ix>=0x40a00000) return two-tiny;/* x < -5 */
+               R=rb0+s*(rb1+s*(rb2+s*(rb3+s*rb4)));
+               S=one+s*(sb1+s*(sb2+s*(sb3+s*sb4)));
            }
-           GET_FLOAT_WORD(ix,x);
-           SET_FLOAT_WORD(z,ix&0xfffff000);
-           r  =  __ieee754_expf(-z*z-(float)0.5625)*
-                       __ieee754_expf((z-x)*(z+x)+R/S);
+           SET_FLOAT_WORD(z,hx&0xffffe000);
+           r  = expf(-z*z-0.5625F)*expf((z-x)*(z+x)+R/S);
            if(hx>0) return r/x; else return two-r/x;
        } else {
            if(hx>0) return tiny*tiny; else return two-tiny;