#include <stdio.h>
#include <math.h>
 
#define N 4
#define EPS 1e-6
#define MAX_ITER 100
 
// 行列2の定義
double A[N][N] = { 
    {2.0,  1.0, 2.0,  5.0}, 
    {3.0, -3.0, 1.0,  6.0}, 
    {0.0,  1.0, 2.0,  7.0}, 
    {1.0,  1.0, 3.0,  1.0} 
};
 
// LU分解を行う関数: (A - lambda_hat * I) = L * U
void lu_decomposition(double B[N][N], double L[N][N], double U[N][N]) {
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            L[i][j] = (i == j) ? 1.0 : 0.0;
            U[i][j] = 0.0;
        }
    }
 
    for (int i = 0; i < N; i++) {
        for (int j = i; j < N; j++) {
            double sum = 0.0;
            for (int k = 0; k < i; k++) sum += L[i][k] * U[k][j];
            U[i][j] = B[i][j] - sum;
        }
        for (int j = i + 1; j < N; j++) {
            double sum = 0.0;
            for (int k = 0; k < i; k++) sum += L[j][k] * U[k][i];
            L[j][i] = (B[j][i] - sum) / U[i][i];
        }
    }
}
 
// 前進代入: L * y = x
void forward_substitution(double L[N][N], double x[N], double y[N]) {
    for (int i = 0; i < N; i++) {
        double sum = 0.0;
        for (int j = 0; j < i; j++) sum += L[i][j] * y[j];
        y[i] = x[i] - sum;
    }
}
 
// 後退代入: U * x_next = y
void backward_substitution(double U[N][N], double y[N], double x[N]) {
    for (int i = N - 1; i >= 0; i--) {
        double sum = 0.0;
        for (int j = i + 1; j < N; j++) sum += U[i][j] * x[j];
        x[i] = (y[i] - sum) / U[i][i];
    }
}
 
// 逆反復法の実行
void find_remaining_eigen(double lambda_hat) {
    double B[N][N], L[N][N], U[N][N];
    double y[N], x_next[N];
 
    // 初期ベクトル x0 = [1, 1, 1, 1]^T
    double x[N] = {1.0, 1.0, 1.0, 1.0};
 
    // シフトした行列 B = A - lambda_hat * I の構成
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            B[i][j] = A[i][j] - (i == j ? lambda_hat : 0.0);
        }
    }
 
    // LU分解（ループの外で1回のみ実行）
    lu_decomposition(B, L, U);
 
    printf("\n==================================================================\n");
    printf(" 行列2 に対する逆反復法を実行中 (近似固有値 lambda_hat = %.3f)\n", lambda_hat);
    printf("==================================================================\n");
 
    int iter = 0;
    double diff;
    double lambda = 0.0;
 
    do {
        // 1. 前進代入により L * y = x^(k) を解く
        forward_substitution(L, x, y);
 
        // 2. 後退代入により U * x^(k+1) = y を解く
        backward_substitution(U, y, x_next);
 
        // 規格化因子の探索 (絶対値が最大の要素)
        double max_val = 0.0;
        for (int i = 0; i < N; i++) {
            if (fabs(x_next[i]) > fabs(max_val)) {
                max_val = x_next[i];
            }
        }
 
        // ベクトル要素の規格化
        for (int i = 0; i < N; i++) {
            x_next[i] /= max_val;
        }
 
        // 収束判定のためのベクトル要素の絶対変化量の計算
        diff = 0.0;
        for (int i = 0; i < N; i++) {
            diff += fabs(x_next[i] - x[i]);
        }
 
        // 既存の作業ベクトル x の更新
        for (int i = 0; i < N; i++) {
            x[i] = x_next[i];
        }
 
        // 現在のステップにおける固有値の計算: lambda = lambda_hat + 1 / max_val
        lambda = lambda_hat + (1.0 / max_val);
 
        // 収束状況の出力レイアウト
        printf("反復 %2d: x = [%.4f, %.4f, %.4f, %.4f], 固有値 = %.6f\n", 
               iter + 1, x[0], x[1], x[2], x[3], lambda);
 
        iter++;
    } while (diff > EPS && iter < MAX_ITER);
 
    // 最終結果出力ブロック
    printf("\n 計算結果:\n");
    printf(" 収束固有値 (lambda)   = %10.6f\n", lambda);
    printf(" 収束固有ベクトル (x)  = [%.4f, %.4f, %.4f, %.4f]\n", x[0], x[1], x[2], x[3]);
    printf("==================================================================\n");
 
    // 固有値・固有ベクトルであることの確認（残差チェック）
    printf("\n 固有対の検証チェック (A*x - lambda*x):\n");
    for (int i = 0; i < N; i++) {
        double ax_i = 0.0;
        for (int j = 0; j < N; j++) {
            ax_i += A[i][j] * x[j];
        }
        double lx_i = lambda * x[i];
        printf(" 第%d行: A*x = %10.6f, lambda*x = %10.6f, 残差絶対値 = %e\n", 
               i + 1, ax_i, lx_i, fabs(ax_i - lx_i));
    }
    printf("==================================================================\n");
}
 
