#include <stdio.h>
 
#define SIZE 4
 
// 行列のトレース（対角成分の和）を計算
double calculate_trace(double mat[SIZE][SIZE]) {
    double trace = 0.0;
    for (int i = 0; i < SIZE; i++) {
        trace += mat[i][i];
    }
    return trace;
}
 
// 行列の積 C = A * B
void matrix_multiply(double A[SIZE][SIZE], double B[SIZE][SIZE], double C[SIZE][SIZE]) {
    for (int i = 0; i < SIZE; i++) {
        for (int j = 0; j < SIZE; j++) {
            C[i][j] = 0.0;
            for (int k = 0; k < SIZE; k++) {
                C[i][j] += A[i][k] * B[k][j];
            }
        }
    }
}
 
int main() {
    // あなたの1つ目の行列 A
    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,  1.0, 3.0,  1.0}
    };
 
    double M[SIZE][SIZE]; // ワーキング行列 AM
    double B[SIZE][SIZE]; // 各ステップの行列 B
    double c[SIZE + 1];   // 多項式の係数 (c0, c1, c2, c3, c4)
 
    // 初期化
    c[0] = 1.0; // λ^4 の係数は必ず 1.0
    for (int i = 0; i < SIZE; i++) {
        for (int j = 0; j < SIZE; j++) {
            B[i][j] = (i == j) ? 1.0 : 0.0; // 単位行列
        }
    }
 
    // Faddeev-LeVerrier アルゴリズムの実行
    for (int m = 1; m <= SIZE; m++) {
        matrix_multiply(A, B, M); // M = A * B_{m-1}
        double trace = calculate_trace(M);
        c[m] = -trace / m;        // 係数の決定
 
        // 次のステップのための B_m を計算: B_m = M + c_m * I
        for (int i = 0; i < SIZE; i++) {
            for (int j = 0; j < SIZE; j++) {
                B[i][j] = M[i][j] + ((i == j) ? c[m] : 0.0);
            }
        }
    }
 
    // === 行列式（特性方程式）の出力 ===
    printf("=== 行列式の展開結果 (特性方程式) ===\n");
    printf("det(A - λI) = 0\n\n");
    printf("f(λ) = λ^4 + (%.1f)λ^3 + (%.1f)λ^2 + (%.1f)λ + (%.1f) = 0\n", 
            c[1], c[2], c[3], c[4]);
 
    return 0;
}

