#include <stdio.h>
#include <math.h>
#include <complex.h>
 
#define N 4
#define EPS 1e-8
#define MAX_ITER 100
 
double complex A[N][N] = {
    {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}
};
 
/* LU分解 */
void lu_decomposition(double complex B[N][N],
                      double complex L[N][N],
                      double complex U[N][N])
{
    for(int i=0;i<N;i++){
        for(int j=0;j<N;j++){
            L[i][j]=(i==j)?1.0+0.0*I:0.0+0.0*I;
            U[i][j]=0.0+0.0*I;
        }
    }
 
    for(int i=0;i<N;i++){
 
        for(int j=i;j<N;j++){
            double complex sum = 0.0 + 0.0*I;
 
            for(int k=0;k<i;k++)
                sum += L[i][k]*U[k][j];
 
            U[i][j] = B[i][j] - sum;
        }
 
        for(int j=i+1;j<N;j++){
            double complex sum = 0.0 + 0.0*I;
 
            for(int k=0;k<i;k++)
                sum += L[j][k]*U[k][i];
 
            L[j][i] = (B[j][i]-sum)/U[i][i];
        }
    }
}
 
/* 前進代入 */
void forward_substitution(double complex L[N][N],
                          double complex x[N],
                          double complex y[N])
{
    for(int i=0;i<N;i++){
 
        double complex sum = 0.0 + 0.0*I;
 
        for(int j=0;j<i;j++)
            sum += L[i][j]*y[j];
 
        y[i] = x[i] - sum;
    }
}
 
/* 後退代入 */
void backward_substitution(double complex U[N][N],
                           double complex y[N],
                           double complex x[N])
{
    for(int i=N-1;i>=0;i--){
 
        double complex sum = 0.0 + 0.0*I;
 
        for(int j=i+1;j<N;j++)
            sum += U[i][j]*x[j];
 
        x[i] = (y[i]-sum)/U[i][i];
    }
}
 
void inverse_iteration(double complex lambda_hat)
{
    double complex B[N][N];
    double complex L[N][N];
    double complex U[N][N];
 
    double complex x[N] =
    {
        1.0+0.0*I,
        1.0+0.0*I,
        1.0+0.0*I,
        1.0+0.0*I
    };
 
    double complex y[N];
    double complex x_next[N];
 
    for(int i=0;i<N;i++){
        for(int j=0;j<N;j++){
 
            B[i][j] = A[i][j];
 
            if(i==j)
                B[i][j] -= lambda_hat;
        }
    }
 
    lu_decomposition(B,L,U);
 
    printf("\n==================================================================\n");
    printf(" 行列1 に対する逆反復法を実行中 (近似固有値 lambda_hat = %.3f %+.3fi)\n",
           creal(lambda_hat), cimag(lambda_hat));
    printf("==================================================================\n");
 
    int iter=0;
    double diff;
    double complex lambda = 0.0 + 0.0*I;
 
    do{
 
        forward_substitution(L,x,y);
        backward_substitution(U,y,x_next);
 
        /* 絶対値が最大の要素を探索 */
        double max_mag = 0.0;
        int idx = 0;
 
        for(int i=0;i<N;i++){
            if(cabs(x_next[i]) > max_mag){
                max_mag = cabs(x_next[i]);
                idx = i;
            }
        }
 
        double complex scale = x_next[idx];
 
        for(int i=0;i<N;i++)
            x_next[i] /= scale;
 
        diff = 0.0;
 
        for(int i=0;i<N;i++)
            diff += cabs(x_next[i]-x[i]);
 
        for(int i=0;i<N;i++)
            x[i]=x_next[i];
 
        lambda = lambda_hat + 1.0/scale;
 
        // 収束状況の出力レイアウト
        printf("反復 %2d: x = [(%.4f%+.4fi), (%.4f%+.4fi), (%.4f%+.4fi), (%.4f%+.4fi)], 固有値 = %.6f%+.6fi\n", 
               iter + 1, 
               creal(x[0]), cimag(x[0]),
               creal(x[1]), cimag(x[1]),
               creal(x[2]), cimag(x[2]),
               creal(x[3]), cimag(x[3]),
               creal(lambda), cimag(lambda));
 
        iter++;
 
    }while(diff>EPS && iter<MAX_ITER);
 
    // 最終結果出力ブロック
    printf("\n 計算結果:\n");
    printf(" 収束固有値 (lambda)   = %10.6f %+.6fi\n", creal(lambda), cimag(lambda));
    printf(" 収束固有ベクトル (x)  = [(%.4f%+.4fi), (%.4f%+.4fi), (%.4f%+.4fi), (%.4f%+.4fi)]\n", 
           creal(x[0]), cimag(x[0]),
           creal(x[1]), cimag(x[1]),
           creal(x[2]), cimag(x[2]),
           creal(x[3]), cimag(x[3]));
    printf("==================================================================\n");
 
    // 固有値・固有ベクトルであることの確認（残差チェック）
    printf("\n 固有対の検証チェック (A*x - lambda*x):\n");
    for (int i = 0; i < N; i++) {
        double complex ax_i = 0.0 + 0.0*I;
        for (int j = 0; j < N; j++) {
            ax_i += A[i][j] * x[j];
        }
        double complex lx_i = lambda * x[i];
        double complex residue = ax_i - lx_i;
        printf(" 第%d行: A*x = (%9.5f%+.5fi), lambda*x = (%9.5f%+.5fi), 残差絶対値 = %e\n", 
               i + 1, 
               creal(ax_i), cimag(ax_i), 
               creal(lx_i), cimag(lx_i), 
               cabs(residue));
    }
    printf("==================================================================\n");
}
 
int main()
{
    printf("=== 行列1: 複素シフトによる残りの固有値スペックトルの探索 ===\n");
 
    // 近似固有値の設定（必要に応じて自由に変更可能）
    inverse_iteration(-5.4 + 0.0*I);
    inverse_iteration(-2.1 + 1.7*I);
    inverse_iteration(-2.1 - 1.7*I);
 
    return 0;
}