#include <stdio.h>
#include <stdlib.h>
#include <math.h>
 
#define SIZE 4
#define EPSILON 1e-9   
#define MAX_ITER 1000  
 
double vector_norm(double *v, int n) {
    double sum = 0.0;
    for (int i = 0; i < n; i++) {
        sum += v[i] * v[i];
    }
    return sqrt(sum);
}
 
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() {
    double A[SIZE][SIZE] = {
        {2.0,  3.0, 2.0,  5.0},
        {3.0, -3.0, 1.0,  6.0},
        {0.0,  1.0, 2.0, 7.0},
        {1.0,  7.0, 3.0,  1.0}
    };
 
    double x[SIZE], y[SIZE], Ax[SIZE];
    double lambda_old = 0.0, lambda_new = 0.0;
 
    for (int i = 0; i < SIZE; i++) {
        x[i] = 1.0;
    }
 
    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)");
    int actual_iters = 0;
    for (int iter = 1; iter <= MAX_ITER; iter++) {
        actual_iters = iter;
        mat_vec_mul(A, x, y, SIZE);
 
        double norm_y = vector_norm(y, SIZE);
        if (norm_y < 1e-12) break;
        for (int i = 0; i < SIZE; i++) {
            y[i] /= norm_y;
        }
 
        mat_vec_mul(A, y, Ax, SIZE);
        double num = 0.0, den = 0.0;
        for (int i = 0; i < SIZE; i++) {
            num += y[i] * Ax[i];
            den += y[i] * y[i];
        }
        lambda_new = num / den;
        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("反復回数 %d 回で収束しました。\n", actual_iters);
    printf("\n最大固有値 (λ1)          : %.10f\n", lambda_new);
    printf("対応する固有ベクトル (x1) : [");
    for (int i = 0; i < SIZE; i++) {
        printf("%.10f%s", y[i], (i == SIZE - 1) ? "]\n" : ", ");
    }
 
    printf("\n   検証: Ax と λx の比較   \n");
    mat_vec_mul(A, y, Ax, SIZE); // Axを再計算
    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;
}
 
 

