#include <stdio.h>
 #include <math.h>
 #include <stdlib.h>
 #define TOL 1e-6
 // 許容誤差
#define MAX_ITER 1000 // 最大反復回数
// 行列とベクトルの積を計算する関数
void mat_vec_mult(int n, double A[n][n], double x[n], double result[n]) {
 for (int i = 0; i < n; i++) {
 result[i] = 0.0;
 for (int j = 0; j < n; j++) {
 result[i] += A[i][j] * x[j];
 }
 }
 }
 // ベクトルのノルムを計算する関数
double norm(int n, double x[n]) {
 double sum = 0.0;
 for (int i = 0; i < n; i++) {
 sum += x[i] * x[i];
 }
 return sqrt(sum);
 }
 // ベクトルを正規化する関数
void normalize(int n, double x[n]) {
 double x_norm = norm(n, x);
 for (int i = 0; i < n; i++) {
 x[i] /= x_norm;
 }
 }
 // 累乗法による最大固有値と対応する固有ベクトルを計算する関数
void power_method(int n, double A[n][n], double x0[n], double *eigenvalue, double eigenvector[n]) {
 double x[n], y[n];
 double lambda_old = 0.0, lambda_new = 0.0;
 // 初期ベクトルをコピー
for (int i = 0; i < n; i++) {
 x[i] = x0[i];
 }
 // 正規化
normalize(n, x);
 for (int iter = 0; iter < MAX_ITER; iter++) {
 // 行列とベクトルの積
mat_vec_mult(n, A, x, y);
 // 新しい固有値の近似値を計算（レイリー商）
lambda_new = 0.0;
 for (int i = 0; i < n; i++) {
 lambda_new += x[i] * y[i];
 }
 // 収束判定
if (fabs(lambda_new- lambda_old) < TOL) {
 *eigenvalue = lambda_new;
 for (int i = 0; i < n; i++) {
 eigenvector[i] = y[i];
 }
 return;
 }
 // 正規化
normalize(n, y);
// 次の反復のために更新
    lambda_old=lambda_new;
 for(int i =0; i<n; i++) {
 x[i] =y[i];
 }
 }
 // 収束しない場合
  fprintf(stderr, "累乗法が収束しませんでした。¥n");
 exit(EXIT_FAILURE);
 }
 // メイン関数
int main() {
 // 行列サイズ
  int n =4;
 // 対象行列
  double A[4][4] = {
 {16, -1, 1, 2},
 {2, 12, 1, -1},
 {1, 3, 24, 2},
 {4, -2, 1, 20},
 };
// 初期ベクトル
double x0[4] = {1, 1, 1, 1};
 // 結果を格納する変数
double eigenvalue;
 double eigenvector[2];
 // 累乗法を実行
power_method(n, A, x0, &eigenvalue, eigenvector);
 // 正規化
normalize(n, eigenvector);
 // 結果を出力
printf("最大固有値: %lf¥n", eigenvalue);
 printf("対応する固有ベクトル: [");
 for (int i = 0; i < n; i++) {
 printf("%lf", eigenvector[i]);
 if (i < n- 1) printf(", ");
 }
 printf("]¥n");
 return 0;
 }