#include <stdio.h>
#include <math.h>
#define N 4
#define EPS 1e-6
#define MAX_ITER 100
 
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} 
};
 
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];
        }
    }
}
 
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;
    }
}
 
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];
	double x[N] = {1.0, 1.0, 1.0, 1.0};
 
    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_decomposition(B, L, U);
    printf("\n               近似固有値 = %.3f)\n", lambda_hat);
 
    int iter = 0;
    double diff;
    double lambda = 0.0;
 
    do {
        forward_substitution(L, x, 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]);
        }
 
        for (int i = 0; i < N; i++) {
            x[i] = x_next[i];
        }
 
        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(" 収束固有値  = %10.6f\n", lambda);
    printf(" 収束固有ベクトル (x)  = [%.4f, %.4f, %.4f, %.4f]\n", x[0], x[1], x[2], x[3]);
    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));
    }
}
 
int main() {
	find_remaining_eigen(-5.0);
    find_remaining_eigen(-2.0);
    find_remaining_eigen(1.5);
 
    return 0;
}

I2luY2x1ZGUgPHN0ZGlvLmg+CiNpbmNsdWRlIDxtYXRoLmg+CiNkZWZpbmUgTiA0CiNkZWZpbmUgRVBTIDFlLTYKI2RlZmluZSBNQVhfSVRFUiAxMDAKIApkb3VibGUgQVtOXVtOXSA9IHsgCiAgICB7Mi4wLCAgMS4wLCAyLjAsICA1LjB9LCAKICAgIHszLjAsIC0zLjAsIDEuMCwgIDYuMH0sIAogICAgezAuMCwgIDEuMCwgMi4wLCAgNy4wfSwgCiAgICB7MS4wLCAgMS4wLCAzLjAsICAxLjB9IAp9OwogCnZvaWQgbHVfZGVjb21wb3NpdGlvbihkb3VibGUgQltOXVtOXSwgZG91YmxlIExbTl1bTl0sIGRvdWJsZSBVW05dW05dKSB7CiAgICBmb3IgKGludCBpID0gMDsgaSA8IE47IGkrKykgewogICAgICAgIGZvciAoaW50IGogPSAwOyBqIDwgTjsgaisrKSB7CiAgICAgICAgICAgIExbaV1bal0gPSAoaSA9PSBqKSA/IDEuMCA6IDAuMDsKICAgICAgICAgICAgVVtpXVtqXSA9IDAuMDsKICAgICAgICB9CiAgICB9CiAKICAgIGZvciAoaW50IGkgPSAwOyBpIDwgTjsgaSsrKSB7CiAgICAgICAgZm9yIChpbnQgaiA9IGk7IGogPCBOOyBqKyspIHsKICAgICAgICAgICAgZG91YmxlIHN1bSA9IDAuMDsKICAgICAgICAgICAgZm9yIChpbnQgayA9IDA7IGsgPCBpOyBrKyspIHN1bSArPSBMW2ldW2tdICogVVtrXVtqXTsKICAgICAgICAgICAgVVtpXVtqXSA9IEJbaV1bal0gLSBzdW07CiAgICAgICAgfQogICAgICAgIGZvciAoaW50IGogPSBpICsgMTsgaiA8IE47IGorKykgewogICAgICAgICAgICBkb3VibGUgc3VtID0gMC4wOwogICAgICAgICAgICBmb3IgKGludCBrID0gMDsgayA8IGk7IGsrKykgc3VtICs9IExbal1ba10gKiBVW2tdW2ldOwogICAgICAgICAgICBMW2pdW2ldID0gKEJbal1baV0gLSBzdW0pIC8gVVtpXVtpXTsKICAgICAgICB9CiAgICB9Cn0KIAp2b2lkIGZvcndhcmRfc3Vic3RpdHV0aW9uKGRvdWJsZSBMW05dW05dLCBkb3VibGUgeFtOXSwgZG91YmxlIHlbTl0pIHsKICAgIGZvciAoaW50IGkgPSAwOyBpIDwgTjsgaSsrKSB7CiAgICAgICAgZG91YmxlIHN1bSA9IDAuMDsKICAgICAgICBmb3IgKGludCBqID0gMDsgaiA8IGk7IGorKykgc3VtICs9IExbaV1bal0gKiB5W2pdOwogICAgICAgIHlbaV0gPSB4W2ldIC0gc3VtOwogICAgfQp9CiAKdm9pZCBiYWNrd2FyZF9zdWJzdGl0dXRpb24oZG91YmxlIFVbTl1bTl0sIGRvdWJsZSB5W05dLCBkb3VibGUgeFtOXSkgewogICAgZm9yIChpbnQgaSA9IE4gLSAxOyBpID49IDA7IGktLSkgewogICAgICAgIGRvdWJsZSBzdW0gPSAwLjA7CiAgICAgICAgZm9yIChpbnQgaiA9IGkgKyAxOyBqIDwgTjsgaisrKSBzdW0gKz0gVVtpXVtqXSAqIHhbal07CiAgICAgICAgeFtpXSA9ICh5W2ldIC0gc3VtKSAvIFVbaV1baV07CiAgICB9Cn0KIAp2b2lkIGZpbmRfcmVtYWluaW5nX2VpZ2VuKGRvdWJsZSBsYW1iZGFfaGF0KSB7CiAgICBkb3VibGUgQltOXVtOXSwgTFtOXVtOXSwgVVtOXVtOXTsKICAgIGRvdWJsZSB5W05dLCB4X25leHRbTl07Cglkb3VibGUgeFtOXSA9IHsxLjAsIDEuMCwgMS4wLCAxLjB9OwogCiAgICBmb3IgKGludCBpID0gMDsgaSA8IE47IGkrKykgewogICAgICAgIGZvciAoaW50IGogPSAwOyBqIDwgTjsgaisrKSB7CiAgICAgICAgICAgIEJbaV1bal0gPSBBW2ldW2pdIC0gKGkgPT0gaiA/IGxhbWJkYV9oYXQgOiAwLjApOwogICAgICAgIH0KICAgIH0KICAgIGx1X2RlY29tcG9zaXRpb24oQiwgTCwgVSk7CiAgICBwcmludGYoIlxuICAgICAgICAgICAgICAg6L+R5Ly85Zu65pyJ5YCkID0gJS4zZilcbiIsIGxhbWJkYV9oYXQpOwogCiAgICBpbnQgaXRlciA9IDA7CiAgICBkb3VibGUgZGlmZjsKICAgIGRvdWJsZSBsYW1iZGEgPSAwLjA7CiAKICAgIGRvIHsKICAgICAgICBmb3J3YXJkX3N1YnN0aXR1dGlvbihMLCB4LCB5KTsKICAgICAgICBiYWNrd2FyZF9zdWJzdGl0dXRpb24oVSwgeSwgeF9uZXh0KTsKICAgICAgICBkb3VibGUgbWF4X3ZhbCA9IDAuMDsKICAgICAgICBmb3IgKGludCBpID0gMDsgaSA8IE47IGkrKykgewogICAgICAgICAgICBpZiAoZmFicyh4X25leHRbaV0pID4gZmFicyhtYXhfdmFsKSkgewogICAgICAgICAgICAgICAgbWF4X3ZhbCA9IHhfbmV4dFtpXTsKICAgICAgICAgICAgfQogICAgICAgIH0KIAogICAgICAgIGZvciAoaW50IGkgPSAwOyBpIDwgTjsgaSsrKSB7CiAgICAgICAgICAgIHhfbmV4dFtpXSAvPSBtYXhfdmFsOwogICAgICAgIH0KIAogICAgICAgIGRpZmYgPSAwLjA7CiAgICAgICAgZm9yIChpbnQgaSA9IDA7IGkgPCBOOyBpKyspIHsKICAgICAgICAgICAgZGlmZiArPSBmYWJzKHhfbmV4dFtpXSAtIHhbaV0pOwogICAgICAgIH0KIAogICAgICAgIGZvciAoaW50IGkgPSAwOyBpIDwgTjsgaSsrKSB7CiAgICAgICAgICAgIHhbaV0gPSB4X25leHRbaV07CiAgICAgICAgfQogCiAgICAgICAgbGFtYmRhID0gbGFtYmRhX2hhdCArICgxLjAgLyBtYXhfdmFsKTsKICAgICAgICBwcmludGYoIuWPjeW+qSAlMmQ6IHggPSBbJS40ZiwgJS40ZiwgJS40ZiwgJS40Zl0sIOWbuuacieWApCA9ICUuNmZcbiIsIAogICAgICAgICAgICAgICBpdGVyICsgMSwgeFswXSwgeFsxXSwgeFsyXSwgeFszXSwgbGFtYmRhKTsKIAogICAgICAgIGl0ZXIrKzsKICAgIH0gd2hpbGUgKGRpZmYgPiBFUFMgJiYgaXRlciA8IE1BWF9JVEVSKTsKCiAgICBwcmludGYoIiDlj47mnZ/lm7rmnInlgKQgID0gJTEwLjZmXG4iLCBsYW1iZGEpOwogICAgcHJpbnRmKCIg5Y+O5p2f5Zu65pyJ44OZ44Kv44OI44OrICh4KSAgPSBbJS40ZiwgJS40ZiwgJS40ZiwgJS40Zl1cbiIsIHhbMF0sIHhbMV0sIHhbMl0sIHhbM10pOwogICAgcHJpbnRmKCJcbiDlm7rmnInlr77jga7mpJzoqLzjg4Hjgqfjg4Pjgq8gKEEqeCAtIGxhbWJkYSp4KTpcbiIpOwogICAgZm9yIChpbnQgaSA9IDA7IGkgPCBOOyBpKyspIHsKICAgICAgICBkb3VibGUgYXhfaSA9IDAuMDsKICAgICAgICBmb3IgKGludCBqID0gMDsgaiA8IE47IGorKykgewogICAgICAgICAgICBheF9pICs9IEFbaV1bal0gKiB4W2pdOwogICAgICAgIH0KICAgICAgICBkb3VibGUgbHhfaSA9IGxhbWJkYSAqIHhbaV07CiAgICAgICAgcHJpbnRmKCIg56ysJWTooYw6IEEqeCA9ICUxMC42ZiwgbGFtYmRhKnggPSAlMTAuNmYsIOaui+W3rue1tuWvvuWApCA9ICVlXG4iLCAKICAgICAgICAgICAgICAgaSArIDEsIGF4X2ksIGx4X2ksIGZhYnMoYXhfaSAtIGx4X2kpKTsKICAgIH0KfQogCmludCBtYWluKCkgewoJZmluZF9yZW1haW5pbmdfZWlnZW4oLTUuMCk7CiAgICBmaW5kX3JlbWFpbmluZ19laWdlbigtMi4wKTsKICAgIGZpbmRfcmVtYWluaW5nX2VpZ2VuKDEuNSk7CiAKICAgIHJldHVybiAwOwp9