#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("近似固有値 lambda_hat = %.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(" 収束固有値 (lambda)   = %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\n", 
               i + 1, ax_i, lx_i, fabs(ax_i - lx_i));
    }
}
 
int main() {
    printf("        残りの固有値スペックトルの探索      \n");
 
    find_remaining_eigen(-5.0);
    find_remaining_eigen(-2.0);
    find_remaining_eigen(1.5);
 
    return 0;
}