Date: Mon, 2 Dec 2013 13:42:40 -0800 From: Steve Kargl <sgk@troutmask.apl.washington.edu> To: freebsd-numerics@freebsd.org Subject: erff, erfl, and erfcl patch Message-ID: <20131202214240.GA44209@troutmask.apl.washington.edu>
next in thread | raw e-mail | index | archive | help
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.
* msun/Makefile:
. Add s_erfl.c.
. Remove nearby extra space.
. Add MLINKS for erfl.3 and erfcl.3.
* msun/Symbol.map:
. Sort erfl and erfcl into proper order.
* msun/ld128/s_erfl.c:
. Implementations for erfl and erfcl in IEEE 128-bit long double format.
* msun/ld80/s_erfl.c:
. Implementations for erfl and erfcl in Intel 80-bit long double format.
* msun/imprecise.c:
. Remove hack that maps erfl to erf and erfcl to erfc.
* msun/src/s_erf.c:
. Fix whitespace issues.
. Add weak references for erfl and erfcl on targets where long double
has 53 bits of precisions.
* msun/src/s_erff.c:
. Consistently use lower case in hex constants.
. Update the rational approximations.
. Fix descriptions of the the domain and range.
. Remove leading zeros in exponents, e.g., 1.28379166e-01 becomes
1.28379166e-1.
--
Steve
Index: Makefile
===================================================================
--- Makefile (revision 258855)
+++ Makefile (working copy)
@@ -99,9 +99,9 @@
e_fmodl.c e_hypotl.c e_remainderl.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_nexttoward.c s_remquol.c s_rintl.c s_roundl.c s_scalbnl.c \
s_sinl.c s_tanl.c s_truncl.c w_cabsl.c
.endif
@@ -162,7 +162,8 @@
MLINKS+=cos.3 cosf.3 cos.3 cosl.3
MLINKS+=cosh.3 coshf.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 258855)
+++ Symbol.map (working copy)
@@ -264,6 +264,8 @@
ctanf;
ctanh;
ctanhf;
+ erfcl;
+ erfl;
expl;
expm1l;
log10l;
@@ -272,8 +274,6 @@
logl;
/* Implemented as weak aliases for imprecise versions */
coshl;
- erfcl;
- erfl;
lgammal;
powl;
sinhl;
Index: ld128/s_erfl.c
===================================================================
--- ld128/s_erfl.c (revision 0)
+++ ld128/s_erfl.c (working copy)
@@ -0,0 +1,347 @@
+/* @(#)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"
+
+/* Sloppy limits. */
+#define LIM0 2.85715 /* ~ 1. / 0.35 */
+#define LIM1 9 /* x**2 = - log(sqrt(pi) * x) - (1 - p) * log(2) */
+#define LIM2 108 /* x**2 = - log(sqrt(pi) * x) - emin */
+
+/* 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 - p(x)/q(x)| < 2**-126
+ */
+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) - p(x)/q(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 - r(x)/s(x)| < 2**-121
+ */
+static const long double
+ra0 = -9.86494292470069026041027804612925719e-3L, /* 0xbff84341, 0x239e8709, 0xeac6f011, 0x288a805d */
+ra1 = -1.13579789657313713897975924115138941L, /* 0xbfff22c3, 0xa6a4a5af, 0xd9327434, 0x30d5c80e */
+ra2 = -4.89737449616464986815099118405280362e+1L, /* 0xc00487ca, 0x3acc6754, 0xf6b0cba0, 0x968544a8 */
+ra3 = -1.10763525821733005976688912399857907e+3L, /* 0xc00914e8, 0xa81214fc, 0xa9daec12, 0x04a27e27 */
+ra4 = -1.49985906868704004330240577932725368e+4L, /* 0xc00cd4b4, 0xb9ba09b4, 0x5f2150b4, 0xe9d43451 */
+ra5 = -1.29799869868909118595022585186399503e+5L, /* 0xc00ffb07, 0xdeafba94, 0x78b8464d, 0x2243a38a */
+ra6 = -7.42786170818211875353524648807020262e+5L, /* 0xc0126ab0, 0x457757c1, 0x321b6122, 0xfbf31ed7 */
+ra7 = -2.85617550342537589290678185465422207e+6L, /* 0xc0145ca7, 0x7c0703e2, 0x2b7e096e, 0x6f5a3f5d */
+ra8 = -7.40614267250515630140845919774792706e+6L, /* 0xc015c408, 0xfab0a531, 0x12d2dc13, 0xa2b8704b */
+ra9 = -1.28641241327410554128542575661432164e+7L, /* 0xc0168894, 0xf843f6a2, 0xb7ab6219, 0x744c300c */
+ra10= -1.47182276924276357440551219889149080e+7L, /* 0xc016c12a, 0x276285e0, 0x04bc1105, 0xe4e10e3c */
+ra11= -1.07799899573272478359709008914120358e+7L, /* 0xc01648fa, 0xabea26cc, 0x0a0ce4a0, 0xb7ec8778 */
+ra12= -4.83480214872222791837505030943085341e+6L, /* 0xc0152717, 0xc8984aa3, 0xc46405f6, 0xa0f6f640 */
+ra13= -1.23938051019472227666411457770036806e+6L, /* 0xc0132e95, 0x4829c1f0, 0xec523a33, 0x776ccfd3 */
+ra14= -1.62262886680556039283486133222303718e+5L, /* 0xc0103ceb, 0x717ebf9b, 0x180a6a70, 0x34417b99 */
+ra15= -8.82734130841920962111981230171892874e+3L, /* 0xc00c13da, 0xbaffe892, 0xd6975b67, 0x05b7a992 */
+ra16= -1.22568436072739969460925172985178212e+2L, /* 0xc005ea46, 0x141b1923, 0xaea7e61c, 0x7207a8df */
+sa1 = 6.44502356908054050802853406380055114e+1L, /* 0x400501cd, 0x0a95be01, 0x360afc70, 0x2d5f256b */
+sa2 = 1.76114826500718135746343166847984357e+3L, /* 0x4009b849, 0x7d2c833f, 0x0c758b06, 0x07d73e83 */
+sa3 = 2.69439831921639091939469234421147772e+4L, /* 0x400da4ff, 0xeec9ed36, 0xb19f5f0d, 0x6e5ff0b1 */
+sa4 = 2.56815637017014469128294345646768041e+5L, /* 0x4010f597, 0xd189c606, 0x11fac05c, 0x9dda4b30 */
+sa5 = 1.60638682704348688266145179445842807e+6L, /* 0x4013882f, 0x2d3b91f3, 0x887e4dcf, 0x8f990fab */
+sa6 = 6.76918282227476395975704195544831096e+6L, /* 0x40159d28, 0x7b4a0265, 0x4e221317, 0x8f19f36b */
+sa7 = 1.94280485423922835890016575512146130e+7L, /* 0x40172872, 0xd08ada38, 0x7f9e15c4, 0x1fbf1deb */
+sa8 = 3.79740375883744863135748410413555680e+7L, /* 0x401821b8, 0x0acb4fda, 0xec4240d5, 0x38ffcf2c */
+sa9 = 5.00608194437408776705107810210078804e+7L, /* 0x40187def, 0x09b8cc80, 0x46bf2089, 0x9ec255ed */
+sa10= 4.36435715828395309369777169953877366e+7L, /* 0x40184cf9, 0x59ca9a7c, 0x5a185bbc, 0x2c480858 */
+sa11= 2.43748298478043040399488401407107134e+7L, /* 0x401773ee, 0x2dd909b3, 0xef427684, 0xb0127a72 */
+sa12= 8.30614655253558095273324523461136470e+6L, /* 0x4015faf7, 0x8a35cbe3, 0x2845fd7e, 0x1c2d314d */
+sa13= 1.60136068063256826643328510848675710e+6L, /* 0x401386f5, 0x0ae3def9, 0xd4bfc49f, 0x9a3fcbe4 */
+sa14= 1.54229628674911158232902134264852520e+5L, /* 0x40102d3a, 0xd0786b63, 0x9ec2ab01, 0xdffa186c */
+sa15= 5.87842946700351712242359416526459677e+3L, /* 0x400b6f66, 0xdf18caed, 0x13140abd, 0x0b881c24 */
+sa16= 4.97176688590085185068333226860548624e+1L; /* 0x40048dbd, 0xc92bb664, 0xce7d48ce, 0x6380778e */
+/*
+ * Domain [2.85715,9], range ~[-7.886e-37,7.918e-37]:
+ * |log(x*erfc(x)) + x**2 + 0.5625 - r(x)/s(x)| < 2**-120
+ */
+static const long double
+rb0 = -9.86494292470008707171371994479162369e-3L, /* 0xbff84341, 0x239e86f4, 0x2f57e561, 0xf4469360 */
+rb1 = -1.57047326624110727986326503729442830L, /* 0xbfff920a, 0x8935bf73, 0x8803b894, 0x4656482d */
+rb2 = -1.03228196364885474342132255440317065e+2L, /* 0xc0059ce9, 0xac4ed0ff, 0x2cff0ff7, 0x5e70d1ab */
+rb3 = -3.74000570653418227179358710865224376e+3L, /* 0xc00ad380, 0x2ebf7835, 0xf6b07ed2, 0x861242f7 */
+rb4 = -8.35435477739098044190860390632813956e+4L, /* 0xc00f4657, 0x8c3ae934, 0x3647d7b3, 0x80e76fb7 */
+rb5 = -1.21398672055223642118716640216747152e+6L, /* 0xc0132862, 0x2b8761c8, 0x27d18c0f, 0x137c9463 */
+rb6 = -1.17669175877248796101665344873273970e+7L, /* 0xc0166719, 0x0b2cea46, 0x81f14174, 0x11602ea5 */
+rb7 = -7.66108006086998253606773064264599615e+7L, /* 0xc019243f, 0x3c26f4f0, 0x1cc05241, 0x3b953728 */
+rb8 = -3.32547117558141845968704725353130804e+8L, /* 0xc01b3d24, 0x42d8ee26, 0x24ef6f3b, 0x604a8c65 */
+rb9 = -9.41561252426350696802167711221739746e+8L, /* 0xc01cc0f8, 0xad23692a, 0x8ddb2310, 0xe9937145 */
+rb10= -1.67157110805390944549427329626281063e+9L, /* 0xc01d8e88, 0x9a903734, 0x09a55fa3, 0xd205c903 */
+rb11= -1.74339631004410841337645931421427373e+9L, /* 0xc01d9fa8, 0x77582d2a, 0xc183b8ab, 0x7e00cb05 */
+rb12= -9.57655233596934915727573141357471703e+8L, /* 0xc01cc8a5, 0x460cc685, 0xd0271fa0, 0x6a70e3da */
+rb13= -2.26320062731339353035254704082495066e+8L, /* 0xc01aafab, 0xd7d76721, 0xc9720e11, 0x6a8bd489 */
+rb14= -1.42777302996263256686002973851837039e+7L, /* 0xc016b3b8, 0xc499689f, 0x2b88d965, 0xc32414f9 */
+sb1 = 1.08512869705594540211033733976348506e+2L, /* 0x4005b20d, 0x2db7528d, 0x00d20dcb, 0x858f6191 */
+sb2 = 5.02757713761390460534494530537572834e+3L, /* 0x400b3a39, 0x3bf4a690, 0x3025d28d, 0xfd40a891 */
+sb3 = 1.31019107205412870059331647078328430e+5L, /* 0x400fffcb, 0x1b71d05e, 0x3b28361d, 0x2a3c3690 */
+sb4 = 2.13021555152296846166736757455018030e+6L, /* 0x40140409, 0x3c6984df, 0xc4491d7c, 0xb04aa08d */
+sb5 = 2.26649105281820861953868568619768286e+7L, /* 0x401759d6, 0xce8736f0, 0xf28ad037, 0x2a901e0c */
+sb6 = 1.61071939490875921812318684143076081e+8L, /* 0x401a3338, 0x686fb541, 0x6bd27d06, 0x4f95c9ac */
+sb7 = 7.66895673844301852676056750497991966e+8L, /* 0x401c6daf, 0x31cec121, 0x54699126, 0x4bd9bf9e */
+sb8 = 2.41884450436101936436023058196042526e+9L, /* 0x401e2059, 0x46b0b8d7, 0x87b64cbf, 0x78bc296d */
+sb9 = 4.92403055884071695093305291535107666e+9L, /* 0x401f257e, 0xbe5ed739, 0x39e17346, 0xcadd2e55 */
+sb10= 6.18627786365587486459633615573786416e+9L, /* 0x401f70bb, 0x1be7a7e7, 0x6a45b5ae, 0x607c70f0 */
+sb11= 4.45898013426501378097430226324743199e+9L, /* 0x401f09c6, 0xa32643d7, 0xf1724620, 0x9ea46c32 */
+sb12= 1.63006115763329848117160344854224975e+9L, /* 0x401d84a3, 0x0996887f, 0x65a4f43b, 0x978c1d74 */
+sb13= 2.39216717012421697446304015847567721e+8L, /* 0x401ac845, 0x09a065c2, 0x30095da7, 0x9d72d6ae */
+sb14= 7.84837329009278694937250358810225609e+6L; /* 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.13742634438922591780742637728666162e+1L, /* 0xc004eafe, 0x7dd51cd8, 0x3c7c5928, 0x751e50cf */
+rc3 = -1.50455835478908280402912854338421517e+3L, /* 0xc0097823, 0xbc15b9ab, 0x3d60745c, 0x523e80a5 */
+rc4 = -2.04415631865861549920184039902945685e+4L, /* 0xc00d3f66, 0x40b3fc04, 0x5388f2ec, 0xb009e1f0 */
+rc5 = -1.57625662981714582753490610560037638e+5L, /* 0xc01033dc, 0xd4dc95b6, 0xfd4da93b, 0xf355b4a9 */
+rc6 = -6.73473451616752528402917538033283794e+5L, /* 0xc01248d8, 0x2e73a4f9, 0xcded49c5, 0xfa3bfeb7 */
+rc7 = -1.47433165421387483167186683764364857e+6L, /* 0xc01367f1, 0xba77a8f7, 0xcfdd0dbb, 0x25d554b3 */
+rc8 = -1.38811981807868828563794929997744139e+6L, /* 0xc01352e5, 0x7d16d9ad, 0xbbdcbf38, 0x38fbc5ea */
+rc9 = -3.59659700530831825640766479698155060e+5L, /* 0xc0115f3a, 0xecd57f45, 0x21f8ad6c, 0x910a5958 */
+sc1 = 7.72730753022908298637508998072635696e+1L, /* 0x40053517, 0xa10d52bc, 0xdabb55b6, 0xbd0328cd */
+sc2 = 2.36825757341694050500333261769082182e+3L, /* 0x400a2808, 0x3e0a9b42, 0x82977842, 0x9c5de29e */
+sc3 = 3.72210540173034735352888847134073099e+4L, /* 0x400e22ca, 0x1ba827ef, 0xac8390d7, 0x1fc39a41 */
+sc4 = 3.24136032646418336712461033591393412e+5L, /* 0x40113c8a, 0x0216e100, 0xc59d1e44, 0xf0e68d9d */
+sc5 = 1.57836135851134393802505823370009175e+6L, /* 0x40138157, 0x95bc7664, 0x17575961, 0xdbe58eeb */
+sc6 = 4.12881981392063738026679089714182355e+6L, /* 0x4014f801, 0x9e82e8d2, 0xb8b3a70e, 0xfd84185d */
+sc7 = 5.24438427289213488410596395361544142e+6L, /* 0x40154017, 0x81177109, 0x2aa6c3b0, 0x1f106625 */
+sc8 = 2.59909544563616121735963429710382149e+6L, /* 0x40143d45, 0xbb90a9b1, 0x12bf9390, 0xa827a700 */
+sc9 = 2.80930665169282501639651995082335693e+5L; /* 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);
+
+ /* erf(+-0) = +-0 */
+ if ((hx == 0 || hx == 0x8000) && (lx | llx) == 0)
+ return (x);
+
+ /* x is INF or NaN */
+ if (hx >= 0x7fff) {
+ if (hx == 0xffff)
+ return (-1);
+ if ((lx | llx) == 0)
+ return (1);
+ return (x * x);
+ }
+
+ ax = fabsl(x);
+ if(ax < 0.84375) {
+ if(ax < 0x1.p-34L) {
+ if (ax < 0x1.p-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 < LIM1) {
+ s = one/(ax*ax);
+ if(ax < LIM0) { /* |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);
+ } else {
+ if(x>=0) return (one-tiny); else return (tiny-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);
+
+ /* x is INF or NaN */
+ if (hx >= 0x7fff) {
+ if (hx == 0xffff)
+ return (2);
+ if ((lx | llx) == 0)
+ return (0);
+ return (x * x);
+ }
+
+ ax = fabsl(x);
+ if(ax < 0.84375) {
+ if(ax < 0x1.p-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.25) { /* x<1/4 */
+ return one-(x+x*y);
+ } else {
+ r = x*y;
+ r += (x-half);
+ return half - r;
+ }
+ }
+
+ 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) {
+ z = one-erx; return z - P/Q;
+ } else {
+ z = erx+P/Q; return one+z;
+ }
+ }
+
+ if (ax < LIM2) { /* |x| < 108 */
+ s = one/(ax*ax);
+ if(ax < LIM0) { /* |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 < LIM1) {
+ 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 < -LIM1) 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: ld80/s_erfl.c
===================================================================
--- ld80/s_erfl.c (revision 0)
+++ ld80/s_erfl.c (working copy)
@@ -0,0 +1,353 @@
+/* @(#)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"
+
+/* Sloppy limits. */
+#define LIM0 2.85715 /* ~ 1. / 0.35 */
+#define LIM1 7 /* x**2 = - log(sqrt(pi) * x) - (1 - p) * log(2) */
+#define LIM2 108 /* x**2 = - log(sqrt(pi) * x) - emin */
+
+/* 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 - p(x)/q(x)| < 2**-73
+ */
+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) - p(x)/q(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 - r(x)/s(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.51682956734308320894e+1L),
+ra3u = LD80C(0x811f905104c4f471, 32775, -1.29123295844711678773e+2L),
+ra4u = LD80C(0x87097aae3b0dd531, 32777, -5.40148112828893618931e+2L),
+ra5u = LD80C(0x8fbad259dd4823ef, 32778, -1.14983817761630154830e+3L),
+ra6u = LD80C(0x97d998d331edd1a6, 32778, -1.21479990539314035369e+3L),
+ra7u = LD80C(0x923a52bb507131c5, 32777, -5.84911299542004396834e+2L),
+ra8u = LD80C(0xd4ee8a37aee580b4, 32774, -1.06465898265933675587e+2L),
+ra9u = LD80C(0x92b266b05a914543, 32770, -4.58427748149964800614L),
+sa1u = LD80C(0xd32e06f57a5beb78, 4, 2.63974742105451605084e+1L),
+sa2u = LD80C(0x838103b3fd116e2e, 8, 2.63007925509041215972e+2L),
+sa3u = LD80C(0x9f5a7c5381e040e7, 10, 1.27482767653814701425e+3L),
+sa4u = LD80C(0xca50d95a53553c99, 11, 3.23705306465675292871e+3L),
+sa5u = LD80C(0x874bbb386365a447, 12, 4.32946641614583907964e+3L),
+sa6u = LD80C(0xb673625fae9ceb9b, 11, 2.91921151703079406192e+3L),
+sa7u = LD80C(0xdf7599ee57dd0b93, 9, 8.93837520204357297959e+2L),
+sa8u = LD80C(0xc76f0ea852bba805, 6, 9.97169087029449912399e+1L),
+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 - r(x)/s(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.92523647707929260794e+1L),
+rb3u = LD80C(0xbf4b6037d8d723c7, 32775, -1.91294436922492550143e+2L),
+rb4u = LD80C(0xdcc0f0c2a30e05bc, 32777, -8.83014694842560938559e+2L),
+rb5u = LD80C(0xe5fb64aebeda2e08, 32778, -1.83985604035640557274e+3L),
+rb6u = LD80C(0xb946be95d55f5988, 32778, -1.48221076480554874966e+3L),
+rb7u = LD80C(0x941fb05d9911fb04, 32776, -2.96247569751492811574e+2L),
+sb1u = LD80C(0x81130861eae89202, 5, 3.22685866641642338307e+1L),
+sb2u = LD80C(0xbdb73fddefb0a7a7, 8, 3.79431636564289246766e+2L),
+sb3u = LD80C(0x80024bbf0b6ea3b3, 11, 2.04814349274125764278e+3L),
+sb4u = LD80C(0xa444b434a0a8549b, 12, 5.25658799100412086380e+3L),
+sb5u = LD80C(0xbb6ff428247dfadc, 12, 5.99799421719083271221e+3L),
+sb6u = LD80C(0x9b76cbd808a6f2e2, 11, 2.48742476657275909036e+3L),
+sb7u = LD80C(0xd3e4fe61373ce96a, 7, 2.11894506526933502549e+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.880e-22,4.903e-22]:
+ * |log(x*erfc(x)) + x**2 + 0.5625 - r(x)/s(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.36903347242128353893e+1L),
+rc3u = LD80C(0xb5cc0cecf47b4b20, 32774, -9.08985361145438643948e+1L),
+rc4u = LD80C(0xd749186e6457c66b, 32775, -2.15285529040811320620e+2L),
+rc5u = LD80C(0xfdd1925f754a9031, 32774, -1.26909319861476689424e+2L),
+sc1u = LD80C(0xc5d40fb999d75384, 4, 2.47285456180642684368e+1L),
+sc2u = LD80C(0xc5c6fe0e300c822b, 7, 1.97777314078074149137e+2L),
+sc3u = LD80C(0x961dc355db85da06, 9, 6.00465047325494566954e+2L),
+sc4u = LD80C(0x997ad30640d6bb4d, 9, 6.13919129908866058554e+2L),
+sc5u = LD80C(0xf7836d4d36940f8a, 6, 1.23756693280136487692e+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)
+{
+ 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);
+
+ /* erf(+-0) = +-0 */
+ if ((hx == 0 || hx == 0x8000) && lx == 0)
+ return (x);
+
+ /* x is INF or NaN */
+ if (hx >= 0x7fff) {
+ if (hx == 0xffff)
+ return (-1);
+ if (lx == 0x8000000000000000)
+ return (1);
+ return (x * x);
+ }
+
+ ENTERI();
+
+ ax = fabsl(x);
+ if(ax < 0.84375) {
+ if(ax < 0x1.p-34L) {
+ if (ax < 0x1.p-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 < LIM1) { /* |x| < 7 */
+ s = one/(ax*ax);
+ if(ax < LIM0) { /* |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);
+ } else {
+ if(x>=0) RETURNI(one-tiny); else RETURNI(tiny-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);
+
+ /* x is INF or NaN */
+ if (hx >= 0x7fff) {
+ if (hx == 0xffff)
+ return (2);
+ if (lx == 0x8000000000000000)
+ return (0);
+ return (x * x);
+ }
+
+ ENTERI();
+
+ ax = fabsl(x);
+ if(ax < 0.84375L) {
+ if(ax < 0x1.p-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 < LIM2) { /* |x|<108 */
+ s = one/(ax*ax);
+ if(ax < LIM0) { /* |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 < LIM1) { /* | |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 < -LIM1) 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: src/imprecise.c
===================================================================
--- src/imprecise.c (revision 258855)
+++ src/imprecise.c (working copy)
@@ -61,8 +61,6 @@
DECLARE_WEAK(f ## l)
DECLARE_IMPRECISE(cosh);
-DECLARE_IMPRECISE(erfc);
-DECLARE_IMPRECISE(erf);
DECLARE_IMPRECISE(lgamma);
DECLARE_IMPRECISE(sinh);
DECLARE_IMPRECISE(tanh);
Index: src/s_erf.c
===================================================================
--- src/s_erf.c (revision 258855)
+++ src/s_erf.c (working copy)
@@ -121,8 +121,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 */
@@ -201,7 +201,7 @@
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 (8*x+efx8*x)/8; /* avoid spurious underflow */
return x + efx*x;
}
z = x*x;
@@ -238,6 +238,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 +303,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 258855)
+++ 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);
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20131202214240.GA44209>