C 语言实现三次样条插值.

/*
System --> Linux & gcc8.1.0
File ----> NaturalCubicSpline.c
Author --> Illusionna
Create --> 2024/2/21 22:16:30
'''
-*- Encoding: UTF-8 -*-
*/

// 生成 Python 调用需要的动态链接库: gcc --share -o spline.dll NaturalCubicSpline.c

# include <stdio.h>
# include <stdlib.h>
# include <string.h>
# include <math.h>
# include <malloc.h>

double ** Spline(double * X, double * Y, int n){
    int i, j;
    n--;
    double stepArray[n], alphaArray[n], lArray[n+1], muArray[n+1], zArray[n+1], cArray[n+1], bArray[n], dArray[n];
    // Step 1.
    for (i=0; i<n; ++i){
        stepArray[i] = X[i+1] - X[i];
    }
    // Step 2.
    for (i=1; i<n; ++i){
        alphaArray[i] = (3 * (Y[i+1] - Y[i]) / stepArray[i]) - (3 * (Y[i] - Y[i-1]) / stepArray[i-1]);
    }
    // Step 3.
    lArray[0] = 1;
    muArray[0] = 0;
    zArray[0] = 0;
    // Step 4.
    for (i=1; i<n; ++i){
        lArray[i] = 2 * (X[i+1] - X[i-1]) - stepArray[i-1] * muArray[i-1];
        muArray[i] = stepArray[i] / lArray[i];
        zArray[i] = (alphaArray[i] - stepArray[i-1] * zArray[i-1]) / lArray[i];
    }
    // Step 5.
    lArray[n] = 1;
    zArray[n] = 0;
    cArray[n] = 0;
    // Step 6.
    for (j=n-1; j>=0; --j){
        cArray[j] = zArray[j] - muArray[j] * cArray[j+1];
        bArray[j] = ((Y[j+1] - Y[j]) / stepArray[j]) - (stepArray[j] * (cArray[j+1] + 2 * cArray[j]) / 3);
        dArray[j] = (cArray[j+1] - cArray[j]) / (3 * stepArray[j]);
    }
    // Information.
    printf("\033[036mThe coefficients of each interval cubic polynomial:\033[0m\n");
    printf("S(x)=a+b(x-xk)+c(x-xk)^2+d(x-xk)^3    x_{k} < x < x_{k+1}\n");
    printf("%5s %9s %8s %8s %8s\n", "xk", "a", "b", "c", "d");
    for (i=0; i<n; ++i){
        printf("%.5f %9.5f %8.5f %9.5f %9.5f\n", X[i], Y[i], bArray[i], cArray[i], dArray[i]);
    }
    double ** result;
    int columns = 5;
    result = (double**)malloc(n * sizeof(double*));
    for (i=0; i<n; ++i){
        result[i] = (double *)malloc(columns * sizeof(double));
        result[i][0] = X[i];
        result[i][1] = Y[i];
        result[i][2] = bArray[i];
        result[i][3] = cArray[i];
        result[i][4] = dArray[i];
    }
    return result;
}


int main(){
    system("");
    printf("\033[H\033[J\n");
    int countX = 1;
    int countY = 1;
    double X[100];
    int i = 0;
    char enter;
    printf("Input X vector: ");
    while (1){
        if (scanf("%lf", & X[i]) == 1){
            ++i;
        }
        else{
            break;
        }
        enter = getchar();
        if (enter == '\n'){
            break;
        }
        else{
            ++countX;
        }
    }
    for (int j=0; j<i; ++j){
        printf("%lf\n", X[j]);
    }
    i = 0;
    double Y[100];
    printf("Input Y vector: ");
    while (1){
        if (scanf("%lf", & Y[i]) == 1){
            ++i;
        }
        else{
            break;
        }
        enter = getchar();
        if (enter == '\n'){
            break;
        }
        else{
            ++countY;
        }
    }
    for (int j=0; j<i; ++j){
        printf("%lf\n", Y[j]);
    }
    if (countX != countY){
        printf("\033[31mVectors X and Y have different lengths.\033[0m");
    }
    else{
        double **cs;
        cs = Spline(X, Y, countX);
        printf("\nDo you want to interpolate? (\033[33mYes:Y, No:N\033[0m)\n");
        char judge[10];
        scanf("%s", judge);
        if ((strcmp(judge, "yes") == 0) || (strcmp(judge, "Yes") == 0) || (strcmp(judge, "Y") == 0) || (strcmp(judge, "y") == 0)){
            printf("You entered: yes\n");
            int count = 1;
            double interpolateX[100];
            i = 0;
            printf("Input interpolation vector or number:\033[33m ");
            while (1){
                if (scanf("%lf", & interpolateX[i]) == 1){
                    ++i;
                }
                else{
                    break;
                }
                enter = getchar();
                if (enter == '\n'){
                    break;
                }
                else{
                    ++count;
                }
            }
            printf("\033[0m");
            for (int j=0; j<i; ++j){
                printf("%lf\n", interpolateX[j]);
            }
            int pos = 0;
            double * predictions = (double*)malloc(count * sizeof(double));
            for (int s=0; s<count; ++s){
                for (int t=0; t<countX-1; ++t){
                    if (pos == 0){
                        if (interpolateX[s] < X[t]){
                            double result = 0;
                            double x = cs[0][0];
                            for (int k=1; k<5; ++k){
                                result = result + cs[0][k]*pow(interpolateX[s]-x, k-1);
                            }
                            predictions[s] = result;
                            pos = 1;
                        }
                        else if(interpolateX[s] >= X[countX-1]){
                            double result = 0;
                            double x = cs[countX-2][0];
                            for (int k=1; k<5; ++k){
                                result = result + cs[countX-2][k]*pow(interpolateX[s]-x, k-1);
                            }
                            predictions[s] = result;
                            pos = 1;
                        }
                        else if((X[t] <= interpolateX[s]) && (interpolateX[s] < X[t+1])){
                            double result = 0;
                            double x = cs[t][0];
                            for (int k=1; k<5; ++k){
                                result = result + cs[t][k]*pow(interpolateX[s]-x, k-1);
                            }
                            predictions[s] = result;
                            pos = 1;
                        }
                        else{

                        }
                    }
                }
                pos = 0;
            }
            printf("Prediction vector: \033[33m");
            for (int idx=0; idx<count; ++idx){
                printf("%f ", predictions[idx]);
            }
        }
        else if ((strcmp(judge, "no") == 0) || (strcmp(judge, "No") == 0) || (strcmp(judge, "N") == 0) || (strcmp(judge, "n") == 0)){

        }
        else{
            printf("Invalid input\n");
        }
        printf("\033[0m\n");
        system("pause");
    }
	return 0;
}