I2luY2x1ZGUgPHN0ZGlvLmg+CiNpbmNsdWRlIDxzdGRsaWIuaD4KI2luY2x1ZGUgPG1hdGguaD4KCiNkZWZpbmUgU0laRSA0CiNkZWZpbmUgRVBTSUxPTiAxZS05ICAgCiNkZWZpbmUgTUFYX0lURVIgMTAwMCAgCgpkb3VibGUgdmVjdG9yX25vcm0oZG91YmxlICp2LCBpbnQgbikgewogICAgZG91YmxlIHN1bSA9IDAuMDsKICAgIGZvciAoaW50IGkgPSAwOyBpIDwgbjsgaSsrKSB7CiAgICAgICAgc3VtICs9IHZbaV0gKiB2W2ldOwogICAgfQogICAgcmV0dXJuIHNxcnQoc3VtKTsKfQoKdm9pZCBtYXRfdmVjX211bChkb3VibGUgQVtTSVpFXVtTSVpFXSwgZG91YmxlICp4LCBkb3VibGUgKnksIGludCBuKSB7CiAgICBmb3IgKGludCBpID0gMDsgaSA8IG47IGkrKykgewogICAgICAgIHlbaV0gPSAwLjA7CiAgICAgICAgZm9yIChpbnQgaiA9IDA7IGogPCBuOyBqKyspIHsKICAgICAgICAgICAgeVtpXSArPSBBW2ldW2pdICogeFtqXTsKICAgICAgICB9CiAgICB9Cn0KCmludCBtYWluKCkgewogICAgZG91YmxlIEFbU0laRV1bU0laRV0gPSB7CiAgICAgICAgezIuMCwgIDMuMCwgMi4wLCAgNS4wfSwKICAgICAgICB7My4wLCAtMy4wLCAxLjAsICA2LjB9LAogICAgICAgIHswLjAsICAxLjAsIDIuMCwgNy4wfSwKICAgICAgICB7MS4wLCAgNy4wLCAzLjAsICAxLjB9CiAgICB9OwoKICAgIGRvdWJsZSB4W1NJWkVdLCB5W1NJWkVdLCBBeFtTSVpFXTsKICAgIGRvdWJsZSBsYW1iZGFfb2xkID0gMC4wLCBsYW1iZGFfbmV3ID0gMC4wOwoKICAgIGZvciAoaW50IGkgPSAwOyBpIDwgU0laRTsgaSsrKSB7CiAgICAgICAgeFtpXSA9IDEuMDsKICAgIH0KCiAgICBwcmludGYoIuWIneacn+ODmeOCr+ODiOODqyB4KDApOiBbJS4xZiwgJS4xZiwgJS4xZiwgJS4xZl1cblxuIiwgeFswXSwgeFsxXSwgeFsyXSwgeFszXSk7CiAgICBwcmludGYoIiUtNXMgfCAlLTEycyB8ICUtMTJzIHwgJS0zMHNcbiIsICJJdGVyIiwgIuWbuuacieWApCjmlrApIiwgIuW3ruWIhuWApCIsICLlm7rmnInjg5njgq/jg4jjg6t4KGspIik7CiAgICBpbnQgYWN0dWFsX2l0ZXJzID0gMDsKICAgIGZvciAoaW50IGl0ZXIgPSAxOyBpdGVyIDw9IE1BWF9JVEVSOyBpdGVyKyspIHsKICAgICAgICBhY3R1YWxfaXRlcnMgPSBpdGVyOwogICAgICAgIG1hdF92ZWNfbXVsKEEsIHgsIHksIFNJWkUpOwoKICAgICAgICBkb3VibGUgbm9ybV95ID0gdmVjdG9yX25vcm0oeSwgU0laRSk7CiAgICAgICAgaWYgKG5vcm1feSA8IDFlLTEyKSBicmVhazsKICAgICAgICBmb3IgKGludCBpID0gMDsgaSA8IFNJWkU7IGkrKykgewogICAgICAgICAgICB5W2ldIC89IG5vcm1feTsKICAgICAgICB9CgogICAgICAgIG1hdF92ZWNfbXVsKEEsIHksIEF4LCBTSVpFKTsKICAgICAgICBkb3VibGUgbnVtID0gMC4wLCBkZW4gPSAwLjA7CiAgICAgICAgZm9yIChpbnQgaSA9IDA7IGkgPCBTSVpFOyBpKyspIHsKICAgICAgICAgICAgbnVtICs9IHlbaV0gKiBBeFtpXTsKICAgICAgICAgICAgZGVuICs9IHlbaV0gKiB5W2ldOwogICAgICAgIH0KICAgICAgICBsYW1iZGFfbmV3ID0gbnVtIC8gZGVuOwogICAgICAgIGRvdWJsZSBkaWZmID0gZmFicyhsYW1iZGFfbmV3IC0gbGFtYmRhX29sZCk7CgogICAKICAgICAgICBpZiAoaXRlciA8PSA1IHx8IGl0ZXIgJSA1ID09IDApIHsKICAgICAgICAgICAgcHJpbnRmKCIlLTVkIHwgJS0xMi44ZiB8ICUtMTIuNGUgfCBbJS41ZiwgJS41ZiwgJS41ZiwgJS41Zl1cbiIsIAogICAgICAgICAgICAgICAgICAgaXRlciwgbGFtYmRhX25ldywgZGlmZiwgeVswXSwgeVsxXSwgeVsyXSwgeVszXSk7CiAgICAgICAgfQoKICAgICAgICBpZiAoZGlmZiA8IEVQU0lMT04pIHsKICAgICAgICAgICAgaWYgKGl0ZXIgJSA1ICE9IDAgJiYgaXRlciA+IDUpIHsgLy8g5pyA57WC44K544OG44OD44OX44KS56K65a6f44Gr6KGo56S644GV44Gb44KLCiAgICAgICAgICAgICAgICBwcmludGYoIiUtNWQgfCAlLTEyLjhmIHwgJS0xMi40ZSB8IFslLjVmLCAlLjVmLCAlLjVmLCAlLjVmXVxuIiwgCiAgICAgICAgICAgICAgICAgICAgICAgaXRlciwgbGFtYmRhX25ldywgZGlmZiwgeVswXSwgeVsxXSwgeVsyXSwgeVszXSk7CiAgICAgICAgICAgIH0KICAgICAgICAgICAgYnJlYWs7CiAgICAgICAgfQoKICAgICAgICAvLyDmrKHjga7jgrnjg4bjg4Pjg5fjgavlkJHjgZHjgaYgeCDjgpLmm7TmlrAKICAgICAgICBmb3IgKGludCBpID0gMDsgaSA8IFNJWkU7IGkrKykgewogICAgICAgICAgICB4W2ldID0geVtpXTsKICAgICAgICB9CiAgICAgICAgbGFtYmRhX29sZCA9IGxhbWJkYV9uZXc7CiAgICB9CgogICAgCiAgICBwcmludGYoIuWPjeW+qeWbnuaVsCAlZCDlm57jgaflj47mnZ/jgZfjgb7jgZfjgZ/jgIJcbiIsIGFjdHVhbF9pdGVycyk7CiAgICBwcmludGYoIlxu5pyA5aSn5Zu65pyJ5YCkICjOuzEpICAgICAgICAgIDogJS4xMGZcbiIsIGxhbWJkYV9uZXcpOwogICAgcHJpbnRmKCLlr77lv5zjgZnjgovlm7rmnInjg5njgq/jg4jjg6sgKHgxKSA6IFsiKTsKICAgIGZvciAoaW50IGkgPSAwOyBpIDwgU0laRTsgaSsrKSB7CiAgICAgICAgcHJpbnRmKCIlLjEwZiVzIiwgeVtpXSwgKGkgPT0gU0laRSAtIDEpID8gIl1cbiIgOiAiLCAiKTsKICAgIH0KCiAgICBwcmludGYoIlxuICAg5qSc6Ki8OiBBeCDjgaggzrt4IOOBruavlOi8gyAgIFxuIik7CiAgICBtYXRfdmVjX211bChBLCB5LCBBeCwgU0laRSk7IC8vIEF444KS5YaN6KiI566XCiAgICBwcmludGYoIkF4ICAgIDogWyIpOwogICAgZm9yIChpbnQgaSA9IDA7IGkgPCBTSVpFOyBpKyspIHsKICAgICAgICBwcmludGYoIiUuNmYlcyIsIEF4W2ldLCAoaSA9PSBTSVpFIC0gMSkgPyAiXVxuIiA6ICIsICIpOwogICAgfQogICAgcHJpbnRmKCLOuyAqIHggOiBbIik7CiAgICBmb3IgKGludCBpID0gMDsgaSA8IFNJWkU7IGkrKykgewogICAgICAgIHByaW50ZigiJS42ZiVzIiwgbGFtYmRhX25ldyAqIHlbaV0sIChpID09IFNJWkUgLSAxKSA/ICJdXG4iIDogIiwgIik7CiAgICB9CgogICAgcmV0dXJuIDA7Cn0KCg==