2

I am attempting to implement the Cublic Spline Extrema Algorithm published in 1996.

I am not an accomplished programmer, and it is clear the author is, and is an order of magnitude or two smarter than I. However when I try and compile this, it seg faults on a simple out-of-bounds array access. I find it hard to believe that this algorithm was published with what seems, to me, to be a trivial mistake. Thus, I conclude the misunderstanding is likely to be mine and not the authors. (e.g perhaps he was using a trick I am unaware of that may have worked in c 20 years ago...).

The first segfault happens here, with an attempt to access diag[-1]:

float *diag;      /* ptr to matrix diagonal array */
for (i = 0; i < num_pnts - 1; i++){
    diag[i-1] = x[i+1] - x[i];
    assert(diag[i-1] > 0);}

And, if I fix the above for loop to start at1, then another segfault will occur at the next loop, for a similir reason, (diag[i-2]-->diag[-1])

for (i = 1; i < num_pnts - 1; i++)
right[i-1] = 6.0 * ((y[i+1]-y[i])/diag[i-1]-(y[i]-y[i-1])/diag[i-2]);

A link to the full source code is here

So my question is, this algorithm is published by a very accomplished individual, not a first year CS student. So while these look like errors to me (and indeed, segfault on a 2019 system), is there something else going on here I am missing? e.g would this OOB error have had a different result with a 1996 c compiler or something?

Question: are these actually mistakes?

Mtl Dev
  • 1,604
  • 20
  • 29
  • 1
    I stopped counting after 6 outright bugs and highly questionable - even for the mid 90's - bits of code in the first listing alone. I wouldn't use this as an example of how to write C. – Shawn Aug 14 '19 at 20:52
  • Thanks,your comments and second eyes appreciated. Strange as the publication seemed so well written and presented! – Mtl Dev Aug 14 '19 at 22:15

1 Answers1

3

Yes, those are actual mistakes. The array diag is allocated using malloc(). The pointer is to the start of the allocated region. By evaluating diag[i - 1] when i equals 0, it access outside the bounds of the allocated array.

It looks like the code meant to have the for loop start from 1, as you correctly tried. Indeed, then diag[i - 2] is also wrong. Since the main_diag[] array is filled starting from index 0, I think that was the intent for the secondary diagonal as well. So I think the second and third loops should have been:

for (i = 0; i < num_pnts - 1; i++) {
  diag[i] = x[i + 1] - x[i];
  assert(diag[i] > 0);
}

/* compute right hand side of equation */
for (i = 1; i < num_pnts - 1; i++) {
  right[i - 1] = 6.0 * ((y[i + 1] - y[i]) / diag[i] - (y[i] - y[i - 1]) / diag[i - 1]);
}

Accessing out of bounds doesn't necessarily result in a segmentation fault, it merely is undefined behaviour. It depends on the compiler and the implementation of malloc() whether memory right outside the bounds of the allocated memory is valid or not. And even if the memory was accessible, diag[i - 2] is on the right hand of a division operator, so if the memory at that location contained a zero, it would have caused a division by zero. So it might have crashed with a compiler from 1999 just as well.

It would be best if you can figure out what formulas were supposed to be implemented, and then check whether the code is correct or not.

G. Sliepen
  • 7,637
  • 1
  • 15
  • 31
  • I think I would have gone with uniform `for (i = 1; i < num_pnts - 1; i++)` for the loop controls. – Jonathan Leffler Aug 14 '19 at 21:07
  • Thanks. It seems rather strange indeed that such an (apparently) well thought, and and well documented, algorithm would be accompanied by such bad coding? – Mtl Dev Aug 14 '19 at 21:41
  • Everyone can make mistakes, even people two orders of magnitude smarter than you or me. The coding style actually looks bad in general now, but was very typical of the late nineties. It's what you would find in widely published books like https://en.wikipedia.org/wiki/Numerical_Recipes. Very condensed. Combined with the fact that compilers were not that strict 20 years ago, it might have actually run and produced a reasonable result on the author's computer. – G. Sliepen Aug 15 '19 at 19:26