Взял для примера подпрограмму IAC1R1 на C и постарался, не меняя алгоритм, упростить её как можно больше. Например, выбросил операторы, обеспечивающие совместимость индексов с фортрановской нумерацией. Так, чтобы можно было легко увидеть, что, собственно, там делается. Вот что получилось.
double iac1r1_c(double *t, double *c, int lx, int j)
{
if (lx==0)
return c[j];
if (lx==1)
{
double
p = t[j+1] - t[j-1],
q = t[j+1] - t[j ],
r = t[j ] - t[j-1];
return (p*p/(q*r) * c[j] - q*q/(r*p) * c[j-1] - r*r/(q*p) * c[j+1]) / 3.;
}
if (lx>1)
return (438*c[j] - 112*(c[j-1] + c[j+1]) + 13*(c[j-2] + c[j+2])) / 240.;
}
Там есть ветвь, которая при LX
вообще никогда не должна выполняться. И такое впечатление, что три ветви, соответствующие LX=0, LX=1 и LX>1, реализуют три совершенно разных подхода.