bug-gsl
[Top][All Lists]
Advanced

[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




reply via email to

[Prev in Thread] Current Thread [Next in Thread]