Now that we have random_normal(), it seems like it would be useful to
add the error functions erf() and erfc(), which I think are
potentially useful to the people who will find random_normal() useful,
and possibly others.
An immediate use for erf() is that it allows us to do a
Kolmogorov-Smirnov test for random_normal(), similar to the one for
random().
Both of these functions are defined in POSIX and C99, so in theory
they should be available on all platforms. If that turns out not to be
the case, then there's a commonly used implementation (e.g., see [1]),
which we could include. I played around with that (replacing the
direct bit manipulation stuff with frexp()/ldexp(), see pg_erf.c
attached), and it appeared to be accurate to +/-1 ULP across the full
range of inputs. Hopefully we won't need that though.
I tested this on a couple of different platforms and found I needed to
reduce extra_float_digits to -1 to get the regression tests to pass
consistently, due to rounding errors. It wouldn't surprise me if that
needs to be reduced further, though perhaps it's not necessary to have
so many tests (I included one test value from each branch, while
testing the hand-rolled implementation).
Regards,
Dean
[1] https://github.com/lsds/musl/blob/master/src/math/erf.c
diff --git a/doc/src/sgml/func.sgml b/doc/src/sgml/func.sgml
new file mode 100644
index 0cbdf63..e30e285
--- a/doc/src/sgml/func.sgml
+++ b/doc/src/sgml/func.sgml
@@ -1289,6 +1289,41 @@ repeat('Pg', 4) <returnvalue>PgPgPgPg</r
<row>
<entry role="func_table_entry"><para role="func_signature">
<indexterm>
+ <primary>erf</primary>
+ </indexterm>
+ <function>erf</function> ( <type>double precision</type> )
+ <returnvalue>double precision</returnvalue>
+ </para>
+ <para>
+ Error function
+ </para>
+ <para>
+ <literal>erf(1.0)</literal>
+ <returnvalue>0.8427007929497149</returnvalue>
+ </para></entry>
+ </row>
+
+ <row>
+ <entry role="func_table_entry"><para role="func_signature">
+ <indexterm>
+ <primary>erfc</primary>
+ </indexterm>
+ <function>erfc</function> ( <type>double precision</type> )
+ <returnvalue>double precision</returnvalue>
+ </para>
+ <para>
+ Complementary error function (<literal>1 - erf(x)</literal>, without
+ loss of precision for large inputs)
+ </para>
+ <para>
+ <literal>erfc(1.0)</literal>
+ <returnvalue>0.15729920705028513</returnvalue>
+ </para></entry>
+ </row>
+
+ <row>
+ <entry role="func_table_entry"><para role="func_signature">
+ <indexterm>
<primary>exp</primary>
</indexterm>
<function>exp</function> ( <type>numeric</type> )
diff --git a/src/backend/utils/adt/float.c b/src/backend/utils/adt/float.c
new file mode 100644
index d290b4c..aa79487
--- a/src/backend/utils/adt/float.c
+++ b/src/backend/utils/adt/float.c
@@ -2742,6 +2742,53 @@ datanh(PG_FUNCTION_ARGS)
}
+/* ========== ERROR FUNCTIONS ========== */
+
+
+/*
+ * derf - returns the error function: erf(arg1)
+ */
+Datum
+derf(PG_FUNCTION_ARGS)
+{
+ float8 arg1 = PG_GETARG_FLOAT8(0);
+ float8 result;
+
+ /*
+ * For erf, we don't need an errno check because it never overflows.
+ */
+ result = erf(arg1);
+
+ if (unlikely(isinf(result)))
+ float_overflow_error();
+
+ PG_RETURN_FLOAT8(result);
+}
+
+/*
+ * derfc - returns the complementary error function: 1 - erf(arg1)
+ */
+Datum
+derfc(PG_FUNCTION_ARGS)
+{
+ float8 arg1 = PG_GETARG_FLOAT8(0);
+ float8 result;
+
+ /*
+ * For erfc, we don't need an errno check because it never overflows.
+ */
+ result = erfc(arg1);
+
+ if (unlikely(isinf(result)))
+ float_overflow_error();
+
+ PG_RETURN_FLOAT8(result);
+}
+
+
+/* ========== RANDOM FUNCTIONS ========== */
+
+
/*
* initialize_drandom_seed - initialize drandom_seed if not yet done
*/
diff --git a/src/include/catalog/pg_proc.dat b/src/include/catalog/pg_proc.dat
new file mode 100644
index e2a7642..4ae8fef
--- a/src/include/catalog/pg_proc.dat
+++ b/src/include/catalog/pg_proc.dat
@@ -3465,6 +3465,13 @@
proname => 'atanh', prorettype => 'float8', proargtypes => 'float8',
prosrc => 'datanh' },
+{ oid => '8788', descr => 'error function',
+ proname => 'erf', prorettype => 'float8', proargtypes => 'float8',
+ prosrc => 'derf' },
+{ oid => '8789', descr => 'complementary error function',
+ proname => 'erfc', prorettype => 'float8', proargtypes => 'float8',
+ prosrc => 'derfc' },
+
{ oid => '1618',
proname => 'interval_mul', prorettype => 'interval',
proargtypes => 'interval float8', prosrc => 'interval_mul' },
diff --git a/src/test/regress/expected/float8.out b/src/test/regress/expected/float8.out
new file mode 100644
index 2b25784..0e89e24
--- a/src/test/regress/expected/float8.out
+++ b/src/test/regress/expected/float8.out
@@ -790,6 +790,45 @@ SELECT atanh(float8 'nan');
NaN
(1 row)
+-- error functions
+-- we run these with extra_float_digits = -1, to get consistently rounded
+-- results on all platforms.
+SET extra_float_digits = -1;
+SELECT x,
+ erf(x),
+ erfc(x)
+FROM (VALUES (float8 '-infinity'),
+ (-28), (-6), (-3.4), (-2.1), (-1.1), (-0.45),
+ (-1.2e-9), (-2.3e-13), (-1.2e-17), (0),
+ (1.2e-17), (2.3e-13), (1.2e-9),
+ (0.45), (1.1), (2.1), (3.4), (6), (28),
+ (float8 'infinity'), (float8 'nan')) AS t(x);
+ x | erf | erfc
+-----------+----------------------+---------------------
+ -Infinity | -1 | 2
+ -28 | -1 | 2
+ -6 | -1 | 2
+ -3.4 | -0.99999847800664 | 1.9999984780066
+ -2.1 | -0.99702053334367 | 1.9970205333437
+ -1.1 | -0.88020506957408 | 1.8802050695741
+ -0.45 | -0.47548171978692 | 1.4754817197869
+ -1.2e-09 | -1.3540550005146e-09 | 1.0000000013541
+ -2.3e-13 | -2.5952720843197e-13 | 1.0000000000003
+ -1.2e-17 | -1.3540550005146e-17 | 1
+ 0 | 0 | 1
+ 1.2e-17 | 1.3540550005146e-17 | 1
+ 2.3e-13 | 2.5952720843197e-13 | 0.99999999999974
+ 1.2e-09 | 1.3540550005146e-09 | 0.99999999864595
+ 0.45 | 0.47548171978692 | 0.52451828021308
+ 1.1 | 0.88020506957408 | 0.11979493042592
+ 2.1 | 0.99702053334367 | 0.002979466656333
+ 3.4 | 0.99999847800664 | 1.5219933628623e-06
+ 6 | 1 | 2.1519736712499e-17
+ 28 | 1 | 0
+ Infinity | 1 | 0
+ NaN | NaN | NaN
+(22 rows)
+
RESET extra_float_digits;
-- test for over- and underflow
INSERT INTO FLOAT8_TBL(f1) VALUES ('10e400');
diff --git a/src/test/regress/expected/random.out b/src/test/regress/expected/random.out
new file mode 100644
index 6dbb43a..2235907
--- a/src/test/regress/expected/random.out
+++ b/src/test/regress/expected/random.out
@@ -88,6 +88,38 @@ GROUP BY r;
-10 | 100
(1 row)
+-- Check standard normal distribution using the Kolmogorov-Smirnov test.
+CREATE FUNCTION ks_test_normal_random()
+RETURNS boolean AS
+$$
+DECLARE
+ n int := 1000; -- Number of samples
+ c float8 := 1.94947; -- Critical value for 99.9% confidence
+ ok boolean;
+BEGIN
+ ok := (
+ WITH samples AS (
+ SELECT random_normal() r FROM generate_series(1, n) ORDER BY 1
+ ), indexed_samples AS (
+ SELECT (row_number() OVER())-1.0 i, r FROM samples
+ )
+ SELECT max(abs((1+erf(r/sqrt(2)))/2 - i/n)) < c / sqrt(n)
+ FROM indexed_samples
+ );
+ RETURN ok;
+END
+$$
+LANGUAGE plpgsql;
+-- As above, ks_test_normal_random() returns true about 99.9%
+-- of the time, so try it 3 times and accept if any test passes.
+SELECT ks_test_normal_random() OR
+ ks_test_normal_random() OR
+ ks_test_normal_random() AS standard_normal;
+ standard_normal
+-----------------
+ t
+(1 row)
+
-- setseed() should produce a reproducible series of random() values.
SELECT setseed(0.5);
setseed
diff --git a/src/test/regress/sql/float8.sql b/src/test/regress/sql/float8.sql
new file mode 100644
index c276a53..de99ef0
--- a/src/test/regress/sql/float8.sql
+++ b/src/test/regress/sql/float8.sql
@@ -229,6 +229,20 @@ SELECT atanh(float8 'infinity');
SELECT atanh(float8 '-infinity');
SELECT atanh(float8 'nan');
+-- error functions
+-- we run these with extra_float_digits = -1, to get consistently rounded
+-- results on all platforms.
+SET extra_float_digits = -1;
+SELECT x,
+ erf(x),
+ erfc(x)
+FROM (VALUES (float8 '-infinity'),
+ (-28), (-6), (-3.4), (-2.1), (-1.1), (-0.45),
+ (-1.2e-9), (-2.3e-13), (-1.2e-17), (0),
+ (1.2e-17), (2.3e-13), (1.2e-9),
+ (0.45), (1.1), (2.1), (3.4), (6), (28),
+ (float8 'infinity'), (float8 'nan')) AS t(x);
+
RESET extra_float_digits;
-- test for over- and underflow
diff --git a/src/test/regress/sql/random.sql b/src/test/regress/sql/random.sql
new file mode 100644
index 6e99e2e..14cc76b
--- a/src/test/regress/sql/random.sql
+++ b/src/test/regress/sql/random.sql
@@ -70,6 +70,36 @@ SELECT r, count(*)
FROM (SELECT random_normal(-10, 0) r FROM generate_series(1, 100)) ss
GROUP BY r;
+-- Check standard normal distribution using the Kolmogorov-Smirnov test.
+
+CREATE FUNCTION ks_test_normal_random()
+RETURNS boolean AS
+$$
+DECLARE
+ n int := 1000; -- Number of samples
+ c float8 := 1.94947; -- Critical value for 99.9% confidence
+ ok boolean;
+BEGIN
+ ok := (
+ WITH samples AS (
+ SELECT random_normal() r FROM generate_series(1, n) ORDER BY 1
+ ), indexed_samples AS (
+ SELECT (row_number() OVER())-1.0 i, r FROM samples
+ )
+ SELECT max(abs((1+erf(r/sqrt(2)))/2 - i/n)) < c / sqrt(n)
+ FROM indexed_samples
+ );
+ RETURN ok;
+END
+$$
+LANGUAGE plpgsql;
+
+-- As above, ks_test_normal_random() returns true about 99.9%
+-- of the time, so try it 3 times and accept if any test passes.
+SELECT ks_test_normal_random() OR
+ ks_test_normal_random() OR
+ ks_test_normal_random() AS standard_normal;
+
-- setseed() should produce a reproducible series of random() values.
SELECT setseed(0.5);
/*-------------------------------------------------------------------------
*
* pg_erf.c
* Portable implementations of the error functions erf() and erfc().
*
* Copyright (c) 2023, PostgreSQL Global Development Group
*
* IDENTIFICATION
* src/port/pg_erf.c
*
*-------------------------------------------------------------------------
*/
#include "c.h"
#include <math.h>
/* ----------
* Original code, including copyright and long comment below, from FreeBSD's
* /usr/src/lib/msun/src/s_erf.c, modified to avoid direct bit manipulation of
* double values.
*
* ====================================================
* 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.
* ====================================================
*
* double erf(double x)
* double erfc(double x)
* x
* 2 |\
* erf(x) = --------- | exp(-t*t)dt
* sqrt(pi) \|
* 0
*
* erfc(x) = 1-erf(x)
* Note that
* erf(-x) = -erf(x)
* erfc(-x) = 2 - erfc(x)
*
* Method:
* 1. For |x| in [0, 0.84375]
* erf(x) = x + x*R(x^2)
* erfc(x) = 1 - erf(x) if x in [-.84375,0.25]
* = 0.5 + ((0.5-x)-x*R) if x in [0.25,0.84375]
* where R = P/Q where P is an odd poly of degree 8 and
* Q is an odd poly of degree 10.
* -57.90
* | R - (erf(x)-x)/x | <= 2
*
*
* Remark. The formula is derived by noting
* erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10 - x^7/42 + ....)
* and that
* 2/sqrt(pi) = 1.128379167095512573896158903121545171688
* is close to one. The interval is chosen because the fix
* point of erf(x) is near 0.6174 (i.e., erf(x)=x when x is
* near 0.6174), and by some experiment, 0.84375 is chosen to
* guarantee the error is less than one ulp for erf.
*
* 2. For |x| in [0.84375,1.25], let s = |x| - 1, and
* c = 0.84506291151 rounded to single (24 bits)
* erf(x) = sign(x) * (c + P1(s)/Q1(s))
* erfc(x) = (1-c) - P1(s)/Q1(s) if x > 0
* 1+(c+P1(s)/Q1(s)) if x < 0
* |P1/Q1 - (erf(|x|)-c)| <= 2**-59.06
* Remark: here we use the taylor series expansion at x=1.
* erf(1+s) = erf(1) + s*Poly(s)
* = 0.845.. + P1(s)/Q1(s)
* That is, we use rational approximation to approximate
* erf(1+s) - (c = (single)0.84506291151)
* Note that |P1/Q1|< 0.078 for x in [0.84375,1.25]
* where
* P1(s) = degree 6 poly in s
* Q1(s) = degree 6 poly in s
*
* 3. For x in [1.25,1/0.35(~2.857143)],
* erfc(x) = (1/x)*exp(-x*x-0.5625+R1/S1)
* erf(x) = 1 - erfc(x)
* where
* R1(z) = degree 7 poly in z, (z=1/x^2)
* S1(z) = degree 8 poly in z
*
* 4. For x in [1/0.35,28]
* erfc(x) = (1/x)*exp(-x*x-0.5625+R2/S2) if x > 0
* = 2.0 - (1/x)*exp(-x*x-0.5625+R2/S2) if -6<x<0
* = 2.0 - tiny (if x <= -6)
* erf(x) = sign(x)*(1.0 - erfc(x)) if x < 6, else
* erf(x) = sign(x)*(1.0 - tiny)
* where
* R2(z) = degree 6 poly in z, (z=1/x^2)
* S2(z) = degree 7 poly in z
*
* Note1:
* To compute exp(-x*x-0.5625+R/S), let s be a single
* precision number and s := x; then
* -x*x = -s*s + (s-x)*(s+x)
* exp(-x*x-0.5626+R/S) =
* exp(-s*s-0.5625)*exp((s-x)*(s+x)+R/S);
* Note2:
* Here 4 and 5 make use of the asymptotic series
* exp(-x*x)
* erfc(x) ~ ---------- * ( 1 + Poly(1/x^2) )
* x*sqrt(pi)
* We use rational approximation to approximate
* g(s)=f(1/x^2) = log(erfc(x)*x) - x*x + 0.5625
* Here is the error bound for R1/S1 and R2/S2
* |R1/S1 - f(x)| < 2**(-62.57)
* |R2/S2 - f(x)| < 2**(-61.52)
*
* 5. For inf > x >= 28
* erf(x) = sign(x) *(1 - tiny) (raise inexact)
* erfc(x) = tiny*tiny (raise underflow) if x > 0
* = 2 - tiny if x<0
*
* 7. Special case:
* erf(0) = 0, erf(inf) = 1, erf(-inf) = -1,
* erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2,
* erfc/erf(NaN) is NaN
*
* ----------
*/
static const double
erx = 8.45062911510467529297e-01, /* 0x3FEB0AC1, 0x60000000 */
/*
* Coefficients for approximation of erf() in the range [0, 0.84375]
*/
efx8 = 1.02703333676410069053e+00, /* 0x3FF06EBA, 0x8214DB69 */
pp0 = 1.28379167095512558561e-01, /* 0x3FC06EBA, 0x8214DB68 */
pp1 = -3.25042107247001499370e-01, /* 0xBFD4CD7D, 0x691CB913 */
pp2 = -2.84817495755985104766e-02, /* 0xBF9D2A51, 0xDBD7194F */
pp3 = -5.77027029648944159157e-03, /* 0xBF77A291, 0x236668E4 */
pp4 = -2.37630166566501626084e-05, /* 0xBEF8EAD6, 0x120016AC */
qq1 = 3.97917223959155352819e-01, /* 0x3FD97779, 0xCDDADC09 */
qq2 = 6.50222499887672944485e-02, /* 0x3FB0A54C, 0x5536CEBA */
qq3 = 5.08130628187576562776e-03, /* 0x3F74D022, 0xC4D36B0F */
qq4 = 1.32494738004321644526e-04, /* 0x3F215DC9, 0x221C1A10 */
qq5 = -3.96022827877536812320e-06, /* 0xBED09C43, 0x42A26120 */
/*
* Coefficients for approximation of erf() in the range [0.84375, 1.25]
*/
pa0 = -2.36211856075265944077e-03, /* 0xBF6359B8, 0xBEF77538 */
pa1 = 4.14856118683748331666e-01, /* 0x3FDA8D00, 0xAD92B34D */
pa2 = -3.72207876035701323847e-01, /* 0xBFD7D240, 0xFBB8C3F1 */
pa3 = 3.18346619901161753674e-01, /* 0x3FD45FCA, 0x805120E4 */
pa4 = -1.10894694282396677476e-01, /* 0xBFBC6398, 0x3D3E28EC */
pa5 = 3.54783043256182359371e-02, /* 0x3FA22A36, 0x599795EB */
pa6 = -2.16637559486879084300e-03, /* 0xBF61BF38, 0x0A96073F */
qa1 = 1.06420880400844228286e-01, /* 0x3FBB3E66, 0x18EEE323 */
qa2 = 5.40397917702171048937e-01, /* 0x3FE14AF0, 0x92EB6F33 */
qa3 = 7.18286544141962662868e-02, /* 0x3FB2635C, 0xD99FE9A7 */
qa4 = 1.26171219808761642112e-01, /* 0x3FC02660, 0xE763351F */
qa5 = 1.36370839120290507362e-02, /* 0x3F8BEDC2, 0x6B51DD1C */
qa6 = 1.19844998467991074170e-02, /* 0x3F888B54, 0x5735151D */
/*
* Coefficients for approximation of erfc() in the range [1.25, 1/0.35]
*/
ra0 = -9.86494403484714822705e-03, /* 0xBF843412, 0x600D6435 */
ra1 = -6.93858572707181764372e-01, /* 0xBFE63416, 0xE4BA7360 */
ra2 = -1.05586262253232909814e+01, /* 0xC0251E04, 0x41B0E726 */
ra3 = -6.23753324503260060396e+01, /* 0xC04F300A, 0xE4CBA38D */
ra4 = -1.62396669462573470355e+02, /* 0xC0644CB1, 0x84282266 */
ra5 = -1.84605092906711035994e+02, /* 0xC067135C, 0xEBCCABB2 */
ra6 = -8.12874355063065934246e+01, /* 0xC0545265, 0x57E4D2F2 */
ra7 = -9.81432934416914548592e+00, /* 0xC023A0EF, 0xC69AC25C */
sa1 = 1.96512716674392571292e+01, /* 0x4033A6B9, 0xBD707687 */
sa2 = 1.37657754143519042600e+02, /* 0x4061350C, 0x526AE721 */
sa3 = 4.34565877475229228821e+02, /* 0x407B290D, 0xD58A1A71 */
sa4 = 6.45387271733267880336e+02, /* 0x40842B19, 0x21EC2868 */
sa5 = 4.29008140027567833386e+02, /* 0x407AD021, 0x57700314 */
sa6 = 1.08635005541779435134e+02, /* 0x405B28A3, 0xEE48AE2C */
sa7 = 6.57024977031928170135e+00, /* 0x401A47EF, 0x8E484A93 */
sa8 = -6.04244152148580987438e-02, /* 0xBFAEEFF2, 0xEE749A62 */
/*
* Coefficients for approximation of erfc() in the range [1/0.35, 28]
*/
rb0 = -9.86494292470009928597e-03, /* 0xBF843412, 0x39E86F4A */
rb1 = -7.99283237680523006574e-01, /* 0xBFE993BA, 0x70C285DE */
rb2 = -1.77579549177547519889e+01, /* 0xC031C209, 0x555F995A */
rb3 = -1.60636384855821916062e+02, /* 0xC064145D, 0x43C5ED98 */
rb4 = -6.37566443368389627722e+02, /* 0xC083EC88, 0x1375F228 */
rb5 = -1.02509513161107724954e+03, /* 0xC0900461, 0x6A2E5992 */
rb6 = -4.83519191608651397019e+02, /* 0xC07E384E, 0x9BDC383F */
sb1 = 3.03380607434824582924e+01, /* 0x403E568B, 0x261D5190 */
sb2 = 3.25792512996573918826e+02, /* 0x40745CAE, 0x221B9F0A */
sb3 = 1.53672958608443695994e+03, /* 0x409802EB, 0x189D5118 */
sb4 = 3.19985821950859553908e+03, /* 0x40A8FFB7, 0x688C246A */
sb5 = 2.55305040643316442583e+03, /* 0x40A3F219, 0xCEDF3BE6 */
sb6 = 4.74528541206955367215e+02, /* 0x407DA874, 0xE79FE763 */
sb7 = -2.24409524465858183362e+01; /* 0xC03670E2, 0x42712D62 */
static double
erfc1(double x) /* 0.84375 <= x < 1.25 */
{
double s,
P,
Q;
s = x - 1;
P = pa0 + s * (pa1 + s * (pa2 + s * (pa3 + s * (pa4 + s * (pa5 + s * pa6)))));
Q = 1 + s * (qa1 + s * (qa2 + s * (qa3 + s * (qa4 + s * (qa5 + s * qa6)))));
return 1 - erx - P / Q;
}
static double
erfc2(double x) /* 0.84375 <= x < 28 */
{
double s,
R,
S,
z;
int xexp;
if (x < 1.25) /* 0.84375 <= x < 1.25 */
return erfc1(x);
s = 1 / (x * x);
if (x < 1 / 0.35) /* 1.25 <= x < 1/0.35 ~ 2.85714 */
{
R = ra0 + s * (ra1 + s * (ra2 + s * (ra3 + s * (ra4 + s * (ra5 + s * (ra6 + s * ra7))))));
S = 1.0 + s * (sa1 + s * (sa2 + s * (sa3 + s * (sa4 + s * (sa5 + s * (sa6 + s * (sa7 + s * sa8)))))));
}
else /* 1/0.35 <= x < 28 */
{
R = rb0 + s * (rb1 + s * (rb2 + s * (rb3 + s * (rb4 + s * (rb5 + s * rb6)))));
S = 1.0 + s * (sb1 + s * (sb2 + s * (sb3 + s * (sb4 + s * (sb5 + s * (sb6 + s * sb7))))));
}
/* ----------
* Original code:
* z = x;
* SET_LOW_WORD(z,0);
* This discards the 32 least significant bits of x, and therefore keeps
* the 20 most significant bits of the mantissa. Use frexp()/ldexp()
* instead, for portability.
*
* Note that frexp() returns a value in the range [0.5, 1), so the most
* significant bit is always set, and doesn't count as one of the 20 bits
* that we need to keep.
* ----------
*/
z = frexp(x, &xexp);
z = ldexp(trunc(ldexp(z, 21)), xexp - 21);
return exp(-z * z - 0.5625) * exp((z - x) * (z + x) + R / S) / x;
}
double
pg_erf(double x)
{
const double two_m28 = 1.0 / (1 << 28);
double r,
s,
z,
y;
int sign;
/* erf(nan)=nan */
if (isnan(x))
return x;
/* erf(+-inf)=+-1 */
sign = x < 0 ? -1 : 1;
if (isinf(x))
return (double) sign;
x = fabs(x);
if (x < 0.84375) /* 0 <= x < 0.84375 */
{
if (x < two_m28) /* 0 <= x < 2^(-28) */
{
/* avoid underflow */
return sign * 0.125 * (8 * x + efx8 * x);
}
z = x * x;
r = pp0 + z * (pp1 + z * (pp2 + z * (pp3 + z * pp4)));
s = 1.0 + z * (qq1 + z * (qq2 + z * (qq3 + z * (qq4 + z * qq5))));
y = r / s;
return sign * (x + x * y);
}
if (x < 6) /* 0.84375 <= x < 6 */
y = 1 - erfc2(x);
else /* x >= 6 */
y = 1;
return sign * y;
}
double
pg_erfc(double x)
{
const double two_m56 = 1.0 / (((int64) 1) << 56);
double r,
s,
z,
y;
int sign;
/* erfc(nan)=nan */
if (isnan(x))
return x;
/* erfc(+-inf)=0,2 */
sign = x < 0 ? -1 : 1;
if (isinf(x))
return (double) (1 - sign);
x = fabs(x);
if (x < 0.84375) /* 0 <= x < 0.84375 */
{
if (x < two_m56) /* 0 <= x < 2^(-56) */
return 1.0 - sign * x;
z = x * x;
r = pp0 + z * (pp1 + z * (pp2 + z * (pp3 + z * pp4)));
s = 1.0 + z * (qq1 + z * (qq2 + z * (qq3 + z * (qq4 + z * qq5))));
y = r / s;
if (x < 0.25) /* 2^(-56) <= x < 1/4 */
return 1.0 - sign * (x + x * y);
return 0.5 - (sign * x - 0.5 + sign * x * y);
}
if (x < 28) /* 0.84375 <= x < 28 */
return sign < 0 ? 2 - erfc2(x) : erfc2(x);
return sign < 0 ? 2 : 0;
}
