Index: head/lib/msun/Makefile
===================================================================
--- head/lib/msun/Makefile	(revision 268592)
+++ head/lib/msun/Makefile	(revision 268593)
@@ -1,220 +1,220 @@
 #  @(#)Makefile 5.1beta 93/09/24
 # $FreeBSD$
 #
 #  ====================================================
 #  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.
 #  ====================================================
 #
 #
 
 .if ${MACHINE_CPUARCH} == "i386"
 ARCH_SUBDIR= i387
 .else
 ARCH_SUBDIR= ${MACHINE_CPUARCH}
 .endif
 
 .include "${ARCH_SUBDIR}/Makefile.inc"
 
 .PATH:	${.CURDIR}/${ARCH_SUBDIR}
 .if ${MACHINE_CPUARCH} == "i386" || ${MACHINE_CPUARCH} == "amd64"
 .PATH:	${.CURDIR}/x86
 CFLAGS+=	-I${.CURDIR}/x86
 .endif
 
 # long double format
 .if ${LDBL_PREC} == 64
 .PATH:  ${.CURDIR}/ld80
 CFLAGS+=	-I${.CURDIR}/ld80
 .elif ${LDBL_PREC} == 113
 .PATH:  ${.CURDIR}/ld128
 CFLAGS+=	-I${.CURDIR}/ld128
 .endif
 
 .PATH:	${.CURDIR}/bsdsrc
 .PATH:	${.CURDIR}/src
 .PATH:	${.CURDIR}/man
 
 LIB=	m
 SHLIBDIR?= /lib
 SHLIB_MAJOR= 5
 WARNS?=	1
 IGNORE_PRAGMA=
 COMMON_SRCS= b_exp.c b_log.c b_tgamma.c \
 	e_acos.c e_acosf.c e_acosh.c e_acoshf.c e_asin.c e_asinf.c \
 	e_atan2.c e_atan2f.c e_atanh.c e_atanhf.c e_cosh.c e_coshf.c e_exp.c \
 	e_expf.c e_fmod.c e_fmodf.c e_gamma.c e_gamma_r.c e_gammaf.c \
 	e_gammaf_r.c e_hypot.c e_hypotf.c e_j0.c e_j0f.c e_j1.c e_j1f.c \
 	e_jn.c e_jnf.c e_lgamma.c e_lgamma_r.c e_lgammaf.c e_lgammaf_r.c \
 	e_log.c e_log10.c e_log10f.c e_log2.c e_log2f.c e_logf.c \
 	e_pow.c e_powf.c e_rem_pio2.c \
 	e_rem_pio2f.c e_remainder.c e_remainderf.c e_scalb.c e_scalbf.c \
 	e_sinh.c e_sinhf.c e_sqrt.c e_sqrtf.c fenv.c \
 	imprecise.c \
 	k_cos.c k_cosf.c k_exp.c k_expf.c k_rem_pio2.c k_sin.c k_sinf.c \
 	k_tan.c k_tanf.c \
 	s_asinh.c s_asinhf.c s_atan.c s_atanf.c s_carg.c s_cargf.c s_cargl.c \
 	s_cbrt.c s_cbrtf.c s_ceil.c s_ceilf.c \
 	s_copysign.c s_copysignf.c s_cos.c s_cosf.c \
 	s_csqrt.c s_csqrtf.c s_erf.c s_erff.c \
 	s_exp2.c s_exp2f.c s_expm1.c s_expm1f.c s_fabsf.c s_fdim.c \
 	s_finite.c s_finitef.c \
 	s_floor.c s_floorf.c s_fma.c s_fmaf.c \
 	s_fmax.c s_fmaxf.c s_fmaxl.c s_fmin.c \
 	s_fminf.c s_fminl.c s_frexp.c s_frexpf.c s_ilogb.c s_ilogbf.c \
 	s_ilogbl.c s_isfinite.c s_isnan.c s_isnormal.c \
 	s_llrint.c s_llrintf.c s_llround.c s_llroundf.c s_llroundl.c \
 	s_log1p.c s_log1pf.c s_logb.c s_logbf.c s_lrint.c s_lrintf.c \
 	s_lround.c s_lroundf.c s_lroundl.c s_modff.c \
 	s_nan.c s_nearbyint.c s_nextafter.c s_nextafterf.c \
 	s_nexttowardf.c s_remquo.c s_remquof.c \
 	s_rint.c s_rintf.c s_round.c s_roundf.c \
 	s_scalbln.c s_scalbn.c s_scalbnf.c s_signbit.c \
 	s_signgam.c s_significand.c s_significandf.c s_sin.c s_sinf.c \
 	s_tan.c s_tanf.c s_tanh.c s_tanhf.c s_tgammaf.c s_trunc.c s_truncf.c \
 	w_cabs.c w_cabsf.c w_drem.c w_dremf.c
 
 # Location of fpmath.h and _fpmath.h
 LIBCDIR=	${.CURDIR}/../libc
 .if exists(${LIBCDIR}/${MACHINE_ARCH})
 LIBC_ARCH=${MACHINE_ARCH}
 .else
 LIBC_ARCH=${MACHINE_CPUARCH}
 .endif
 CFLAGS+=	-I${.CURDIR}/src -I${LIBCDIR}/include \
 	-I${LIBCDIR}/${LIBC_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_llrintl.c s_lrintl.c s_modfl.c
 .if ${LDBL_PREC} != 53
 # If long double != double use these; otherwise, we alias the double versions.
 COMMON_SRCS+=	e_acoshl.c e_acosl.c e_asinl.c e_atan2l.c e_atanhl.c \
 	e_coshl.c e_fmodl.c e_hypotl.c \
 	e_remainderl.c e_sinhl.c e_sqrtl.c \
 	invtrig.c k_cosl.c k_sinl.c k_tanl.c \
 	s_asinhl.c s_atanl.c s_cbrtl.c s_ceill.c s_cosl.c s_cprojl.c \
-	s_csqrtl.c s_exp2l.c s_expl.c s_floorl.c s_fmal.c \
+	s_csqrtl.c s_erfl.c s_exp2l.c s_expl.c s_floorl.c s_fmal.c \
 	s_frexpl.c s_logbl.c s_logl.c s_nanl.c s_nextafterl.c \
 	s_nexttoward.c s_remquol.c s_rintl.c s_roundl.c s_scalbnl.c \
 	s_sinl.c s_tanhl.c s_tanl.c s_truncl.c w_cabsl.c
 .endif
 
 # C99 complex functions
 COMMON_SRCS+=	catrig.c catrigf.c \
 	s_ccosh.c s_ccoshf.c s_cexp.c s_cexpf.c \
 	s_cimag.c s_cimagf.c s_cimagl.c \
 	s_conj.c s_conjf.c s_conjl.c \
 	s_cproj.c s_cprojf.c s_creal.c s_crealf.c s_creall.c \
 	s_csinh.c s_csinhf.c s_ctanh.c s_ctanhf.c
 
 # FreeBSD's C library supplies these functions:
 #COMMON_SRCS+=	s_fabs.c s_frexp.c s_isnan.c s_ldexp.c s_modf.c
 
 # Exclude the generic versions of what we provide in the MD area.
 .if defined(ARCH_SRCS)
 .for i in ${ARCH_SRCS}
 COMMON_SRCS:=  ${COMMON_SRCS:N${i:R}.c}
 .endfor
 .endif
 
 SRCS=	${COMMON_SRCS} ${ARCH_SRCS}
 
 INCS+=	fenv.h math.h
 
 MAN=	acos.3 acosh.3 asin.3 asinh.3 atan.3 atan2.3 atanh.3 \
 	ceil.3 cacos.3 ccos.3 ccosh.3 cexp.3 \
 	cimag.3 copysign.3 cos.3 cosh.3 csqrt.3 erf.3 exp.3 fabs.3 fdim.3 \
 	feclearexcept.3 feenableexcept.3 fegetenv.3 \
 	fegetround.3 fenv.3 floor.3 \
 	fma.3 fmax.3 fmod.3 hypot.3 ieee.3 ieee_test.3 ilogb.3 j0.3 \
 	lgamma.3 log.3 lrint.3 lround.3 math.3 nan.3 \
 	nextafter.3 remainder.3 rint.3 \
 	round.3 scalbn.3 signbit.3 sin.3 sinh.3 sqrt.3 tan.3 tanh.3 trunc.3 \
 	complex.3
 
 MLINKS+=acos.3 acosf.3 acos.3 acosl.3
 MLINKS+=acosh.3 acoshf.3 acosh.3 acoshl.3
 MLINKS+=asin.3 asinf.3 asin.3 asinl.3
 MLINKS+=asinh.3 asinhf.3 asinh.3 asinhl.3
 MLINKS+=atan.3 atanf.3 atan.3 atanl.3
 MLINKS+=atanh.3 atanhf.3 atanh.3 atanhl.3
 MLINKS+=atan2.3 atan2f.3 atan2.3 atan2l.3 \
 	atan2.3 carg.3 atan2.3 cargf.3 atan2.3 cargl.3
 MLINKS+=cacos.3 cacosf.3 cacos.3 cacosh.3 cacos.3 cacoshf.3 \
 	cacos.3 casin.3 cacos.3 casinf.3 cacos.3 casinh.3 cacos.3 casinhf.3 \
 	cacos.3 catan.3 cacos.3 catanf.3 cacos.3 catanh.3 cacos.3 catanhf.3
 MLINKS+=ccos.3 ccosf.3 ccos.3 csin.3 ccos.3 csinf.3 ccos.3 ctan.3 ccos.3 ctanf.3
 MLINKS+=ccosh.3 ccoshf.3 ccosh.3 csinh.3 ccosh.3 csinhf.3 \
 	ccosh.3 ctanh.3 ccosh.3 ctanhf.3
 MLINKS+=ceil.3 ceilf.3 ceil.3 ceill.3
 MLINKS+=cexp.3 cexpf.3
 MLINKS+=cimag.3 cimagf.3 cimag.3 cimagl.3 \
 	cimag.3 conj.3 cimag.3 conjf.3 cimag.3 conjl.3 \
 	cimag.3 cproj.3 cimag.3 cprojf.3 cimag.3 cprojl.3 \
 	cimag.3 creal.3 cimag.3 crealf.3 cimag.3 creall.3
 MLINKS+=copysign.3 copysignf.3 copysign.3 copysignl.3
 MLINKS+=cos.3 cosf.3 cos.3 cosl.3
 MLINKS+=cosh.3 coshf.3 cosh.3 coshl.3
 MLINKS+=csqrt.3 csqrtf.3 csqrt.3 csqrtl.3
-MLINKS+=erf.3 erfc.3 erf.3 erff.3 erf.3 erfcf.3
+MLINKS+=erf.3 erfc.3 erf.3 erff.3 erf.3 erfcf.3 erf.3 erfl.3 erf.3 erfcl.3
 MLINKS+=exp.3 expm1.3 exp.3 expm1f.3 exp.3 expm1l.3 exp.3 pow.3 exp.3 powf.3 \
 	exp.3 exp2.3 exp.3 exp2f.3 exp.3 exp2l.3 exp.3 expf.3 exp.3 expl.3
 MLINKS+=fabs.3 fabsf.3 fabs.3 fabsl.3
 MLINKS+=fdim.3 fdimf.3 fdim.3 fdiml.3
 MLINKS+=feclearexcept.3 fegetexceptflag.3 feclearexcept.3 feraiseexcept.3 \
 	feclearexcept.3 fesetexceptflag.3 feclearexcept.3 fetestexcept.3
 MLINKS+=feenableexcept.3 fedisableexcept.3 feenableexcept.3 fegetexcept.3
 MLINKS+=fegetenv.3 feholdexcept.3 fegetenv.3 fesetenv.3 \
 	fegetenv.3 feupdateenv.3
 MLINKS+=fegetround.3 fesetround.3
 MLINKS+=floor.3 floorf.3 floor.3 floorl.3
 MLINKS+=fma.3 fmaf.3 fma.3 fmal.3
 MLINKS+=fmax.3 fmaxf.3 fmax.3 fmaxl.3 \
 	fmax.3 fmin.3 fmax.3 fminf.3 fmax.3 fminl.3
 MLINKS+=fmod.3 fmodf.3 fmod.3 fmodl.3
 MLINKS+=hypot.3 cabs.3 hypot.3 cabsf.3 hypot.3 cabsl.3 \
 	hypot.3 hypotf.3 hypot.3 hypotl.3
 MLINKS+=ieee_test.3 scalb.3 ieee_test.3 scalbf.3
 MLINKS+=ieee_test.3 significand.3 ieee_test.3 significandf.3
 MLINKS+=ilogb.3 ilogbf.3 ilogb.3 ilogbl.3 \
 	ilogb.3 logb.3 ilogb.3 logbf.3 ilogb.3 logbl.3
 MLINKS+=j0.3 j1.3 j0.3 jn.3 j0.3 y0.3 j0.3 y1.3 j0.3 y1f.3 j0.3 yn.3
 MLINKS+=j0.3 j0f.3 j0.3 j1f.3 j0.3 jnf.3 j0.3 y0f.3 j0.3 ynf.3
 MLINKS+=lgamma.3 gamma.3 lgamma.3 gammaf.3 lgamma.3 lgammaf.3 \
 	lgamma.3 tgamma.3 lgamma.3 tgammaf.3
 MLINKS+=log.3 log10.3 log.3 log10f.3 log.3 log10l.3 \
 	log.3 log1p.3 log.3 log1pf.3 log.3 log1pl.3 \
 	log.3 logf.3 log.3 logl.3 \
 	log.3 log2.3 log.3 log2f.3 log.3 log2l.3
 MLINKS+=lrint.3 llrint.3 lrint.3 llrintf.3 lrint.3 llrintl.3 \
 	lrint.3 lrintf.3 lrint.3 lrintl.3
 MLINKS+=lround.3 llround.3 lround.3 llroundf.3 lround.3 llroundl.3 \
 	lround.3 lroundf.3 lround.3 lroundl.3
 MLINKS+=nan.3 nanf.3 nan.3 nanl.3
 MLINKS+=nextafter.3 nextafterf.3 nextafter.3 nextafterl.3
 MLINKS+=nextafter.3 nexttoward.3 nextafter.3 nexttowardf.3
 MLINKS+=nextafter.3 nexttowardl.3
 MLINKS+=remainder.3 remainderf.3 remainder.3 remainderl.3 \
 	remainder.3 remquo.3 remainder.3 remquof.3 remainder.3 remquol.3
 MLINKS+=rint.3 rintf.3 rint.3 rintl.3 \
 	rint.3 nearbyint.3 rint.3 nearbyintf.3 rint.3 nearbyintl.3
 MLINKS+=round.3 roundf.3 round.3 roundl.3
 MLINKS+=scalbn.3 scalbln.3 scalbn.3 scalblnf.3 scalbn.3 scalblnl.3
 MLINKS+=scalbn.3 scalbnf.3 scalbn.3 scalbnl.3
 MLINKS+=sin.3 sinf.3 sin.3 sinl.3
 MLINKS+=sinh.3 sinhf.3 sinh.3 sinhl.3
 MLINKS+=sqrt.3 cbrt.3 sqrt.3 cbrtf.3 sqrt.3 cbrtl.3 sqrt.3 sqrtf.3 \
 	sqrt.3 sqrtl.3
 MLINKS+=tan.3 tanf.3 tan.3 tanl.3
 MLINKS+=tanh.3 tanhf.3 tanh.3 tanhl.3
 MLINKS+=trunc.3 truncf.3 trunc.3 truncl.3
 
 .include <bsd.lib.mk>
Index: head/lib/msun/Symbol.map
===================================================================
--- head/lib/msun/Symbol.map	(revision 268592)
+++ head/lib/msun/Symbol.map	(revision 268593)
@@ -1,282 +1,282 @@
 /*
  * $FreeBSD$
  */
 
 /* 7.0-CURRENT */
 FBSD_1.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;
 };
 
 /* First added in 8.0-CURRENT */
 FBSD_1.1 {
 	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;
 };
 
 /* First added in 9.0-CURRENT */
 FBSD_1.2 {
 	__isnanf;
 	cbrtl;
 	cexp;
 	cexpf;
 	log2;
 	log2f;
 };
 
 /* First added in 10.0-CURRENT */
 FBSD_1.3 {
 	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;
 	coshl;
 	ctan;
 	ctanf;
 	ctanh;
 	ctanhf;
+	erfcl;
+	erfl;
 	expl;
 	expm1l;
 	log10l;
 	log1pl;
 	log2l;
 	logl;
 	sinhl;
 	tanhl;
 	/* Implemented as weak aliases for imprecise versions */
-	erfcl;
-	erfl;
 	lgammal;
 	powl;
 	tgammal;
 };
