erff, erfl, and erfcl patch
Steve Kargl
sgk at troutmask.apl.washington.edu
Sun Mar 9 00:57:27 UTC 2014
On Mon, Dec 02, 2013 at 08:01:27PM -0800, Steve Kargl wrote:
> On Tue, Dec 03, 2013 at 11:56:49AM +1100, Bruce Evans wrote:
> > On Mon, 2 Dec 2013, Steve Kargl wrote:
> >
> > > I intend to commit the patch that follows by sig in the next day or
> > > two. In particular, I want to get the ld80 and ld128 versions of
> > > erfl and erfcl into tree as the current hack provided by imprecise.c
> > > is ugly.
> >
> > Looks good, and a lot of work. I don' have time to look at it closely
> > for a week or 2.
>
> For not having time to look over the diff closely, you've
> pointed out a few things I can fix/investigate.
>
Reviving an old. I finally had time to lookover the email
threads that I could find in my archive and I hoped to have
addressed your comments adequately. I plan to commit the
patch that follows my ./sig tomorrow. Whitespace changes
will be done first in a separate commit followed by the
addition of functional changes.
* lib/msun/Makefile:
. Hook s_erfl.c to the build.
. Add MLINKS for erfl.3 and erfcl.3
* lib/msun/Symbol.map:
. Move erfl and erfcl to their rightful place.
* lib/msun/man/erf.3:
. Document erfl and erfcl.
. Remove a nonsense sentence.
. Add ISO C99 conformance.
* lib/msun/src/s_erf.c:
. Adjust whitespace in declarations of efx and efx8.
. Add weak references for erfl and erfcl for architectures where
long double is the same as double.
* lib/msun/src/s_erff.c:
. Consistently use lower case in hexidecimal constants.
. Adjust the number of terms in some of the rational approximations.
Fix the domain and range statements for new coefficients.
* lib/msun/ld80/s_erfl.c:
. Intel 80-bit format long double implementations of erfl and erfcl.
* lib/msun/ld128/s_erfl.c:
. IEEE 128-bit format long double implementations of erfl and erfcl.
* lib/msun/src/imprecise.c:
. Remove the undocumented and questionable mapping of erf[c]l to erf[c].
--
Steve
Index: Makefile
===================================================================
--- Makefile (revision 262934)
+++ Makefile (working copy)
@@ -101,7 +101,7 @@
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
@@ -164,7 +164,7 @@
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
Index: Symbol.map
===================================================================
--- Symbol.map (revision 262934)
+++ Symbol.map (working copy)
@@ -265,6 +265,8 @@
ctanf;
ctanh;
ctanhf;
+ erfcl;
+ erfl;
expl;
expm1l;
log10l;
@@ -274,8 +276,6 @@
sinhl;
tanhl;
/* Implemented as weak aliases for imprecise versions */
- erfcl;
- erfl;
lgammal;
powl;
tgammal;
Index: man/erf.3
===================================================================
--- man/erf.3 (revision 262934)
+++ man/erf.3 (working copy)
@@ -28,14 +28,16 @@
.\" from: @(#)erf.3 6.4 (Berkeley) 4/20/91
.\" $FreeBSD$
.\"
-.Dd April 20, 1991
+.Dd March 5, 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
@@ -45,18 +47,23 @@
.Fn erf "double x"
.Ft float
.Fn erff "float x"
+.Ft "double long"
+.Fn erfl "long double x"
.Ft double
.Fn erfc "double x"
.Ft float
.Fn erfcf "float x"
+.Ft "double long"
+.Fn erfcl "long double x"
.Sh DESCRIPTION
These functions calculate the error function of
.Fa x .
.Pp
The
-.Fn erf
-and the
-.Fn erff
+.Fn erf ,
+.Fn erff ,
+and
+.Fn erfl
functions calculate the error function of x; where
.Bd -ragged -offset indent
.if n \{\
@@ -69,9 +76,10 @@
.Ed
.Pp
The
-.Fn erfc
-and the
-.Fn erfcf
+.Fn erfc ,
+.Fn erfcf ,
+and
+.Fn erfcl
functions calculate the complementary error function of
.Fa x ;
that is
@@ -79,9 +87,6 @@
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
@@ -91,3 +96,6 @@
.Fn erfc
functions appeared in
.Bx 4.3 .
+.Sh STANDARDS
+These functions conform to
+.St -isoC-99 .
Index: src/s_erf.c
===================================================================
--- src/s_erf.c (revision 262934)
+++ src/s_erf.c (working copy)
@@ -107,6 +107,7 @@
* erfc/erf(NaN) is NaN
*/
+#include <float.h>
#include "math.h"
#include "math_private.h"
@@ -121,8 +122,8 @@
/*
* Coefficients for approximation to erf on [0,0.84375]
*/
-efx = 1.28379167095512586316e-01, /* 0x3FC06EBA, 0x8214DB69 */
-efx8= 1.02703333676410069053e+00, /* 0x3FF06EBA, 0x8214DB69 */
+efx = 1.28379167095512586316e-01, /* 0x3FC06EBA, 0x8214DB69 */
+efx8= 1.02703333676410069053e+00, /* 0x3FF06EBA, 0x8214DB69 */
pp0 = 1.28379167095512558561e-01, /* 0x3FC06EBA, 0x8214DB68 */
pp1 = -3.25042107247001499370e-01, /* 0xBFD4CD7D, 0x691CB913 */
pp2 = -2.84817495755985104766e-02, /* 0xBF9D2A51, 0xDBD7194F */
@@ -238,6 +239,10 @@
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)
{
@@ -299,3 +304,7 @@
if(hx>0) return tiny*tiny; else return two-tiny;
}
}
+
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(erfc, erfcl);
+#endif
Index: src/s_erff.c
===================================================================
--- src/s_erff.c (revision 262934)
+++ src/s_erff.c (working copy)
@@ -21,62 +21,54 @@
static const float
tiny = 1e-30,
-half= 5.0000000000e-01, /* 0x3F000000 */
-one = 1.0000000000e+00, /* 0x3F800000 */
-two = 2.0000000000e+00, /* 0x40000000 */
+half= 5.0000000000e-1, /* 0x3f000000 */
+one = 1.0000000000e+0, /* 0x3f800000 */
+two = 2.0000000000e+0, /* 0x40000000 */
/*
- * Coefficients for approximation to erf on [0,0.84375]
+ * Domain [0, 0.84375], range ~[-5.4419e-10,5.5179e-10]:
+ * |(erf(x) - x)/x - p(x)/q(x)| < 2**-31
*/
-efx = 1.2837916613e-01, /* 0x3e0375d4 */
-efx8= 1.0270333290e+00, /* 0x3f8375d4 */
+efx = 1.28379166e-1, /* 0x3e0375d4 */
+efx8 = 1.02703333, /* 0x3f8375d4 */
+pp0 = 1.28379166e-1, /* 0x3e0375d4 */
+pp1 = -3.36030394e-1, /* 0xbeac0c2d */
+pp2 = -1.86261395e-3, /* 0xbaf422f4 */
+qq1 = 3.12324315e-1, /* 0x3e9fe8f9 */
+qq2 = 2.16070414e-2, /* 0x3cb10140 */
+qq3 = -1.98859372e-3, /* 0xbb025311 */
/*
- * Domain [0, 0.84375], range ~[-5.4446e-10,5.5197e-10]:
- * |(erf(x) - x)/x - p(x)/q(x)| < 2**-31.
+ * Domain [0.84375, 1.25], range ~[-1.023e-9,1.023e-9]:
+ * |(erf(x) - erx) - 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 */
+erx = 8.42697144e-1, /* 0x3f57bb00 */
+pa0 = 3.65041046e-6, /* 0x3674f993 */
+pa1 = 4.15109307e-1, /* 0x3ed48935 */
+pa2 = -2.09395722e-1, /* 0xbe566bd5 */
+pa3 = 8.67677554e-2, /* 0x3db1b34b */
+qa1 = 4.95560974e-1, /* 0x3efdba2b */
+qa2 = 3.71248513e-1, /* 0x3ebe1449 */
+qa3 = 3.92478965e-2, /* 0x3d20c267 */
/*
- * Domain [0.84375, 1.25], range ~[-1.953e-11,1.940e-11]:
- * |(erf(x) - erx) - p(x)/q(x)| < 2**-36.
+ * Domain [1.25,1/0.35], range ~[-4.821e-9,4.927e-9]:
+ * |log(x*erfc(x)) + x**2 + 0.5625 - r(x)/s(x)| < 2**-28
*/
-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 */
+ra0 = -9.88156721e-3, /* 0xbc21e64c */
+ra1 = -5.43658376e-1, /* 0xbf0b2d32 */
+ra2 = -1.66828310, /* 0xbfd58a4d */
+ra3 = -6.91554189e-1, /* 0xbf3109b2 */
+sa1 = 4.48581553, /* 0x408f8bcd */
+sa2 = 4.10799170, /* 0x408374ab */
+sa3 = 5.53855181e-1, /* 0x3f0dc974 */
/*
- * Domain [1.25,1/0.35], range ~[-7.043e-10,7.457e-10]:
+ * Domain [2.85715, 11], range ~[-1.484e-9,1.505e-9]:
* |log(x*erfc(x)) + x**2 + 0.5625 - r(x)/s(x)| < 2**-30
*/
-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 */
-/*
- * 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.86494310e-03F, /* -0x1.434124p-7 */
-rb1 = -6.25171244e-01F, /* -0x1.401672p-1 */
-rb2 = -6.16498327e+00F, /* -0x1.8a8f16p+2 */
-rb3 = -1.66696873e+01F, /* -0x1.0ab70ap+4 */
-rb4 = -9.53764343e+00F, /* -0x1.313460p+3 */
-sb1 = 1.26884899e+01F, /* 0x1.96081cp+3 */
-sb2 = 4.51839523e+01F, /* 0x1.6978bcp+5 */
-sb3 = 4.72810211e+01F, /* 0x1.7a3f88p+5 */
-sb4 = 8.93033314e+00F; /* 0x1.1dc54ap+3 */
+rb0 = -9.86496918e-3, /* 0xbc21a0ae */
+rb1 = -5.48049808e-1, /* 0xbf0c4cfe */
+rb2 = -1.84115684, /* 0xbfebab07 */
+sb1 = 4.87132740, /* 0x409be1ea */
+sb2 = 3.04982710, /* 0x4043305e */
+sb3 = -7.61900663e-1; /* 0xbf430bec */
float
erff(float x)
@@ -105,7 +97,7 @@
if(ix < 0x3fa00000) { /* 0.84375 <= |x| < 1.25 */
s = fabsf(x)-one;
P = pa0+s*(pa1+s*(pa2+s*pa3));
- Q = one+s*(qa1+s*(qa2+s*(qa3+s*qa4)));
+ Q = one+s*(qa1+s*(qa2+s*qa3));
if(hx>=0) return erx + P/Q; else return -erx - P/Q;
}
if (ix >= 0x40800000) { /* inf>|x|>=4 */
@@ -113,12 +105,12 @@
}
x = fabsf(x);
s = one/(x*x);
- if(ix< 0x4036DB6E) { /* |x| < 1/0.35 */
+ if(ix< 0x4036db8c) { /* |x| < 2.85715 ~ 1/0.35 */
R=ra0+s*(ra1+s*(ra2+s*ra3));
- S=one+s*(sa1+s*(sa2+s*(sa3+s*sa4)));
- } else { /* |x| >= 1/0.35 */
- R=rb0+s*(rb1+s*(rb2+s*(rb3+s*rb4)));
- S=one+s*(sb1+s*(sb2+s*(sb3+s*sb4)));
+ S=one+s*(sa1+s*(sa2+s*sa3));
+ } else { /* |x| >= 2.85715 ~ 1/0.35 */
+ R=rb0+s*(rb1+s*rb2);
+ S=one+s*(sb1+s*(sb2+s*sb3));
}
SET_FLOAT_WORD(z,hx&0xffffe000);
r = expf(-z*z-0.5625F)*expf((z-x)*(z+x)+R/S);
@@ -155,7 +147,7 @@
if(ix < 0x3fa00000) { /* 0.84375 <= |x| < 1.25 */
s = fabsf(x)-one;
P = pa0+s*(pa1+s*(pa2+s*pa3));
- Q = one+s*(qa1+s*(qa2+s*(qa3+s*qa4)));
+ Q = one+s*(qa1+s*(qa2+s*qa3));
if(hx>=0) {
z = one-erx; return z - P/Q;
} else {
@@ -165,13 +157,13 @@
if (ix < 0x41300000) { /* |x|<11 */
x = fabsf(x);
s = one/(x*x);
- if(ix< 0x4036DB6D) { /* |x| < 1/.35 ~ 2.857143*/
+ if(ix< 0x4036db8c) { /* |x| < 2.85715 ~ 1/.35 */
R=ra0+s*(ra1+s*(ra2+s*ra3));
- S=one+s*(sa1+s*(sa2+s*(sa3+s*sa4)));
- } else { /* |x| >= 1/.35 ~ 2.857143 */
+ S=one+s*(sa1+s*(sa2+s*sa3));
+ } else { /* |x| >= 2.85715 ~ 1/.35 */
if(hx<0&&ix>=0x40a00000) return two-tiny;/* x < -5 */
- R=rb0+s*(rb1+s*(rb2+s*(rb3+s*rb4)));
- S=one+s*(sb1+s*(sb2+s*(sb3+s*sb4)));
+ R=rb0+s*(rb1+s*rb2);
+ S=one+s*(sb1+s*(sb2+s*sb3));
}
SET_FLOAT_WORD(z,hx&0xffffe000);
r = expf(-z*z-0.5625F)*expf((z-x)*(z+x)+R/S);
Index: ld80/s_erfl.c
===================================================================
--- ld80/s_erfl.c (revision 0)
+++ ld80/s_erfl.c (working copy)
@@ -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 these: */
+static const volatile long double tiny = 0x1p-10000L;
+
+static const double
+half = 0.5,
+one = 1,
+two = 2;
+
+/*
+ * Domain [0, 0.84375], range ~[-1.573e-22, 1.524e-22]:
+ * |(erf(x) - x)/x - pp(x)/qq(x)| < 2**-73
+ *
+ * Note the domain [0, 0.84375] contains the subdomain [0,0x1p-34]
+ * where simple polynomial approximations are used. The error estimate
+ * does not apply to the subdomain.
+ */
+static const union IEEEl2bits
+efxu = LD80C(0x8375d410a6db446c, -3, 1.28379167095512573902e-1L),
+efx8u = LD80C(0x8375d410a6db446c, 0, 1.02703333676410059122L),
+pp0u = LD80C(0x8375d410a6db446c, -3, 1.28379167095512573902e-1L),
+pp1u = LD80C(0xa46c3b45612b6a1e, 32766, -3.21138241023240692798e-1L),
+pp2u = LD80C(0x9b3719d3c4ac4895, 32763, -3.78943451762934832069e-2L),
+pp3u = LD80C(0x804ca11f3c7b884e, 32761, -7.83076986858308357277e-3L),
+pp4u = LD80C(0x9f4ea36c9af744be, 32756, -3.03854334175109160953e-4L),
+pp5u = LD80C(0x9ec5abb84f9e646d, 32752, -1.89271153794195656738e-5L),
+qq1u = LD80C(0xdb4d8f02145058e4, -2, 4.28326100353076006451e-1L),
+qq2u = LD80C(0xa5783a50a8d83c35, -4, 8.07957225671139745698e-2L),
+qq3u = LD80C(0x8b8a0942511cd01a, -7, 8.51679710582846710136e-3L),
+qq4u = LD80C(0x873879d9502010a8, -11, 5.15825688424454091100e-4L),
+qq5u = LD80C(0x834af3bb90a69cc0, -16, 1.56513192135679610005e-5L),
+qq6u = LD80C(0xf59f11a3404946fb, -24, 1.14376360028953741942e-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)
+/*
+ * Domain [0.84375, 1.25], range ~[-8.396e-22, 8.417e-22]:
+ * |(erf(x) - erx) - pax)/qa(x)| < 2**-71
+ */
+static const union IEEEl2bits
+erxu = LD80C(0xd7bb3d0000000000, -1, 8.42700779438018798828e-1L),
+pa0u = LD80C(0xe8211158d9ff1d17, -27, 1.35116960705129358047e-8L),
+pa1u = LD80C(0xd488f89f5a81a191, -2, 4.15107507183395851624e-1L),
+pa2u = LD80C(0xebe046a54d22eaa5, 32764, -1.15173866195293202610e-1L),
+pa3u = LD80C(0xc7e7f3955bcda99b, -4, 9.76103811221723127828e-2L),
+pa4u = LD80C(0x99ef8f412193ce21, -5, 3.75819774509169158152e-2L),
+pa5u = LD80C(0xaa117d1eca526d19, 32760, -5.19007310284054296686e-3L),
+pa6u = LD80C(0x872c6dc079ba8e90, -8, 4.12516936280645110337e-3L),
+pa7u = LD80C(0xdb1858d6cf405001, -13, 2.08945375039182583868e-4L),
+qa1u = LD80C(0xb8f8ac7beb11603e, -1, 7.22544460538363523291e-1L),
+qa2u = LD80C(0x9fd5ca8e1b934583, -1, 6.24355945295588098283e-1L),
+qa3u = LD80C(0x9d6078f0c864ecd2, -2, 3.07376651184492511034e-1L),
+qa4u = LD80C(0x883087d5c1cfe6f2, -3, 1.32997629567810717548e-1L),
+qa5u = LD80C(0xbb4f5c04cb521d1a, -5, 4.57299799424686600875e-2L),
+qa6u = LD80C(0xa706ad813a9b5ab9, -7, 1.01944631941621374419e-2L),
+qa7u = LD80C(0x90985f65f597dbcf, -9, 2.20634774264849368479e-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)
+/*
+ * Domain [1.25,2.85715], range ~[-2.469e-22,2.475e-22]:
+ * |log(x*erfc(x)) + x**2 + 0.5625 - ra(x)/sa(x)| < 2**-71
+ */
+static const union IEEEl2bits
+ra0u = LD80C(0xa1a091e11c233978, 32761, -9.86494298962466012389e-3L),
+ra1u = LD80C(0xc2aa33978c3ef213, 32767, -7.60409569262488622607e-1L),
+ra2u = LD80C(0xf2b156cdd7197a15, 32771, -1.51682956734308320894e1L),
+ra3u = LD80C(0x811f905104c4f471, 32775, -1.29123295844711678773e2L),
+ra4u = LD80C(0x87097aae3b0dd531, 32777, -5.40148112828893618931e2L),
+ra5u = LD80C(0x8fbad259dd4823ef, 32778, -1.14983817761630154830e3L),
+ra6u = LD80C(0x97d998d331edd1a6, 32778, -1.21479990539314035369e3L),
+ra7u = LD80C(0x923a52bb507131c5, 32777, -5.84911299542004396834e2L),
+ra8u = LD80C(0xd4ee8a37aee580b4, 32774, -1.06465898265933675587e2L),
+ra9u = LD80C(0x92b266b05a914543, 32770, -4.58427748149964800614L),
+sa1u = LD80C(0xd32e06f57a5beb78, 4, 2.63974742105451605084e1L),
+sa2u = LD80C(0x838103b3fd116e2e, 8, 2.63007925509041215972e2L),
+sa3u = LD80C(0x9f5a7c5381e040e7, 10, 1.27482767653814701425e3L),
+sa4u = LD80C(0xca50d95a53553c99, 11, 3.23705306465675292871e3L),
+sa5u = LD80C(0x874bbb386365a447, 12, 4.32946641614583907964e3L),
+sa6u = LD80C(0xb673625fae9ceb9b, 11, 2.91921151703079406192e3L),
+sa7u = LD80C(0xdf7599ee57dd0b93, 9, 8.93837520204357297959e2L),
+sa8u = LD80C(0xc76f0ea852bba805, 6, 9.97169087029449912399e1L),
+sa9u = LD80C(0x92523c9ccba4ec38, 1, 2.28626933395442112455L);
+#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,11], range ~[-8.323e-22,8.390e-22]:
+ * |log(x*erfc(x)) + x**2 + 0.5625 - rb(x)/sb(x)| < 2**-70
+ */
+static const union IEEEl2bits
+rb0u = LD80C(0xa1a091cf43aed35d, 32761, -9.86494292470301435071e-3L),
+rb1u = LD80C(0xd17dedaf18eda123, 32767, -8.18327765701659009211e-1L),
+rb2u = LD80C(0x9a04d7d229bec743, 32772, -1.92523647707929260794e1L),
+rb3u = LD80C(0xbf4b6037d8d723c7, 32775, -1.91294436922492550143e2L),
+rb4u = LD80C(0xdcc0f0c2a30e05bc, 32777, -8.83014694842560938559e2L),
+rb5u = LD80C(0xe5fb64aebeda2e08, 32778, -1.83985604035640557274e3L),
+rb6u = LD80C(0xb946be95d55f5988, 32778, -1.48221076480554874966e3L),
+rb7u = LD80C(0x941fb05d9911fb04, 32776, -2.96247569751492811574e2L),
+sb1u = LD80C(0x81130861eae89202, 5, 3.22685866641642338307e1L),
+sb2u = LD80C(0xbdb73fddefb0a7a7, 8, 3.79431636564289246766e2L),
+sb3u = LD80C(0x80024bbf0b6ea3b3, 11, 2.04814349274125764278e3L),
+sb4u = LD80C(0xa444b434a0a8549b, 12, 5.25658799100412086380e3L),
+sb5u = LD80C(0xbb6ff428247dfadc, 12, 5.99799421719083271221e3L),
+sb6u = LD80C(0x9b76cbd808a6f2e2, 11, 2.48742476657275909036e3L),
+sb7u = LD80C(0xd3e4fe61373ce96a, 7, 2.11894506526933502549e2L);
+#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.880e-22,4.903e-22]:
+ * |log(x*erfc(x)) + x**2 + 0.5625 - rc(x)/sc(x)| < 2**-71
+ */
+static const union IEEEl2bits
+rc0u = LD80C(0xa1a091cf437a17ad, 32761, -9.86494292470008707260e-3L),
+rc1u = LD80C(0xbe73398d6aa320bf, 32767, -7.43945691133046437622e-1L),
+rc2u = LD80C(0xdb0b9c6c7c98cf95, 32771, -1.36903347242128353893e1L),
+rc3u = LD80C(0xb5cc0cecf47b4b20, 32774, -9.08985361145438643948e1L),
+rc4u = LD80C(0xd749186e6457c66b, 32775, -2.15285529040811320620e2L),
+rc5u = LD80C(0xfdd1925f754a9031, 32774, -1.26909319861476689424e2L),
+sc1u = LD80C(0xc5d40fb999d75384, 4, 2.47285456180642684368e1L),
+sc2u = LD80C(0xc5c6fe0e300c822b, 7, 1.97777314078074149137e2L),
+sc3u = LD80C(0x961dc355db85da06, 9, 6.00465047325494566954e2L),
+sc4u = LD80C(0x997ad30640d6bb4d, 9, 6.13919129908866058554e2L),
+sc5u = LD80C(0xf7836d4d36940f8a, 6, 1.23756693280136487692e2L);
+#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)
+{
+ int32_t i;
+ 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) { /* x is +-Inf or NaN */
+ i = (hx>>15)<<1;
+ return (1-i)+one/x;
+ }
+
+ 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) { /* |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 {
+ 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) { /* erfc(nan)=nan */
+ /* erfc(+-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: ld80/s_erfl.c
___________________________________________________________________
Added: svn:eol-style
## -0,0 +1 ##
+native
\ No newline at end of property
Added: svn:mime-type
## -0,0 +1 ##
+text/plain
\ No newline at end of property
Added: svn:keywords
## -0,0 +1 ##
+FreeBSD=%H
\ No newline at end of property
Index: ld128/s_erfl.c
===================================================================
--- ld128/s_erfl.c (revision 0)
+++ ld128/s_erfl.c (working copy)
@@ -0,0 +1,327 @@
+/* @(#)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 long double
+half = 0.5L,
+one = 1,
+two = 2;
+/*
+ * Domain [0, 0.84375], range ~[-2.076e-38, 2.074e-38]:
+ * |(erf(x) - x)/x - pp(x)/qq(x)| < 2**-126
+ *
+ * Note the domain [0, 0.84375] contains the subdomain [0,0x1p-34]
+ * where simple polynomial approximations are used. The error estimate
+ * does not apply to the subdomain.
+ */
+static const long double
+efx = 1.28379167095512573896158903121545167e-1L, /* 0xecbff6a7, 0x481dd788, 0xb64d21a8, 0xeb06fc3f */
+efx8 = 1.02703333676410059116927122497236133L, /* 0xecbff6a7, 0x481dd788, 0xb64d21a8, 0xeb06ff3f */
+pp0 = 1.28379167095512573896158903121545167e-1L, /* 0x3ffc06eb, 0xa8214db6, 0x88d71d48, 0xa7f6bfec */
+pp1 = -3.14930453779220199897729762882897733e-1L, /* 0xbffd427d, 0x20fdfc19, 0x5395ba59, 0x26b231dc */
+pp2 = -5.27541911766414675560926940708709654e-2L, /* 0xbffab029, 0x8eba9454, 0xd7a92b48, 0x80023341 */
+pp3 = -1.13207381741395209251907161443436107e-2L, /* 0xbff872f5, 0x3c1cbcad, 0x90c8a6f3, 0x2b227613 */
+pp4 = -9.18720900388670035443046669196822569e-4L, /* 0xbff4e1ac, 0xa1c543f3, 0x522d6f94, 0x4ebd756e */
+pp5 = -7.87594031651048602426379053816416599e-5L, /* 0xbff14a57, 0x43e4e927, 0xd8815b94, 0xc2cd92da */
+pp6 = -3.42417294394664810489445209450592648e-6L, /* 0xbfeccb95, 0xafbd3483, 0x5b58d3e7, 0x712e6f9f */
+pp7 = -1.37340867367364131528690800390865353e-7L, /* 0xbfe826ef, 0xf0b8c986, 0xe445805b, 0x60654151 */
+pp8 = -2.71188885190546726082043702148966261e-9L, /* 0xbfe274b8, 0x1b3eee47, 0x1739ceb3, 0x2217b257 */
+pp9 = -3.38021083561888762328182153549042469e-11L, /* 0xbfdc2953, 0x94cffcc1, 0x249ea8a2, 0x2152a05a */
+qq1 = 4.76681198648751759882307858251781949e-1L, /* 0x3ffde81f, 0x1dbb4203, 0xf3a919ed, 0x9d4721ed */
+qq2 = 1.06717237624212542329851335448543918e-1L, /* 0x3ffbb51d, 0x22583f4c, 0x480ed405, 0x2b999881 */
+qq3 = 1.47756708414471579803258620664775360e-2L, /* 0x3ff8e42b, 0x4f852e54, 0x30ab6e4d, 0x5ffcefe3 */
+qq4 = 1.39951717291319836305643111144723543e-3L, /* 0x3ff56ee0, 0x01f38f9f, 0x5ee018dc, 0x355c6b48 */
+qq5 = 9.44415276581511130343426049763432505e-5L, /* 0x3ff18c1d, 0xd1738632, 0xf03254be, 0x6a6e2071 */
+qq6 = 4.56270013341289983962750459103111817e-6L, /* 0x3fed3232, 0x97672910, 0xc334df35, 0xfe01507b */
+qq7 = 1.53049464020884224321556705278110797e-7L, /* 0x3fe848ab, 0xd52911d5, 0xabe051cc, 0x8abf1cd2 */
+qq8 = 3.25623515051818985116653695841960050e-9L, /* 0x3fe2bf88, 0x96d22da3, 0x0207a0ee, 0x17463bb2 */
+qq9 = 3.37510289528866336760710014443284362e-11L; /* 0x3fdc28e0, 0x8f8a8bcc, 0x6cc0a227, 0x6fce7fe7 */
+/*
+ * Domain [0.84375, 1.25], range ~[-2.372e-36, 2.373e-36]:
+ * |(erf(x) - erx) - pa(x)/qa(x)| < 2**-120
+ */
+static const long double
+erx = 8.42700792949714894142232424201210961e-1L, /* 0x3ffeaf76, 0x7a741088, 0xb0000000, 0x00000000 */
+pa0 = -2.48010117891186017025877408964183854e-17L, /* 0xbfc7c97f, 0x77812279, 0x6c87b01b, 0xec38d68e */
+pa1 = 4.15107497420594680851927951552407784e-1L, /* 0x3ffda911, 0xf096fbc2, 0x5562ff1b, 0x3b544e2b */
+pa2 = -3.87083915464088350525375137436107590e-2L, /* 0xbffa3d19, 0x6178b90c, 0x2c86df33, 0x5319973e */
+pa3 = 4.44954959590143486203056389208894541e-2L, /* 0x3ffa6c81, 0xd17ed314, 0xdcacac2f, 0xc58a3516 */
+pa4 = 8.05166892893185510859535256028511853e-2L, /* 0x3ffb49cb, 0xde347a21, 0x08774372, 0xabfb5f6a */
+pa5 = -1.03051186773465597646984141058502267e-2L, /* 0xbff851ad, 0x99d9ad55, 0xa4289298, 0xc77d3e5a */
+pa6 = 5.72604838819243779197220502236739196e-3L, /* 0x3ff77743, 0x2690068e, 0x45e8dce5, 0x9ca78cd1 */
+pa7 = 1.22625964697528485542108776707051120e-3L, /* 0x3ff54174, 0xe4521657, 0x9b318890, 0x918d9b04 */
+pa8 = 5.39388406655659100369616014543606682e-4L, /* 0x3ff41acb, 0x7c880b32, 0x296affe7, 0xed6f9a01 */
+pa9 = -1.98335880680615762490364057251014162e-4L, /* 0xbff29ff0, 0xc3e06ed7, 0x3227ca50, 0x5db44d09 */
+pa10 = 6.20889395425360375062078472544482827e-5L, /* 0x3ff1046b, 0x7dbeee60, 0xd2792787, 0x19c9aa7b */
+pa11 = -5.40647780776317096102550495545612307e-6L, /* 0xbfed6ad2, 0x94dc1b75, 0xedf4970c, 0x9f2fcd17 */
+qa1 = 9.06750921660206127928124875031300980e-1L, /* 0x3ffed041, 0xa8244c00, 0xc4182e08, 0x27bd05f1 */
+qa2 = 6.80607885383663092673769524543423470e-1L, /* 0x3ffe5c78, 0xa3023ebe, 0xcb5d7d74, 0x6465c592 */
+qa3 = 4.05656787299304601494351680677926251e-1L, /* 0x3ffd9f64, 0x7e2b675d, 0x307e6db8, 0x612204de */
+qa4 = 1.69503824422159751714827838146178306e-1L, /* 0x3ffc5b24, 0xd2338554, 0x5192bb8c, 0x0e186eff */
+qa5 = 7.46584241233528910281045656042597495e-2L, /* 0x3ffb31cd, 0x081fb0f0, 0x878f1e11, 0xce245503 */
+qa6 = 2.03534629098115442535903549876232764e-2L, /* 0x3ff94d78, 0x9c63b619, 0xed805e2c, 0x9467d788 */
+qa7 = 6.96756479371300019285814024610131574e-3L, /* 0x3ff7c8a0, 0x56ebf85a, 0x78d461d0, 0x5ac7c4fe */
+qa8 = 1.13424099180080891902336864678392922e-3L, /* 0x3ff52955, 0x9fdcbd9d, 0xc035851a, 0xfa944ed7 */
+qa9 = 3.15259605619527194057660981542654259e-4L, /* 0x3ff34a92, 0xdb225912, 0x7b9010de, 0x3fe39023 */
+qa10 = 1.18798605497472051761873112336132064e-5L, /* 0x3fee8e9f, 0x399f21c2, 0x4e53cc38, 0x5e6e1c31 */
+qa11 = 4.65215947978206806357776624485104779e-6L, /* 0x3fed3833, 0x7dc4e756, 0x4c73bce3, 0x43a896cc */
+qa12 = -1.02073578608782110388686097137679971e-6L; /* 0xbfeb1200, 0x6dd9de6c, 0x647bcf5b, 0xec577463 */
+/*
+ * Domain [1.25,2.85715], range ~[-2.932e-37,2.932e-37]:
+ * |log(x*erfc(x)) + x**2 + 0.5625 - ra(x)/sa(x)| < 2**-121
+ */
+static const long double
+ra0 = -9.86494292470069026041027804612925719e-3L, /* 0xbff84341, 0x239e8709, 0xeac6f011, 0x288a805d */
+ra1 = -1.13579789657313713897975924115138941L, /* 0xbfff22c3, 0xa6a4a5af, 0xd9327434, 0x30d5c80e */
+ra2 = -4.89737449616464986815099118405280362e1L, /* 0xc00487ca, 0x3acc6754, 0xf6b0cba0, 0x968544a8 */
+ra3 = -1.10763525821733005976688912399857907e3L, /* 0xc00914e8, 0xa81214fc, 0xa9daec12, 0x04a27e27 */
+ra4 = -1.49985906868704004330240577932725368e4L, /* 0xc00cd4b4, 0xb9ba09b4, 0x5f2150b4, 0xe9d43451 */
+ra5 = -1.29799869868909118595022585186399503e5L, /* 0xc00ffb07, 0xdeafba94, 0x78b8464d, 0x2243a38a */
+ra6 = -7.42786170818211875353524648807020262e5L, /* 0xc0126ab0, 0x457757c1, 0x321b6122, 0xfbf31ed7 */
+ra7 = -2.85617550342537589290678185465422207e6L, /* 0xc0145ca7, 0x7c0703e2, 0x2b7e096e, 0x6f5a3f5d */
+ra8 = -7.40614267250515630140845919774792706e6L, /* 0xc015c408, 0xfab0a531, 0x12d2dc13, 0xa2b8704b */
+ra9 = -1.28641241327410554128542575661432164e7L, /* 0xc0168894, 0xf843f6a2, 0xb7ab6219, 0x744c300c */
+ra10 = -1.47182276924276357440551219889149080e7L, /* 0xc016c12a, 0x276285e0, 0x04bc1105, 0xe4e10e3c */
+ra11 = -1.07799899573272478359709008914120358e7L, /* 0xc01648fa, 0xabea26cc, 0x0a0ce4a0, 0xb7ec8778 */
+ra12 = -4.83480214872222791837505030943085341e6L, /* 0xc0152717, 0xc8984aa3, 0xc46405f6, 0xa0f6f640 */
+ra13 = -1.23938051019472227666411457770036806e6L, /* 0xc0132e95, 0x4829c1f0, 0xec523a33, 0x776ccfd3 */
+ra14 = -1.62262886680556039283486133222303718e5L, /* 0xc0103ceb, 0x717ebf9b, 0x180a6a70, 0x34417b99 */
+ra15 = -8.82734130841920962111981230171892874e3L, /* 0xc00c13da, 0xbaffe892, 0xd6975b67, 0x05b7a992 */
+ra16 = -1.22568436072739969460925172985178212e2L, /* 0xc005ea46, 0x141b1923, 0xaea7e61c, 0x7207a8df */
+sa1 = 6.44502356908054050802853406380055114e1L, /* 0x400501cd, 0x0a95be01, 0x360afc70, 0x2d5f256b */
+sa2 = 1.76114826500718135746343166847984357e3L, /* 0x4009b849, 0x7d2c833f, 0x0c758b06, 0x07d73e83 */
+sa3 = 2.69439831921639091939469234421147772e4L, /* 0x400da4ff, 0xeec9ed36, 0xb19f5f0d, 0x6e5ff0b1 */
+sa4 = 2.56815637017014469128294345646768041e5L, /* 0x4010f597, 0xd189c606, 0x11fac05c, 0x9dda4b30 */
+sa5 = 1.60638682704348688266145179445842807e6L, /* 0x4013882f, 0x2d3b91f3, 0x887e4dcf, 0x8f990fab */
+sa6 = 6.76918282227476395975704195544831096e6L, /* 0x40159d28, 0x7b4a0265, 0x4e221317, 0x8f19f36b */
+sa7 = 1.94280485423922835890016575512146130e7L, /* 0x40172872, 0xd08ada38, 0x7f9e15c4, 0x1fbf1deb */
+sa8 = 3.79740375883744863135748410413555680e7L, /* 0x401821b8, 0x0acb4fda, 0xec4240d5, 0x38ffcf2c */
+sa9 = 5.00608194437408776705107810210078804e7L, /* 0x40187def, 0x09b8cc80, 0x46bf2089, 0x9ec255ed */
+sa10 = 4.36435715828395309369777169953877366e7L, /* 0x40184cf9, 0x59ca9a7c, 0x5a185bbc, 0x2c480858 */
+sa11 = 2.43748298478043040399488401407107134e7L, /* 0x401773ee, 0x2dd909b3, 0xef427684, 0xb0127a72 */
+sa12 = 8.30614655253558095273324523461136470e6L, /* 0x4015faf7, 0x8a35cbe3, 0x2845fd7e, 0x1c2d314d */
+sa13 = 1.60136068063256826643328510848675710e6L, /* 0x401386f5, 0x0ae3def9, 0xd4bfc49f, 0x9a3fcbe4 */
+sa14 = 1.54229628674911158232902134264852520e5L, /* 0x40102d3a, 0xd0786b63, 0x9ec2ab01, 0xdffa186c */
+sa15 = 5.87842946700351712242359416526459677e3L, /* 0x400b6f66, 0xdf18caed, 0x13140abd, 0x0b881c24 */
+sa16 = 4.97176688590085185068333226860548624e1L; /* 0x40048dbd, 0xc92bb664, 0xce7d48ce, 0x6380778e */
+/*
+ * 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)
+{
+ int32_t i;
+ 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) { /* x is +-Inf or NaN */
+ i = (hx>>15)<<1;
+ return (1-i)+one/x;
+ }
+
+ ax = fabsl(x);
+ if(ax < 0.84375) {
+ if(ax < 0x1p-34L) {
+ 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) {
+ 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) { /* erfc(nan)=nan */
+ /* erfc(+-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: ld128/s_erfl.c
___________________________________________________________________
Added: svn:eol-style
## -0,0 +1 ##
+native
\ No newline at end of property
Added: svn:mime-type
## -0,0 +1 ##
+text/plain
\ No newline at end of property
Added: svn:keywords
## -0,0 +1 ##
+FreeBSD=%H
\ No newline at end of property
Index: src/imprecise.c
===================================================================
--- src/imprecise.c (revision 262934)
+++ src/imprecise.c (working copy)
@@ -60,7 +60,5 @@
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);
More information about the freebsd-numerics
mailing list