#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("回数  残差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%.2e\t%.4f\t\t%.4f\t%.4ft\t%.4f\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);
    }
}
 

I2luY2x1ZGUgPHN0ZGlvLmg+CiNpbmNsdWRlIDxtYXRoLmg+CgojZGVmaW5lIE4gNAojZGVmaW5lIE1BWF9JVEVSIDE1MAojZGVmaW5lIEVQUyAxZS03Cgp2b2lkIHByaW50X3N5c3RlbShkb3VibGUgQVtOXVtOXSwgZG91YmxlIGJbTl0pOwp2b2lkIGphY29iaV9zb2x2ZXIoZG91YmxlIEFbTl1bTl0sIGRvdWJsZSBiW05dKTsKdm9pZCBnYXVzc19zZWlkZWxfc29sdmVyKGRvdWJsZSBBW05dW05dLCBkb3VibGUgYltOXSk7CnZvaWQgdmVyaWZ5X3NvbHV0aW9uKGRvdWJsZSBBW05dW05dLCBkb3VibGUgYltOXSwgZG91YmxlIHhbTl0pOwoKaW50IG1haW4oKSB7CiAgICBkb3VibGUgQVtOXVtOXSA9IHsKICAgICAgICB7IDUuMCwgLTEuMCwgLTIuMCwgIDAuMH0sCiAgICAgICAgeyAwLjAsICA4LjAsICAxLjAsICA3LjB9LAogICAgICAgIHstMi4wLCAtNy4wLCAtNC4wLCAgMC4wfSwKICAgICAgICB7IDQuMCwgIDAuMCwgLTMuMCwgLTI1LjB9CiAgICB9OwogICAgZG91YmxlIGJbTl0gPSB7LTEuMzIsIDE0LjQxLCAwLjcwLCAtMC4wM307CgogICAgcHJpbnRfc3lzdGVtKEEsIGIpOwogICAgCiAgICBwcmludGYoIlxuPT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09XG4iKTsKICAgIHByaW50ZigiICAgICAgICAgICAgICAgMS4g44Ok44Kz44OT44Gu5Y+N5b6p5rOVXG4iKTsKICAgIHByaW50ZigiPT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09XG4iKTsKICAgIGphY29iaV9zb2x2ZXIoQSwgYik7CgogICAgcHJpbnRmKCJcbj09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PVxuIik7CiAgICBwcmludGYoIiAgICAgICAgICAgICAgIDIuIOOCrOOCpuOCueODu+OCtuOCpOODh+ODq+azlVxuIik7CiAgICBwcmludGYoIj09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PVxuIik7CiAgICBnYXVzc19zZWlkZWxfc29sdmVyKEEsIGIpOwoKICAgIHJldHVybiAwOwp9Cgp2b2lkIHByaW50X3N5c3RlbShkb3VibGUgQVtOXVtOXSwgZG91YmxlIGJbTl0pIHsKICAgIHByaW50Zigi5a++6KeS5YSq5L2N5YyW44GX44Gf6YCj56uL5LiA5qyh5pa556iL5byPOlxuIik7CiAgICBmb3IgKGludCBpID0gMDsgaSA8IE47IGkrKykgewogICAgICAgIHByaW50ZigiWyAiKTsKICAgICAgICBmb3IgKGludCBqID0gMDsgaiA8IE47IGorKykgewogICAgICAgICAgICBwcmludGYoIiU2LjFmICIsIEFbaV1bal0pOwogICAgICAgIH0KICAgICAgICBwcmludGYoIl0gW3glZF0gPSBbICU2LjJmIF1cbiIsIGkgKyAxLCBiW2ldKTsKICAgIH0KfQoKdm9pZCBqYWNvYmlfc29sdmVyKGRvdWJsZSBBW05dW05dLCBkb3VibGUgYltOXSkgewogICAgZG91YmxlIHhbTl0gPSB7MC4wLCAwLjAsIDAuMCwgMC4wfTsKICAgIGRvdWJsZSB4X25ld1tOXTsKICAgIGludCBjb252ZXJnZWQgPSAwOwoKICAgIHByaW50Zigi5Zue5pWwICDmrovlt64y5LmX5ZKMXHRcdHgxXHRcdHgyXHRcdHgzXHRcdHg0XG4iKTsKICAgIHByaW50ZigiLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tXG4iKTsKCiAgICBmb3IgKGludCBrID0gMTsgayA8PSBNQVhfSVRFUjsgaysrKSB7CiAgICAgICAgZm9yIChpbnQgaSA9IDA7IGkgPCBOOyBpKyspIHsKICAgICAgICAgICAgZG91YmxlIHN1bSA9IDAuMDsKICAgICAgICAgICAgZm9yIChpbnQgaiA9IDA7IGogPCBOOyBqKyspIHsKICAgICAgICAgICAgICAgIGlmIChqICE9IGkpIHsKICAgICAgICAgICAgICAgICAgICBzdW0gKz0gQVtpXVtqXSAqIHhbal07CiAgICAgICAgICAgICAgICB9CiAgICAgICAgICAgIH0KICAgICAgICAgICAgeF9uZXdbaV0gPSAoYltpXSAtIHN1bSkgLyBBW2ldW2ldOwogICAgICAgIH0KCiAgICAgICAgZG91YmxlIHJzcyA9IDAuMDsKICAgICAgICBmb3IgKGludCBpID0gMDsgaSA8IE47IGkrKykgewogICAgICAgICAgICByc3MgKz0gcG93KHhfbmV3W2ldIC0geFtpXSwgMik7CiAgICAgICAgfQoKICAgICAgICBwcmludGYoIiVkICBcdCUuMmVcdCUuNGZcdFx0JS40Zlx0JS40ZnRcdCUuNGZcbiIsIGssIHJzcywgeF9uZXdbMF0sIHhfbmV3WzFdLCB4X25ld1syXSwgeF9uZXdbM10pOwoKICAgICAgICBmb3IgKGludCBpID0gMDsgaSA8IE47IGkrKykgewogICAgICAgICAgICB4W2ldID0geF9uZXdbaV07CiAgICAgICAgfQoKICAgICAgICBpZiAocnNzIDwgRVBTKSB7CiAgICAgICAgICAgIHByaW50ZigiXG4+PiDlj47mnZ/liKTlrprmnaHku7bjgpLmuoDotrPjgZfjgb7jgZfjgZ/vvIjoqIjnrpflm57mlbA6ICVk5Zue77yJXG4iLCBrKTsKICAgICAgICAgICAgY29udmVyZ2VkID0gMTsKICAgICAgICAgICAgYnJlYWs7CiAgICAgICAgfQogICAgfQoKICAgIGlmICghY29udmVyZ2VkKSB7CiAgICAgICAgcHJpbnRmKCJcbuitpuWRijog5oyH5a6a44GV44KM44Gf6KiI566X5Zue5pWw5Lul5YaF44Gr5Y+O5p2f44GX44G+44Gb44KT44Gn44GX44Gf44CCXG4iKTsKICAgIH0KCiAgICBwcmludGYoIlxuPT09IOacgOe1guinoyAo44Ok44Kz44OT44Gu5Y+N5b6p5rOVKSA9PT1cbiIpOwogICAgZm9yIChpbnQgaSA9IDA7IGkgPCBOOyBpKyspIHsKICAgICAgICBwcmludGYoInhbJWRdID0gJTEwLjZmXG4iLCBpICsgMSwgeFtpXSk7CiAgICB9CgogICAgdmVyaWZ5X3NvbHV0aW9uKEEsIGIsIHgpOwp9Cgp2b2lkIGdhdXNzX3NlaWRlbF9zb2x2ZXIoZG91YmxlIEFbTl1bTl0sIGRvdWJsZSBiW05dKSB7CiAgICBkb3VibGUgeFtOXSA9IHswLjAsIDAuMCwgMC4wLCAwLjB9OwogICAgZG91YmxlIHhfb2xkW05dOwogICAgaW50IGNvbnZlcmdlZCA9IDA7CgogICAgcHJpbnRmKCLoqIjnrpflm57mlbBcdOaui+W3rjLkuZflkoxcdFx0eDFcdFx0eDJcdFx0eDNcdFx0eDRcbiIpOwogICAgcHJpbnRmKCItLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS1cbiIpOwoKICAgIGZvciAoaW50IGsgPSAxOyBrIDw9IE1BWF9JVEVSOyBrKyspIHsKICAgICAgICBmb3IgKGludCBpID0gMDsgaSA8IE47IGkrKykgewogICAgICAgICAgICB4X29sZFtpXSA9IHhbaV07CiAgICAgICAgfQoKICAgICAgICBmb3IgKGludCBpID0gMDsgaSA8IE47IGkrKykgewogICAgICAgICAgICBkb3VibGUgc3VtID0gMC4wOwogICAgICAgICAgICBmb3IgKGludCBqID0gMDsgaiA8IE47IGorKykgewogICAgICAgICAgICAgICAgaWYgKGogIT0gaSkgewogICAgICAgICAgICAgICAgICAgIHN1bSArPSBBW2ldW2pdICogeFtqXTsgCiAgICAgICAgICAgICAgICB9CiAgICAgICAgICAgIH0KICAgICAgICAgICAgeFtpXSA9IChiW2ldIC0gc3VtKSAvIEFbaV1baV07CiAgICAgICAgfQoKICAgICAgICBkb3VibGUgcnNzID0gMC4wOwogICAgICAgIGZvciAoaW50IGkgPSAwOyBpIDwgTjsgaSsrKSB7CiAgICAgICAgICAgIHJzcyArPSBwb3coeFtpXSAtIHhfb2xkW2ldLCAyKTsKICAgICAgICB9CgogICAgICAgIHByaW50ZigiJWRcdCUuMmVcdFx0JS40Zlx0JS40Zlx0JS40Zlx0JS40ZlxuIiwgaywgcnNzLCB4WzBdLCB4WzFdLCB4WzJdLCB4WzNdKTsKCiAgICAgICAgaWYgKHJzcyA8IEVQUykgewogICAgICAgICAgICBwcmludGYoIlxuPj4g5Y+O5p2f5Yik5a6a5p2h5Lu244KS5rqA6Laz44GX44G+44GX44Gf77yI6KiI566XOiAlZOWbnu+8iVxuIiwgayk7CiAgICAgICAgICAgIGNvbnZlcmdlZCA9IDE7CiAgICAgICAgICAgIGJyZWFrOwogICAgICAgIH0KICAgIH0KCiAgICBpZiAoIWNvbnZlcmdlZCkgewogICAgICAgIHByaW50ZigiXG7orablkYo6IOaMh+WumuOBleOCjOOBn+ioiOeul+WbnuaVsOS7peWGheOBq+WPjuadn+OBl+OBvuOBm+OCk+OBp+OBl+OBn+OAglxuIik7CiAgICB9CgogICAgcHJpbnRmKCJcbj09PSDmnIDntYLop6MgKOOCrOOCpuOCueODu+OCtuOCpOODh+ODq+azlSkgPT09XG4iKTsKICAgIGZvciAoaW50IGkgPSAwOyBpIDwgTjsgaSsrKSB7CiAgICAgICAgcHJpbnRmKCJ4WyVkXSA9ICUxMC41ZlxuIiwgaSArIDEsIHhbaV0pOwogICAgfQoKICAgIHZlcmlmeV9zb2x1dGlvbihBLCBiLCB4KTsKfQoKdm9pZCB2ZXJpZnlfc29sdXRpb24oZG91YmxlIEFbTl1bTl0sIGRvdWJsZSBiW05dLCBkb3VibGUgeFtOXSkgewogICAgcHJpbnRmKCJcbi0tLSDjgrPjg7zjg4nlhoXjgafjga7op6Pjga7norroqo3vvIjmpJznrpc6IEF4ID09IGLvvIkgLS0tXG4iKTsKICAgIGZvciAoaW50IGkgPSAwOyBpIDwgTjsgaSsrKSB7CiAgICAgICAgZG91YmxlIGxoc19jYWxjID0gMC4wOwogICAgICAgIGZvciAoaW50IGogPSAwOyBqIDwgTjsgaisrKSB7CiAgICAgICAgICAgIGxoc19jYWxjICs9IEFbaV1bal0gKiB4W2pdOwogICAgICAgIH0KICAgICAgICBkb3VibGUgZXJyb3IgPSBsaHNfY2FsYyAtIGJbaV07CiAgICAgICAgcHJpbnRmKCIlZOihjOebrjog6KiI566X5YCkIEF4ID0gJTEwLjVmIHwg6Kit5a6a5YCkIGIgPSAlMTAuNWYgfCDoqqTlt64gPSAlZVxuIiwgCiAgICAgICAgICAgICAgIGkgKyAxLCBsaHNfY2FsYywgYltpXSwgZXJyb3IpOwogICAgfQp9Cg==