#include <stdio.h>
#include <stdlib.h>

// 计算n次拉格朗日插值多项式的系数
// 参数：n - 多项式次数，X - 插值节点数组，Y - 函数值数组，a - 输出系数数组
// 返回值：0表示成功，-1表示错误
int lagPolynomial(int n, double* X, double* Y, double* a) {
    if (n < 0 || X == NULL || Y == NULL || a == NULL) {
        return -1; // 参数检查
    }
    
    // 初始化系数数组为0
    for (int i = 0; i <= n; i++) {
        a[i] = 0.0;
    }
    
    // 临时数组，用于存储每个拉格朗日基函数的系数
    double* temp_coeff = (double*)malloc((n + 1) * sizeof(double));
    if (temp_coeff == NULL) {
        return -1; // 内存分配失败
    }
    
    // 对于每个基函数 l_k(x)
    for (int k = 0; k <= n; k++) {
        // 初始化当前基函数的系数
        for (int i = 0; i <= n; i++) {
            temp_coeff[i] = 0.0;
        }
        temp_coeff[0] = 1.0; // 常数项初始化为1
        
        double denominator = 1.0; // 分母
        
        // 构造基函数 l_k(x) = ∏(x - x_j) / (x_k - x_j), j≠k
        for (int j = 0; j <= n; j++) {
            if (j != k) {
                denominator *= (X[k] - X[j]);
                
                // 乘以 (x - x_j)，相当于多项式乘以一次项
                // 使用临时数组存储新的系数
                double* new_coeff = (double*)malloc((n + 1) * sizeof(double));
                for (int i = 0; i <= n; i++) {
                    new_coeff[i] = 0.0;
                }
                
                // 多项式乘法：temp_coeff * (x - x_j)
                for (int i = 0; i <= n; i++) {
                    if (i < n) {
                        new_coeff[i + 1] += temp_coeff[i]; // 乘以x
                    }
                    new_coeff[i] -= X[j] * temp_coeff[i]; // 乘以 -x_j
                }
                
                // 更新系数
                for (int i = 0; i <= n; i++) {
                    temp_coeff[i] = new_coeff[i];
                }
                
                free(new_coeff);
            }
        }
        
        // 将基函数乘以 y_k/denominator 并加到最终多项式
        for (int i = 0; i <= n; i++) {
            a[i] += Y[k] * temp_coeff[i] / denominator;
        }
    }
    
    free(temp_coeff);
    return 0;
}

// 测试函数
void testLagPolynomial() {
    // 测试用例1：线性插值 (n=1)
    printf("测试用例1：线性插值\n");
    int n1 = 1;
    double X1[] = {0.0, 1.0};
    double Y1[] = {1.0, 2.0};
    double a1[2];
    
    if (lagPolynomial(n1, X1, Y1, a1) == 0) {
        printf("系数: ");
        for (int i = 0; i <= n1; i++) {
            printf("a[%d]=%.6f ", i, a1[i]);
        }
        printf("\n");
        printf("多项式: L(x) = %.6f + %.6f*x\n", a1[0], a1[1]);
    }
    
    // 测试用例2：抛物线插值 (n=2)
    printf("\n测试用例2：抛物线插值\n");
    int n2 = 2;
    double X2[] = {100.0, 121.0, 144.0};
    double Y2[] = {10.0, 11.0, 12.0};
    double a2[3];
    
    if (lagPolynomial(n2, X2, Y2, a2) == 0) {
        printf("系数: ");
        for (int i = 0; i <= n2; i++) {
            printf("a[%d]=%.6f ", i, a2[i]);
        }
        printf("\n");
        printf("多项式: L(x) = %.6f + %.6f*x + %.6f*x^2\n", a2[0], a2[1], a2[2]);
        
        // 计算 sqrt(119) 的近似值
        double x = 119.0;
        double result = a2[0] + a2[1] * x + a2[2] * x * x;
        printf("L(119) = %.6f, 真实值 sqrt(119) ≈ 10.9087\n", result);
    }
    
    // 测试用例3：三次插值 (n=3)
    printf("\n测试用例3：三次插值\n");
    int n3 = 3;
    double X3[] = {0.0, 1.0, 2.0, 3.0};
    double Y3[] = {1.0, 2.0, 1.0, 4.0};
    double a3[4];
    
    if (lagPolynomial(n3, X3, Y3, a3) == 0) {
        printf("系数: ");
        for (int i = 0; i <= n3; i++) {
            printf("a[%d]=%.6f ", i, a3[i]);
        }
        printf("\n");
        printf("多项式: L(x) = %.6f + %.6f*x + %.6f*x^2 + %.6f*x^3\n", 
               a3[0], a3[1], a3[2], a3[3]);
    }
}

int main() {
    printf("拉格朗日插值多项式系数计算测试\n");
    printf("==============================\n");
    
    testLagPolynomial();
    
    return 0;
}