[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Bug-gsl] [PATCH] handle large values correctly in (log_)erf(c) function
From: |
Alex Henrie |
Subject: |
[Bug-gsl] [PATCH] handle large values correctly in (log_)erf(c) functions |
Date: |
Mon, 12 Feb 2018 22:18:15 -0700 |
specfunc/erfc.c | 31 ++++++++++++++++++++++++++-----
specfunc/test_sf.c | 6 ++++++
2 files changed, 32 insertions(+), 5 deletions(-)
diff --git a/specfunc/erfc.c b/specfunc/erfc.c
index 819dbb02..4a7925ca 100644
--- a/specfunc/erfc.c
+++ b/specfunc/erfc.c
@@ -285,6 +285,12 @@ int gsl_sf_erfc_e(double x, gsl_sf_result * result)
}
else {
e_val = erfc8(ax);
+ if(isnan(e_val) && !isnan(x)) {
+ /* overflow */
+ result->val = (x > 0) ? 0 : 2;
+ result->err = 0;
+ return GSL_SUCCESS;
+ }
e_err = (x*x + 1.0) * GSL_DBL_EPSILON * fabs(e_val);
}
@@ -339,16 +345,31 @@ int gsl_sf_log_erfc_e(double x, gsl_sf_result * result)
}
*/
else if(x > 8.0) {
- result->val = log_erfc8(x);
- result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
+ double e_val = log_erfc8(x);
+ if (isnan(e_val) && !isnan(x)) {
+ /* positive overflow */
+ result->val = -HUGE_VAL;
+ result->err = 0;
+ }
+ else {
+ result->val = e_val;
+ result->err = 2.0 * GSL_DBL_EPSILON * fabs(e_val);
+ }
return GSL_SUCCESS;
}
else {
gsl_sf_result result_erfc;
gsl_sf_erfc_e(x, &result_erfc);
- result->val = log(result_erfc.val);
- result->err = fabs(result_erfc.err / result_erfc.val);
- result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
+ if (isnan(result_erfc.val) && !isnan(x)) {
+ /* negative overflow */
+ result->val = M_LN2;
+ result->err = 0;
+ }
+ else {
+ result->val = log(result_erfc.val);
+ result->err = fabs(result_erfc.err / result_erfc.val);
+ result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
+ }
return GSL_SUCCESS;
}
}
diff --git a/specfunc/test_sf.c b/specfunc/test_sf.c
index f3eff9d4..a55dbc9f 100644
--- a/specfunc/test_sf.c
+++ b/specfunc/test_sf.c
@@ -889,6 +889,7 @@ int test_erf(void)
gsl_sf_result r;
int s = 0;
+ TEST_SF(s, gsl_sf_erfc_e, (-HUGE_VAL, &r), 2.0, TEST_TOL0, GSL_SUCCESS);
TEST_SF(s, gsl_sf_erfc_e, (-10.0, &r), 2.0, TEST_TOL0, GSL_SUCCESS);
TEST_SF(s, gsl_sf_erfc_e, (-5.0000002, &r), 1.9999999999984625433,
TEST_TOL0, GSL_SUCCESS);
TEST_SF(s, gsl_sf_erfc_e, (-5.0, &r), 1.9999999999984625402, TEST_TOL0,
GSL_SUCCESS);
@@ -898,7 +899,9 @@ int test_erf(void)
TEST_SF(s, gsl_sf_erfc_e, (3.0, &r), 0.000022090496998585441373, TEST_TOL1,
GSL_SUCCESS);
TEST_SF(s, gsl_sf_erfc_e, (7.0, &r), 4.183825607779414399e-23, TEST_TOL2,
GSL_SUCCESS);
TEST_SF(s, gsl_sf_erfc_e, (10.0, &r), 2.0884875837625447570e-45, TEST_TOL2,
GSL_SUCCESS);
+ TEST_SF(s, gsl_sf_erfc_e, (HUGE_VAL, &r), 0.0, TEST_TOL2, GSL_SUCCESS);
+ TEST_SF(s, gsl_sf_log_erfc_e, (-HUGE_VAL, &r), M_LN2, TEST_TOL0,
GSL_SUCCESS);
TEST_SF(s, gsl_sf_log_erfc_e, (-1.0, &r), log(1.842700792949714869),
TEST_TOL0, GSL_SUCCESS);
TEST_SF(s, gsl_sf_log_erfc_e, (-0.1, &r), 0.106576400586522485015,
TEST_TOL0, GSL_SUCCESS);
TEST_SF(s, gsl_sf_log_erfc_e, (-1e-10, &r), 1.1283791670318505967e-10,
TEST_TOL0, GSL_SUCCESS);
@@ -908,11 +911,14 @@ int test_erf(void)
TEST_SF(s, gsl_sf_log_erfc_e, (0.1, &r), -0.119304973737395598329,
TEST_TOL0, GSL_SUCCESS);
TEST_SF(s, gsl_sf_log_erfc_e, (1.0, &r), log(0.15729920705028513066),
TEST_TOL0, GSL_SUCCESS);
TEST_SF(s, gsl_sf_log_erfc_e, (10.0, &r), log(2.0884875837625447570e-45),
TEST_TOL0, GSL_SUCCESS);
+ TEST_SF(s, gsl_sf_log_erfc_e, (HUGE_VAL, &r), -HUGE_VAL, TEST_TOL0,
GSL_SUCCESS);
+ TEST_SF(s, gsl_sf_erf_e, (-HUGE_VAL, &r), -1.0000000000000000000, TEST_TOL0,
GSL_SUCCESS);
TEST_SF(s, gsl_sf_erf_e, (-10.0, &r), -1.0000000000000000000, TEST_TOL0,
GSL_SUCCESS);
TEST_SF(s, gsl_sf_erf_e, (0.5, &r), 0.5204998778130465377, TEST_TOL0,
GSL_SUCCESS);
TEST_SF(s, gsl_sf_erf_e, (1.0, &r), 0.8427007929497148693, TEST_TOL0,
GSL_SUCCESS);
TEST_SF(s, gsl_sf_erf_e, (10.0, &r), 1.0000000000000000000, TEST_TOL0,
GSL_SUCCESS);
+ TEST_SF(s, gsl_sf_erf_e, (HUGE_VAL, &r), 1.0000000000000000000, TEST_TOL0,
GSL_SUCCESS);
TEST_SF(s, gsl_sf_erf_Z_e, (1.0, &r), 0.24197072451914334980, TEST_TOL0,
GSL_SUCCESS);
TEST_SF(s, gsl_sf_erf_Q_e, (10.0, &r), 7.619853024160526066e-24, TEST_TOL2,
GSL_SUCCESS);
--
2.16.1
[Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Bug-gsl] [PATCH] handle large values correctly in (log_)erf(c) functions,
Alex Henrie <=