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

#define N 4         
#define EPS 1e-6    // Convergence threshold
#define MAX_ITER 100

void print_matrix(double mat[N][N]) {
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            printf("%7.4f ", mat[i][j]);
        }
        printf("\n");
    }
}

int main() {
    // Initial symmetric matrix A [cite: 13]
    double A[N][N] = {
        {5.0, 4.0, 1.0, 1.0},
        {4.0, 5.0, 1.0, 1.0},
        {1.0, 1.0, 4.0, 2.0},
        {1.0, 1.0, 2.0, 4.0}
    };

    // Backup the original matrix for verification (Ax = lambda * x)
    double A_orig[N][N];
    for (int i = 0; i < N; i++)
        for (int j = 0; j < N; j++)
            A_orig[i][j] = A[i][j];

    // Initialize P as an Identity Matrix. Columns will store eigenvectors[cite: 39].
    double P[N][N] = {0};
    for (int i = 0; i < N; i++) P[i][i] = 1.0;

    int iter = 0;
    
    // Tidy, aligned table header for convergence history
    printf("=== CONVERGENCE HISTORY ===\n");
    printf("%-5s   %-16s   %-8s\n", "Iter", "Max Off-Diagonal", "Position");
    printf("---------------------------------------\n");

    while (iter < MAX_ITER) {
        // 1. Find the maximum off-diagonal element A[p][q] [cite: 76]
        int p = 0, q = 1;
        double max_val = fabs(A[0][1]);
        for (int i = 0; i < N; i++) {
            for (int j = i + 1; j < N; j++) {
                if (fabs(A[i][j]) > max_val) {
                    max_val = fabs(A[i][j]);
                    p = i;
                    q = j;
                }
            }
        }

        // Tabled log row output
        printf("%-5d   %-16.6f   A[%d][%d]\n", iter, max_val, p, q);

        // Check convergence status [cite: 76]
        if (max_val < EPS) {
            printf("---------------------------------------\n");
            printf("Status: Successfully converged.\n\n");
            break;
        }

        // 2. Calculate rotation angle using acos(-1.0) dynamically [cite: 54]
        double phi, cos_t, sin_t;
        if (fabs(A[p][p] - A[q][q]) < 1e-12) {
            phi = acos(-1.0) / 4.0; 
        } else {
            phi = 0.5 * atan2(2.0 * A[p][q], A[p][p] - A[q][q]);
        }
        cos_t = cos(phi);
        sin_t = sin(phi);

        // 3. Update Matrix A [cite: 51]
        double Ap_old = A[p][p];
        double Aq_old = A[q][q];
        
        A[p][p] = Ap_old * cos_t * cos_t + Aq_old * sin_t * sin_t + 2.0 * A[p][q] * sin_t * cos_t;
        A[q][q] = Ap_old * sin_t * sin_t + Aq_old * cos_t * cos_t - 2.0 * A[p][q] * sin_t * cos_t;
        A[p][q] = A[q][p] = 0.0; 

        for (int i = 0; i < N; i++) {
            if (i != p && i != q) {
                double a_ip = A[i][p];
                double a_iq = A[i][q];
                A[i][p] = A[p][i] = a_ip * cos_t + a_iq * sin_t;
                A[i][q] = A[q][i] = -a_ip * sin_t + a_iq * cos_t;
            }
        }

        // 4. Update Cumulative Eigenvector Matrix P [cite: 77]
        for (int i = 0; i < N; i++) {
            double p_ip = P[i][p];
            double p_iq = P[i][q];
            P[i][p] = p_ip * cos_t + p_iq * sin_t;
            P[i][q] = -p_ip * sin_t + p_iq * cos_t;
        }

        iter++;
    }

    // === Output Final Results ===
    printf("=== EIGENVALUES & EIGENVECTORS ===\n");
    for (int j = 0; j < N; j++) {
        printf("Eigenvalue %d (固有値): %f\n", j + 1, A[j][j]); // [cite: 24]
        printf("Eigenvector %d (固有ベクトル):\n[ ", j + 1);
        for (int i = 0; i < N; i++) {
            printf("%f ", P[i][j]); // Prints column vectors horizontally [cite: 39]
        }
        printf("]\n\n");
    }

    // === Automated Verification Step (Ax = lambda * x) ===
    printf("=== 固有対の検証チェック (A*x - lambda*x) ===\n");
    
    for (int j = 0; j < N; j++) {
        double lambda = A[j][j];
        printf("Eigenvalue %d (%f):\n", j + 1, lambda);
        
        for (int i = 0; i < N; i++) {
            // Compute elements of Ax
            double Ax_i = 0.0;
            for (int k = 0; k < N; k++) {
                Ax_i += A_orig[i][k] * P[k][j];
            }
            // Compute elements of lambda * x
            double lambda_x_i = lambda * P[i][j];
            
            // Calculate absolute residual error
            double residual = fabs(Ax_i - lambda_x_i);
            
            // Print out the numerical matching comparison with scientific notation residuals
            printf("  第 %d 行: A*x = %9.6f, lambda*x = %9.6f, 残差絶対値 = %e\n", 
                   i + 1, Ax_i, lambda_x_i, residual);
        }
        printf("\n");
    }

    return 0;
}