#include <stdio.h>  
#include <stdlib.h>  
  
#define N 100  
  
// 声明函数原型  
void spline(double x[], double y[], double b[], double c[], double d[], int n);  
double cubicSplineEval(double x[], double y[], double b[], double c[], double d[], int n, double xval);  
  
int main() {  
    double x[N], y[N], b[N], c[N], d[N];  
    int i;  
  
    // 初始化数据点  
    for (i = 0; i < N; i++) {  
        x[i] = i;  
        y[i] = x[i] * x[i] - x[i] + 3;  // 示例函数：f(x) = x^2 - x + 3  
    }  
  
    // 计算三次样条插值的系数  
    spline(x, y, b, c, d, N);  
  
    // 输出插值结果  
    printf("Cubic Spline Interpolation results:\n");  
    for (i = 0; i < N; i += 1) {  
        printf("f(%5.2f) = %8.4f\n", x[i], cubicSplineEval(x, y, b, c, d, N, x[i]));  
    }  
  
    return 0;  
}  
  
// 计算三次样条插值的系数  
void spline(double x[], double y[], double b[], double c[], double d[], int n) {  
    double h;  
    int i;  
    double *a = malloc((n - 1) * sizeof(double));  
  
    // 计算斜率差  
    for (i = 0; i < n - 1; i++) {  
        h = x[i + 1] - x[i];  
        a[i] = (3 / h) * (y[i + 1] - y[i]) - (3 / h) * (y[i] - y[i - 1]);  
    }  
  
    // 求解三对角线性方程组  
    c[0] = 0.0;  
    for (i = 1; i < n - 1; i++) {  
        double l = 2 * (x[i + 1] - x[i - 1]) - (x[i] - x[i - 1]) * c[i - 1];  
        c[i] = (a[i] - (x[i] - x[i - 1]) * c[i - 1]) / l;  
        b[i] = a[i + 1] - (x[i + 2] - x[i + 1]) * c[i + 1];  
        b[i] /= x[i + 2] - x[i];  
        d[i] = (c[i + 1] - c[i]) / (3 * (x[i + 2] - x[i]));  
    }  
  
    // 边界条件  
    b[n - 1] = 0.0;  
    c[n - 1] = 0.0;  
    d[n - 1] = 0.0;  
  
    free(a);  
}  
  
// 使用三次样条插值计算给定点的函数值  
double cubicSplineEval(double x[], double y[], double b[], double c[], double d[], int n, double xval) {  
    int i = 0, j = n - 1;  
    while (i < j) {  
        int m = (i + j) / 2;  
        if (xval < x[m])  
            j = m;  
        else  
            i = m + 1;  
    }  
    double dx = xval - x[i];  
    return y[i] + b[i] * dx + c[i] * dx * dx + d[i] * dx * dx * dx;  
}

/*
三次样条插值应该保证插值曲线穿过所有的数据节点。如果在计算插值系数时没有得到这样的结果，那可能是由于算法实现中的错误。下面我将提供一个优化后的计算插值系数的算法，确保插值曲线穿过所有节点。

首先，我们需要明确三次样条插值的基本形式。对于每个区间 [x 
i

 ,x 
i+1

 ]，插值多项式 S 
i

 (x) 可以表示为：

S 
i

 (x)=a 
i

 +b 
i

 (x−x 
i

 )+c 
i

 (x−x 
i

 ) 
2
 +d 
i

 (x−x 
i

 ) 
3
 

其中，a 
i

 =y 
i

 ，以确保曲线穿过节点 (x 
i

 ,y 
i

 )。

接下来，我们需要确定系数 b 
i

 ,c 
i

 ,d 
i

 。这可以通过求解以下方程组来实现：

S 
i

 (x 
i

 )=y 
i

 （已经满足）
S 
i

 (x 
i+1

 )=y 
i+1

 
S 
i
′

 (x 
i+1

 )=S 
i+1
′

 (x 
i+1

 )（一阶导数连续）
S 
i
′′

 (x 
i+1

 )=S 
i+1
′′

 (x 
i+1

 )（二阶导数连续）
对于第一个和第二个方程，我们已经知道 a 
i

 =y 
i

  和 a 
i+1

 =y 
i+1

 。对于第三个和第四个方程，我们需要计算一阶和二阶导数，并设置它们相等。
*/