#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;
}

I2luY2x1ZGUgPHN0ZGlvLmg+CiNpbmNsdWRlIDxtYXRoLmg+CgojZGVmaW5lIE4gNAojZGVmaW5lIEVQUyAxZS02CiNkZWZpbmUgTUFYX0lURVIgMTAwCgovLyBNQVRSSVggMiBkZWZpbml0aW9uCmRvdWJsZSBBW05dW05dID0geyAKICAgIHsyLjAsICAxLjAsIDIuMCwgIDUuMH0sIAogICAgezMuMCwgLTMuMCwgMS4wLCAgNi4wfSwgCiAgICB7MC4wLCAgMS4wLCAyLjAsICA3LjB9LCAKICAgIHsxLjAsICAxLjAsIDMuMCwgIDEuMH0gCn07CgovLyBGdW5jdGlvbiB0byBwZXJmb3JtIExVIGRlY29tcG9zaXRpb246IChBIC0gbGFtYmRhX2hhdCAqIEkpID0gTCAqIFUKdm9pZCBsdV9kZWNvbXBvc2l0aW9uKGRvdWJsZSBCW05dW05dLCBkb3VibGUgTFtOXVtOXSwgZG91YmxlIFVbTl1bTl0pIHsKICAgIGZvciAoaW50IGkgPSAwOyBpIDwgTjsgaSsrKSB7CiAgICAgICAgZm9yIChpbnQgaiA9IDA7IGogPCBOOyBqKyspIHsKICAgICAgICAgICAgTFtpXVtqXSA9IChpID09IGopID8gMS4wIDogMC4wOwogICAgICAgICAgICBVW2ldW2pdID0gMC4wOwogICAgICAgIH0KICAgIH0KCiAgICBmb3IgKGludCBpID0gMDsgaSA8IE47IGkrKykgewogICAgICAgIGZvciAoaW50IGogPSBpOyBqIDwgTjsgaisrKSB7CiAgICAgICAgICAgIGRvdWJsZSBzdW0gPSAwLjA7CiAgICAgICAgICAgIGZvciAoaW50IGsgPSAwOyBrIDwgaTsgaysrKSBzdW0gKz0gTFtpXVtrXSAqIFVba11bal07CiAgICAgICAgICAgIFVbaV1bal0gPSBCW2ldW2pdIC0gc3VtOwogICAgICAgIH0KICAgICAgICBmb3IgKGludCBqID0gaSArIDE7IGogPCBOOyBqKyspIHsKICAgICAgICAgICAgZG91YmxlIHN1bSA9IDAuMDsKICAgICAgICAgICAgZm9yIChpbnQgayA9IDA7IGsgPCBpOyBrKyspIHN1bSArPSBMW2pdW2tdICogVVtrXVtpXTsKICAgICAgICAgICAgTFtqXVtpXSA9IChCW2pdW2ldIC0gc3VtKSAvIFVbaV1baV07CiAgICAgICAgfQogICAgfQp9CgovLyBGb3J3YXJkIFN1YnN0aXR1dGlvbjogTCAqIHkgPSB4CnZvaWQgZm9yd2FyZF9zdWJzdGl0dXRpb24oZG91YmxlIExbTl1bTl0sIGRvdWJsZSB4W05dLCBkb3VibGUgeVtOXSkgewogICAgZm9yIChpbnQgaSA9IDA7IGkgPCBOOyBpKyspIHsKICAgICAgICBkb3VibGUgc3VtID0gMC4wOwogICAgICAgIGZvciAoaW50IGogPSAwOyBqIDwgaTsgaisrKSBzdW0gKz0gTFtpXVtqXSAqIHlbal07CiAgICAgICAgeVtpXSA9IHhbaV0gLSBzdW07CiAgICB9Cn0KCi8vIEJhY2t3YXJkIFN1YnN0aXR1dGlvbjogVSAqIHhfbmV4dCA9IHkKdm9pZCBiYWNrd2FyZF9zdWJzdGl0dXRpb24oZG91YmxlIFVbTl1bTl0sIGRvdWJsZSB5W05dLCBkb3VibGUgeFtOXSkgewogICAgZm9yIChpbnQgaSA9IE4gLSAxOyBpID49IDA7IGktLSkgewogICAgICAgIGRvdWJsZSBzdW0gPSAwLjA7CiAgICAgICAgZm9yIChpbnQgaiA9IGkgKyAxOyBqIDwgTjsgaisrKSBzdW0gKz0gVVtpXVtqXSAqIHhbal07CiAgICAgICAgeFtpXSA9ICh5W2ldIC0gc3VtKSAvIFVbaV1baV07CiAgICB9Cn0KCi8vIEludmVyc2UgSXRlcmF0aW9uIEV4ZWN1dGlvbgp2b2lkIGZpbmRfcmVtYWluaW5nX2VpZ2VuKGRvdWJsZSBsYW1iZGFfaGF0KSB7CiAgICBkb3VibGUgQltOXVtOXSwgTFtOXVtOXSwgVVtOXVtOXTsKICAgIGRvdWJsZSB5W05dLCB4X25leHRbTl07CiAgICAKICAgIC8vIEluaXRpYWwgc3RhcnRpbmcgdmVjdG9yIHgwID0gWzEsIDEsIDEsIDFdXlQKICAgIGRvdWJsZSB4W05dID0gezEuMCwgMS4wLCAxLjAsIDEuMH07CgogICAgLy8gQ29uc3RydWN0IG1hdHJpeCBCID0gQSAtIGxhbWJkYV9oYXQgKiBJCiAgICBmb3IgKGludCBpID0gMDsgaSA8IE47IGkrKykgewogICAgICAgIGZvciAoaW50IGogPSAwOyBqIDwgTjsgaisrKSB7CiAgICAgICAgICAgIEJbaV1bal0gPSBBW2ldW2pdIC0gKGkgPT0gaiA/IGxhbWJkYV9oYXQgOiAwLjApOwogICAgICAgIH0KICAgIH0KCiAgICAvLyBTaW5nbGUgTFUgZGVjb21wb3NpdGlvbiBzdGVwCiAgICBsdV9kZWNvbXBvc2l0aW9uKEIsIEwsIFUpOwoKICAgIHByaW50ZigiXG49PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT1cbiIpOwogICAgcHJpbnRmKCIgUnVubmluZyBJbnZlcnNlIEl0ZXJhdGlvbiBmb3IgTWF0cml4IDIgKGxhbWJkYV9oYXQgPSAlLjNmKVxuIiwgbGFtYmRhX2hhdCk7CiAgICBwcmludGYoIj09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PVxuIik7CiAgICAKICAgIGludCBpdGVyID0gMDsKICAgIGRvdWJsZSBkaWZmOwogICAgZG91YmxlIGxhbWJkYSA9IDAuMDsKCiAgICBkbyB7CiAgICAgICAgLy8gMS4gU29sdmUgTCAqIHkgPSB4XihrKQogICAgICAgIGZvcndhcmRfc3Vic3RpdHV0aW9uKEwsIHgsIHkpOwogICAgICAgIAogICAgICAgIC8vIDIuIFNvbHZlIFUgKiB4XihrKzEpID0geQogICAgICAgIGJhY2t3YXJkX3N1YnN0aXR1dGlvbihVLCB5LCB4X25leHQpOwoKICAgICAgICAvLyBGaW5kIHNjYWxpbmcgZmFjdG9yIChtYXggZWxlbWVudCBieSBhYnNvbHV0ZSB2YWx1ZSkKICAgICAgICBkb3VibGUgbWF4X3ZhbCA9IDAuMDsKICAgICAgICBmb3IgKGludCBpID0gMDsgaSA8IE47IGkrKykgewogICAgICAgICAgICBpZiAoZmFicyh4X25leHRbaV0pID4gZmFicyhtYXhfdmFsKSkgewogICAgICAgICAgICAgICAgbWF4X3ZhbCA9IHhfbmV4dFtpXTsKICAgICAgICAgICAgfQogICAgICAgIH0KCiAgICAgICAgLy8gTm9ybWFsaXplIHZlY3RvciBjb21wb25lbnRzCiAgICAgICAgZm9yIChpbnQgaSA9IDA7IGkgPCBOOyBpKyspIHsKICAgICAgICAgICAgeF9uZXh0W2ldIC89IG1heF92YWw7CiAgICAgICAgfQoKICAgICAgICAvLyBDYWxjdWxhdGUgYWJzb2x1dGUgY2hhbmdlIGluIHZlY3RvciBlbGVtZW50cyB0byB0cmFjZSBjb252ZXJnZW5jZQogICAgICAgIGRpZmYgPSAwLjA7CiAgICAgICAgZm9yIChpbnQgaSA9IDA7IGkgPCBOOyBpKyspIHsKICAgICAgICAgICAgZGlmZiArPSBmYWJzKHhfbmV4dFtpXSAtIHhbaV0pOwogICAgICAgIH0KCiAgICAgICAgLy8gVXBkYXRlIHdvcmtpbmcgdmVjdG9yIHgKICAgICAgICBmb3IgKGludCBpID0gMDsgaSA8IE47IGkrKykgewogICAgICAgICAgICB4W2ldID0geF9uZXh0W2ldOwogICAgICAgIH0KCiAgICAgICAgLy8gQ29tcHV0ZSBjdXJyZW50IHRydWUgdGFyZ2V0IGVpZ2VudmFsdWUgc3RlcDogbGFtYmRhID0gbGFtYmRhX2hhdCArIDEgLyBtYXhfdmFsCiAgICAgICAgbGFtYmRhID0gbGFtYmRhX2hhdCArICgxLjAgLyBtYXhfdmFsKTsKICAgICAgICAKICAgICAgICAvLyBUcmFjayBjb250aW51b3VzIHN0ZXAgYmVoYXZpb3IKICAgICAgICBwcmludGYoIkl0ZXIgJTJkOiB4ID0gWyUuNGYsICUuNGYsICUuNGYsICUuNGZdLCBsYW1iZGEgPSAlLjZmXG4iLCAKICAgICAgICAgICAgICAgaXRlciArIDEsIHhbMF0sIHhbMV0sIHhbMl0sIHhbM10sIGxhbWJkYSk7CgogICAgICAgIGl0ZXIrKzsKICAgIH0gd2hpbGUgKGRpZmYgPiBFUFMgJiYgaXRlciA8IE1BWF9JVEVSKTsKCiAgICAvLyBGaW5hbCBhbnN3ZXIgb3V0cHV0IGJsb2NrCiAgICBwcmludGYoIlxuIERpc2NvdmVyZWQgQW5zd2VyOlxuIik7CiAgICBwcmludGYoIiBFaWdlbnZhbHVlIChsYW1iZGEpICAgPSAlMTAuNmZcbiIsIGxhbWJkYSk7CiAgICBwcmludGYoIiBFaWdlbnZlY3RvciAoeCkgICAgICAgID0gWyUuNGYsICUuNGYsICUuNGYsICUuNGZdXG4iLCB4WzBdLCB4WzFdLCB4WzJdLCB4WzNdKTsKICAgIHByaW50ZigiPT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09XG4iKTsKfQoKaW50IG1haW4oKSB7CiAgICBwcmludGYoIj09PSBNYXRyaXggMjogRmluZGluZyBNaXNzaW5nIFNwZWN0cnVtIEVsZW1lbnRzIHZpYSBSZWFsLVNoaWZ0cyA9PT1cbiIpOwogICAgCiAgICAvLyBDaGFuZ2luZyB0aGUgZ3Vlc3MgdmFsdWVzIHRvIGNhcHR1cmUgdGhlIHRocmVlIHJlbWFpbmluZyBzcGVjdHJ1bSBib3VuZHMKICAgIC8vIChBdm9pZGluZyB0aGUgZG9taW5hbnQgZWlnZW52YWx1ZSA3LjgxNyBmb3VuZCBpbiBhc3NpZ25tZW50IDkpCiAgICBmaW5kX3JlbWFpbmluZ19laWdlbigtNS4wKTsgLy8gVGFyZ2V0IGxhcmdlIG5lZ2F0aXZlIGJvdW5kCiAgICBmaW5kX3JlbWFpbmluZ19laWdlbigtMi4wKTsgLy8gVGFyZ2V0IGludGVybWVkaWF0ZSBuZWdhdGl2ZSBib3VuZAogICAgZmluZF9yZW1haW5pbmdfZWlnZW4oMS41KTsgIC8vIFRhcmdldCBpbnRlcmlvciBwb3NpdGl2ZSBib3VuZAoKICAgIHJldHVybiAwOwp9