Index: head/lib/msun/ld128/s_erfl.c
===================================================================
--- head/lib/msun/ld128/s_erfl.c	(nonexistent)
+++ head/lib/msun/ld128/s_erfl.c	(revision 268593)
@@ -0,0 +1,329 @@
+/* @(#)s_erf.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/*
+ * See s_erf.c for complete comments.
+ *
+ * Converted to long double by Steven G. Kargl.
+ */
+#include <float.h>
+
+#include "fpmath.h"
+#include "math.h"
+#include "math_private.h"
+
+/* XXX Prevent compilers from erroneously constant folding these: */
+static const volatile long double tiny = 0x1p-10000L;
+
+static const double
+half= 0.5,
+one = 1,
+two = 2;
+/*
+ * In the domain [0, 2**-40], only the first term in the power series
+ * expansion of erf(x) is used.  The magnitude of the first neglected
+ * terms is less than 2**-120.
+ */
+static const long double
+efx  =  1.28379167095512573896158903121545167e-01L,	/* 0xecbff6a7, 0x481dd788, 0xb64d21a8, 0xeb06fc3f */
+efx8 =  1.02703333676410059116927122497236133e+00L,	/* 0xecbff6a7, 0x481dd788, 0xb64d21a8, 0xeb06ff3f */
+/*
+ * Domain [0, 0.84375], range ~[-1.919e-38, 1.919e-38]:
+ * |(erf(x) - x)/x - pp(x)/qq(x)| < 2**-125.29
+ */
+pp0  =  1.28379167095512573896158903121545167e-01L,	/* 0x3ffc06eb, 0xa8214db6, 0x88d71d48, 0xa7f6bfec */
+pp1  = -3.14931554396568573802046931159683404e-01L,	/* 0xbffd427d, 0x6ada7263, 0x547eb096, 0x95f37463 */
+pp2  = -5.27514920282183487103576956956725309e-02L,	/* 0xbffab023, 0xe5a271e3, 0xb0e79b01, 0x2f7ac962 */
+pp3  = -1.13202828509005281355609495523452713e-02L,	/* 0xbff872f1, 0x6a5023a1, 0xe08b3884, 0x326af20f */
+pp4  = -9.18626155872522453865998391206048506e-04L,	/* 0xbff4e19f, 0xea5fb024, 0x43247a37, 0xe430b06c */
+pp5  = -7.87518862406176274922506447157284230e-05L,	/* 0xbff14a4f, 0x31a85fe0, 0x7fff2204, 0x09c49b37 */
+pp6  = -3.42357944472240436548115331090560881e-06L,	/* 0xbfeccb81, 0x4b43c336, 0xcd2eb6c2, 0x903f2d87 */
+pp7  = -1.37317432573890412634717890726745428e-07L,	/* 0xbfe826e3, 0x0e915eb6, 0x42aee414, 0xf7e36805 */
+pp8  = -2.71115170113861755855049008732113726e-09L,	/* 0xbfe2749e, 0x2b94fd00, 0xecb4d166, 0x0efb91f8 */
+pp9  = -3.37925756196555959454018189718117864e-11L,	/* 0xbfdc293e, 0x1d9060cb, 0xd043204a, 0x314cd7f0 */
+qq1  =  4.76672625471551170489978555182449450e-01L,	/* 0x3ffde81c, 0xde6531f0, 0x76803bee, 0x526e29e9 */
+qq2  =  1.06713144672281502058807525850732240e-01L,	/* 0x3ffbb518, 0xd7a6bb74, 0xcd9bdd33, 0x7601eee5 */
+qq3  =  1.47747613127513761102189201923147490e-02L,	/* 0x3ff8e423, 0xae527e18, 0xf12cb447, 0x723b4749 */
+qq4  =  1.39939377672028671891148770908874816e-03L,	/* 0x3ff56ed7, 0xba055d84, 0xc21b45c4, 0x388d1812 */
+qq5  =  9.44302939359455241271983309378738276e-05L,	/* 0x3ff18c11, 0xc18c99a4, 0x86d0fe09, 0x46387b4c */
+qq6  =  4.56199342312522842161301671745365650e-06L,	/* 0x3fed3226, 0x73421d05, 0x08875300, 0x32fa1432 */
+qq7  =  1.53019260483764773845294600092361197e-07L,	/* 0x3fe8489b, 0x3a63f627, 0x2b9ad2ce, 0x26516e57 */
+qq8  =  3.25542691121324805094777901250005508e-09L,	/* 0x3fe2bf6c, 0x26d93a29, 0x9142be7c, 0x9f1dd043 */
+qq9  =  3.37405581964478060434410167262684979e-11L;	/* 0x3fdc28c8, 0xfb8fa1be, 0x10e57eec, 0xaa19e49f */
+
+static const long double
+erx  =  8.42700792949714894142232424201210961e-01L,	/* 0x3ffeaf76, 0x7a741088, 0xb0000000, 0x00000000 */
+/*
+ * Domain [0.84375, 1.25], range ~[-2.521e-36, 2.523e-36]:
+ * |(erf(x) - erx) - pa(x)/qa(x)| < 2**-120.15
+ */
+pa0  = -2.48010117891186017024438233323795897e-17L,	/* 0xbfc7c97f, 0x77812279, 0x6c877f22, 0xef4bfb2e */
+pa1  =  4.15107497420594680894327969504526489e-01L,	/* 0x3ffda911, 0xf096fbc2, 0x55662005, 0x2337fa64 */
+pa2  = -3.94180628087084846724448515851892609e-02L,	/* 0xbffa42e9, 0xab54528c, 0xad529da1, 0x6efc2af3 */
+pa3  =  4.48897599625192107295954790681677462e-02L,	/* 0x3ffa6fbc, 0xa65edba1, 0x0e4cbcea, 0x73ef9a31 */
+pa4  =  8.02069252143016600110972019232995528e-02L,	/* 0x3ffb4887, 0x0e8b548e, 0x3230b417, 0x11b553b3 */
+pa5  = -1.02729816533435279443621120242391295e-02L,	/* 0xbff850a0, 0x041de3ee, 0xd5bca6c9, 0x4ef5f9f2 */
+pa6  =  5.70777694530755634864821094419982095e-03L,	/* 0x3ff77610, 0x9b501e10, 0x4c978382, 0x742df68f */
+pa7  =  1.22635150233075521018231779267077071e-03L,	/* 0x3ff5417b, 0x0e623682, 0x60327da0, 0x96b9219e */
+pa8  =  5.36100234820204569428412542856666503e-04L,	/* 0x3ff41912, 0x27ceb4c1, 0x1d3298ec, 0x84ced627 */
+pa9  = -1.97753571846365167177187858667583165e-04L,	/* 0xbff29eb8, 0x23f5bcf3, 0x15c83c46, 0xe4fda98b */
+pa10 =  6.19333039900846970674794789568415105e-05L,	/* 0x3ff103c4, 0x60f88e46, 0xc0c9fb02, 0x13cc7fc1 */
+pa11 = -5.40531400436645861492290270311751349e-06L,	/* 0xbfed6abe, 0x9665f8a8, 0xdd0ad3ba, 0xe5dc0ee3 */
+qa1  =  9.05041313265490487793231810291907851e-01L,	/* 0x3ffecf61, 0x93340222, 0xe9930620, 0xc4e61168 */
+qa2  =  6.79848064708886864767240880834868092e-01L,	/* 0x3ffe5c15, 0x0ba858dc, 0xf7900ae9, 0xfea1e09a */
+qa3  =  4.04720609926471677581066689316516445e-01L,	/* 0x3ffd9e6f, 0x145e9b00, 0x6d8c1749, 0xd2928623 */
+qa4  =  1.69183273898369996364661075664302225e-01L,	/* 0x3ffc5a7c, 0xc2a363c1, 0xd6c19097, 0xef9b4063 */
+qa5  =  7.44476185988067992342479750486764248e-02L,	/* 0x3ffb30ef, 0xfc7259ef, 0x1bcbb089, 0x686dd62d */
+qa6  =  2.02981172725892407200420389604788573e-02L,	/* 0x3ff94c90, 0x7976cb0e, 0x21e1d36b, 0x0f09ca2b */
+qa7  =  6.94281866271607668268269403102277234e-03L,	/* 0x3ff7c701, 0x2b193250, 0xc5d46ecc, 0x374843d8 */
+qa8  =  1.12952275469171559611651594706820034e-03L,	/* 0x3ff52818, 0xfd2a7c06, 0xd13e38fd, 0xda4b34f5 */
+qa9  =  3.13736683241992737197226578597710179e-04L,	/* 0x3ff348fa, 0x0cb48d18, 0x051f849b, 0x135ccf74 */
+qa10 =  1.17037675204033225470121134087771410e-05L,	/* 0x3fee88b6, 0x98f47704, 0xa5d8f8f2, 0xc6422e11 */
+qa11 =  4.61312518293853991439362806880973592e-06L,	/* 0x3fed3594, 0xe31db94f, 0x3592b693, 0xed4386b4 */
+qa12 = -1.02158572037456893687737553657431771e-06L;	/* 0xbfeb123a, 0xd60d9b1e, 0x1f6fdeb9, 0x7dc8410a */
+/*
+ * Domain [1.25,2.85715], range ~[-2.922e-37,2.922e-37]:
+ * |log(x*erfc(x)) + x**2 + 0.5625 - ra(x)/sa(x)| < 2**-121.36
+ */
+static const long double
+ra0  = -9.86494292470069009555706994426014461e-03L,	/* 0xbff84341, 0x239e8709, 0xe941b06a, 0xcb4b6ec5 */
+ra1  = -1.13580436992565640457579040117568870e+00L,	/* 0xbfff22c4, 0x133f7c0d, 0x72d5e231, 0x2eb1ee3f */
+ra2  = -4.89744330295291950661185707066921755e+01L,	/* 0xc00487cb, 0xa38b4fc2, 0xc136695b, 0xc1df8047 */
+ra3  = -1.10766149300215937173768072715352140e+03L,	/* 0xc00914ea, 0x55e6beb3, 0xabc50e07, 0xb6e5664d */
+ra4  = -1.49991031232170934967642795601952100e+04L,	/* 0xc00cd4b8, 0xd33243e6, 0xffbf6545, 0x3c57ef6e */
+ra5  = -1.29805749738318462882524181556996692e+05L,	/* 0xc00ffb0d, 0xbfeed9b6, 0x5b2a3ff4, 0xe245bd3c */
+ra6  = -7.42828497044940065828871976644647850e+05L,	/* 0xc0126ab5, 0x8fe7caca, 0x473352d9, 0xcd4e0c90 */
+ra7  = -2.85637299581890734287995171242421106e+06L,	/* 0xc0145cad, 0xa7f76fe7, 0x3e358051, 0x1799f927 */
+ra8  = -7.40674797129824999383748865571026084e+06L,	/* 0xc015c412, 0x6fe29c02, 0x298ad158, 0x7d24e45c */
+ra9  = -1.28653420911930973914078724204151759e+07L,	/* 0xc016889e, 0x7c2eb0dc, 0x95d5863b, 0x0aa34dc3 */
+ra10 = -1.47198163599330179552932489109452638e+07L,	/* 0xc016c136, 0x90b84923, 0xf9bcb497, 0x19bbd0f5 */
+ra11 = -1.07812992258382800318665248311522624e+07L,	/* 0xc0164904, 0xe673a113, 0x35d7f079, 0xe13701f3 */
+ra12 = -4.83545565681708642630419905537756076e+06L,	/* 0xc0152721, 0xfea094a8, 0x869eb39d, 0x413d6f13 */
+ra13 = -1.23956521201673964822976917356685286e+06L,	/* 0xc0132ea0, 0xd3646baa, 0x2fe62b0d, 0xbae5ce85 */
+ra14 = -1.62289333553652417591275333240371812e+05L,	/* 0xc0103cf8, 0xaab1e2d6, 0x4c25e014, 0x248d76ab */
+ra15 = -8.82890392601176969729168894389833110e+03L,	/* 0xc00c13e7, 0x3b3d8f94, 0x6fbda6f6, 0xe7049a82 */
+ra16 = -1.22591866337261720023681535568334619e+02L,	/* 0xc005ea5e, 0x12358891, 0xcfa712c5, 0x77f050d4 */
+sa1  =  6.44508918884710829371852723353794047e+01L,	/* 0x400501cd, 0xb69a6c0f, 0x5716de14, 0x47161af6 */
+sa2  =  1.76118475473171481523704824327358534e+03L,	/* 0x4009b84b, 0xd305829f, 0xc4c771b0, 0xbf1f7f9b */
+sa3  =  2.69448346969488374857087646131950188e+04L,	/* 0x400da503, 0x56bacc05, 0x4fdba68d, 0x2cca27e6 */
+sa4  =  2.56826633369941456778326497384543763e+05L,	/* 0x4010f59d, 0x51124428, 0x69c41de6, 0xbd0d5753 */
+sa5  =  1.60647413092257206847700054645905859e+06L,	/* 0x40138834, 0xa2184244, 0x557a1bed, 0x68c9d556 */
+sa6  =  6.76963075165099718574753447122393797e+06L,	/* 0x40159d2f, 0x7b01b0cc, 0x8bac9e95, 0x5d35d56e */
+sa7  =  1.94295690905361884290986932493647741e+07L,	/* 0x40172878, 0xc1172d61, 0x3068501e, 0x2f3c71da */
+sa8  =  3.79774781017759149060839255547073541e+07L,	/* 0x401821be, 0xc30d06fe, 0x410563d7, 0x032111fd */
+sa9  =  5.00659831846029484248302236457727397e+07L,	/* 0x40187df9, 0x1f97a111, 0xc51d6ac2, 0x4b389793 */
+sa10 =  4.36486287620506484276130525941972541e+07L,	/* 0x40184d03, 0x3a618ae0, 0x2a723357, 0xfa45c60a */
+sa11 =  2.43779678791333894255510508253951934e+07L,	/* 0x401773fa, 0x6fe10ee2, 0xc467850d, 0xc6b7ff30 */
+sa12 =  8.30732360384443202039372372212966542e+06L,	/* 0x4015fb09, 0xee6a5631, 0xdd98de7e, 0x8b00461a */
+sa13 =  1.60160846942050515734192397495105693e+06L,	/* 0x40138704, 0x8782bf13, 0x5b8fb315, 0xa898abe5 */
+sa14 =  1.54255505242533291014555153757001825e+05L,	/* 0x40102d47, 0xc0abc98e, 0x843c9490, 0xb4352440 */
+sa15 =  5.87949220002375547561467275493888824e+03L,	/* 0x400b6f77, 0xe00d21d1, 0xec4d41e8, 0x2f8e1673 */
+sa16 =  4.97272976346793193860385983372237710e+01L;	/* 0x40048dd1, 0x816c1b3f, 0x24f540a6, 0x4cfe03cc */
+/*
+ * Domain [2.85715,9], range ~[-7.886e-37,7.918e-37]:
+ * |log(x*erfc(x)) + x**2 + 0.5625 - rb(x)/sb(x)| < 2**-120
+ */
+static const long double
+rb0  = -9.86494292470008707171371994479162369e-3L, /* 0xbff84341, 0x239e86f4, 0x2f57e561, 0xf4469360 */
+rb1  = -1.57047326624110727986326503729442830L,    /* 0xbfff920a, 0x8935bf73, 0x8803b894, 0x4656482d */
+rb2  = -1.03228196364885474342132255440317065e2L,  /* 0xc0059ce9, 0xac4ed0ff, 0x2cff0ff7, 0x5e70d1ab */
+rb3  = -3.74000570653418227179358710865224376e3L,  /* 0xc00ad380, 0x2ebf7835, 0xf6b07ed2, 0x861242f7 */
+rb4  = -8.35435477739098044190860390632813956e4L,  /* 0xc00f4657, 0x8c3ae934, 0x3647d7b3, 0x80e76fb7 */
+rb5  = -1.21398672055223642118716640216747152e6L,  /* 0xc0132862, 0x2b8761c8, 0x27d18c0f, 0x137c9463 */
+rb6  = -1.17669175877248796101665344873273970e7L,  /* 0xc0166719, 0x0b2cea46, 0x81f14174, 0x11602ea5 */
+rb7  = -7.66108006086998253606773064264599615e7L,  /* 0xc019243f, 0x3c26f4f0, 0x1cc05241, 0x3b953728 */
+rb8  = -3.32547117558141845968704725353130804e8L,  /* 0xc01b3d24, 0x42d8ee26, 0x24ef6f3b, 0x604a8c65 */
+rb9  = -9.41561252426350696802167711221739746e8L,  /* 0xc01cc0f8, 0xad23692a, 0x8ddb2310, 0xe9937145 */
+rb10 = -1.67157110805390944549427329626281063e9L,  /* 0xc01d8e88, 0x9a903734, 0x09a55fa3, 0xd205c903 */
+rb11 = -1.74339631004410841337645931421427373e9L,  /* 0xc01d9fa8, 0x77582d2a, 0xc183b8ab, 0x7e00cb05 */
+rb12 = -9.57655233596934915727573141357471703e8L,  /* 0xc01cc8a5, 0x460cc685, 0xd0271fa0, 0x6a70e3da */
+rb13 = -2.26320062731339353035254704082495066e8L,  /* 0xc01aafab, 0xd7d76721, 0xc9720e11, 0x6a8bd489 */
+rb14 = -1.42777302996263256686002973851837039e7L,  /* 0xc016b3b8, 0xc499689f, 0x2b88d965, 0xc32414f9 */
+sb1  =  1.08512869705594540211033733976348506e2L,  /* 0x4005b20d, 0x2db7528d, 0x00d20dcb, 0x858f6191 */
+sb2  =  5.02757713761390460534494530537572834e3L,  /* 0x400b3a39, 0x3bf4a690, 0x3025d28d, 0xfd40a891 */
+sb3  =  1.31019107205412870059331647078328430e5L,  /* 0x400fffcb, 0x1b71d05e, 0x3b28361d, 0x2a3c3690 */
+sb4  =  2.13021555152296846166736757455018030e6L,  /* 0x40140409, 0x3c6984df, 0xc4491d7c, 0xb04aa08d */
+sb5  =  2.26649105281820861953868568619768286e7L,  /* 0x401759d6, 0xce8736f0, 0xf28ad037, 0x2a901e0c */
+sb6  =  1.61071939490875921812318684143076081e8L,  /* 0x401a3338, 0x686fb541, 0x6bd27d06, 0x4f95c9ac */
+sb7  =  7.66895673844301852676056750497991966e8L,  /* 0x401c6daf, 0x31cec121, 0x54699126, 0x4bd9bf9e */
+sb8  =  2.41884450436101936436023058196042526e9L,  /* 0x401e2059, 0x46b0b8d7, 0x87b64cbf, 0x78bc296d */
+sb9  =  4.92403055884071695093305291535107666e9L,  /* 0x401f257e, 0xbe5ed739, 0x39e17346, 0xcadd2e55 */
+sb10 =  6.18627786365587486459633615573786416e9L,  /* 0x401f70bb, 0x1be7a7e7, 0x6a45b5ae, 0x607c70f0 */
+sb11 =  4.45898013426501378097430226324743199e9L,  /* 0x401f09c6, 0xa32643d7, 0xf1724620, 0x9ea46c32 */
+sb12 =  1.63006115763329848117160344854224975e9L,  /* 0x401d84a3, 0x0996887f, 0x65a4f43b, 0x978c1d74 */
+sb13 =  2.39216717012421697446304015847567721e8L,  /* 0x401ac845, 0x09a065c2, 0x30095da7, 0x9d72d6ae */
+sb14 =  7.84837329009278694937250358810225609e6L;  /* 0x4015df06, 0xd5290e15, 0x63031fac, 0x4d9c894c */
+/*
+ * Domain [9,108], range ~[-5.324e-38,5.340e-38]:
+ * |log(x*erfc(x)) + x**2 + 0.5625 - r(x)/s(x)| < 2**-124
+ */
+static const long double
+rc0  = -9.86494292470008707171367567652935673e-3L, /* 0xbff84341, 0x239e86f4, 0x2f57e55b, 0x1aa10fd3 */
+rc1  = -1.26229447747315096406518846411562266L,    /* 0xbfff4325, 0xbb1aab28, 0xda395cd9, 0xfb861c15 */
+rc2  = -6.13742634438922591780742637728666162e1L,  /* 0xc004eafe, 0x7dd51cd8, 0x3c7c5928, 0x751e50cf */
+rc3  = -1.50455835478908280402912854338421517e3L,  /* 0xc0097823, 0xbc15b9ab, 0x3d60745c, 0x523e80a5 */
+rc4  = -2.04415631865861549920184039902945685e4L,  /* 0xc00d3f66, 0x40b3fc04, 0x5388f2ec, 0xb009e1f0 */
+rc5  = -1.57625662981714582753490610560037638e5L,  /* 0xc01033dc, 0xd4dc95b6, 0xfd4da93b, 0xf355b4a9 */
+rc6  = -6.73473451616752528402917538033283794e5L,  /* 0xc01248d8, 0x2e73a4f9, 0xcded49c5, 0xfa3bfeb7 */
+rc7  = -1.47433165421387483167186683764364857e6L,  /* 0xc01367f1, 0xba77a8f7, 0xcfdd0dbb, 0x25d554b3 */
+rc8  = -1.38811981807868828563794929997744139e6L,  /* 0xc01352e5, 0x7d16d9ad, 0xbbdcbf38, 0x38fbc5ea */
+rc9  = -3.59659700530831825640766479698155060e5L,  /* 0xc0115f3a, 0xecd57f45, 0x21f8ad6c, 0x910a5958 */
+sc1  =  7.72730753022908298637508998072635696e1L,  /* 0x40053517, 0xa10d52bc, 0xdabb55b6, 0xbd0328cd */
+sc2  =  2.36825757341694050500333261769082182e3L,  /* 0x400a2808, 0x3e0a9b42, 0x82977842, 0x9c5de29e */
+sc3  =  3.72210540173034735352888847134073099e4L,  /* 0x400e22ca, 0x1ba827ef, 0xac8390d7, 0x1fc39a41 */
+sc4  =  3.24136032646418336712461033591393412e5L,  /* 0x40113c8a, 0x0216e100, 0xc59d1e44, 0xf0e68d9d */
+sc5  =  1.57836135851134393802505823370009175e6L,  /* 0x40138157, 0x95bc7664, 0x17575961, 0xdbe58eeb */
+sc6  =  4.12881981392063738026679089714182355e6L,  /* 0x4014f801, 0x9e82e8d2, 0xb8b3a70e, 0xfd84185d */
+sc7  =  5.24438427289213488410596395361544142e6L,  /* 0x40154017, 0x81177109, 0x2aa6c3b0, 0x1f106625 */
+sc8  =  2.59909544563616121735963429710382149e6L,  /* 0x40143d45, 0xbb90a9b1, 0x12bf9390, 0xa827a700 */
+sc9  =  2.80930665169282501639651995082335693e5L;  /* 0x40111258, 0xaa92222e, 0xa97e3216, 0xa237fa6c */
+
+long double
+erfl(long double x)
+{
+	long double ax,R,S,P,Q,s,y,z,r;
+	uint64_t lx, llx;
+	int32_t i;
+	uint16_t hx;
+
+	EXTRACT_LDBL128_WORDS(hx, lx, llx, x);
+
+	if((hx & 0x7fff) == 0x7fff) {	/* erfl(nan)=nan */
+		i = (hx>>15)<<1;
+		return (1-i)+one/x;	/* erfl(+-inf)=+-1 */
+	}
+
+	ax = fabsl(x);
+	if(ax < 0.84375) {
+	    if(ax < 0x1p-40L) {
+	        if(ax < 0x1p-16373L)	
+		    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+z*(pp5+z*(pp6+z*(pp7+
+		z*(pp8+z*pp9))))))));
+	    s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*(qq5+z*(qq6+z*(qq7+
+		z*(qq8+z*qq9))))))));
+	    y = r/s;
+	    return x + x*y;
+	}
+	if(ax < 1.25) {
+	    s = ax-one;
+	    P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*(pa6+s*(pa7+
+		s*(pa8+s*(pa9+s*(pa10+s*pa11))))))))));
+	    Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*(qa6+s*(qa7+
+		s*(qa8+s*(qa9+s*(qa10+s*(qa11+s*qa12)))))))))));
+	    if(x>=0) return (erx + P/Q); else return (-erx - P/Q);
+	}
+	if (ax >= 9) {			/* inf>|x|>= 9 */
+	    if(x>=0) return (one-tiny); else return (tiny-one);
+	}
+	s = one/(ax*ax);
+	if(ax < 2.85715) {	/* |x| < 2.85715 */
+	    R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(ra5+s*(ra6+s*(ra7+
+		s*(ra8+s*(ra9+s*(ra10+s*(ra11+s*(ra12+s*(ra13+s*(ra14+
+		s*(ra15+s*ra16)))))))))))))));
+	    S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(sa5+s*(sa6+s*(sa7+
+		s*(sa8+s*(sa9+s*(sa10+s*(sa11+s*(sa12+s*(sa13+s*(sa14+
+		s*(sa15+s*sa16)))))))))))))));
+	} else {	/* |x| >= 2.85715 */
+	    R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(rb5+s*(rb6+s*(rb7+
+		s*(rb8+s*(rb9+s*(rb10+s*(rb11+s*(rb12+s*(rb13+
+		s*rb14)))))))))))));
+	    S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(sb5+s*(sb6+s*(sb7+
+		s*(sb8+s*(sb9+s*(sb10+s*(sb11+s*(sb12+s*(sb13+
+		s*sb14)))))))))))));
+	}
+	z = (float)ax;
+	r = expl(-z*z-0.5625)*expl((z-ax)*(z+ax)+R/S);
+	if(x>=0) return (one-r/ax); else return (r/ax-one);
+}
+
+long double
+erfcl(long double x)
+{
+	long double ax,R,S,P,Q,s,y,z,r;
+	uint64_t lx, llx;
+	uint16_t hx;
+
+	EXTRACT_LDBL128_WORDS(hx, lx, llx, x);
+
+	if((hx & 0x7fff) == 0x7fff) {	/* erfcl(nan)=nan */
+					/* erfcl(+-inf)=0,2 */
+	    return ((hx>>15)<<1)+one/x;
+	}
+
+	ax = fabsl(x);
+	if(ax < 0.84375L) {
+	    if(ax < 0x1p-34L)
+		return one-x;
+	    z = x*x;
+	    r = pp0+z*(pp1+z*(pp2+z*(pp3+z*(pp4+z*(pp5+z*(pp6+z*(pp7+
+		z*(pp8+z*pp9))))))));
+	    s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*(qq5+z*(qq6+z*(qq7+
+		z*(qq8+z*qq9))))))));
+	    y = r/s;
+	    if(ax < 0.25L) {  	/* x<1/4 */
+		return one-(x+x*y);
+	    } else {
+		r = x*y;
+		r += (x-half);
+	       return half - r;
+	    }
+	}
+	if(ax < 1.25L) {
+	    s = ax-one;
+	    P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*(pa6+s*(pa7+
+		    s*(pa8+s*(pa9+s*(pa10+s*pa11))))))))));
+	    Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*(qa6+s*(qa7+
+		    s*(qa8+s*(qa9+s*(qa10+s*(qa11+s*qa12)))))))))));
+	    if(x>=0) {
+	        z  = one-erx; return z - P/Q;
+	    } else {
+		z = erx+P/Q; return one+z;
+	    }
+	}
+
+	if(ax < 108) {			/* |x| < 108 */
+ 	    s = one/(ax*ax);
+	    if(ax < 2.85715) {		/* |x| < 2.85715 */
+	        R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(ra5+s*(ra6+s*(ra7+
+		    s*(ra8+s*(ra9+s*(ra10+s*(ra11+s*(ra12+s*(ra13+s*(ra14+
+		    s*(ra15+s*ra16)))))))))))))));
+	        S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(sa5+s*(sa6+s*(sa7+
+		    s*(sa8+s*(sa9+s*(sa10+s*(sa11+s*(sa12+s*(sa13+s*(sa14+
+		    s*(sa15+s*sa16)))))))))))))));
+	    } else if(ax < 9) {
+		R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(rb5+s*(rb6+s*(rb7+
+		    s*(rb8+s*(rb9+s*(rb10+s*(rb11+s*(rb12+s*(rb13+
+		    s*rb14)))))))))))));
+		S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(sb5+s*(sb6+s*(sb7+
+		    s*(sb8+s*(sb9+s*(sb10+s*(sb11+s*(sb12+s*(sb13+
+		    s*sb14)))))))))))));
+	    } else {
+		if(x < -9) return two-tiny;	/* x < -9 */
+		R=rc0+s*(rc1+s*(rc2+s*(rc3+s*(rc4+s*(rc5+s*(rc6+s*(rc7+
+		    s*(rc8+s*rc9))))))));
+		S=one+s*(sc1+s*(sc2+s*(sc3+s*(sc4+s*(sc5+s*(sc6+s*(sc7+
+		    s*(sc8+s*sc9))))))));
+	    }
+	    z = (float)ax;
+	    r = expl(-z*z-0.5625)*expl((z-ax)*(z+ax)+R/S);
+	    if(x>0) return r/ax; else return two-r/ax;
+	} else {
+	    if(x>0) return tiny*tiny; else return two-tiny;
+	}
+}

