#include <stdio.h>
#include <stdlib.h>
#include <math.h>
 
#define SIZE 4
#define EPSILON 1e-9   // 収束判定条件
#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];
    double lambda_old = 0.0, lambda_new = 0.0;
 
    // 1. 初期値の設定 [1, 1, 1, 1]^T
    for (int i = 0; i < SIZE; i++) {
        x[i] = 1.0;
    }
 
    printf("=== スライドに準拠したべき乗法（第1成分を1に正規化） ===\n");
    printf("初期ベクトル x(0): [%.1f, %.1f, %.1f, %.1f]\n\n", x[0], x[1], x[2], x[3]);
    printf("%-5s | %-12s | %-12s | %-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になるように全体を割る (くくり出す)
        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);
 
        // 収束状況の出力表示 (スライドのように推移がわかる工夫)
        if (iter <= 5 || iter % 5 == 0) {
            printf("%-5d | %-12.8f | %-12.4e | [%.5f, %.5f, %.5f, %.5f]\n", 
                   iter, lambda_new, diff, y[0], y[1], y[2], y[3]);
        }
 
        if (diff < EPSILON) {
            if (iter % 5 != 0 && iter > 5) {
                printf("%-5d | %-12.8f | %-12.4e | [%.5f, %.5f, %.5f, %.5f]\n", 
                       iter, lambda_new, diff, y[0], y[1], y[2], y[3]);
            }
            break;
        }
 
        // 次のステップに向けて x を更新
        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" : ", ");
    }
 
    return 0;
}

