#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,  1.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;
}