Property changes on: head/lib/msun/ld128/s_erfl.c
___________________________________________________________________
Added: svn:eol-style
## -0,0 +1 ##
+native
\ No newline at end of property
Added: svn:keywords
## -0,0 +1 ##
+FreeBSD=%H
\ No newline at end of property
Added: svn:mime-type
## -0,0 +1 ##
+text/plain
\ No newline at end of property
Index: head/lib/msun/ld80/s_erfl.c
===================================================================
--- head/lib/msun/ld80/s_erfl.c	(nonexistent)
+++ head/lib/msun/ld80/s_erfl.c	(revision 268593)
@@ -0,0 +1,337 @@
+/* @(#)s_erf.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/*
+ * See s_erf.c for complete comments.
+ *
+ * Converted to long double by Steven G. Kargl.
+ */
+#include <float.h>
+#ifdef __i386__
+#include <ieeefp.h>
+#endif
+
+#include "fpmath.h"
+#include "math.h"
+#include "math_private.h"
+
+/* XXX Prevent compilers from erroneously constant folding: */
+static const volatile long double tiny = 0x1p-10000L;
+
+static const double
+half= 0.5,
+one = 1,
+two = 2;
+/*
+ * In the domain [0, 2**-34], only the first term in the power series
+ * expansion of erf(x) is used.  The magnitude of the first neglected
+ * terms is less than 2**-102.
+ */
+static const union IEEEl2bits
+efxu  = LD80C(0x8375d410a6db446c, -3,  1.28379167095512573902e-1L),
+efx8u = LD80C(0x8375d410a6db446c,  0,  1.02703333676410059122e+0L),
+/*
+ * Domain [0, 0.84375], range ~[-1.423e-22, 1.423e-22]:
+ * |(erf(x) - x)/x - pp(x)/qq(x)| < 2**-72.573
+ */
+pp0u  = LD80C(0x8375d410a6db446c, -3,   1.28379167095512573902e-1L),
+pp1u  = LD80C(0xa46c7d09ec3d0cec, -2,  -3.21140201054840180596e-1L),
+pp2u  = LD80C(0x9b31e66325576f86, -5,  -3.78893851760347812082e-2L),
+pp3u  = LD80C(0x804ac72c9a0b97dd, -7,  -7.83032847030604679616e-3L),
+pp4u  = LD80C(0x9f42bcbc3d5a601d, -12, -3.03765663857082048459e-4L),
+pp5u  = LD80C(0x9ec4ad6193470693, -16, -1.89266527398167917502e-5L),
+qq1u  = LD80C(0xdb4b8eb713188d6b, -2,   4.28310832832310510579e-1L),
+qq2u  = LD80C(0xa5750835b2459bd1, -4,   8.07896272074540216658e-2L),
+qq3u  = LD80C(0x8b85d6bd6a90b51c, -7,   8.51579638189385354266e-3L),
+qq4u  = LD80C(0x87332f82cff4ff96, -11,  5.15746855583604912827e-4L),
+qq5u  = LD80C(0x83466cb6bf9dca00, -16,  1.56492109706256700009e-5L),
+qq6u  = LD80C(0xf5bf98c2f996bf63, -24,  1.14435527803073879724e-7L);
+#define	efx	(efxu.e)
+#define	efx8	(efx8u.e)
+#define	pp0	(pp0u.e)
+#define	pp1	(pp1u.e)
+#define	pp2	(pp2u.e)
+#define	pp3	(pp3u.e)
+#define	pp4	(pp4u.e)
+#define	pp5	(pp5u.e)
+#define	qq1	(qq1u.e)
+#define	qq2	(qq2u.e)
+#define	qq3	(qq3u.e)
+#define	qq4	(qq4u.e)
+#define	qq5	(qq5u.e)
+#define	qq6	(qq6u.e)
+static const union IEEEl2bits
+erxu  = LD80C(0xd7bb3d0000000000, -1,  8.42700779438018798828e-1L),
+/*
+ * Domain [0.84375, 1.25], range ~[-8.132e-22, 8.113e-22]:
+ * |(erf(x) - erx) - pa(x)/qa(x)| < 2**-71.762
+ */
+pa0u  = LD80C(0xe8211158da02c692, -27,  1.35116960705131296711e-8L),
+pa1u  = LD80C(0xd488f89f36988618, -2,   4.15107507167065612570e-1L),
+pa2u  = LD80C(0xece74f8c63fa3942, -4,  -1.15675565215949226989e-1L),
+pa3u  = LD80C(0xc8d31e020727c006, -4,   9.80589241379624665791e-2L),
+pa4u  = LD80C(0x985d5d5fafb0551f, -5,   3.71984145558422368847e-2L),
+pa5u  = LD80C(0xa5b6c4854d2f5452, -8,  -5.05718799340957673661e-3L),
+pa6u  = LD80C(0x85c8d58fe3993a47, -8,   4.08277919612202243721e-3L),
+pa7u  = LD80C(0xddbfbc23677b35cf, -13,  2.11476292145347530794e-4L),
+qa1u  = LD80C(0xb8a977896f5eff3f, -1,   7.21335860303380361298e-1L),
+qa2u  = LD80C(0x9fcd662c3d4eac86, -1,   6.24227891731886593333e-1L),
+qa3u  = LD80C(0x9d0b618eac67ba07, -2,   3.06727455774491855801e-1L),
+qa4u  = LD80C(0x881a4293f6d6c92d, -3,   1.32912674218195890535e-1L),
+qa5u  = LD80C(0xbab144f07dea45bf, -5,   4.55792134233613027584e-2L),
+qa6u  = LD80C(0xa6c34ba438bdc900, -7,   1.01783980070527682680e-2L),
+qa7u  = LD80C(0x8fa866dc20717a91, -9,   2.19204436518951438183e-3L);
+#define erx	(erxu.e)
+#define pa0	(pa0u.e)
+#define pa1	(pa1u.e)
+#define pa2	(pa2u.e)
+#define pa3	(pa3u.e)
+#define pa4	(pa4u.e)
+#define pa5	(pa5u.e)
+#define pa6	(pa6u.e)
+#define pa7	(pa7u.e)
+#define qa1	(qa1u.e)
+#define qa2	(qa2u.e)
+#define qa3	(qa3u.e)
+#define qa4	(qa4u.e)
+#define qa5	(qa5u.e)
+#define qa6	(qa6u.e)
+#define qa7	(qa7u.e)
+static const union IEEEl2bits
+/*
+ * Domain [1.25,2.85715], range ~[-2.334e-22,2.334e-22]:
+ * |log(x*erfc(x)) + x**2 + 0.5625 - ra(x)/sa(x)| < 2**-71.860
+ */
+ra0u  = LD80C(0xa1a091e0fb4f335a, -7, -9.86494298915814308249e-3L),
+ra1u  = LD80C(0xc2b0d045ae37df6b, -1, -7.60510460864878271275e-1L),
+ra2u  = LD80C(0xf2cec3ee7da636c5, 3,  -1.51754798236892278250e+1L),
+ra3u  = LD80C(0x813cc205395adc7d, 7,  -1.29237335516455333420e+2L),
+ra4u  = LD80C(0x8737c8b7b4062c2f, 9,  -5.40871625829510494776e+2L),
+ra5u  = LD80C(0x8ffe5383c08d4943, 10, -1.15194769466026108551e+3L),
+ra6u  = LD80C(0x983573e64d5015a9, 10, -1.21767039790249025544e+3L),
+ra7u  = LD80C(0x92a794e763a6d4db, 9,  -5.86618463370624636688e+2L),
+ra8u  = LD80C(0xd5ad1fae77c3d9a3, 6,  -1.06838132335777049840e+2L),
+ra9u  = LD80C(0x934c1a247807bb9c, 2,  -4.60303980944467334806e+0L),
+sa1u  = LD80C(0xd342f90012bb1189, 4,   2.64077014928547064865e+1L),
+sa2u  = LD80C(0x839be13d9d5da883, 8,   2.63217811300123973067e+2L),
+sa3u  = LD80C(0x9f8cba6d1ae1b24b, 10,  1.27639775710344617587e+3L),
+sa4u  = LD80C(0xcaa83f403713e33e, 11,  3.24251544209971162003e+3L),
+sa5u  = LD80C(0x8796aff2f3c47968, 12,  4.33883591261332837874e+3L),
+sa6u  = LD80C(0xb6ef97f9c753157b, 11,  2.92697460344182158454e+3L),
+sa7u  = LD80C(0xe02aee5f83773d1c, 9,   8.96670799139389559818e+2L),
+sa8u  = LD80C(0xc82b83855b88e07e, 6,   1.00084987800048510018e+2L),
+sa9u  = LD80C(0x92f030aefadf28ad, 1,   2.29591004455459083843e+0L);
+#define ra0	(ra0u.e)
+#define ra1	(ra1u.e)
+#define ra2	(ra2u.e)
+#define ra3	(ra3u.e)
+#define ra4	(ra4u.e)
+#define ra5	(ra5u.e)
+#define ra6	(ra6u.e)
+#define ra7	(ra7u.e)
+#define ra8	(ra8u.e)
+#define ra9	(ra9u.e)
+#define sa1	(sa1u.e)
+#define sa2	(sa2u.e)
+#define sa3	(sa3u.e)
+#define sa4	(sa4u.e)
+#define sa5	(sa5u.e)
+#define sa6	(sa6u.e)
+#define sa7	(sa7u.e)
+#define sa8	(sa8u.e)
+#define sa9	(sa9u.e)
+/*
+ * Domain [2.85715,7], range ~[-8.323e-22,8.390e-22]:
+ * |log(x*erfc(x)) + x**2 + 0.5625 - rb(x)/sb(x)| < 2**-70.326
+ */
+static const union IEEEl2bits
+rb0u = LD80C(0xa1a091cf43abcd26, -7, -9.86494292470284646962e-3L),
+rb1u = LD80C(0xd19d2df1cbb8da0a, -1, -8.18804618389296662837e-1L),
+rb2u = LD80C(0x9a4dd1383e5daf5b, 4,  -1.92879967111618594779e+1L),
+rb3u = LD80C(0xbff0ae9fc0751de6, 7,  -1.91940164551245394969e+2L),
+rb4u = LD80C(0xdde08465310b472b, 9,  -8.87508080766577324539e+2L),
+rb5u = LD80C(0xe796e1d38c8c70a9, 10, -1.85271506669474503781e+3L),
+rb6u = LD80C(0xbaf655a76e0ab3b5, 10, -1.49569795581333675349e+3L),
+rb7u = LD80C(0x95d21e3e75503c21, 8,  -2.99641547972948019157e+2L),
+sb1u = LD80C(0x814487ed823c8cbd, 5,   3.23169247732868256569e+1L),
+sb2u = LD80C(0xbe4bfbb1301304be, 8,   3.80593618534539961773e+2L),
+sb3u = LD80C(0x809c4ade46b927c7, 11,  2.05776827838541292848e+3L),
+sb4u = LD80C(0xa55284359f3395a8, 12,  5.29031455540062116327e+3L),
+sb5u = LD80C(0xbcfa72da9b820874, 12,  6.04730608102312640462e+3L),
+sb6u = LD80C(0x9d09a35988934631, 11,  2.51260238030767176221e+3L),
+sb7u = LD80C(0xd675bbe542c159fa, 7,   2.14459898308561015684e+2L);
+#define rb0	(rb0u.e)
+#define rb1	(rb1u.e)
+#define rb2	(rb2u.e)
+#define rb3	(rb3u.e)
+#define rb4	(rb4u.e)
+#define rb5	(rb5u.e)
+#define rb6	(rb6u.e)
+#define rb7	(rb7u.e)
+#define sb1	(sb1u.e)
+#define sb2	(sb2u.e)
+#define sb3	(sb3u.e)
+#define sb4	(sb4u.e)
+#define sb5	(sb5u.e)
+#define sb6	(sb6u.e)
+#define sb7	(sb7u.e)
+/*
+ * Domain [7,108], range ~[-4.422e-22,4.422e-22]:
+ * |log(x*erfc(x)) + x**2 + 0.5625 - rc(x)/sc(x)| < 2**-70.938
+ */
+static const union IEEEl2bits
+/* err = -4.422092275318925082e-22 -70.937689 */
+rc0u = LD80C(0xa1a091cf437a17ad, -7, -9.86494292470008707260e-3L),
+rc1u = LD80C(0xbe79c5a978122b00, -1, -7.44045595049165939261e-1L),
+rc2u = LD80C(0xdb26f9bbe31a2794, 3,  -1.36970155085888424425e+1L),
+rc3u = LD80C(0xb5f69a38f5747ac8, 6,  -9.09816453742625888546e+1L),
+rc4u = LD80C(0xd79676d970d0a21a, 7,  -2.15587750997584074147e+2L),
+rc5u = LD80C(0xfe528153c45ec97c, 6,  -1.27161142938347796666e+2L),
+sc1u = LD80C(0xc5e8cd46d5604a96, 4,   2.47386727842204312937e+1L),
+sc2u = LD80C(0xc5f0f5a5484520eb, 7,   1.97941248254913378865e+2L),
+sc3u = LD80C(0x964e3c7b34db9170, 9,   6.01222441484087787522e+2L),
+sc4u = LD80C(0x99be1b89faa0596a, 9,   6.14970430845978077827e+2L),
+sc5u = LD80C(0xf80dfcbf37ffc5ea, 6,   1.24027318931184605891e+2L);
+#define rc0	(rc0u.e)
+#define rc1	(rc1u.e)
+#define rc2	(rc2u.e)
+#define rc3	(rc3u.e)
+#define rc4	(rc4u.e)
+#define rc5	(rc5u.e)
+#define sc1	(sc1u.e)
+#define sc2	(sc2u.e)
+#define sc3	(sc3u.e)
+#define sc4	(sc4u.e)
+#define sc5	(sc5u.e)
+
+long double
+erfl(long double x)
+{
+	long double ax,R,S,P,Q,s,y,z,r;
+	uint64_t lx;
+	int32_t i;
+	uint16_t hx;
+
+	EXTRACT_LDBL80_WORDS(hx, lx, x);
+
+	if((hx & 0x7fff) == 0x7fff) {	/* erfl(nan)=nan */
+		i = (hx>>15)<<1;
+		return (1-i)+one/x;	/* erfl(+-inf)=+-1 */
+	}
+
+	ENTERI();
+
+	ax = fabsl(x);
+	if(ax < 0.84375) {
+	    if(ax < 0x1p-34L) {
+	        if(ax < 0x1p-16373L)	
+		    RETURNI((8*x+efx8*x)/8);	/* avoid spurious underflow */
+		RETURNI(x + efx*x);
+	    }
+	    z = x*x;
+	    r = pp0+z*(pp1+z*(pp2+z*(pp3+z*(pp4+z*pp5))));
+	    s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*(qq5+z*qq6)))));
+	    y = r/s;
+	    RETURNI(x + x*y);
+	}
+	if(ax < 1.25) {
+	    s = ax-one;
+	    P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*(pa6+s*pa7))))));
+	    Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*(qa6+s*qa7))))));
+	    if(x>=0) RETURNI(erx + P/Q); else RETURNI(-erx - P/Q);
+	}
+	if(ax >= 7) {			/* inf>|x|>= 7 */
+	    if(x>=0) RETURNI(one-tiny); else RETURNI(tiny-one);
+	}
+	s = one/(ax*ax);
+	if(ax < 2.85715) {	/* |x| < 2.85715 */
+	    R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(ra5+s*(ra6+s*(ra7+
+		s*(ra8+s*ra9))))))));
+	    S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(sa5+s*(sa6+s*(sa7+
+		s*(sa8+s*sa9))))))));
+	} else {	/* |x| >= 2.85715 */
+	    R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(rb5+s*(rb6+s*rb7))))));
+	    S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(sb5+s*(sb6+s*sb7))))));
+	}
+	z=(float)ax;
+	r=expl(-z*z-0.5625)*expl((z-ax)*(z+ax)+R/S);
+	if(x>=0) RETURNI(one-r/ax); else RETURNI(r/ax-one);
+}
+
+long double
+erfcl(long double x)
+{
+	long double ax,R,S,P,Q,s,y,z,r;
+	uint64_t lx;
+	uint16_t hx;
+
+	EXTRACT_LDBL80_WORDS(hx, lx, x);
+
+	if((hx & 0x7fff) == 0x7fff) {	/* erfcl(nan)=nan */
+					/* erfcl(+-inf)=0,2 */
+	    return ((hx>>15)<<1)+one/x;
+	}
+
+	ENTERI();
+
+	ax = fabsl(x);
+	if(ax < 0.84375L) {
+	    if(ax < 0x1p-34L)
+		RETURNI(one-x);
+	    z = x*x;
+	    r = pp0+z*(pp1+z*(pp2+z*(pp3+z*(pp4+z*pp5))));
+	    s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*(qq5+z*qq6)))));
+	    y = r/s;
+	    if(ax < 0.25L) {  	/* x<1/4 */
+		RETURNI(one-(x+x*y));
+	    } else {
+		r = x*y;
+		r += (x-half);
+	       RETURNI(half - r);
+	    }
+	}
+	if(ax < 1.25L) {
+	    s = ax-one;
+	    P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*(pa6+s*pa7))))));
+	    Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*(qa6+s*qa7))))));
+	    if(x>=0) {
+	        z  = one-erx; RETURNI(z - P/Q);
+	    } else {
+		z = (erx+P/Q); RETURNI(one+z);
+	    }
+	}
+
+	if(ax < 108) {			/* |x| < 108 */
+ 	    s = one/(ax*ax);
+	    if(ax < 2.85715) {		/* |x| < 2.85715 */
+		R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(ra5+s*(ra6+s*(ra7+
+		    s*(ra8+s*ra9))))))));
+		S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(sa5+s*(sa6+s*(sa7+
+		    s*(sa8+s*sa9))))))));
+	    } else if(ax < 7) {		/* | |x| < 7 */
+		R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(rb5+s*(rb6+s*rb7))))));
+		S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(sb5+s*(sb6+s*sb7))))));
+	    } else {
+		if(x < -7) RETURNI(two-tiny);/* x < -7 */
+		R=rc0+s*(rc1+s*(rc2+s*(rc3+s*(rc4+s*rc5))));
+		S=one+s*(sc1+s*(sc2+s*(sc3+s*(sc4+s*sc5))));
+	    }
+	    z = (float)ax;
+	    r = expl(-z*z-0.5625)*expl((z-ax)*(z+ax)+R/S);
+	    if(x>0) RETURNI(r/ax); else RETURNI(two-r/ax);
+	} else {
+	    if(x>0) RETURNI(tiny*tiny); else RETURNI(two-tiny);
+	}
+}

