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

#define SIZE 4
#define EPSILON 1e-7   // 収束判定条件
#define MAX_ITER 1000  // 最大反復回数

// 行列とベクトルの積 y = A * x
void mat_vec_mul(double A[SIZE][SIZE], double *x, double *y, int n) {
    for (int i = 0; i < n; i++) {
        y[i] = 0.0;
        for (int j = 0; j < n; j++) {
            y[i] += A[i][j] * x[j];
        }
    }
}

int main() {
    // 自身で設定した行列A
    double A[SIZE][SIZE] = {
        {1.0,  2.0, 3.0,  5.0},
        {3.0, -5.0, 1.0,  4.0},
        {5.0,  9.0, 2.0, -6.0},
        {1.0,  7.0, 4.0,  1.0}
    };

    double x[SIZE], y[SIZE], Ax[SIZE];
    double lambda_old = 0.0, lambda_new = 0.0;

    // 初期値の設定 [1, 1, 1, 1]^T
    for (int i = 0; i < SIZE; i++) {
        x[i] = 1.0;
    }

    printf("=== べき乗法の計算および収束状況 ===\n");
    printf("初期ベクトル x(0): [%.1f, %.1f, %.1f, %.1f]\n\n", x[0], x[1], x[2], x[3]);
    printf("%-5s | %-14s | %-14s | %-30s\n", "Iter", "固有値の近似値", "変化量(差分)", "固有ベクトル x(k)");
    printf("----------------------------------------------------------------------------\n");

    int actual_iters = 0;
    for (int iter = 1; iter <= MAX_ITER; iter++) {
        actual_iters = iter;
        
        // y = A * x の計算
        mat_vec_mul(A, x, y, SIZE);

        // スライドのアルゴリズムに従い、第1成分を基準（係数）として抽出
        lambda_new = y[0]; 

        // 第1成分が 1.0 になるように全体を正規化
        if (fabs(lambda_new) < 1e-12) {
            break;
        }
        for (int i = 0; i < SIZE; i++) {
            y[i] /= lambda_new;
        }

        // 前のステップからの変化量を算出
        double diff = fabs(lambda_new - lambda_old);

        // 収束過程の表示 (最初の5回、5回ごと、および最終ステップ)
        if (iter <= 5 || iter % 5 == 0) {
            printf("%-5d | %-14.10f | %-14.4e | [%.6f, %.6f, %.6f, %.6f]\n", 
                   iter, lambda_new, diff, y[0], y[1], y[2], y[3]);
        }

        // 収束判定条件のチェック
        if (diff < EPSILON) {
            if (iter % 5 != 0 && iter > 5) {
                printf("%-5d | %-14.10f | %-14.4e | [%.6f, %.6f, %.6f, %.6f]\n", 
                       iter, lambda_new, diff, y[0], y[1], y[2], y[3]);
            }
            break;
        }

        // 次の反復のためにベクトルと値を更新
        for (int i = 0; i < SIZE; i++) {
            x[i] = y[i];
        }
        lambda_old = lambda_new;
    }

    printf("----------------------------------------------------------------------------\n");
    printf("反復回数 %d 回で収束判定条件を満たしました。\n", actual_iters);

    // === 最終計算結果の出力 ===
    printf("\n=== 最終結果 ===\n");
    printf("最大固有値 (λ1)          : %.10f\n", lambda_new);
    printf("対応する固有ベクトル (x1) : [");
    for (int i = 0; i < SIZE; i++) {
        printf("%.10f%s", y[i], (i == SIZE - 1) ? "]\n" : ", ");
    }

    // === Ax = λx の検証出力 ===
    printf("\n=== 固有値・固有ベクトルであることの検証 (Ax = λx) ===\n");
    mat_vec_mul(A, y, Ax, SIZE); 
    printf("Ax    : [");
    for (int i = 0; i < SIZE; i++) {
        printf("%.6f%s", Ax[i], (i == SIZE - 1) ? "]\n" : ", ");
    }
    printf("λ * x : [");
    for (int i = 0; i < SIZE; i++) {
        printf("%.6f%s", lambda_new * y[i], (i == SIZE - 1) ? "]\n" : ", ");
    }

    return 0;
}