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

#define EPSILON 1e-7   // べき乗法と合わせた収束判定条件
#define MAX_ITER 100

// 1つ目の行列の特性方程式 f(λ)
double f(double lambda) {
    return pow(lambda, 4) + pow(lambda, 3) - 88.0 * pow(lambda, 2) - 132.0 * lambda + 603.0;
}

// 導関数 f'(λ)
double df(double lambda) {
    return 4.0 * pow(lambda, 3) + 3.0 * pow(lambda, 2) - 176.0 * lambda - 132.0;
}

void solve_newton(int attempt, double lambda0) {
    double lambda = lambda0;
    int iter = 0;

    printf("\n--- 試行 #%d: (初期値 lambda0 = %.2f) ---\n", attempt, lambda0);
    printf("%-5s %-15s\n", "反復", "lambda");
    printf("-------------------------\n");

    while (iter <= MAX_ITER) {
        double f_val = f(lambda);
        double df_val = df(lambda);

        // lambda の値のみを綺麗に出力
        printf("%-5d %-15.10f\n", iter, lambda);

        // 収束判定
        if (fabs(f_val) < EPSILON) {
            printf(">> 収束しました（反復回数: %d回）\n", iter);
            printf(">> 得られた固有値 (lambda): %.10f\n", lambda);
            return;
        }

        // 導関数がゼロに近すぎる場合の安全対策
        if (fabs(df_val) < 1e-12) {
            printf(">> 失敗: 導関数が0に近づきすぎました。\n");
            return;
        }

        // ニュートン・ラフソン法の更新式: λ = λ - f(λ)/f'(λ)
        lambda -= f_val / df_val;
        iter++;
    }
    printf(">> 失敗: %d回以内に収束しませんでした。\n", MAX_ITER);
}

int main() {
    // すべての実数解（固有値）を見つけられるように複数の初期値を設定
    // 8.0 からスタートすると、べき乗法で求めた最大固有値に収束します
    double initial_lambda[] = {8.0, 5.0, -1.0, -9.0};

    for (int i = 0; i < 4; i++) {
        solve_newton(i + 1, initial_lambda[i]);
    }

    return 0;
}