Re: [Bug-gsl] Re: Interpolation

[Евгений Курбатов]

Hi Brian!

Here the patches, the hardcoded range checkup and fix for error codes returning.

interp.c.patch - performs a range checkup in gsl_interp_eval_... functions

akima.c.patch
cspline.c.patch - return GSL_EINVAL (as in linear.c), instead of GSL_FAILURE, 
when ordering violation of x_array.

Regards,
Evgeny Kurbatov
diff --git a/interpolation/akima.c b/interpolation/akima.c
index 551bb8b..4244e0f 100644
--- a/interpolation/akima.c
+++ b/interpolation/akima.c
@@ -350,7 +350,7 @@ akima_eval_integ (const void * vstate,
     }
     else {
       *result = 0.0;
-      return GSL_FAILURE;
+      return GSL_EINVAL;
     }
   }
   

diff --git a/interpolation/cspline.c b/interpolation/cspline.c
index e6a1abc..37c73d4 100644
--- a/interpolation/cspline.c
+++ b/interpolation/cspline.c
@@ -334,7 +334,7 @@ cspline_eval_deriv (const void * vstate,
   else
     {
       *dydx = 0.0;
-      return GSL_FAILURE;
+      return GSL_EINVAL;
     }
 }
 
@@ -380,7 +380,7 @@ cspline_eval_deriv2 (const void * vstate,
   else
     {
       *y_pp = 0.0;
-      return GSL_FAILURE;
+      return GSL_EINVAL;
     }
 }
 
@@ -435,7 +435,7 @@ cspline_eval_integ (const void * vstate,
     }
     else {
       *result = 0.0;
-      return GSL_FAILURE;
+      return GSL_EINVAL;
     }
   }
   

diff --git a/interpolation/interp.c b/interpolation/interp.c
index bb75e10..275735f 100644
--- a/interpolation/interp.c
+++ b/interpolation/interp.c
@@ -142,7 +142,18 @@ gsl_interp_eval (const gsl_interp * interp,
                  gsl_interp_accel * a)
 {
   double y;
-  int status = interp->type->eval (interp->state, xa, ya, interp->size, x, a, &y);
+  int status;
+
+  if (x < interp->xmin)
+    {
+      GSL_ERROR_VAL("interpolation error", GSL_EDOM, ya[0]);
+    }
+  else if (x > interp->xmax)
+    {
+      GSL_ERROR_VAL("interpolation error", GSL_EDOM, ya[interp->size - 1]);
+    }
+
+  status = interp->type->eval (interp->state, xa, ya, interp->size, x, a, &y);
 
   DISCARD_STATUS(status);
 
@@ -176,7 +187,18 @@ gsl_interp_eval_deriv (const gsl_interp * interp,
                        gsl_interp_accel * a)
 {
   double dydx;
-  int status = interp->type->eval_deriv (interp->state, xa, ya, interp->size, x, a, &dydx);
+  int status;
+
+  if (x < interp->xmin)
+    {
+      GSL_ERROR_VAL("interpolation error", GSL_EDOM, 0.0);
+    }
+  else if (x > interp->xmax)
+    {
+      GSL_ERROR_VAL("interpolation error", GSL_EDOM, 0.0);
+    }
+
+  status = interp->type->eval_deriv (interp->state, xa, ya, interp->size, x, a, &dydx);
 
   DISCARD_STATUS(status);
 
@@ -210,7 +232,18 @@ gsl_interp_eval_deriv2 (const gsl_interp * interp,
                         gsl_interp_accel * a)
 {
   double d2;
-  int status = interp->type->eval_deriv2 (interp->state, xa, ya, interp->size, x, a, &d2);
+  int status;
+
+  if (x < interp->xmin)
+    {
+      GSL_ERROR_VAL("interpolation error", GSL_EDOM, 0.0);
+    }
+  else if (x > interp->xmax)
+    {
+      GSL_ERROR_VAL("interpolation error", GSL_EDOM, 0.0);
+    }
+
+  status = interp->type->eval_deriv2 (interp->state, xa, ya, interp->size, x, a, &d2);
 
   DISCARD_STATUS(status);
 
@@ -247,7 +280,18 @@ gsl_interp_eval_integ (const gsl_interp * interp,
                        gsl_interp_accel * acc)
 {
   double result;
-  int status = interp->type->eval_integ (interp->state, xa, ya, interp->size, acc, a, b, &result);
+  int status;
+
+  if (a > b || a < interp->xmin || b > interp->xmax)
+    {
+      GSL_ERROR_VAL("interpolation error", GSL_EDOM, 0.0);
+    }
+  else if(a == b)
+    {
+      return 0.0;
+    }
+
+  status = interp->type->eval_integ (interp->state, xa, ya, interp->size, acc, a, b, &result);
 
   DISCARD_STATUS(status);
 

_______________________________________________
Bug-gsl mailing list
Bug-gsl@gnu.org
http://lists.gnu.org/mailman/listinfo/bug-gsl