Bug in cubic spline interpolator?

[swharden]

@Marco-Sanguinetti recently commented the following on a blog post of mine swharden/SWHarden.com#27

Now arrays A, B, and r are defined for index "n-1" but array C is not. However, C[n-1] is required in the indicated for loop.

The referenced code is also in ScottPlot, here

ScottPlot/src/ScottPlot4/ScottPlot/Statistics/Interpolation/Cubic.cs

Lines 99 to 106 in 093994e

	A[n - 1] = 1.0f / dx1;
	B[n - 1] = 2.0f * A[n - 1];
	r[n - 1] = 3 * (dy1 / (dx1 * dx1));
	double[] cPrime = new double[n];
	cPrime[0] = C[0] / B[0];
	for (int i = 1; i < n; i++)
	cPrime[i] = C[i] / (B[i] - cPrime[i - 1] * A[i]);

This issue will track better understanding this issue and coming up with a fix. I'll ping @drolevar in case they're interested in providing feedback, but they recently contributed some great related work in #3623 #3629

[swharden]

Adding more context, I think the issue may be that C[n-1] never gets set in the loop above, so it's always 0

ScottPlot/src/ScottPlot4/ScottPlot/Statistics/Interpolation/Cubic.cs

Lines 85 to 95 in 093994e

	for (int i = 1; i < n - 1; i++)
	{
	dx1 = x[i] - x[i - 1];
	dx2 = x[i + 1] - x[i];
	A[i] = 1.0f / dx1;
	C[i] = 1.0f / dx2;
	B[i] = 2.0f * (A[i] + C[i]);
	dy1 = y[i] - y[i - 1];
	dy2 = y[i + 1] - y[i];
	r[i] = 3 * (dy1 / (dx1 * dx1) + dy2 / (dx2 * dx2));
	}

Perhaps the fix is to add:

C[n - 1] = 1.0f / (x[n] - x[n-1]);