I2luY2x1ZGUgPHN0ZGlvLmg+CgojZGVmaW5lIFNJWkUgNAoKLy8g6KGM5YiX44Gu44OI44Os44O844K577yI5a++6KeS5oiQ5YiG44Gu5ZKM77yJ44KS6KiI566XCmRvdWJsZSBjYWxjdWxhdGVfdHJhY2UoZG91YmxlIG1hdFtTSVpFXVtTSVpFXSkgewogICAgZG91YmxlIHRyYWNlID0gMC4wOwogICAgZm9yIChpbnQgaSA9IDA7IGkgPCBTSVpFOyBpKyspIHsKICAgICAgICB0cmFjZSArPSBtYXRbaV1baV07CiAgICB9CiAgICByZXR1cm4gdHJhY2U7Cn0KCi8vIOihjOWIl+OBruepjSBDID0gQSAqIEIKdm9pZCBtYXRyaXhfbXVsdGlwbHkoZG91YmxlIEFbU0laRV1bU0laRV0sIGRvdWJsZSBCW1NJWkVdW1NJWkVdLCBkb3VibGUgQ1tTSVpFXVtTSVpFXSkgewogICAgZm9yIChpbnQgaSA9IDA7IGkgPCBTSVpFOyBpKyspIHsKICAgICAgICBmb3IgKGludCBqID0gMDsgaiA8IFNJWkU7IGorKykgewogICAgICAgICAgICBDW2ldW2pdID0gMC4wOwogICAgICAgICAgICBmb3IgKGludCBrID0gMDsgayA8IFNJWkU7IGsrKykgewogICAgICAgICAgICAgICAgQ1tpXVtqXSArPSBBW2ldW2tdICogQltrXVtqXTsKICAgICAgICAgICAgfQogICAgICAgIH0KICAgIH0KfQoKaW50IG1haW4oKSB7CiAgICAvLyDjgYLjgarjgZ/jga4x44Gk55uu44Gu6KGM5YiXIEEKICAgIGRvdWJsZSBBW1NJWkVdW1NJWkVdID0gewogICAgICAgIHsyLjAsICAzLjAsIDIuMCwgIDUuMH0sCiAgICAgICAgezMuMCwgLTMuMCwgMS4wLCAgNi4wfSwKICAgICAgICB7MC4wLCAgMS4wLCAyLjAsIDcuMH0sCiAgICAgICAgezEuMCwgIDEuMCwgMy4wLCAgMS4wfQogICAgfTsKCiAgICBkb3VibGUgTVtTSVpFXVtTSVpFXTsgLy8g44Ov44O844Kt44Oz44Kw6KGM5YiXIEFNCiAgICBkb3VibGUgQltTSVpFXVtTSVpFXTsgLy8g5ZCE44K544OG44OD44OX44Gu6KGM5YiXIEIKICAgIGRvdWJsZSBjW1NJWkUgKyAxXTsgICAvLyDlpJrpoIXlvI/jga7kv4LmlbAgKGMwLCBjMSwgYzIsIGMzLCBjNCkKCiAgICAvLyDliJ3mnJ/ljJYKICAgIGNbMF0gPSAxLjA7IC8vIM67XjQg44Gu5L+C5pWw44Gv5b+F44GaIDEuMAogICAgZm9yIChpbnQgaSA9IDA7IGkgPCBTSVpFOyBpKyspIHsKICAgICAgICBmb3IgKGludCBqID0gMDsgaiA8IFNJWkU7IGorKykgewogICAgICAgICAgICBCW2ldW2pdID0gKGkgPT0gaikgPyAxLjAgOiAwLjA7IC8vIOWNmOS9jeihjOWIlwogICAgICAgIH0KICAgIH0KCiAgICAvLyBGYWRkZWV2LUxlVmVycmllciDjgqLjg6vjgrTjg6rjgrrjg6Djga7lrp/ooYwKICAgIGZvciAoaW50IG0gPSAxOyBtIDw9IFNJWkU7IG0rKykgewogICAgICAgIG1hdHJpeF9tdWx0aXBseShBLCBCLCBNKTsgLy8gTSA9IEEgKiBCX3ttLTF9CiAgICAgICAgZG91YmxlIHRyYWNlID0gY2FsY3VsYXRlX3RyYWNlKE0pOwogICAgICAgIGNbbV0gPSAtdHJhY2UgLyBtOyAgICAgICAgLy8g5L+C5pWw44Gu5rG65a6aCgogICAgICAgIC8vIOasoeOBruOCueODhuODg+ODl+OBruOBn+OCgeOBriBCX20g44KS6KiI566XOiBCX20gPSBNICsgY19tICogSQogICAgICAgIGZvciAoaW50IGkgPSAwOyBpIDwgU0laRTsgaSsrKSB7CiAgICAgICAgICAgIGZvciAoaW50IGogPSAwOyBqIDwgU0laRTsgaisrKSB7CiAgICAgICAgICAgICAgICBCW2ldW2pdID0gTVtpXVtqXSArICgoaSA9PSBqKSA/IGNbbV0gOiAwLjApOwogICAgICAgICAgICB9CiAgICAgICAgfQogICAgfQoKICAgIC8vID09PSDooYzliJflvI/vvIjnibnmgKfmlrnnqIvlvI/vvInjga7lh7rlipsgPT09CiAgICBwcmludGYoIj09PSDooYzliJflvI/jga7lsZXplovntZDmnpwgKOeJueaAp+aWueeoi+W8jykgPT09XG4iKTsKICAgIHByaW50ZigiZGV0KEEgLSDOu0kpID0gMFxuXG4iKTsKICAgIHByaW50ZigiZijOuykgPSDOu140ICsgKCUuMWYpzrteMyArICglLjFmKc67XjIgKyAoJS4xZinOuyArICglLjFmKSA9IDBcbiIsIAogICAgICAgICAgICBjWzFdLCBjWzJdLCBjWzNdLCBjWzRdKTsKCiAgICByZXR1cm4gMDsKfQ==