[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Bug-gsl] Bug in ./interpolation/integ_eval.h (patch supplied)
From: |
Patricio Rojo |
Subject: |
[Bug-gsl] Bug in ./interpolation/integ_eval.h (patch supplied) |
Date: |
Thu, 24 Jun 2004 20:28:52 -0400 |
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
[Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Bug-gsl] Bug in ./interpolation/integ_eval.h (patch supplied),
Patricio Rojo <=