#include <stdio.h>
#include <math.h>

#define N 4
#define MAX_ITER 150
#define EPS 1e-7

void print_system(double A[N][N], double b[N]);
void jacobi_solver(double A[N][N], double b[N]);
void gauss_seidel_solver(double A[N][N], double b[N]);
void verify_solution(double A[N][N], double b[N], double x[N]);

int main() {
    double A[N][N] = {
        { 5.0, -1.0, -2.0,  0.0},
        { 0.0,  8.0,  1.0,  7.0},
        {-2.0, -7.0, -4.0,  0.0},
        { 4.0,  0.0, -3.0, -25.0}
    };
    double b[N] = {-1.32, 14.41, 0.70, -0.03};

    print_system(A, b);
    
    printf("\n======================================================\n");
    printf("               1. ヤコビの反復法\n");
    printf("======================================================\n");
    jacobi_solver(A, b);

    printf("\n======================================================\n");
    printf("               2. ガウス・ザイデル法\n");
    printf("======================================================\n");
    gauss_seidel_solver(A, b);

    return 0;
}

void print_system(double A[N][N], double b[N]) {
    printf("対角優位化した連立一次方程式:\n");
    for (int i = 0; i < N; i++) {
        printf("[ ");
        for (int j = 0; j < N; j++) {
            printf("%6.1f ", A[i][j]);
        }
        printf("] [x%d] = [ %6.2f ]\n", i + 1, b[i]);
    }
}

void jacobi_solver(double A[N][N], double b[N]) {
    double x[N] = {0.0, 0.0, 0.0, 0.0};
    double x_new[N];
    int converged = 0;

    printf("計算回数\t残差2乗和\t\tx1\t\tx2\t\tx3\t\tx4\n");
    printf("------------------------------------------------------------------------\n");

    for (int k = 1; k <= MAX_ITER; k++) {
        for (int i = 0; i < N; i++) {
            double sum = 0.0;
            for (int j = 0; j < N; j++) {
                if (j != i) {
                    sum += A[i][j] * x[j];
                }
            }
            x_new[i] = (b[i] - sum) / A[i][i];
        }

        double rss = 0.0;
        for (int i = 0; i < N; i++) {
            rss += pow(x_new[i] - x[i], 2);
        }

        printf("%d\t\t%.2e\t\t%.5f\t%.5f\t%.5f\t%.5f\n", k, rss, x_new[0], x_new[1], x_new[2], x_new[3]);

        for (int i = 0; i < N; i++) {
            x[i] = x_new[i];
        }

        if (rss < EPS) {
            printf("\n>> 収束判定条件を満足しました（計算回数: %d回）\n", k);
            converged = 1;
            break;
        }
    }

    if (!converged) {
        printf("\n警告: 指定された計算回数以内に収束しませんでした。\n");
    }

    printf("\n=== 最終解 (ヤコビの反復法) ===\n");
    for (int i = 0; i < N; i++) {
        printf("x[%d] = %10.6f\n", i + 1, x[i]);
    }

    verify_solution(A, b, x);
}

void gauss_seidel_solver(double A[N][N], double b[N]) {
    double x[N] = {0.0, 0.0, 0.0, 0.0};
    double x_old[N];
    int converged = 0;

    printf("計算回数\t残差2乗和\t\tx1\t\tx2\t\tx3\t\tx4\n");
    printf("------------------------------------------------------------------------\n");

    for (int k = 1; k <= MAX_ITER; k++) {
        for (int i = 0; i < N; i++) {
            x_old[i] = x[i];
        }

        for (int i = 0; i < N; i++) {
            double sum = 0.0;
            for (int j = 0; j < N; j++) {
                if (j != i) {
                    sum += A[i][j] * x[j]; 
                }
            }
            x[i] = (b[i] - sum) / A[i][i];
        }

        double rss = 0.0;
        for (int i = 0; i < N; i++) {
            rss += pow(x[i] - x_old[i], 2);
        }

        printf("%d\t%.2e\t\t%.4f\t%.4f\t%.4f\t%.4f\n", k, rss, x[0], x[1], x[2], x[3]);

        if (rss < EPS) {
            printf("\n>> 収束判定条件を満足しました（計算: %d回）\n", k);
            converged = 1;
            break;
        }
    }

    if (!converged) {
        printf("\n警告: 指定された計算回数以内に収束しませんでした。\n");
    }

    printf("\n=== 最終解 (ガウス・ザイデル法) ===\n");
    for (int i = 0; i < N; i++) {
        printf("x[%d] = %10.5f\n", i + 1, x[i]);
    }

    verify_solution(A, b, x);
}

void verify_solution(double A[N][N], double b[N], double x[N]) {
    printf("\n--- コード内での解の確認（検算: Ax == b） ---\n");
    for (int i = 0; i < N; i++) {
        double lhs_calc = 0.0;
        for (int j = 0; j < N; j++) {
            lhs_calc += A[i][j] * x[j];
        }
        double error = lhs_calc - b[i];
        printf("%d行目: 計算値 Ax = %10.5f | 設定値 b = %10.5f | 誤差 = %e\n", 
               i + 1, lhs_calc, b[i], error);
    }
}