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