#include <stdio.h>
#include <math.h>

#define N 4
#define EPS 1e-6
#define MAX_ITER 100

// MATRIX 2 definition
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} 
};

// Function to perform LU decomposition: (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];
        }
    }
}

// Forward Substitution: 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;
    }
}

// Backward Substitution: 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];
    }
}

// Inverse Iteration Execution
void find_remaining_eigen(double lambda_hat) {
    double B[N][N], L[N][N], U[N][N];
    double y[N], x_next[N];
    
    // Initial starting vector x0 = [1, 1, 1, 1]^T
    double x[N] = {1.0, 1.0, 1.0, 1.0};

    // Construct matrix 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);
        }
    }

    // Single LU decomposition step
    lu_decomposition(B, L, U);

    printf("\n==================================================================\n");
    printf(" Running Inverse Iteration for Matrix 2 (lambda_hat = %.3f)\n", lambda_hat);
    printf("==================================================================\n");
    
    int iter = 0;
    double diff;
    double lambda = 0.0;

    do {
        // 1. Solve L * y = x^(k)
        forward_substitution(L, x, y);
        
        // 2. Solve U * x^(k+1) = y
        backward_substitution(U, y, x_next);

        // Find scaling factor (max element by absolute value)
        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];
            }
        }

        // Normalize vector components
        for (int i = 0; i < N; i++) {
            x_next[i] /= max_val;
        }

        // Calculate absolute change in vector elements to trace convergence
        diff = 0.0;
        for (int i = 0; i < N; i++) {
            diff += fabs(x_next[i] - x[i]);
        }

        // Update working vector x
        for (int i = 0; i < N; i++) {
            x[i] = x_next[i];
        }

        // Compute current true target eigenvalue step: lambda = lambda_hat + 1 / max_val
        lambda = lambda_hat + (1.0 / max_val);
        
        // Track continuous step behavior
        printf("Iter %2d: x = [%.4f, %.4f, %.4f, %.4f], lambda = %.6f\n", 
               iter + 1, x[0], x[1], x[2], x[3], lambda);

        iter++;
    } while (diff > EPS && iter < MAX_ITER);

    // Final answer output block
    printf("\n Discovered Answer:\n");
    printf(" Eigenvalue (lambda)   = %10.6f\n", lambda);
    printf(" Eigenvector (x)        = [%.4f, %.4f, %.4f, %.4f]\n", x[0], x[1], x[2], x[3]);
    printf("==================================================================\n");
}

int main() {
    printf("=== Matrix 2: Finding Missing Spectrum Elements via Real-Shifts ===\n");
    
    // Changing the guess values to capture the three remaining spectrum bounds
    // (Avoiding the dominant eigenvalue 7.817 found in assignment 9)
    find_remaining_eigen(-5.0); // Target large negative bound
    find_remaining_eigen(-2.0); // Target intermediate negative bound
    find_remaining_eigen(1.5);  // Target interior positive bound

    return 0;
}