#include <stdio.h>  
#include <stdlib.h>  
#include <math.h>  
  
#define MAX_SIZE 20  
#define EPS 1e-12  
  
// 高斯消去法（列主元）  
int gauss_elimination(int n, double A[MAX_SIZE][MAX_SIZE], double b[MAX_SIZE], double x[MAX_SIZE]) {  
    int i, j, k;  
    double max_val, temp, factor;  
    int max_row;  
  
    // 创建矩阵和向量的副本，避免修改原始数据  
    double A_copy[MAX_SIZE][MAX_SIZE];  
    double b_copy[MAX_SIZE];  
  
    for (i = 0; i < n; i++) {  
        for (j = 0; j < n; j++) {  
            A_copy[i][j] = A[i][j];  
        }  
        b_copy[i] = b[i];  
    }  
  
    // 消元过程  
    for (k = 0; k < n - 1; k++) {  
        // 列主元选取  
        max_val = fabs(A_copy[k][k]);  
        max_row = k;  
        for (i = k + 1; i < n; i++) {  
            if (fabs(A_copy[i][k]) > max_val) {  
                max_val = fabs(A_copy[i][k]);  
                max_row = i;  
            }  
        }  
  
        // 检查主元是否接近0  
        if (max_val < EPS) {  
            printf("矩阵奇异，无法求解\n");  
            return -1;  
        }  
  
        // 交换行  
        if (max_row != k) {  
            // 交换系数矩阵行  
            for (j = k; j < n; j++) {  
                temp = A_copy[k][j];  
                A_copy[k][j] = A_copy[max_row][j];  
                A_copy[max_row][j] = temp;  
            }  
            // 交换常数向量元素  
            temp = b_copy[k];  
            b_copy[k] = b_copy[max_row];  
            b_copy[max_row] = temp;  
        }  
  
        // 消元  
        for (i = k + 1; i < n; i++) {  
            factor = A_copy[i][k] / A_copy[k][k];  
            for (j = k; j < n; j++) {  
                A_copy[i][j] -= factor * A_copy[k][j];  
            }  
            b_copy[i] -= factor * b_copy[k];  
        }  
    }  
  
    // 检查最后一个主元  
    if (fabs(A_copy[n-1][n-1]) < EPS) {  
        printf("矩阵奇异，无法求解\n");  
        return -1;  
    }  
  
    // 回代过程  
    x[n-1] = b_copy[n-1] / A_copy[n-1][n-1];  
    for (i = n - 2; i >= 0; i--) {  
        temp = b_copy[i];  
        for (j = i + 1; j < n; j++) {  
            temp -= A_copy[i][j] * x[j];  
        }  
        x[i] = temp / A_copy[i][i];  
    }  
  
    return 0;  
}  
  
// 生成Hilbert矩阵  
void generate_hilbert_matrix(int n, double A[MAX_SIZE][MAX_SIZE]) {  
    int i, j;  
    for (i = 0; i < n; i++) {  
        for (j = 0; j < n; j++) {  
            A[i][j] = 1.0 / (i + j + 1.0);  
        }  
    }  
}  
  
// 生成全1向量  
void generate_ones_vector(int n, double b[MAX_SIZE]) {  
    int i;  
    for (i = 0; i < n; i++) {  
        b[i] = 1.0;  
    }  
}  
  
// 计算残差（Ax - b的最大绝对值）  
double calculate_residual(int n, double A[MAX_SIZE][MAX_SIZE], double b[MAX_SIZE], double x[MAX_SIZE]) {  
    double max_residual = 0.0;  
    int i, j;  
  
    for (i = 0; i < n; i++) {  
        double residual = 0.0;  
        for (j = 0; j < n; j++) {  
            residual += A[i][j] * x[j];  
        }  
        residual = fabs(residual - b[i]);  
        if (residual > max_residual) {  
            max_residual = residual;  
        }  
    }  
    return max_residual;  
}  
  
// 打印矩阵  
void print_matrix(int n, double A[MAX_SIZE][MAX_SIZE]) {  
    int i, j;  
    for (i = 0; i < n; i++) {  
        for (j = 0; j < n; j++) {  
            printf("%10.6f ", A[i][j]);  
        }  
        printf("\n");  
    }  
}  
  
// 打印向量  
void print_vector(int n, double v[MAX_SIZE]) {  
    int i;  
    for (i = 0; i < n; i++) {  
        printf("x[%d] = %15.10f\n", i, v[i]);  
    }  
}  
  
int main() {  
    int n, i;  
    double A[MAX_SIZE][MAX_SIZE];  
    double b[MAX_SIZE];  
    double x[MAX_SIZE];  
    double residual;  
  
    printf("高斯消去法测试 - Hilbert矩阵\n");  
    printf("==============================\n\n");  
  
    // 测试5阶Hilbert矩阵  
    n = 5;  
    printf("1. %d阶Hilbert矩阵测试:\n", n);  
    generate_hilbert_matrix(n, A);  
    generate_ones_vector(n, b);  
  
    printf("Hilbert矩阵:\n");  
    print_matrix(n, A);  
  
    if (gauss_elimination(n, A, b, x) == 0) {  
        printf("\n解向量:\n");  
        print_vector(n, x);  
        residual = calculate_residual(n, A, b, x);  
        printf("最大残差: %e\n\n", residual);  
    }  
  
    // 测试10阶Hilbert矩阵  
    n = 10;  
    printf("2. %d阶Hilbert矩阵测试:\n", n);  
    generate_hilbert_matrix(n, A);  
    generate_ones_vector(n, b);  
  
    if (gauss_elimination(n, A, b, x) == 0) {  
        printf("\n解向量:\n");  
        print_vector(n, x);  
        residual = calculate_residual(n, A, b, x);  
        printf("最大残差: %e\n\n", residual);  
    }  
  
    // 测试15阶Hilbert矩阵  
    n = 15;  
    printf("3. %d阶Hilbert矩阵测试:\n", n);  
    generate_hilbert_matrix(n, A);  
    generate_ones_vector(n, b);  
  
    if (gauss_elimination(n, A, b, x) == 0) {  
        printf("\n解向量:\n");  
        print_vector(n, x);  
        residual = calculate_residual(n, A, b, x);  
        printf("最大残差: %e\n\n", residual);  
    } else {  
        printf("矩阵奇异，无法求解\n\n");  
    }  
  
    // 验证PDF中的例子  
    printf("4. 验证PDF中的例子:\n");  
    n = 3;  
    double test_A[3][3] = {  
            {2, -1, 3},  
            {4, 2, 5},  
            {1, 2, 0}  
    };  
    double test_b[3] = {1, 4, 7};  
  
    for (i = 0; i < n; i++) {  
        for (int j = 0; j < n; j++) {  
            A[i][j] = test_A[i][j];  
        }  
        b[i] = test_b[i];  
    }  
  
    printf("系数矩阵:\n");  
    print_matrix(n, A);  
    printf("常数向量: ");  
    for (i = 0; i < n; i++) {  
        printf("%.1f ", b[i]);  
    }  
    printf("\n");  
  
    if (gauss_elimination(n, A, b, x) == 0) {  
        printf("解向量:\n");  
        print_vector(n, x);  
        residual = calculate_residual(n, A, b, x);  
        printf("最大残差: %e\n", residual);  
        printf("预期解: x1=9, x2=-1, x3=-6\n");  
    }  
  
    return 0;  
}