Dear Patrick and help-gsl,
I had found, a few days ago, a test failure in interpolation/test.c.
Then just yesterday (May 10) there was a revert of interpolation/test.c to a
previous version, which yanked out the test failure I had been investigating.
But the problem still exists, so let me outline here what I found and make some
suggestions.
The test_cspline4() in git hash e647c6484702b638062d6e2154055f6535790cd7 has a
failure.
The problem has two parts: what brings it out, and what
1. test_cspline4() calls make_xy_table() with an argument of 2, which is not OK
- I guess you need at least 3 to do csplines, and things go wrong down the line
when you start doing linear algebra.
The test should be updated to have an argument of 3, or the cspline subsystem
should be updated to warn about that.
2. the real problem is that all the intervening function calls do not check for
that size properly. At one point two gets subtracted from it, and we end up
with calls to malloc(0). Note that if by chance we had called with 1 or 0, we
might end up with malloc(-1) or malloc(-2).
The behavior of malloc(0) surprised me: I had thought it would just return
NULL, but the man page states:
"""The memory is not initialized. If size is 0, then malloc() returns either
NULL, or a unique pointer value that can later be successfully passed to
free()."""
which means that this code in tridiag.c:
static
int
solve_tridiag(
const double diag[], size_t d_stride,
const double offdiag[], size_t o_stride,
const double b[], size_t b_stride,
double x[], size_t x_stride,
size_t N)
{
int status = GSL_SUCCESS;
double *gamma = (double *) malloc (N * sizeof (double));
double *alpha = (double *) malloc (N * sizeof (double));
double *c = (double *) malloc (N * sizeof (double));
double *z = (double *) malloc (N * sizeof (double));
if (gamma == 0 || alpha == 0 || c == 0 || z == 0)
{
GSL_ERROR("failed to allocate working space", GSL_ENOMEM);
}
else
does not capture this kind of error, and we end up with bad access, since the
interpolation index arithmetic is not up to handling a size of 2 (and hence N
of 0).
I would suggest some things:
a. gsl vectors allow a size of zero, so we should not trap at the level of the
gsl_vector subsystem
b. tridiag.c needs to require N > 0. Right now it looks at the return values
of malloc():
double *gamma = (double *) malloc (N * sizeof (double));
double *alpha = (double *) malloc (N * sizeof (double));
double *c = (double *) malloc (N * sizeof (double));
double *z = (double *) malloc (N * sizeof (double));
if (gamma == 0 || alpha == 0 || c == 0 || z == 0)
{
GSL_ERROR("failed to allocate working space", GSL_ENOMEM);
}
else
and then it does unsigned arithmetic here:
for (i = 1; i < N - 1; i++)
{
alpha[i] = diag[d_stride * i] - offdiag[o_stride*(i - 1)] * gamma[i -
1];
where N-1 is:
(gdb) print N-1
$50 = 18446744073709551615
(gdb)
so clearly a bug in allowing N to be zero. Maybe something like:
if (N <= 0)
{
GSL_ERROR("matrix size must be positive", GSL_EDOM);
}
else if (gamma == 0 || alpha == 0 || c == 0 || z == 0)
{
GSL_ERROR("failed to allocate working space", GSL_ENOMEM);
}
else
and probably also do that kind of checking in the higher level function
gsl_linalg_solve_symm_tridiag(), and probably many other parts of linear
algebra. This is a real project to examine the linear algebra package.
c. take cspline_init() and give it a check on size, insisting that it be >= 3.
