#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");
 
    // Verification block integrated safely inside the function scope
    printf("\n Verification Check (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(" Row %d: A*x = %10.6f, lambda*x = %10.6f, Diff = %e\n", 
               i + 1, ax_i, lx_i, fabs(ax_i - lx_i));
    }
    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;
}

I2luY2x1ZGUgPHN0ZGlvLmg+CiNpbmNsdWRlIDxtYXRoLmg+CiAKI2RlZmluZSBOIDQKI2RlZmluZSBFUFMgMWUtNgojZGVmaW5lIE1BWF9JVEVSIDEwMAogCi8vIE1BVFJJWCAyIGRlZmluaXRpb24KZG91YmxlIEFbTl1bTl0gPSB7IAogICAgezIuMCwgIDEuMCwgMi4wLCAgNS4wfSwgCiAgICB7My4wLCAtMy4wLCAxLjAsICA2LjB9LCAKICAgIHswLjAsICAxLjAsIDIuMCwgIDcuMH0sIAogICAgezEuMCwgIDEuMCwgMy4wLCAgMS4wfSAKfTsKIAovLyBGdW5jdGlvbiB0byBwZXJmb3JtIExVIGRlY29tcG9zaXRpb246IChBIC0gbGFtYmRhX2hhdCAqIEkpID0gTCAqIFUKdm9pZCBsdV9kZWNvbXBvc2l0aW9uKGRvdWJsZSBCW05dW05dLCBkb3VibGUgTFtOXVtOXSwgZG91YmxlIFVbTl1bTl0pIHsKICAgIGZvciAoaW50IGkgPSAwOyBpIDwgTjsgaSsrKSB7CiAgICAgICAgZm9yIChpbnQgaiA9IDA7IGogPCBOOyBqKyspIHsKICAgICAgICAgICAgTFtpXVtqXSA9IChpID09IGopID8gMS4wIDogMC4wOwogICAgICAgICAgICBVW2ldW2pdID0gMC4wOwogICAgICAgIH0KICAgIH0KIAogICAgZm9yIChpbnQgaSA9IDA7IGkgPCBOOyBpKyspIHsKICAgICAgICBmb3IgKGludCBqID0gaTsgaiA8IE47IGorKykgewogICAgICAgICAgICBkb3VibGUgc3VtID0gMC4wOwogICAgICAgICAgICBmb3IgKGludCBrID0gMDsgayA8IGk7IGsrKykgc3VtICs9IExbaV1ba10gKiBVW2tdW2pdOwogICAgICAgICAgICBVW2ldW2pdID0gQltpXVtqXSAtIHN1bTsKICAgICAgICB9CiAgICAgICAgZm9yIChpbnQgaiA9IGkgKyAxOyBqIDwgTjsgaisrKSB7CiAgICAgICAgICAgIGRvdWJsZSBzdW0gPSAwLjA7CiAgICAgICAgICAgIGZvciAoaW50IGsgPSAwOyBrIDwgaTsgaysrKSBzdW0gKz0gTFtqXVtrXSAqIFVba11baV07CiAgICAgICAgICAgIExbal1baV0gPSAoQltqXVtpXSAtIHN1bSkgLyBVW2ldW2ldOwogICAgICAgIH0KICAgIH0KfQogCi8vIEZvcndhcmQgU3Vic3RpdHV0aW9uOiBMICogeSA9IHgKdm9pZCBmb3J3YXJkX3N1YnN0aXR1dGlvbihkb3VibGUgTFtOXVtOXSwgZG91YmxlIHhbTl0sIGRvdWJsZSB5W05dKSB7CiAgICBmb3IgKGludCBpID0gMDsgaSA8IE47IGkrKykgewogICAgICAgIGRvdWJsZSBzdW0gPSAwLjA7CiAgICAgICAgZm9yIChpbnQgaiA9IDA7IGogPCBpOyBqKyspIHN1bSArPSBMW2ldW2pdICogeVtqXTsKICAgICAgICB5W2ldID0geFtpXSAtIHN1bTsKICAgIH0KfQogCi8vIEJhY2t3YXJkIFN1YnN0aXR1dGlvbjogVSAqIHhfbmV4dCA9IHkKdm9pZCBiYWNrd2FyZF9zdWJzdGl0dXRpb24oZG91YmxlIFVbTl1bTl0sIGRvdWJsZSB5W05dLCBkb3VibGUgeFtOXSkgewogICAgZm9yIChpbnQgaSA9IE4gLSAxOyBpID49IDA7IGktLSkgewogICAgICAgIGRvdWJsZSBzdW0gPSAwLjA7CiAgICAgICAgZm9yIChpbnQgaiA9IGkgKyAxOyBqIDwgTjsgaisrKSBzdW0gKz0gVVtpXVtqXSAqIHhbal07CiAgICAgICAgeFtpXSA9ICh5W2ldIC0gc3VtKSAvIFVbaV1baV07CiAgICB9Cn0KIAovLyBJbnZlcnNlIEl0ZXJhdGlvbiBFeGVjdXRpb24Kdm9pZCBmaW5kX3JlbWFpbmluZ19laWdlbihkb3VibGUgbGFtYmRhX2hhdCkgewogICAgZG91YmxlIEJbTl1bTl0sIExbTl1bTl0sIFVbTl1bTl07CiAgICBkb3VibGUgeVtOXSwgeF9uZXh0W05dOwogCiAgICAvLyBJbml0aWFsIHN0YXJ0aW5nIHZlY3RvciB4MCA9IFsxLCAxLCAxLCAxXV5UCiAgICBkb3VibGUgeFtOXSA9IHsxLjAsIDEuMCwgMS4wLCAxLjB9OwogCiAgICAvLyBDb25zdHJ1Y3QgbWF0cml4IEIgPSBBIC0gbGFtYmRhX2hhdCAqIEkKICAgIGZvciAoaW50IGkgPSAwOyBpIDwgTjsgaSsrKSB7CiAgICAgICAgZm9yIChpbnQgaiA9IDA7IGogPCBOOyBqKyspIHsKICAgICAgICAgICAgQltpXVtqXSA9IEFbaV1bal0gLSAoaSA9PSBqID8gbGFtYmRhX2hhdCA6IDAuMCk7CiAgICAgICAgfQogICAgfQogCiAgICAvLyBTaW5nbGUgTFUgZGVjb21wb3NpdGlvbiBzdGVwCiAgICBsdV9kZWNvbXBvc2l0aW9uKEIsIEwsIFUpOwogCiAgICBwcmludGYoIlxuPT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09XG4iKTsKICAgIHByaW50ZigiIFJ1bm5pbmcgSW52ZXJzZSBJdGVyYXRpb24gZm9yIE1hdHJpeCAyIChsYW1iZGFfaGF0ID0gJS4zZilcbiIsIGxhbWJkYV9oYXQpOwogICAgcHJpbnRmKCI9PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT1cbiIpOwogCiAgICBpbnQgaXRlciA9IDA7CiAgICBkb3VibGUgZGlmZjsKICAgIGRvdWJsZSBsYW1iZGEgPSAwLjA7CiAKICAgIGRvIHsKICAgICAgICAvLyAxLiBTb2x2ZSBMICogeSA9IHheKGspCiAgICAgICAgZm9yd2FyZF9zdWJzdGl0dXRpb24oTCwgeCwgeSk7CiAKICAgICAgICAvLyAyLiBTb2x2ZSBVICogeF4oaysxKSA9IHkKICAgICAgICBiYWNrd2FyZF9zdWJzdGl0dXRpb24oVSwgeSwgeF9uZXh0KTsKIAogICAgICAgIC8vIEZpbmQgc2NhbGluZyBmYWN0b3IgKG1heCBlbGVtZW50IGJ5IGFic29sdXRlIHZhbHVlKQogICAgICAgIGRvdWJsZSBtYXhfdmFsID0gMC4wOwogICAgICAgIGZvciAoaW50IGkgPSAwOyBpIDwgTjsgaSsrKSB7CiAgICAgICAgICAgIGlmIChmYWJzKHhfbmV4dFtpXSkgPiBmYWJzKG1heF92YWwpKSB7CiAgICAgICAgICAgICAgICBtYXhfdmFsID0geF9uZXh0W2ldOwogICAgICAgICAgICB9CiAgICAgICAgfQogCiAgICAgICAgLy8gTm9ybWFsaXplIHZlY3RvciBjb21wb25lbnRzCiAgICAgICAgZm9yIChpbnQgaSA9IDA7IGkgPCBOOyBpKyspIHsKICAgICAgICAgICAgeF9uZXh0W2ldIC89IG1heF92YWw7CiAgICAgICAgfQogCiAgICAgICAgLy8gQ2FsY3VsYXRlIGFic29sdXRlIGNoYW5nZSBpbiB2ZWN0b3IgZWxlbWVudHMgdG8gdHJhY2UgY29udmVyZ2VuY2UKICAgICAgICBkaWZmID0gMC4wOwogICAgICAgIGZvciAoaW50IGkgPSAwOyBpIDwgTjsgaSsrKSB7CiAgICAgICAgICAgIGRpZmYgKz0gZmFicyh4X25leHRbaV0gLSB4W2ldKTsKICAgICAgICB9CiAKICAgICAgICAvLyBVcGRhdGUgd29ya2luZyB2ZWN0b3IgeAogICAgICAgIGZvciAoaW50IGkgPSAwOyBpIDwgTjsgaSsrKSB7CiAgICAgICAgICAgIHhbaV0gPSB4X25leHRbaV07CiAgICAgICAgfQogCiAgICAgICAgLy8gQ29tcHV0ZSBjdXJyZW50IHRydWUgdGFyZ2V0IGVpZ2VudmFsdWUgc3RlcDogbGFtYmRhID0gbGFtYmRhX2hhdCArIDEgLyBtYXhfdmFsCiAgICAgICAgbGFtYmRhID0gbGFtYmRhX2hhdCArICgxLjAgLyBtYXhfdmFsKTsKIAogICAgICAgIC8vIFRyYWNrIGNvbnRpbnVvdXMgc3RlcCBiZWhhdmlvcgogICAgICAgIHByaW50ZigiSXRlciAlMmQ6IHggPSBbJS40ZiwgJS40ZiwgJS40ZiwgJS40Zl0sIGxhbWJkYSA9ICUuNmZcbiIsIAogICAgICAgICAgICAgICBpdGVyICsgMSwgeFswXSwgeFsxXSwgeFsyXSwgeFszXSwgbGFtYmRhKTsKIAogICAgICAgIGl0ZXIrKzsKICAgIH0gd2hpbGUgKGRpZmYgPiBFUFMgJiYgaXRlciA8IE1BWF9JVEVSKTsKIAogICAgLy8gRmluYWwgYW5zd2VyIG91dHB1dCBibG9jawogICAgcHJpbnRmKCJcbiBEaXNjb3ZlcmVkIEFuc3dlcjpcbiIpOwogICAgcHJpbnRmKCIgRWlnZW52YWx1ZSAobGFtYmRhKSAgID0gJTEwLjZmXG4iLCBsYW1iZGEpOwogICAgcHJpbnRmKCIgRWlnZW52ZWN0b3IgKHgpICAgICAgICA9IFslLjRmLCAlLjRmLCAlLjRmLCAlLjRmXVxuIiwgeFswXSwgeFsxXSwgeFsyXSwgeFszXSk7CiAgICBwcmludGYoIj09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PVxuIik7CiAKICAgIC8vIFZlcmlmaWNhdGlvbiBibG9jayBpbnRlZ3JhdGVkIHNhZmVseSBpbnNpZGUgdGhlIGZ1bmN0aW9uIHNjb3BlCiAgICBwcmludGYoIlxuIFZlcmlmaWNhdGlvbiBDaGVjayAoQSp4IC0gbGFtYmRhKngpOlxuIik7CiAgICBmb3IgKGludCBpID0gMDsgaSA8IE47IGkrKykgewogICAgICAgIGRvdWJsZSBheF9pID0gMC4wOwogICAgICAgIGZvciAoaW50IGogPSAwOyBqIDwgTjsgaisrKSB7CiAgICAgICAgICAgIGF4X2kgKz0gQVtpXVtqXSAqIHhbal07CiAgICAgICAgfQogICAgICAgIGRvdWJsZSBseF9pID0gbGFtYmRhICogeFtpXTsKICAgICAgICBwcmludGYoIiBSb3cgJWQ6IEEqeCA9ICUxMC42ZiwgbGFtYmRhKnggPSAlMTAuNmYsIERpZmYgPSAlZVxuIiwgCiAgICAgICAgICAgICAgIGkgKyAxLCBheF9pLCBseF9pLCBmYWJzKGF4X2kgLSBseF9pKSk7CiAgICB9CiAgICBwcmludGYoIj09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PVxuIik7Cn0KIAppbnQgbWFpbigpIHsKICAgIHByaW50ZigiPT09IE1hdHJpeCAyOiBGaW5kaW5nIE1pc3NpbmcgU3BlY3RydW0gRWxlbWVudHMgdmlhIFJlYWwtU2hpZnRzID09PVxuIik7CiAKICAgIC8vIENoYW5naW5nIHRoZSBndWVzcyB2YWx1ZXMgdG8gY2FwdHVyZSB0aGUgdGhyZWUgcmVtYWluaW5nIHNwZWN0cnVtIGJvdW5kcwogICAgLy8gKEF2b2lkaW5nIHRoZSBkb21pbmFudCBlaWdlbnZhbHVlIDcuODE3IGZvdW5kIGluIGFzc2lnbm1lbnQgOSkKICAgIGZpbmRfcmVtYWluaW5nX2VpZ2VuKC01LjApOyAvLyBUYXJnZXQgbGFyZ2UgbmVnYXRpdmUgYm91bmQKICAgIGZpbmRfcmVtYWluaW5nX2VpZ2VuKC0yLjApOyAvLyBUYXJnZXQgaW50ZXJtZWRpYXRlIG5lZ2F0aXZlIGJvdW5kCiAgICBmaW5kX3JlbWFpbmluZ19laWdlbigxLjUpOyAgLy8gVGFyZ2V0IGludGVyaW9yIHBvc2l0aXZlIGJvdW5kCiAKICAgIHJldHVybiAwOwp9