int main() {
    printf("=== 行列2: 実数シフトによる残りの固有値スペックトルの探索 ===\n");
 
    // 近似固有値の設定（課題9で求めた主固有値 7.817 を避ける初期推測値）
    find_remaining_eigen(-5.0); // 負の大きな固有値領域
    find_remaining_eigen(-2.0); // 負の中間固有値領域
    find_remaining_eigen(1.5);  // 正の内部固有値領域
 
    return 0;
}

I2luY2x1ZGUgPHN0ZGlvLmg+CiNpbmNsdWRlIDxtYXRoLmg+CiAKI2RlZmluZSBOIDQKI2RlZmluZSBFUFMgMWUtNgojZGVmaW5lIE1BWF9JVEVSIDEwMAogCi8vIOihjOWIlzLjga7lrprnvqkKZG91YmxlIEFbTl1bTl0gPSB7IAogICAgezIuMCwgIDEuMCwgMi4wLCAgNS4wfSwgCiAgICB7My4wLCAtMy4wLCAxLjAsICA2LjB9LCAKICAgIHswLjAsICAxLjAsIDIuMCwgIDcuMH0sIAogICAgezEuMCwgIDEuMCwgMy4wLCAgMS4wfSAKfTsKIAovLyBMVeWIhuino+OCkuihjOOBhumWouaVsDogKEEgLSBsYW1iZGFfaGF0ICogSSkgPSBMICogVQp2b2lkIGx1X2RlY29tcG9zaXRpb24oZG91YmxlIEJbTl1bTl0sIGRvdWJsZSBMW05dW05dLCBkb3VibGUgVVtOXVtOXSkgewogICAgZm9yIChpbnQgaSA9IDA7IGkgPCBOOyBpKyspIHsKICAgICAgICBmb3IgKGludCBqID0gMDsgaiA8IE47IGorKykgewogICAgICAgICAgICBMW2ldW2pdID0gKGkgPT0gaikgPyAxLjAgOiAwLjA7CiAgICAgICAgICAgIFVbaV1bal0gPSAwLjA7CiAgICAgICAgfQogICAgfQogCiAgICBmb3IgKGludCBpID0gMDsgaSA8IE47IGkrKykgewogICAgICAgIGZvciAoaW50IGogPSBpOyBqIDwgTjsgaisrKSB7CiAgICAgICAgICAgIGRvdWJsZSBzdW0gPSAwLjA7CiAgICAgICAgICAgIGZvciAoaW50IGsgPSAwOyBrIDwgaTsgaysrKSBzdW0gKz0gTFtpXVtrXSAqIFVba11bal07CiAgICAgICAgICAgIFVbaV1bal0gPSBCW2ldW2pdIC0gc3VtOwogICAgICAgIH0KICAgICAgICBmb3IgKGludCBqID0gaSArIDE7IGogPCBOOyBqKyspIHsKICAgICAgICAgICAgZG91YmxlIHN1bSA9IDAuMDsKICAgICAgICAgICAgZm9yIChpbnQgayA9IDA7IGsgPCBpOyBrKyspIHN1bSArPSBMW2pdW2tdICogVVtrXVtpXTsKICAgICAgICAgICAgTFtqXVtpXSA9IChCW2pdW2ldIC0gc3VtKSAvIFVbaV1baV07CiAgICAgICAgfQogICAgfQp9CiAKLy8g5YmN6YCy5Luj5YWlOiBMICogeSA9IHgKdm9pZCBmb3J3YXJkX3N1YnN0aXR1dGlvbihkb3VibGUgTFtOXVtOXSwgZG91YmxlIHhbTl0sIGRvdWJsZSB5W05dKSB7CiAgICBmb3IgKGludCBpID0gMDsgaSA8IE47IGkrKykgewogICAgICAgIGRvdWJsZSBzdW0gPSAwLjA7CiAgICAgICAgZm9yIChpbnQgaiA9IDA7IGogPCBpOyBqKyspIHN1bSArPSBMW2ldW2pdICogeVtqXTsKICAgICAgICB5W2ldID0geFtpXSAtIHN1bTsKICAgIH0KfQogCi8vIOW+jOmAgOS7o+WFpTogVSAqIHhfbmV4dCA9IHkKdm9pZCBiYWNrd2FyZF9zdWJzdGl0dXRpb24oZG91YmxlIFVbTl1bTl0sIGRvdWJsZSB5W05dLCBkb3VibGUgeFtOXSkgewogICAgZm9yIChpbnQgaSA9IE4gLSAxOyBpID49IDA7IGktLSkgewogICAgICAgIGRvdWJsZSBzdW0gPSAwLjA7CiAgICAgICAgZm9yIChpbnQgaiA9IGkgKyAxOyBqIDwgTjsgaisrKSBzdW0gKz0gVVtpXVtqXSAqIHhbal07CiAgICAgICAgeFtpXSA9ICh5W2ldIC0gc3VtKSAvIFVbaV1baV07CiAgICB9Cn0KIAovLyDpgIblj43lvqnms5Xjga7lrp/ooYwKdm9pZCBmaW5kX3JlbWFpbmluZ19laWdlbihkb3VibGUgbGFtYmRhX2hhdCkgewogICAgZG91YmxlIEJbTl1bTl0sIExbTl1bTl0sIFVbTl1bTl07CiAgICBkb3VibGUgeVtOXSwgeF9uZXh0W05dOwogCiAgICAvLyDliJ3mnJ/jg5njgq/jg4jjg6sgeDAgPSBbMSwgMSwgMSwgMV1eVAogICAgZG91YmxlIHhbTl0gPSB7MS4wLCAxLjAsIDEuMCwgMS4wfTsKIAogICAgLy8g44K344OV44OI44GX44Gf6KGM5YiXIEIgPSBBIC0gbGFtYmRhX2hhdCAqIEkg44Gu5qeL5oiQCiAgICBmb3IgKGludCBpID0gMDsgaSA8IE47IGkrKykgewogICAgICAgIGZvciAoaW50IGogPSAwOyBqIDwgTjsgaisrKSB7CiAgICAgICAgICAgIEJbaV1bal0gPSBBW2ldW2pdIC0gKGkgPT0gaiA/IGxhbWJkYV9oYXQgOiAwLjApOwogICAgICAgIH0KICAgIH0KIAogICAgLy8gTFXliIbop6PvvIjjg6vjg7zjg5fjga7lpJbjgacx5Zue44Gu44G/5a6f6KGM77yJCiAgICBsdV9kZWNvbXBvc2l0aW9uKEIsIEwsIFUpOwogCiAgICBwcmludGYoIlxuPT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09XG4iKTsKICAgIHByaW50ZigiIOihjOWIlzIg44Gr5a++44GZ44KL6YCG5Y+N5b6p5rOV44KS5a6f6KGM5LitICjov5HkvLzlm7rmnInlgKQgbGFtYmRhX2hhdCA9ICUuM2YpXG4iLCBsYW1iZGFfaGF0KTsKICAgIHByaW50ZigiPT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09XG4iKTsKIAogICAgaW50IGl0ZXIgPSAwOwogICAgZG91YmxlIGRpZmY7CiAgICBkb3VibGUgbGFtYmRhID0gMC4wOwogCiAgICBkbyB7CiAgICAgICAgLy8gMS4g5YmN6YCy5Luj5YWl44Gr44KI44KKIEwgKiB5ID0geF4oaykg44KS6Kej44GPCiAgICAgICAgZm9yd2FyZF9zdWJzdGl0dXRpb24oTCwgeCwgeSk7CiAKICAgICAgICAvLyAyLiDlvozpgIDku6PlhaXjgavjgojjgoogVSAqIHheKGsrMSkgPSB5IOOCkuino+OBjwogICAgICAgIGJhY2t3YXJkX3N1YnN0aXR1dGlvbihVLCB5LCB4X25leHQpOwogCiAgICAgICAgLy8g6KaP5qC85YyW5Zug5a2Q44Gu5o6i57SiICjntbblr77lgKTjgYzmnIDlpKfjga7opoHntKApCiAgICAgICAgZG91YmxlIG1heF92YWwgPSAwLjA7CiAgICAgICAgZm9yIChpbnQgaSA9IDA7IGkgPCBOOyBpKyspIHsKICAgICAgICAgICAgaWYgKGZhYnMoeF9uZXh0W2ldKSA+IGZhYnMobWF4X3ZhbCkpIHsKICAgICAgICAgICAgICAgIG1heF92YWwgPSB4X25leHRbaV07CiAgICAgICAgICAgIH0KICAgICAgICB9CiAKICAgICAgICAvLyDjg5njgq/jg4jjg6vopoHntKDjga7opo/moLzljJYKICAgICAgICBmb3IgKGludCBpID0gMDsgaSA8IE47IGkrKykgewogICAgICAgICAgICB4X25leHRbaV0gLz0gbWF4X3ZhbDsKICAgICAgICB9CiAKICAgICAgICAvLyDlj47mnZ/liKTlrprjga7jgZ/jgoHjga7jg5njgq/jg4jjg6vopoHntKDjga7ntbblr77lpInljJbph4/jga7oqIjnrpcKICAgICAgICBkaWZmID0gMC4wOwogICAgICAgIGZvciAoaW50IGkgPSAwOyBpIDwgTjsgaSsrKSB7CiAgICAgICAgICAgIGRpZmYgKz0gZmFicyh4X25leHRbaV0gLSB4W2ldKTsKICAgICAgICB9CiAKICAgICAgICAvLyDml6LlrZjjga7kvZzmpa3jg5njgq/jg4jjg6sgeCDjga7mm7TmlrAKICAgICAgICBmb3IgKGludCBpID0gMDsgaSA8IE47IGkrKykgewogICAgICAgICAgICB4W2ldID0geF9uZXh0W2ldOwogICAgICAgIH0KIAogICAgICAgIC8vIOePvuWcqOOBruOCueODhuODg+ODl+OBq+OBiuOBkeOCi+WbuuacieWApOOBruioiOeulzogbGFtYmRhID0gbGFtYmRhX2hhdCArIDEgLyBtYXhfdmFsCiAgICAgICAgbGFtYmRhID0gbGFtYmRhX2hhdCArICgxLjAgLyBtYXhfdmFsKTsKIAogICAgICAgIC8vIOWPjuadn+eKtuazgeOBruWHuuWKm+ODrOOCpOOCouOCpuODiAogICAgICAgIHByaW50Zigi5Y+N5b6pICUyZDogeCA9IFslLjRmLCAlLjRmLCAlLjRmLCAlLjRmXSwg5Zu65pyJ5YCkID0gJS42ZlxuIiwgCiAgICAgICAgICAgICAgIGl0ZXIgKyAxLCB4WzBdLCB4WzFdLCB4WzJdLCB4WzNdLCBsYW1iZGEpOwogCiAgICAgICAgaXRlcisrOwogICAgfSB3aGlsZSAoZGlmZiA+IEVQUyAmJiBpdGVyIDwgTUFYX0lURVIpOwogCiAgICAvLyDmnIDntYLntZDmnpzlh7rlipvjg5bjg63jg4Pjgq8KICAgIHByaW50ZigiXG4g6KiI566X57WQ5p6cOlxuIik7CiAgICBwcmludGYoIiDlj47mnZ/lm7rmnInlgKQgKGxhbWJkYSkgICA9ICUxMC42ZlxuIiwgbGFtYmRhKTsKICAgIHByaW50ZigiIOWPjuadn+WbuuacieODmeOCr+ODiOODqyAoeCkgID0gWyUuNGYsICUuNGYsICUuNGYsICUuNGZdXG4iLCB4WzBdLCB4WzFdLCB4WzJdLCB4WzNdKTsKICAgIHByaW50ZigiPT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09XG4iKTsKIAogICAgLy8g5Zu65pyJ5YCk44O75Zu65pyJ44OZ44Kv44OI44Or44Gn44GC44KL44GT44Go44Gu56K66KqN77yI5q6L5beu44OB44Kn44OD44Kv77yJCiAgICBwcmludGYoIlxuIOWbuuacieWvvuOBruaknOiovOODgeOCp+ODg+OCryAoQSp4IC0gbGFtYmRhKngpOlxuIik7CiAgICBmb3IgKGludCBpID0gMDsgaSA8IE47IGkrKykgewogICAgICAgIGRvdWJsZSBheF9pID0gMC4wOwogICAgICAgIGZvciAoaW50IGogPSAwOyBqIDwgTjsgaisrKSB7CiAgICAgICAgICAgIGF4X2kgKz0gQVtpXVtqXSAqIHhbal07CiAgICAgICAgfQogICAgICAgIGRvdWJsZSBseF9pID0gbGFtYmRhICogeFtpXTsKICAgICAgICBwcmludGYoIiDnrKwlZOihjDogQSp4ID0gJTEwLjZmLCBsYW1iZGEqeCA9ICUxMC42Ziwg5q6L5beu57W25a++5YCkID0gJWVcbiIsIAogICAgICAgICAgICAgICBpICsgMSwgYXhfaSwgbHhfaSwgZmFicyhheF9pIC0gbHhfaSkpOwogICAgfQogICAgcHJpbnRmKCI9PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT1cbiIpOwp9CiAKaW50IG1haW4oKSB7CiAgICBwcmludGYoIj09PSDooYzliJcyOiDlrp/mlbDjgrfjg5Xjg4jjgavjgojjgovmrovjgorjga7lm7rmnInlgKTjgrnjg5rjg4Pjgq/jg4jjg6vjga7mjqLntKIgPT09XG4iKTsKIAogICAgLy8g6L+R5Ly85Zu65pyJ5YCk44Gu6Kit5a6a77yI6Kqy6aGMOeOBp+axguOCgeOBn+S4u+WbuuacieWApCA3LjgxNyDjgpLpgb/jgZHjgovliJ3mnJ/mjqjmuKzlgKTvvIkKICAgIGZpbmRfcmVtYWluaW5nX2VpZ2VuKC01LjApOyAvLyDosqDjga7lpKfjgY3jgarlm7rmnInlgKTpoJjln58KICAgIGZpbmRfcmVtYWluaW5nX2VpZ2VuKC0yLjApOyAvLyDosqDjga7kuK3plpPlm7rmnInlgKTpoJjln58KICAgIGZpbmRfcmVtYWluaW5nX2VpZ2VuKDEuNSk7ICAvLyDmraPjga7lhoXpg6jlm7rmnInlgKTpoJjln58KIAogICAgcmV0dXJuIDA7Cn0=