I2luY2x1ZGUgPHN0ZGlvLmg+CiNpbmNsdWRlIDxzdGRsaWIuaD4KI2luY2x1ZGUgPG1hdGguaD4KCiNkZWZpbmUgU0laRSA0CiNkZWZpbmUgRVBTSUxPTiAxZS05ICAgLy8g5Y+O5p2f5Yik5a6a5p2h5Lu2CiNkZWZpbmUgTUFYX0lURVIgMTAwMCAgLy8g5pyA5aSn5Y+N5b6p5Zue5pWwCgovLyDooYzliJfjgajjg5njgq/jg4jjg6vjga7nqY0geSA9IEEgKiB4CnZvaWQgbWF0X3ZlY19tdWwoZG91YmxlIEFbU0laRV1bU0laRV0sIGRvdWJsZSAqeCwgZG91YmxlICp5LCBpbnQgbikgewogICAgZm9yIChpbnQgaSA9IDA7IGkgPCBuOyBpKyspIHsKICAgICAgICB5W2ldID0gMC4wOwogICAgICAgIGZvciAoaW50IGogPSAwOyBqIDwgbjsgaisrKSB7CiAgICAgICAgICAgIHlbaV0gKz0gQVtpXVtqXSAqIHhbal07CiAgICAgICAgfQogICAgfQp9CgppbnQgbWFpbigpIHsKICAgIC8vIOiHqui6q+OBp+ioreWumuOBl+OBn+ihjOWIl0EKICAgIGRvdWJsZSBBW1NJWkVdW1NJWkVdID0gewogICAgICAgIHsxLjAsICAyLjAsIDMuMCwgIDUuMH0sCiAgICAgICAgezMuMCwgLTUuMCwgMS4wLCAgNC4wfSwKICAgICAgICB7NS4wLCAgOS4wLCAyLjAsIC02LjB9LAogICAgICAgIHsxLjAsICA3LjAsIDQuMCwgIDEuMH0KICAgIH07CgogICAgZG91YmxlIHhbU0laRV0sIHlbU0laRV07CiAgICBkb3VibGUgbGFtYmRhX29sZCA9IDAuMCwgbGFtYmRhX25ldyA9IDAuMDsKCiAgICAvLyAxLiDliJ3mnJ/lgKTjga7oqK3lrpogWzEsIDEsIDEsIDFdXlQKICAgIGZvciAoaW50IGkgPSAwOyBpIDwgU0laRTsgaSsrKSB7CiAgICAgICAgeFtpXSA9IDEuMDsKICAgIH0KCiAgICBwcmludGYoIj09PSDjgrnjg6njgqTjg4njgavmupbmi6DjgZfjgZ/jgbnjgY3kuZfms5XvvIjnrKwx5oiQ5YiG44KSMeOBq+ato+imj+WMlu+8iSA9PT1cbiIpOwogICAgcHJpbnRmKCLliJ3mnJ/jg5njgq/jg4jjg6sgeCgwKTogWyUuMWYsICUuMWYsICUuMWYsICUuMWZdXG5cbiIsIHhbMF0sIHhbMV0sIHhbMl0sIHhbM10pOwogICAgcHJpbnRmKCIlLTVzIHwgJS0xMnMgfCAlLTEycyB8ICUtMzBzXG4iLCAiSXRlciIsICLlm7rmnInlgKTjga7ov5HkvLwiLCAi5beu5YiG5YCkIiwgIuWbuuacieODmeOCr+ODiOODq3goaykiKTsKICAgIHByaW50ZigiLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLVxuIik7CgogICAgaW50IGFjdHVhbF9pdGVycyA9IDA7CiAgICBmb3IgKGludCBpdGVyID0gMTsgaXRlciA8PSBNQVhfSVRFUjsgaXRlcisrKSB7CiAgICAgICAgYWN0dWFsX2l0ZXJzID0gaXRlcjsKICAgICAgICAKICAgICAgICAvLyB5ID0gQSAqIHgKICAgICAgICBtYXRfdmVjX211bChBLCB4LCB5LCBTSVpFKTsKCiAgICAgICAgLy8g44CQ44K544Op44Kk44OJ5rqW5oug44CRMeihjOebruOBruaIkOWIhuOCkuWbuuacieWApOOBruaOqOWumuWApO+8iOS/guaVsO+8ieOBqOOBl+OBpuaKveWHugogICAgICAgIGxhbWJkYV9uZXcgPSB5WzBdOyAKCiAgICAgICAgLy8g44CQ44K544Op44Kk44OJ5rqW5oug44CR56ysMeaIkOWIhuOBjDHjgavjgarjgovjgojjgYbjgavlhajkvZPjgpLlibLjgosgKOOBj+OBj+OCiuWHuuOBmSkKICAgICAgICBpZiAoZmFicyhsYW1iZGFfbmV3KSA8IDFlLTEyKSBicmVhazsKICAgICAgICBmb3IgKGludCBpID0gMDsgaSA8IFNJWkU7IGkrKykgewogICAgICAgICAgICB5W2ldIC89IGxhbWJkYV9uZXc7CiAgICAgICAgfQoKICAgICAgICAvLyDlj47mnZ/liKTlrpogKOWbuuacieWApOOBruWkieWMlumHj+OCkui/vei3oSkKICAgICAgICBkb3VibGUgZGlmZiA9IGZhYnMobGFtYmRhX25ldyAtIGxhbWJkYV9vbGQpOwoKICAgICAgICAvLyDlj47mnZ/nirbms4Hjga7lh7rlipvooajnpLogKOOCueODqeOCpOODieOBruOCiOOBhuOBq+aOqOenu+OBjOOCj+OBi+OCi+W3peWkqykKICAgICAgICBpZiAoaXRlciA8PSA1IHx8IGl0ZXIgJSA1ID09IDApIHsKICAgICAgICAgICAgcHJpbnRmKCIlLTVkIHwgJS0xMi44ZiB8ICUtMTIuNGUgfCBbJS41ZiwgJS41ZiwgJS41ZiwgJS41Zl1cbiIsIAogICAgICAgICAgICAgICAgICAgaXRlciwgbGFtYmRhX25ldywgZGlmZiwgeVswXSwgeVsxXSwgeVsyXSwgeVszXSk7CiAgICAgICAgfQoKICAgICAgICBpZiAoZGlmZiA8IEVQU0lMT04pIHsKICAgICAgICAgICAgaWYgKGl0ZXIgJSA1ICE9IDAgJiYgaXRlciA+IDUpIHsKICAgICAgICAgICAgICAgIHByaW50ZigiJS01ZCB8ICUtMTIuOGYgfCAlLTEyLjRlIHwgWyUuNWYsICUuNWYsICUuNWYsICUuNWZdXG4iLCAKICAgICAgICAgICAgICAgICAgICAgICBpdGVyLCBsYW1iZGFfbmV3LCBkaWZmLCB5WzBdLCB5WzFdLCB5WzJdLCB5WzNdKTsKICAgICAgICAgICAgfQogICAgICAgICAgICBicmVhazsKICAgICAgICB9CgogICAgICAgIC8vIOasoeOBruOCueODhuODg+ODl+OBq+WQkeOBkeOBpiB4IOOCkuabtOaWsAogICAgICAgIGZvciAoaW50IGkgPSAwOyBpIDwgU0laRTsgaSsrKSB7CiAgICAgICAgICAgIHhbaV0gPSB5W2ldOwogICAgICAgIH0KICAgICAgICBsYW1iZGFfb2xkID0gbGFtYmRhX25ldzsKICAgIH0KCiAgICBwcmludGYoIi0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS1cbiIpOwogICAgcHJpbnRmKCLlj43lvqnlm57mlbAgJWQg5Zue44Gn5Y+O5p2f44GX44G+44GX44Gf44CCXG4iLCBhY3R1YWxfaXRlcnMpOwoKICAgIC8vID09PSDoqIjnrpfntZDmnpzjga7ooajnpLogPT09CiAgICBwcmludGYoIlxuPT09IOioiOeul+e1kOaenCA9PT1cbiIpOwogICAgcHJpbnRmKCLmnIDlpKflm7rmnInlgKQgKM67MSkgICAgICAgICAgOiAlLjEwZlxuIiwgbGFtYmRhX25ldyk7CiAgICBwcmludGYoIuWvvuW/nOOBmeOCi+WbuuacieODmeOCr+ODiOODqyAoeDEpIDogWyIpOwogICAgZm9yIChpbnQgaSA9IDA7IGkgPCBTSVpFOyBpKyspIHsKICAgICAgICBwcmludGYoIiUuMTBmJXMiLCB5W2ldLCAoaSA9PSBTSVpFIC0gMSkgPyAiXVxuIiA6ICIsICIpOwogICAgfQoKICAgIHJldHVybiAwOwp9