#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("計算回数k\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("計算回数k\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\t%.2e\t\t%.5f\t%.5f\t%.5f\t%.5f\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.6f\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+CgojZGVmaW5lIE4gNAojZGVmaW5lIE1BWF9JVEVSIDE1MAojZGVmaW5lIEVQUyAxZS03Cgp2b2lkIHByaW50X3N5c3RlbShkb3VibGUgQVtOXVtOXSwgZG91YmxlIGJbTl0pOwp2b2lkIGphY29iaV9zb2x2ZXIoZG91YmxlIEFbTl1bTl0sIGRvdWJsZSBiW05dKTsKdm9pZCBnYXVzc19zZWlkZWxfc29sdmVyKGRvdWJsZSBBW05dW05dLCBkb3VibGUgYltOXSk7CnZvaWQgdmVyaWZ5X3NvbHV0aW9uKGRvdWJsZSBBW05dW05dLCBkb3VibGUgYltOXSwgZG91YmxlIHhbTl0pOwoKaW50IG1haW4oKSB7CiAgICBkb3VibGUgQVtOXVtOXSA9IHsKICAgICAgICB7IDUuMCwgLTEuMCwgLTIuMCwgIDAuMH0sCiAgICAgICAgeyAwLjAsICA4LjAsICAxLjAsICA3LjB9LAogICAgICAgIHstMi4wLCAtNy4wLCAtNC4wLCAgMC4wfSwKICAgICAgICB7IDQuMCwgIDAuMCwgLTMuMCwgLTI1LjB9CiAgICB9OwogICAgZG91YmxlIGJbTl0gPSB7LTEuMzIsIDE0LjQxLCAwLjcwLCAtMC4wM307CgogICAgcHJpbnRfc3lzdGVtKEEsIGIpOwogICAgCiAgICBwcmludGYoIlxuPT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09XG4iKTsKICAgIHByaW50ZigiIDEuIOODpOOCs+ODk+OBruWPjeW+qeazlVxuIik7CiAgICBwcmludGYoIj09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PVxuIik7CiAgICBqYWNvYmlfc29sdmVyKEEsIGIpOwoKICAgIHByaW50ZigiXG49PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT1cbiIpOwogICAgcHJpbnRmKCIgMi4g44Ks44Km44K544O744K244Kk44OH44Or5rOVXG4iKTsKICAgIHByaW50ZigiPT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09XG4iKTsKICAgIGdhdXNzX3NlaWRlbF9zb2x2ZXIoQSwgYik7CgogICAgcmV0dXJuIDA7Cn0KCnZvaWQgcHJpbnRfc3lzdGVtKGRvdWJsZSBBW05dW05dLCBkb3VibGUgYltOXSkgewogICAgcHJpbnRmKCLlr77op5LlhKrkvY3ljJblvozjga7mlrnnqIvlvI86XG4iKTsKICAgIGZvciAoaW50IGkgPSAwOyBpIDwgTjsgaSsrKSB7CiAgICAgICAgcHJpbnRmKCJbICIpOwogICAgICAgIGZvciAoaW50IGogPSAwOyBqIDwgTjsgaisrKSB7CiAgICAgICAgICAgIHByaW50ZigiJTYuMWYgIiwgQVtpXVtqXSk7CiAgICAgICAgfQogICAgICAgIHByaW50ZigiXSBbeCVkXSA9IFsgJTYuMmYgXVxuIiwgaSArIDEsIGJbaV0pOwogICAgfQp9Cgp2b2lkIGphY29iaV9zb2x2ZXIoZG91YmxlIEFbTl1bTl0sIGRvdWJsZSBiW05dKSB7CiAgICBkb3VibGUgeFtOXSA9IHswLjAsIDAuMCwgMC4wLCAwLjB9OwogICAgZG91YmxlIHhfbmV3W05dOwogICAgaW50IGNvbnZlcmdlZCA9IDA7CgogICAgcHJpbnRmKCLoqIjnrpflm57mlbBrXHTmrovlt64y5LmX5ZKMXHRcdHgxXHRcdHgyXHRcdHgzXHRcdHg0XG4iKTsKICAgIHByaW50ZigiLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tXG4iKTsKCiAgICBmb3IgKGludCBrID0gMTsgayA8PSBNQVhfSVRFUjsgaysrKSB7CiAgICAgICAgZm9yIChpbnQgaSA9IDA7IGkgPCBOOyBpKyspIHsKICAgICAgICAgICAgZG91YmxlIHN1bSA9IDAuMDsKICAgICAgICAgICAgZm9yIChpbnQgaiA9IDA7IGogPCBOOyBqKyspIHsKICAgICAgICAgICAgICAgIGlmIChqICE9IGkpIHsKICAgICAgICAgICAgICAgICAgICBzdW0gKz0gQVtpXVtqXSAqIHhbal07CiAgICAgICAgICAgICAgICB9CiAgICAgICAgICAgIH0KICAgICAgICAgICAgeF9uZXdbaV0gPSAoYltpXSAtIHN1bSkgLyBBW2ldW2ldOwogICAgICAgIH0KCiAgICAgICAgZG91YmxlIHJzcyA9IDAuMDsKICAgICAgICBmb3IgKGludCBpID0gMDsgaSA8IE47IGkrKykgewogICAgICAgICAgICByc3MgKz0gcG93KHhfbmV3W2ldIC0geFtpXSwgMik7CiAgICAgICAgfQoKICAgICAgICBwcmludGYoIiVkXHRcdCUuMmVcdFx0JS41Zlx0JS41Zlx0JS41Zlx0JS41ZlxuIiwgaywgcnNzLCB4X25ld1swXSwgeF9uZXdbMV0sIHhfbmV3WzJdLCB4X25ld1szXSk7CgogICAgICAgIGZvciAoaW50IGkgPSAwOyBpIDwgTjsgaSsrKSB7CiAgICAgICAgICAgIHhbaV0gPSB4X25ld1tpXTsKICAgICAgICB9CgogICAgICAgIGlmIChyc3MgPCBFUFMpIHsKICAgICAgICAgICAgcHJpbnRmKCJcbj4+IOWPjuadn+WIpOWumuadoeS7tuOCkua6gOi2s+OBl+OBvuOBl+OBn++8iOioiOeul+WbnuaVsDogJWTlm57vvIlcbiIsIGspOwogICAgICAgICAgICBjb252ZXJnZWQgPSAxOwogICAgICAgICAgICBicmVhazsKICAgICAgICB9CiAgICB9CgogICAgaWYgKCFjb252ZXJnZWQpIHsKICAgICAgICBwcmludGYoIlxu6K2m5ZGKOiDmjIflrprjgZXjgozjgZ/oqIjnrpflm57mlbDku6XlhoXjgavlj47mnZ/jgZfjgb7jgZvjgpPjgafjgZfjgZ/jgIJcbiIpOwogICAgfQoKICAgIHByaW50ZigiXG49PT0g6YCj56uL5pa556iL5byP44Gu5pyA57WC6KejICjjg6TjgrPjg5Pjga7lj43lvqnms5UpID09PVxuIik7CiAgICBmb3IgKGludCBpID0gMDsgaSA8IE47IGkrKykgewogICAgICAgIHByaW50ZigieFslZF0gPSAlMTAuNmZcbiIsIGkgKyAxLCB4W2ldKTsKICAgIH0KCiAgICB2ZXJpZnlfc29sdXRpb24oQSwgYiwgeCk7Cn0KCnZvaWQgZ2F1c3Nfc2VpZGVsX3NvbHZlcihkb3VibGUgQVtOXVtOXSwgZG91YmxlIGJbTl0pIHsKICAgIGRvdWJsZSB4W05dID0gezAuMCwgMC4wLCAwLjAsIDAuMH07CiAgICBkb3VibGUgeF9vbGRbTl07CiAgICBpbnQgY29udmVyZ2VkID0gMDsKCiAgICBwcmludGYoIuioiOeul+WbnuaVsGtcdOaui+W3rjLkuZflkoxcdFx0eDFcdFx0eDJcdFx0eDNcdFx0eDRcbiIpOwogICAgcHJpbnRmKCItLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS1cbiIpOwoKICAgIGZvciAoaW50IGsgPSAxOyBrIDw9IE1BWF9JVEVSOyBrKyspIHsKICAgICAgICBmb3IgKGludCBpID0gMDsgaSA8IE47IGkrKykgewogICAgICAgICAgICB4X29sZFtpXSA9IHhbaV07CiAgICAgICAgfQoKICAgICAgICBmb3IgKGludCBpID0gMDsgaSA8IE47IGkrKykgewogICAgICAgICAgICBkb3VibGUgc3VtID0gMC4wOwogICAgICAgICAgICBmb3IgKGludCBqID0gMDsgaiA8IE47IGorKykgewogICAgICAgICAgICAgICAgaWYgKGogIT0gaSkgewogICAgICAgICAgICAgICAgICAgIHN1bSArPSBBW2ldW2pdICogeFtqXTsgCiAgICAgICAgICAgICAgICB9CiAgICAgICAgICAgIH0KICAgICAgICAgICAgeFtpXSA9IChiW2ldIC0gc3VtKSAvIEFbaV1baV07CiAgICAgICAgfQoKICAgICAgICBkb3VibGUgcnNzID0gMC4wOwogICAgICAgIGZvciAoaW50IGkgPSAwOyBpIDwgTjsgaSsrKSB7CiAgICAgICAgICAgIHJzcyArPSBwb3coeFtpXSAtIHhfb2xkW2ldLCAyKTsKICAgICAgICB9CgogICAgICAgIHByaW50ZigiJWRcdFx0JS4yZVx0XHQlLjVmXHQlLjVmXHQlLjVmXHQlLjVmXG4iLCBrLCByc3MsIHhbMF0sIHhbMV0sIHhbMl0sIHhbM10pOwoKICAgICAgICBpZiAocnNzIDwgRVBTKSB7CiAgICAgICAgICAgIHByaW50ZigiXG4+PiDlj47mnZ/liKTlrprmnaHku7bjgpLmuoDotrPjgZfjgb7jgZfjgZ/vvIjoqIjnrpflm57mlbA6ICVk5Zue77yJXG4iLCBrKTsKICAgICAgICAgICAgY29udmVyZ2VkID0gMTsKICAgICAgICAgICAgYnJlYWs7CiAgICAgICAgfQogICAgfQoKICAgIGlmICghY29udmVyZ2VkKSB7CiAgICAgICAgcHJpbnRmKCJcbuitpuWRijog5oyH5a6a44GV44KM44Gf6KiI566X5Zue5pWw5Lul5YaF44Gr5Y+O5p2f44GX44G+44Gb44KT44Gn44GX44Gf44CCXG4iKTsKICAgIH0KCiAgICBwcmludGYoIlxuPT09IOmAo+eri+aWueeoi+W8j+OBruacgOe1guinoyAo44Ks44Km44K544O744K244Kk44OH44Or5rOVKSA9PT1cbiIpOwogICAgZm9yIChpbnQgaSA9IDA7IGkgPCBOOyBpKyspIHsKICAgICAgICBwcmludGYoInhbJWRdID0gJTEwLjZmXG4iLCBpICsgMSwgeFtpXSk7CiAgICB9CgogICAgdmVyaWZ5X3NvbHV0aW9uKEEsIGIsIHgpOwp9Cgp2b2lkIHZlcmlmeV9zb2x1dGlvbihkb3VibGUgQVtOXVtOXSwgZG91YmxlIGJbTl0sIGRvdWJsZSB4W05dKSB7CiAgICBwcmludGYoIlxuLS0tIOaxguOBvuOBo+OBn+OCguOBruOBjOino+OBp+OBguOCi+eiuuiqje+8iOOCs+ODvOODieWGheOBp+OBruaknOeulzogQXggPT0gYu+8iSAtLS1cbiIpOwogICAgZm9yIChpbnQgaSA9IDA7IGkgPCBOOyBpKyspIHsKICAgICAgICBkb3VibGUgbGhzX2NhbGMgPSAwLjA7CiAgICAgICAgZm9yIChpbnQgaiA9IDA7IGogPCBOOyBqKyspIHsKICAgICAgICAgICAgbGhzX2NhbGMgKz0gQVtpXVtqXSAqIHhbal07CiAgICAgICAgfQogICAgICAgIGRvdWJsZSBlcnJvciA9IGxoc19jYWxjIC0gYltpXTsKICAgICAgICBwcmludGYoIiVk6KGM55uuOiDoqIjnrpflgKQgQXggPSAlMTAuNWYgfCDoqK3lrprlgKQgYiA9ICUxMC41ZiB8IOiqpOW3riA9ICVlXG4iLCAKICAgICAgICAgICAgICAgaSArIDEsIGxoc19jYWxjLCBiW2ldLCBlcnJvcik7CiAgICB9Cn0=