Property changes on: head/lib/msun/ld80/s_erfl.c
___________________________________________________________________
Added: svn:eol-style
## -0,0 +1 ##
+native
\ No newline at end of property
Added: svn:keywords
## -0,0 +1 ##
+FreeBSD=%H
\ No newline at end of property
Added: svn:mime-type
## -0,0 +1 ##
+text/plain
\ No newline at end of property
Index: head/lib/msun/man/erf.3
===================================================================
--- head/lib/msun/man/erf.3	(revision 268592)
+++ head/lib/msun/man/erf.3	(revision 268593)
@@ -1,93 +1,98 @@
 .\" Copyright (c) 1985, 1991 Regents of the University of California.
 .\" 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.
 .\" 4. Neither the name of the University nor the names of its contributors
 .\"    may be used to endorse or promote products derived from this software
 .\"    without specific prior written permission.
 .\"
 .\" THIS SOFTWARE IS PROVIDED BY THE REGENTS 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 REGENTS 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.
 .\"
 .\"     from: @(#)erf.3	6.4 (Berkeley) 4/20/91
 .\" $FreeBSD$
 .\"
-.Dd April 20, 1991
+.Dd July 13, 2014
 .Dt ERF 3
 .Os
 .Sh NAME
 .Nm erf ,
 .Nm erff ,
