C语言实现三次样条插值

算法推导过程参考：
https://wenku.baidu.com/view/88c5af9f185f312b3169a45177232f60dccce772.html?rec_flag=default&fr=pc_newview_relate-1001_1-3-wk_rec_doc2-1001_1-3-88c5af9f185f312b3169a45177232f60dccce772&sxts=1609484456691

#include <stdio.h>


# define MAX_N 20             // 定义(x_i,y_i)的最大的维数  

typedef struct tagPOINT         // 插值点的结构体成员有x, y
{
	double x;
	double y;
} POINT;

POINT points[MAX_N + 1];
//需要拟合的六个点（x0 y0 x1 y1...）
double passData[12] = { 0, 0, 1, 2.15, 2, 12.1, 3, 32.9, 4, 42.87, 5, 45 };

int i, j, k, n, num;
double h[MAX_N + 1], b[MAX_N + 1], c[MAX_N + 1], d[MAX_N + 1], M[MAX_N + 1];
double u[MAX_N + 1], v[MAX_N + 1], y[MAX_N + 1];

//计算M矩阵
void CaculateM(void) {
	num = sizeof(passData) / sizeof(passData[0]);	//获取数组元素个数
	n = num / 2;	//获取x个数
	M[0] = M[n] = 0;	//自然边界

	for (i = 0; i < num; i += 2) {
		points[j].x = passData[i];
		points[j].y = passData[i + 1];
		j++;
	}
	// 计算M关系式中各参数的值  
	h[0] = points[1].x - points[0].x;
	for (i = 1; i<n; i++) {
		h[i] = points[i + 1].x - points[i].x;
		b[i] = h[i] / (h[i] + h[i - 1]);
		c[i] = 1 - b[i];
		d[i] = 6 * ((points[i + 1].y - points[i].y) / h[i] - (points[i].y - points[i - 1].y) / h[i - 1]) / (h[i] + h[i - 1]);
	}
	// 用追赶法计算Mi,i=1,…,n-1  
	d[1] -= c[1] * M[0];
	d[n - 1] -= b[n - 1] * M[n];
	b[n - 1] = 0; c[1] = 0; v[0] = 0;
	for (i = 1; i<n; i++) {
		u[i] = 2 - c[i] * v[i - 1];
		v[i] = b[i] / u[i];
		y[i] = (d[i] - c[i] * y[i - 1]) / u[i];
	}
	for (i = 1; i<n; i++) {
		M[n - i] = y[n - i] - v[n - i] * M[n - i + 1];
	}
}

//计算要插值的x对应的y值
double Spline(double spline_x) {
	double p, q, S;
	M[0] = M[n] = 0;

	while (spline_x >= points[k].x) k++;
	k = k - 1;
	p = points[k + 1].x - spline_x;
	q = spline_x - points[k].x;
	S = (p* p* p* M[k] + q* q* q* M[k + 1]) / (6 * h[k]) + (p* points[k].y + q*points[k + 1].y) / h[k] - h[k] * (p* M[k] + q* M[k + 1]) / 6;

	return S;
}

int main() {
	double ret;

	CaculateM();
	ret = Spline(3);
	printf("result = %f\n", ret);
	// 输出 

	system("pause");
}