[Bug-gsl] Bug in ./interpolation/integ_eval.h (patch supplied)

[Patricio Rojo]

Hi!,

  More than an important error, this is a problem that appears only when
using a very "unlucky" set of values (which were just those I was
using:D).

  In the interpolation/integ_eval.h file, to compute dterm:
 const double dterm =
  di / 4.0 * (t2 - 2.0 * xi * (2.0 * t1 - xi * (3.0 * t0 - 2.0 * xi)));

  it can happen with some values that t2 is very similar to what it is
subtracted from it. Similarity that can go beyond precision, but which
is non-zero, nevertheless.
  I couldn't find any clear mathematical circumstance under which this
problem happen. It seems it just was an incredible coincidence that this
happened for the values I was using (xi=3.3255576, a=3.3255576,
b=3.3255578).

   At the end, the result was that this function was returning a value
which was more than 2 order of magnitude different from what it should
have been.

   The solution I propose in the patch below not only gets rid of this
problem by having only sums (no subtractions), but it is also more
efficient by having less floating point operations.

   I discover this problem and its solution by using cspline.c. Because
my solution is only an algebraic rearrangements of terms, I expect it to
work with Akima interpolation as well, but I didn't test that:).

   I'm not sure if I edited the copyright in the right way, feel free to
point out my collaboration in a different way if you decide to accept
it...

   I'll be happy to provide more info, just let me know. Thank you!...

                              Patricio  

---------------------------------
GSL version: 1.4 (From the Debian distribution)
-----
#### gsl_1.4.pato.diff below ####
#####################################################
--- interpolation/integ_eval.h.old      2004-06-24 10:46:55.000000000 -0400
+++ interpolation/integ_eval.h  2004-06-24 16:21:08.000000000 -0400
@@ -1,6 +1,7 @@
 /* interpolation/integ_eval_macro.h
  *
  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
+ *               2004 modified by Patricio Rojo
  *
  * This program is free software; you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -25,13 +26,12 @@
 integ_eval (double ai, double bi, double ci, double di, double xi, double a,
             double b)
 {
-  const double t0 = b + a;
-  const double t1 = a * a + a * b + b * b;
-  const double t2 = a * a * a + a * a * b + b * b * a + b * b * b;
-  const double bterm = 0.5 * bi * (t0 - 2.0 * xi);
-  const double cterm = ci / 3.0 * (t1 - 3.0 * xi * (t0 - xi));
-  const double dterm =
-    di / 4.0 * (t2 - 2.0 * xi * (2.0 * t1 - xi * (3.0 * t0 - 2.0 * xi)));
+  const double r1 = a - xi;
+  const double r2 = b - xi;
+  const double rpr = r1 + r2;
+  const double bterm = bi * 0.5  * rpr;
+  const double cterm = ci / 3.0  * ( r1 * r1 + r2 * r2 + r1 * r2 );
+  const double dterm = di * 0.25 * rpr * ( r1 * r1 + r2 * r2 );

   return (b - a) * (ai + bterm + cterm + dterm);
 }
######################################################

-- 
 ----------------------------------------------------------------
      . .         /.
     .       *   /`'\  .                          Patricio Rojo
        .       /    ./ \..         516 Space Sciences Building
             ../   .'      \.                Cornell University
          ../    /    ../    \                 Ithaca, NY 14853
        ./     ./   ./        \                   (607)255-6438
      ./     ./            \   \.        address@hidden
 ----------------------------------------------------------------
 
 http://www.gnu.org/philosophy/no-word-attachments.html