Skip site navigation (1)Skip section navigation (2)
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: <http://docs.FreeBSD.org/cgi/mid.cgi?20131202214240.GA44209>