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

#define SIZE 4
#define EPSILON_POWER 1e-9  
#define EPSILON_NEWTON 1e-7  
#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];
        }
    }
}

double f(double lambda) {
    return pow(lambda, 4) + pow(lambda, 3) - 54.0 * pow(lambda, 2) - 224.0 * lambda - 345.0;
}
 
double df(double lambda) {
    return 4.0 * pow(lambda, 3) + 3.0 * pow(lambda, 2) - 108.0 * lambda - 224.0;
}
 

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

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

        if (fabs(df_val) < 1e-12) return; 

        double delta = f_val / df_val;
        lambda -= delta;
        iter++;

        if (fabs(delta) < EPSILON_NEWTON) {
            printf("試行 #%d (初期値 %5.1f) -> 収束成功 (%2d回) | 固有値解 λ = %.10f\n", 
                   attempt, lambda0, iter, lambda);
            return;
        }
    }
}


int main() {
    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], Ax[SIZE];
    double lambda_old = 0.0, lambda_new = 0.0;
	 printf("                べき乗法              \n");
    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_POWER) {
            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;
        }

        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);
    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" : ", ");
    }

    printf("\n\n");
    printf("   固有値が最大固有値であることの確認:ニュートン・ラフソン法     \n");
    printf("対象の特性方程式: f(λ) = λ^4 - 2.0λ^3 - 43.0λ^2 - 36.0λ + 130.0 = 0\n\n");

    double initial_lambda[] = {10.0, 2.0, -2.5, -6.0};
    for (int i = 0; i < 4; i++) {
        solve_newton(i + 1, initial_lambda[i]);
    }
    return 0;}