这个发布的 Cubic Spline Extrema 代码有错误吗?

Does this published Cubic Spline Extrema code have errors?

我正在尝试实施 1996 年发布的 Cublic Spline Extrema Algorithm

我不是一个有成就的程序员,很明显作者是,而且比我聪明一两个数量级。但是当我尝试编译它时,它在一个简单的错误上出现了错误-边界数组访问。我发现很难相信这个算法是在我看来是一个微不足道的错误的情况下发布的。因此,我得出结论,误解很可能是我的,而不是作者的。 (例如,也许他使用了一个我不知道的技巧,它可能在 20 年前用 c 语言工作过......)。

第一个段错误发生在这里,试图访问 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);}

并且,如果我将上面的 for 循环修复为从 1 开始,那么出于类似的原因,下一个循环将发生另一个段错误,(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

所以我的问题是,这个算法是由一个非常有成就的人发布的,而不是一年级 CS 学生。因此,虽然这些对我来说看起来像是错误(实际上,2019 系统上的段错误),但我还缺少其他东西吗?例如,这个 OOB 错误在 1996 年的 c 编译器或其他东西上会产生不同的结果吗?

提问:这些真的是错误吗?

是的,这些都是实际的错误。数组 diag 是使用 malloc() 分配的。指针指向分配区域的开始。通过在 i 等于 0 时评估 diag[i - 1],它访问分配数组的边界之外。

看起来该代码旨在让 for 循环从 1 开始,正如您正确尝试的那样。确实,那么diag[i - 2]也是错误的。由于 main_diag[] 数组是从索引 0 开始填充的,我认为这也是辅助对角线的意图。所以我认为第二个和第三个循环应该是:

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]);
}

越界访问不一定会导致分段错误,它只是未定义的行为。这取决于编译器和 malloc() 的实现,分配内存范围之外的内存是否有效。即使内存是可访问的,diag[i - 2] 位于除法运算符的右侧,因此如果该位置的内存包含零,它会导致被零除。所以它也可能与 1999 年的编译器一起崩溃。

如果你能搞清楚应该实现什么公式,然后检查代码是否正确,那就最好了。