+.Nm erfl ,
 .Nm erfc ,
-.Nm erfcf
+.Nm erfcf ,
+.Nm erfcl
 .Nd error function operators
 .Sh LIBRARY
 .Lb libm
 .Sh SYNOPSIS
 .In math.h
 .Ft double
 .Fn erf "double x"
 .Ft float
 .Fn erff "float x"
+.Ft "long double"
+.Fn erfl "long double x"
 .Ft double
 .Fn erfc "double x"
 .Ft float
 .Fn erfcf "float x"
+.Ft "long double"
+.Fn erfcl "long double x"
 .Sh DESCRIPTION
 These functions calculate the error function of
 .Fa x .
 .Pp
 The
-.Fn erf
-and the
-.Fn erff
+.Fn erf ,
+.Fn erff ,
+and
+.Fn erfl
 functions calculate the error function of x; where
 .Bd -ragged -offset indent
 .if n \{\
 erf(x) = 2/sqrt(pi)\(**\|integral from 0 to x of exp(\-t\(**t) dt.
 \}
 .if t \{\
 erf\|(x) :=
 (2/\(sr\(*p)\|\(is\d\s8\z0\s10\u\u\s8x\s10\d\|exp(\-t\u\s82\s10\d)\|dt.
 \}
 .Ed
 .Pp
 The
-.Fn erfc
-and the
-.Fn erfcf
+.Fn erfc ,
+.Fn erfcf ,
+and
+.Fn erfcl
 functions calculate the complementary error function of
 .Fa x ;
 that is
 .Fn erfc
 subtracts the result of the error function
 .Fn erf x
 from 1.0.
-This is useful, since for large
-.Fa x
-places disappear.
 .Sh SEE ALSO
 .Xr math 3
 .Sh HISTORY
 The
 .Fn erf
 and
 .Fn erfc
 functions appeared in
 .Bx 4.3 .
Index: head/lib/msun/src/imprecise.c
===================================================================
--- head/lib/msun/src/imprecise.c	(revision 268592)
+++ head/lib/msun/src/imprecise.c	(revision 268593)
@@ -1,66 +1,64 @@
 /*-
  * 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$
  */
 
 #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(erfc);
-DECLARE_IMPRECISE(erf);
 DECLARE_IMPRECISE(lgamma);
 DECLARE_IMPRECISE(tgamma);
Index: head/lib/msun/src/math.h
===================================================================
--- head/lib/msun/src/math.h	(revision 268592)
+++ head/lib/msun/src/math.h	(revision 268593)
@@ -1,522 +1,522 @@
 /*
  * ====================================================
  * 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.
  * ====================================================
  */
 
 /*
  * from: @(#)fdlibm.h 5.1 93/09/24
  * $FreeBSD$
  */
 
 #ifndef _MATH_H_
 #define	_MATH_H_
 
 #include <sys/cdefs.h>
 #include <sys/_types.h>
 #include <machine/_limits.h>
 
 /*
  * ANSI/POSIX
  */
 extern const union __infinity_un {
 	unsigned char	__uc[8];
 	double		__ud;
 } __infinity;
 
 extern const union __nan_un {
 	unsigned char	__uc[sizeof(float)];
 	float		__uf;
 } __nan;
 
 #if __GNUC_PREREQ__(3, 3) || (defined(__INTEL_COMPILER) && __INTEL_COMPILER >= 800)
 #define	__MATH_BUILTIN_CONSTANTS
 #endif
 
 #if __GNUC_PREREQ__(3, 0) && !defined(__INTEL_COMPILER)
 #define	__MATH_BUILTIN_RELOPS
 #endif
 
 #ifdef __MATH_BUILTIN_CONSTANTS
 #define	HUGE_VAL	__builtin_huge_val()
 #else
 #define	HUGE_VAL	(__infinity.__ud)
 #endif
 
 #if __ISO_C_VISIBLE >= 1999
 #define	FP_ILOGB0	(-__INT_MAX)
 #define	FP_ILOGBNAN	__INT_MAX
 
 #ifdef __MATH_BUILTIN_CONSTANTS
 #define	HUGE_VALF	__builtin_huge_valf()
 #define	HUGE_VALL	__builtin_huge_vall()
 #define	INFINITY	__builtin_inff()
 #define	NAN		__builtin_nanf("")
 #else
 #define	HUGE_VALF	(float)HUGE_VAL
 #define	HUGE_VALL	(long double)HUGE_VAL
 #define	INFINITY	HUGE_VALF
 #define	NAN		(__nan.__uf)
 #endif /* __MATH_BUILTIN_CONSTANTS */
 
 #define	MATH_ERRNO	1
 #define	MATH_ERREXCEPT	2
 #define	math_errhandling	MATH_ERREXCEPT
 
 #define	FP_FAST_FMAF	1
 
 /* Symbolic constants to classify floating point numbers. */
 #define	FP_INFINITE	0x01
 #define	FP_NAN		0x02
 #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) \
 	__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))
 #define	isgreaterequal(x, y)	__builtin_isgreaterequal((x), (y))
 #define	isless(x, y)		__builtin_isless((x), (y))
 #define	islessequal(x, y)	__builtin_islessequal((x), (y))
 #define	islessgreater(x, y)	__builtin_islessgreater((x), (y))
 #define	isunordered(x, y)	__builtin_isunordered((x), (y))
 #else
 #define	isgreater(x, y)		(!isunordered((x), (y)) && (x) > (y))
 #define	isgreaterequal(x, y)	(!isunordered((x), (y)) && (x) >= (y))
 #define	isless(x, y)		(!isunordered((x), (y)) && (x) < (y))
 #define	islessequal(x, y)	(!isunordered((x), (y)) && (x) <= (y))
 #define	islessgreater(x, y)	(!isunordered((x), (y)) && \
 					((x) > (y) || (y) > (x)))
 #define	isunordered(x, y)	(isnan(x) || isnan(y))
 #endif /* __MATH_BUILTIN_RELOPS */
 
 #define	signbit(x) __fp_type_select(x, __signbitf, __signbit, __signbitl)
 
 typedef	__double_t	double_t;
 typedef	__float_t	float_t;
 #endif /* __ISO_C_VISIBLE >= 1999 */
 
 /*
  * XOPEN/SVID
  */
 #if __BSD_VISIBLE || __XSI_VISIBLE
 #define	M_E		2.7182818284590452354	/* e */
 #define	M_LOG2E		1.4426950408889634074	/* log 2e */
 #define	M_LOG10E	0.43429448190325182765	/* log 10e */
 #define	M_LN2		0.69314718055994530942	/* log e2 */
 #define	M_LN10		2.30258509299404568402	/* log e10 */
 #define	M_PI		3.14159265358979323846	/* pi */
 #define	M_PI_2		1.57079632679489661923	/* pi/2 */
 #define	M_PI_4		0.78539816339744830962	/* pi/4 */
 #define	M_1_PI		0.31830988618379067154	/* 1/pi */
 #define	M_2_PI		0.63661977236758134308	/* 2/pi */
 #define	M_2_SQRTPI	1.12837916709551257390	/* 2/sqrt(pi) */
 #define	M_SQRT2		1.41421356237309504880	/* sqrt(2) */
 #define	M_SQRT1_2	0.70710678118654752440	/* 1/sqrt(2) */
 
 #define	MAXFLOAT	((float)3.40282346638528860e+38)
 extern int signgam;
 #endif /* __BSD_VISIBLE || __XSI_VISIBLE */
 
 #if __BSD_VISIBLE
 #if 0
 /* Old value from 4.4BSD-Lite math.h; this is probably better. */
 #define	HUGE		HUGE_VAL
 #else
 #define	HUGE		MAXFLOAT
 #endif
 #endif /* __BSD_VISIBLE */
 
 /*
  * Most of these functions depend on the rounding mode and have the side
  * effect of raising floating-point exceptions, so they are not declared
  * as __pure2.  In C99, FENV_ACCESS affects the purity of these functions.
  */
 __BEGIN_DECLS
 /*
  * ANSI/POSIX
  */
 int	__fpclassifyd(double) __pure2;
 int	__fpclassifyf(float) __pure2;
 int	__fpclassifyl(long double) __pure2;
 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	__isnormalf(float) __pure2;
 int	__isnormal(double) __pure2;
 int	__isnormall(long double) __pure2;
 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);
 double	atan2(double, double);
 double	cos(double);
 double	sin(double);
 double	tan(double);
 
 double	cosh(double);
 double	sinh(double);
 double	tanh(double);
 
 double	exp(double);
 double	frexp(double, int *);	/* fundamentally !__pure2 */
 double	ldexp(double, int);
 double	log(double);
 double	log10(double);
 double	modf(double, double *);	/* fundamentally !__pure2 */
 
 double	pow(double, double);
 double	sqrt(double);
 
 double	ceil(double);
 double	fabs(double) __pure2;
 double	floor(double);
 double	fmod(double, double);
 
 /*
  * These functions are not in C90.
  */
 #if __BSD_VISIBLE || __ISO_C_VISIBLE >= 1999 || __XSI_VISIBLE
 double	acosh(double);
 double	asinh(double);
 double	atanh(double);
 double	cbrt(double);
 double	erf(double);
 double	erfc(double);
 double	exp2(double);
 double	expm1(double);
 double	fma(double, double, double);
 double	hypot(double, double);
 int	ilogb(double) __pure2;
 double	lgamma(double);
 long long llrint(double);
 long long llround(double);
 double	log1p(double);
 double	log2(double);
 double	logb(double);
 long	lrint(double);
 long	lround(double);
 double	nan(const char *) __pure2;
 double	nextafter(double, double);
 double	remainder(double, double);
 double	remquo(double, double, int *);
 double	rint(double);
 #endif /* __BSD_VISIBLE || __ISO_C_VISIBLE >= 1999 || __XSI_VISIBLE */
 
 #if __BSD_VISIBLE || __XSI_VISIBLE
 double	j0(double);
 double	j1(double);
 double	jn(int, double);
 double	y0(double);
 double	y1(double);
 double	yn(int, double);
 
 #if __XSI_VISIBLE <= 500 || __BSD_VISIBLE
 double	gamma(double);
 #endif
 
 #if __XSI_VISIBLE <= 600 || __BSD_VISIBLE
 double	scalb(double, double);
 #endif
 #endif /* __BSD_VISIBLE || __XSI_VISIBLE */
 
 #if __BSD_VISIBLE || __ISO_C_VISIBLE >= 1999
 double	copysign(double, double) __pure2;
 double	fdim(double, double);
 double	fmax(double, double) __pure2;
 double	fmin(double, double) __pure2;
 double	nearbyint(double);
 double	round(double);
 double	scalbln(double, long);
 double	scalbn(double, int);
 double	tgamma(double);
 double	trunc(double);
 #endif
 
 /*
  * BSD math library entry points
  */
 #if __BSD_VISIBLE
 double	drem(double, double);
 int	finite(double) __pure2;
 int	isnanf(float) __pure2;
 
 /*
  * Reentrant version of gamma & lgamma; passes signgam back by reference
  * as the second argument; user must allocate space for signgam.
  */
 double	gamma_r(double, int *);
 double	lgamma_r(double, int *);
 
 /*
  * IEEE Test Vector
  */
 double	significand(double);
 #endif /* __BSD_VISIBLE */
 
 /* float versions of ANSI/POSIX functions */
 #if __ISO_C_VISIBLE >= 1999
 float	acosf(float);
 float	asinf(float);
 float	atanf(float);
 float	atan2f(float, float);
 float	cosf(float);
 float	sinf(float);
 float	tanf(float);
 
 float	coshf(float);
 float	sinhf(float);
 float	tanhf(float);
 
 float	exp2f(float);
 float	expf(float);
 float	expm1f(float);
 float	frexpf(float, int *);	/* fundamentally !__pure2 */
 int	ilogbf(float) __pure2;
 float	ldexpf(float, int);
 float	log10f(float);
 float	log1pf(float);
 float	log2f(float);
 float	logf(float);
 float	modff(float, float *);	/* fundamentally !__pure2 */
 
 float	powf(float, float);
 float	sqrtf(float);
 
 float	ceilf(float);
 float	fabsf(float) __pure2;
 float	floorf(float);
 float	fmodf(float, float);
 float	roundf(float);
 
 float	erff(float);
 float	erfcf(float);
 float	hypotf(float, float);
 float	lgammaf(float);
 float	tgammaf(float);
 
 float	acoshf(float);
 float	asinhf(float);
 float	atanhf(float);
 float	cbrtf(float);
 float	logbf(float);
 float	copysignf(float, float) __pure2;
 long long llrintf(float);
 long long llroundf(float);
 long	lrintf(float);
 long	lroundf(float);
 float	nanf(const char *) __pure2;
 float	nearbyintf(float);
 float	nextafterf(float, float);
 float	remainderf(float, float);
 float	remquof(float, float, int *);
 float	rintf(float);
 float	scalblnf(float, long);
 float	scalbnf(float, int);
 float	truncf(float);
 
 float	fdimf(float, float);
 float	fmaf(float, float, float);
 float	fmaxf(float, float) __pure2;
 float	fminf(float, float) __pure2;
 #endif
 
 /*
  * float versions of BSD math library entry points
  */
 #if __BSD_VISIBLE
 float	dremf(float, float);
 int	finitef(float) __pure2;
 float	gammaf(float);
 float	j0f(float);
 float	j1f(float);
 float	jnf(int, float);
 float	scalbf(float, float);
 float	y0f(float);
 float	y1f(float);
 float	ynf(int, float);
 
 /*
  * Float versions of reentrant version of gamma & lgamma; passes
  * signgam back by reference as the second argument; user must
  * allocate space for signgam.
  */
 float	gammaf_r(float, int *);
 float	lgammaf_r(float, int *);
 
 /*
  * float version of IEEE Test Vector
  */
 float	significandf(float);
 #endif	/* __BSD_VISIBLE */
 
 /*
  * 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);
 long double	copysignl(long double, long double) __pure2;
 long double	coshl(long double);
 long double	cosl(long double);
+long double	erfcl(long double);
+long double	erfl(long double);
 long double	exp2l(long double);
 long double	expl(long double);
 long double	expm1l(long double);
 long double	fabsl(long double) __pure2;
 long double	fdiml(long double, long double);
 long double	floorl(long double);
 long double	fmal(long double, long double, long double);
 long double	fmaxl(long double, long double) __pure2;
 long double	fminl(long double, long double) __pure2;
 long double	fmodl(long double, long double);
 long double	frexpl(long double value, int *); /* fundamentally !__pure2 */
 long double	hypotl(long double, long double);
 int		ilogbl(long double) __pure2;
 long double	ldexpl(long double, int);
 long 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 */
 long double	nanl(const char *) __pure2;
 long double	nearbyintl(long double);
 long double	nextafterl(long double, long double);
 double		nexttoward(double, long double);
 float		nexttowardf(float, long double);
 long double	nexttowardl(long double, long double);
 long double	remainderl(long double, long double);
 long double	remquol(long double, long double, int *);
 long double	rintl(long double);
 long double	roundl(long double);
 long double	scalblnl(long double, long);
 long double	scalbnl(long double, int);
 long double	sinhl(long double);
 long double	sinl(long double);
 long double	sqrtl(long double);
 long double	tanhl(long double);
 long double	tanl(long double);
 long double	truncl(long double);
 
 #endif /* __ISO_C_VISIBLE >= 1999 */
 __END_DECLS
 
 #endif /* !_MATH_H_ */
 
 /* separate header for cmath */
 #ifndef _MATH_EXTRA_H_
 #if __ISO_C_VISIBLE >= 1999
 #if _DECLARE_C99_LDBL_MATH
 
 #define _MATH_EXTRA_H_
 
 /*
  * extra long double versions of math functions for C99 and cmath
  */
 __BEGIN_DECLS
 
-long double	erfcl(long double);
-long double	erfl(long double);
 long double	lgammal(long double);
 long double	powl(long double, long double);
 long double	tgammal(long double);
 
 __END_DECLS
 
 #endif /* !_DECLARE_C99_LDBL_MATH */
 #endif /* __ISO_C_VISIBLE >= 1999 */
 #endif /* !_MATH_EXTRA_H_ */
Index: head/lib/msun/src/s_erf.c
===================================================================
--- head/lib/msun/src/s_erf.c	(revision 268592)
+++ head/lib/msun/src/s_erf.c	(revision 268593)
@@ -1,301 +1,309 @@
 /* @(#)s_erf.c 5.1 93/09/24 */
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
  * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
 
 #include <sys/cdefs.h>
 __FBSDID("$FreeBSD$");
 
 /* double erf(double x)
  * double erfc(double x)
  *			     x
  *		      2      |\
  *     erf(x)  =  ---------  | exp(-t*t)dt
  *	 	   sqrt(pi) \|
  *			     0
  *
  *     erfc(x) =  1-erf(x)
  *  Note that
  *		erf(-x) = -erf(x)
  *		erfc(-x) = 2 - erfc(x)
  *
  * Method:
  *	1. For |x| in [0, 0.84375]
  *	    erf(x)  = x + x*R(x^2)
  *          erfc(x) = 1 - erf(x)           if x in [-.84375,0.25]
  *                  = 0.5 + ((0.5-x)-x*R)  if x in [0.25,0.84375]
  *	   where R = P/Q where P is an odd poly of degree 8 and
  *	   Q is an odd poly of degree 10.
  *						 -57.90
  *			| R - (erf(x)-x)/x | <= 2
  *
  *
  *	   Remark. The formula is derived by noting
  *          erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10 - x^7/42 + ....)
  *	   and that
  *          2/sqrt(pi) = 1.128379167095512573896158903121545171688
  *	   is close to one. The interval is chosen because the fix
  *	   point of erf(x) is near 0.6174 (i.e., erf(x)=x when x is
  *	   near 0.6174), and by some experiment, 0.84375 is chosen to
  * 	   guarantee the error is less than one ulp for erf.
  *
  *      2. For |x| in [0.84375,1.25], let s = |x| - 1, and
  *         c = 0.84506291151 rounded to single (24 bits)
  *         	erf(x)  = sign(x) * (c  + P1(s)/Q1(s))
  *         	erfc(x) = (1-c)  - P1(s)/Q1(s) if x > 0
  *			  1+(c+P1(s)/Q1(s))    if x < 0
  *         	|P1/Q1 - (erf(|x|)-c)| <= 2**-59.06
  *	   Remark: here we use the taylor series expansion at x=1.
  *		erf(1+s) = erf(1) + s*Poly(s)
  *			 = 0.845.. + P1(s)/Q1(s)
  *	   That is, we use rational approximation to approximate
  *			erf(1+s) - (c = (single)0.84506291151)
  *	   Note that |P1/Q1|< 0.078 for x in [0.84375,1.25]
  *	   where
  *		P1(s) = degree 6 poly in s
  *		Q1(s) = degree 6 poly in s
  *
  *      3. For x in [1.25,1/0.35(~2.857143)],
  *         	erfc(x) = (1/x)*exp(-x*x-0.5625+R1/S1)
  *         	erf(x)  = 1 - erfc(x)
  *	   where
  *		R1(z) = degree 7 poly in z, (z=1/x^2)
  *		S1(z) = degree 8 poly in z
  *
  *      4. For x in [1/0.35,28]
  *         	erfc(x) = (1/x)*exp(-x*x-0.5625+R2/S2) if x > 0
  *			= 2.0 - (1/x)*exp(-x*x-0.5625+R2/S2) if -6<x<0
  *			= 2.0 - tiny		(if x <= -6)
  *         	erf(x)  = sign(x)*(1.0 - erfc(x)) if x < 6, else
  *         	erf(x)  = sign(x)*(1.0 - tiny)
  *	   where
  *		R2(z) = degree 6 poly in z, (z=1/x^2)
  *		S2(z) = degree 7 poly in z
  *
  *      Note1:
  *	   To compute exp(-x*x-0.5625+R/S), let s be a single
  *	   precision number and s := x; then
  *		-x*x = -s*s + (s-x)*(s+x)
  *	        exp(-x*x-0.5626+R/S) =
  *			exp(-s*s-0.5625)*exp((s-x)*(s+x)+R/S);
  *      Note2:
  *	   Here 4 and 5 make use of the asymptotic series
  *			  exp(-x*x)
  *		erfc(x) ~ ---------- * ( 1 + Poly(1/x^2) )
  *			  x*sqrt(pi)
  *	   We use rational approximation to approximate
  *      	g(s)=f(1/x^2) = log(erfc(x)*x) - x*x + 0.5625
  *	   Here is the error bound for R1/S1 and R2/S2
  *      	|R1/S1 - f(x)|  < 2**(-62.57)
  *      	|R2/S2 - f(x)|  < 2**(-61.52)
  *
  *      5. For inf > x >= 28
  *         	erf(x)  = sign(x) *(1 - tiny)  (raise inexact)
  *         	erfc(x) = tiny*tiny (raise underflow) if x > 0
  *			= 2 - tiny if x<0
  *
  *      7. Special case:
  *         	erf(0)  = 0, erf(inf)  = 1, erf(-inf) = -1,
  *         	erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2,
  *	   	erfc/erf(NaN) is NaN
  */
 
 
 #include "math.h"
 #include "math_private.h"
 
 /* XXX Prevent compilers from erroneously constant folding: */
 static const volatile double tiny= 1e-300;
 
 static const double
 half= 0.5,
 one = 1,
 two = 2,
 /* c = (float)0.84506291151 */
 erx =  8.45062911510467529297e-01, /* 0x3FEB0AC1, 0x60000000 */
 /*
  * In the domain [0, 2**-28], only the first term in the power series
  * expansion of erf(x) is used.  The magnitude of the first neglected
  * terms is less than 2**-84.
  */
 efx =  1.28379167095512586316e-01, /* 0x3FC06EBA, 0x8214DB69 */
 efx8=  1.02703333676410069053e+00, /* 0x3FF06EBA, 0x8214DB69 */
 /*
  * Coefficients for approximation to erf on [0,0.84375]
  */
 pp0  =  1.28379167095512558561e-01, /* 0x3FC06EBA, 0x8214DB68 */
 pp1  = -3.25042107247001499370e-01, /* 0xBFD4CD7D, 0x691CB913 */
 pp2  = -2.84817495755985104766e-02, /* 0xBF9D2A51, 0xDBD7194F */
 pp3  = -5.77027029648944159157e-03, /* 0xBF77A291, 0x236668E4 */
 pp4  = -2.37630166566501626084e-05, /* 0xBEF8EAD6, 0x120016AC */
 qq1  =  3.97917223959155352819e-01, /* 0x3FD97779, 0xCDDADC09 */
 qq2  =  6.50222499887672944485e-02, /* 0x3FB0A54C, 0x5536CEBA */
 qq3  =  5.08130628187576562776e-03, /* 0x3F74D022, 0xC4D36B0F */
 qq4  =  1.32494738004321644526e-04, /* 0x3F215DC9, 0x221C1A10 */
 qq5  = -3.96022827877536812320e-06, /* 0xBED09C43, 0x42A26120 */
 /*
  * Coefficients for approximation to erf in [0.84375,1.25]
  */
 pa0  = -2.36211856075265944077e-03, /* 0xBF6359B8, 0xBEF77538 */
 pa1  =  4.14856118683748331666e-01, /* 0x3FDA8D00, 0xAD92B34D */
 pa2  = -3.72207876035701323847e-01, /* 0xBFD7D240, 0xFBB8C3F1 */
 pa3  =  3.18346619901161753674e-01, /* 0x3FD45FCA, 0x805120E4 */
 pa4  = -1.10894694282396677476e-01, /* 0xBFBC6398, 0x3D3E28EC */
 pa5  =  3.54783043256182359371e-02, /* 0x3FA22A36, 0x599795EB */
 pa6  = -2.16637559486879084300e-03, /* 0xBF61BF38, 0x0A96073F */
 qa1  =  1.06420880400844228286e-01, /* 0x3FBB3E66, 0x18EEE323 */
 qa2  =  5.40397917702171048937e-01, /* 0x3FE14AF0, 0x92EB6F33 */
 qa3  =  7.18286544141962662868e-02, /* 0x3FB2635C, 0xD99FE9A7 */
 qa4  =  1.26171219808761642112e-01, /* 0x3FC02660, 0xE763351F */
 qa5  =  1.36370839120290507362e-02, /* 0x3F8BEDC2, 0x6B51DD1C */
 qa6  =  1.19844998467991074170e-02, /* 0x3F888B54, 0x5735151D */
 /*
  * Coefficients for approximation to erfc in [1.25,1/0.35]
  */
 ra0  = -9.86494403484714822705e-03, /* 0xBF843412, 0x600D6435 */
 ra1  = -6.93858572707181764372e-01, /* 0xBFE63416, 0xE4BA7360 */
 ra2  = -1.05586262253232909814e+01, /* 0xC0251E04, 0x41B0E726 */
 ra3  = -6.23753324503260060396e+01, /* 0xC04F300A, 0xE4CBA38D */
 ra4  = -1.62396669462573470355e+02, /* 0xC0644CB1, 0x84282266 */
 ra5  = -1.84605092906711035994e+02, /* 0xC067135C, 0xEBCCABB2 */
 ra6  = -8.12874355063065934246e+01, /* 0xC0545265, 0x57E4D2F2 */
 ra7  = -9.81432934416914548592e+00, /* 0xC023A0EF, 0xC69AC25C */
 sa1  =  1.96512716674392571292e+01, /* 0x4033A6B9, 0xBD707687 */
 sa2  =  1.37657754143519042600e+02, /* 0x4061350C, 0x526AE721 */
 sa3  =  4.34565877475229228821e+02, /* 0x407B290D, 0xD58A1A71 */
 sa4  =  6.45387271733267880336e+02, /* 0x40842B19, 0x21EC2868 */
 sa5  =  4.29008140027567833386e+02, /* 0x407AD021, 0x57700314 */
 sa6  =  1.08635005541779435134e+02, /* 0x405B28A3, 0xEE48AE2C */
 sa7  =  6.57024977031928170135e+00, /* 0x401A47EF, 0x8E484A93 */
 sa8  = -6.04244152148580987438e-02, /* 0xBFAEEFF2, 0xEE749A62 */
 /*
  * Coefficients for approximation to erfc in [1/.35,28]
  */
 rb0  = -9.86494292470009928597e-03, /* 0xBF843412, 0x39E86F4A */
 rb1  = -7.99283237680523006574e-01, /* 0xBFE993BA, 0x70C285DE */
 rb2  = -1.77579549177547519889e+01, /* 0xC031C209, 0x555F995A */
 rb3  = -1.60636384855821916062e+02, /* 0xC064145D, 0x43C5ED98 */
 rb4  = -6.37566443368389627722e+02, /* 0xC083EC88, 0x1375F228 */
 rb5  = -1.02509513161107724954e+03, /* 0xC0900461, 0x6A2E5992 */
 rb6  = -4.83519191608651397019e+02, /* 0xC07E384E, 0x9BDC383F */
 sb1  =  3.03380607434824582924e+01, /* 0x403E568B, 0x261D5190 */
 sb2  =  3.25792512996573918826e+02, /* 0x40745CAE, 0x221B9F0A */
 sb3  =  1.53672958608443695994e+03, /* 0x409802EB, 0x189D5118 */
 sb4  =  3.19985821950859553908e+03, /* 0x40A8FFB7, 0x688C246A */
 sb5  =  2.55305040643316442583e+03, /* 0x40A3F219, 0xCEDF3BE6 */
 sb6  =  4.74528541206955367215e+02, /* 0x407DA874, 0xE79FE763 */
 sb7  = -2.24409524465858183362e+01; /* 0xC03670E2, 0x42712D62 */
 
 double
 erf(double x)
 {
 	int32_t hx,ix,i;
 	double R,S,P,Q,s,y,z,r;
 	GET_HIGH_WORD(hx,x);
 	ix = hx&0x7fffffff;
 	if(ix>=0x7ff00000) {		/* erf(nan)=nan */
 	    i = ((u_int32_t)hx>>31)<<1;
 	    return (double)(1-i)+one/x;	/* erf(+-inf)=+-1 */
 	}
 
 	if(ix < 0x3feb0000) {		/* |x|<0.84375 */
 	    if(ix < 0x3e300000) { 	/* |x|<2**-28 */
 	        if (ix < 0x00800000)
 		    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))));
 	    y = r/s;
 	    return x + x*y;
 	}
 	if(ix < 0x3ff40000) {		/* 0.84375 <= |x| < 1.25 */
 	    s = fabs(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)))));
 	    if(hx>=0) return erx + P/Q; else return -erx - P/Q;
 	}
 	if (ix >= 0x40180000) {		/* inf>|x|>=6 */
 	    if(hx>=0) return one-tiny; else return tiny-one;
 	}
 	x = fabs(x);
  	s = one/(x*x);
 	if(ix< 0x4006DB6E) {	/* |x| < 1/0.35 */
 	    R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(ra5+s*(ra6+s*ra7))))));
 	    S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(sa5+s*(sa6+s*(sa7+
 		s*sa8)))))));
 	} 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))))));
 	}
 	z  = x;
 	SET_LOW_WORD(z,0);
 	r  =  __ieee754_exp(-z*z-0.5625)*__ieee754_exp((z-x)*(z+x)+R/S);
 	if(hx>=0) return one-r/x; else return  r/x-one;
 }
 
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(erf, erfl);
+#endif
+
 double
 erfc(double x)
 {
 	int32_t hx,ix;
 	double R,S,P,Q,s,y,z,r;
 	GET_HIGH_WORD(hx,x);
 	ix = hx&0x7fffffff;
 	if(ix>=0x7ff00000) {			/* erfc(nan)=nan */
 						/* erfc(+-inf)=0,2 */
 	    return (double)(((u_int32_t)hx>>31)<<1)+one/x;
 	}
 
 	if(ix < 0x3feb0000) {		/* |x|<0.84375 */
 	    if(ix < 0x3c700000)  	/* |x|<2**-56 */
 		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))));
 	    y = r/s;
 	    if(hx < 0x3fd00000) {  	/* x<1/4 */
 		return one-(x+x*y);
 	    } else {
 		r = x*y;
 		r += (x-half);
 	        return half - r ;
 	    }
 	}
 	if(ix < 0x3ff40000) {		/* 0.84375 <= |x| < 1.25 */
 	    s = fabs(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)))));
 	    if(hx>=0) {
 	        z  = one-erx; return z - P/Q;
 	    } else {
 		z = erx+P/Q; return one+z;
 	    }
 	}
 	if (ix < 0x403c0000) {		/* |x|<28 */
 	    x = fabs(x);
  	    s = one/(x*x);
 	    if(ix< 0x4006DB6D) {	/* |x| < 1/.35 ~ 2.857143*/
 		R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(ra5+s*(ra6+s*ra7))))));
 		S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(sa5+s*(sa6+s*(sa7+
 		    s*sa8)))))));
 	    } else {			/* |x| >= 1/.35 ~ 2.857143 */
 		if(hx<0&&ix>=0x40180000) return two-tiny;/* x < -6 */
 		R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(rb5+s*rb6)))));
 		S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(sb5+s*(sb6+s*sb7))))));
 	    }
 	    z  = x;
 	    SET_LOW_WORD(z,0);
 	    r  =  __ieee754_exp(-z*z-0.5625)*__ieee754_exp((z-x)*(z+x)+R/S);
 	    if(hx>0) return r/x; else return two-r/x;
 	} else {
 	    if(hx>0) return tiny*tiny; else return two-tiny;
 	}
 }
+
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(erfc, erfcl);
+#endif