#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             近似固有値 lambda_hat = %.3f %+.3fi)\n",
           creal(lambda_hat), cimag(lambda_hat));
 
 
    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(" 収束固有値   = %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 固有対の検証チェック (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()
{
    inverse_iteration(-5.4 + 0.0*I);
    inverse_iteration(-2.1 + 1.7*I);
    inverse_iteration(-2.1 - 1.7*I);
 
    return 0;
}
 

I2luY2x1ZGUgPHN0ZGlvLmg+CiNpbmNsdWRlIDxtYXRoLmg+CiNpbmNsdWRlIDxjb21wbGV4Lmg+CiNkZWZpbmUgTiA0CiNkZWZpbmUgRVBTIDFlLTgKI2RlZmluZSBNQVhfSVRFUiAxMDAKIApkb3VibGUgY29tcGxleCBBW05dW05dID0gewogICAgezEuMCwgMi4wLCAzLjAsIDUuMH0sCiAgICB7My4wLCAtNS4wLCAxLjAsIDQuMH0sCiAgICB7NS4wLCA5LjAsIDIuMCwgLTYuMH0sCiAgICB7MS4wLCA3LjAsIDQuMCwgMS4wfQp9OwogCi8qIExV5YiG6KejICovCnZvaWQgbHVfZGVjb21wb3NpdGlvbihkb3VibGUgY29tcGxleCBCW05dW05dLAogICAgICAgICAgICAgICAgICAgICAgZG91YmxlIGNvbXBsZXggTFtOXVtOXSwKICAgICAgICAgICAgICAgICAgICAgIGRvdWJsZSBjb21wbGV4IFVbTl1bTl0pCnsKICAgIGZvcihpbnQgaT0wO2k8TjtpKyspewogICAgICAgIGZvcihpbnQgaj0wO2o8TjtqKyspewogICAgICAgICAgICBMW2ldW2pdPShpPT1qKT8xLjArMC4wKkk6MC4wKzAuMCpJOwogICAgICAgICAgICBVW2ldW2pdPTAuMCswLjAqSTsKICAgICAgICB9CiAgICB9CiAKICAgIGZvcihpbnQgaT0wO2k8TjtpKyspewogCiAgICAgICAgZm9yKGludCBqPWk7ajxOO2orKyl7CiAgICAgICAgICAgIGRvdWJsZSBjb21wbGV4IHN1bSA9IDAuMCArIDAuMCpJOwogCiAgICAgICAgICAgIGZvcihpbnQgaz0wO2s8aTtrKyspCiAgICAgICAgICAgICAgICBzdW0gKz0gTFtpXVtrXSpVW2tdW2pdOwogCiAgICAgICAgICAgIFVbaV1bal0gPSBCW2ldW2pdIC0gc3VtOwogICAgICAgIH0KIAogICAgICAgIGZvcihpbnQgaj1pKzE7ajxOO2orKyl7CiAgICAgICAgICAgIGRvdWJsZSBjb21wbGV4IHN1bSA9IDAuMCArIDAuMCpJOwogCiAgICAgICAgICAgIGZvcihpbnQgaz0wO2s8aTtrKyspCiAgICAgICAgICAgICAgICBzdW0gKz0gTFtqXVtrXSpVW2tdW2ldOwogCiAgICAgICAgICAgIExbal1baV0gPSAoQltqXVtpXS1zdW0pL1VbaV1baV07CiAgICAgICAgfQogICAgfQp9CiAKLyog5YmN6YCy5Luj5YWlICovCnZvaWQgZm9yd2FyZF9zdWJzdGl0dXRpb24oZG91YmxlIGNvbXBsZXggTFtOXVtOXSwKICAgICAgICAgICAgICAgICAgICAgICAgICBkb3VibGUgY29tcGxleCB4W05dLAogICAgICAgICAgICAgICAgICAgICAgICAgIGRvdWJsZSBjb21wbGV4IHlbTl0pCnsKICAgIGZvcihpbnQgaT0wO2k8TjtpKyspewogCiAgICAgICAgZG91YmxlIGNvbXBsZXggc3VtID0gMC4wICsgMC4wKkk7CiAKICAgICAgICBmb3IoaW50IGo9MDtqPGk7aisrKQogICAgICAgICAgICBzdW0gKz0gTFtpXVtqXSp5W2pdOwogCiAgICAgICAgeVtpXSA9IHhbaV0gLSBzdW07CiAgICB9Cn0KIAovKiDlvozpgIDku6PlhaUgKi8Kdm9pZCBiYWNrd2FyZF9zdWJzdGl0dXRpb24oZG91YmxlIGNvbXBsZXggVVtOXVtOXSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgZG91YmxlIGNvbXBsZXggeVtOXSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgZG91YmxlIGNvbXBsZXggeFtOXSkKewogICAgZm9yKGludCBpPU4tMTtpPj0wO2ktLSl7CiAKICAgICAgICBkb3VibGUgY29tcGxleCBzdW0gPSAwLjAgKyAwLjAqSTsKIAogICAgICAgIGZvcihpbnQgaj1pKzE7ajxOO2orKykKICAgICAgICAgICAgc3VtICs9IFVbaV1bal0qeFtqXTsKIAogICAgICAgIHhbaV0gPSAoeVtpXS1zdW0pL1VbaV1baV07CiAgICB9Cn0KIAp2b2lkIGludmVyc2VfaXRlcmF0aW9uKGRvdWJsZSBjb21wbGV4IGxhbWJkYV9oYXQpCnsKICAgIGRvdWJsZSBjb21wbGV4IEJbTl1bTl07CiAgICBkb3VibGUgY29tcGxleCBMW05dW05dOwogICAgZG91YmxlIGNvbXBsZXggVVtOXVtOXTsKIAogICAgZG91YmxlIGNvbXBsZXggeFtOXSA9CiAgICB7CiAgICAgICAgMS4wKzAuMCpJLAogICAgICAgIDEuMCswLjAqSSwKICAgICAgICAxLjArMC4wKkksCiAgICAgICAgMS4wKzAuMCpJCiAgICB9OwogCiAgICBkb3VibGUgY29tcGxleCB5W05dOwogICAgZG91YmxlIGNvbXBsZXggeF9uZXh0W05dOwogCiAgICBmb3IoaW50IGk9MDtpPE47aSsrKXsKICAgICAgICBmb3IoaW50IGo9MDtqPE47aisrKXsKIAogICAgICAgICAgICBCW2ldW2pdID0gQVtpXVtqXTsKIAogICAgICAgICAgICBpZihpPT1qKQogICAgICAgICAgICAgICAgQltpXVtqXSAtPSBsYW1iZGFfaGF0OwogICAgICAgIH0KICAgIH0KIAogICAgbHVfZGVjb21wb3NpdGlvbihCLEwsVSk7CiAKICAgIHByaW50ZigiXG4gICAgICAgICAgICAg6L+R5Ly85Zu65pyJ5YCkIGxhbWJkYV9oYXQgPSAlLjNmICUrLjNmaSlcbiIsCiAgICAgICAgICAgY3JlYWwobGFtYmRhX2hhdCksIGNpbWFnKGxhbWJkYV9oYXQpKTsKCiAKICAgIGludCBpdGVyPTA7CiAgICBkb3VibGUgZGlmZjsKICAgIGRvdWJsZSBjb21wbGV4IGxhbWJkYSA9IDAuMCArIDAuMCpJOwogCiAgICBkb3sKIAogICAgICAgIGZvcndhcmRfc3Vic3RpdHV0aW9uKEwseCx5KTsKICAgICAgICBiYWNrd2FyZF9zdWJzdGl0dXRpb24oVSx5LHhfbmV4dCk7CiAKICAgICAgICBkb3VibGUgbWF4X21hZyA9IDAuMDsKICAgICAgICBpbnQgaWR4ID0gMDsKIAogICAgICAgIGZvcihpbnQgaT0wO2k8TjtpKyspewogICAgICAgICAgICBpZihjYWJzKHhfbmV4dFtpXSkgPiBtYXhfbWFnKXsKICAgICAgICAgICAgICAgIG1heF9tYWcgPSBjYWJzKHhfbmV4dFtpXSk7CiAgICAgICAgICAgICAgICBpZHggPSBpOwogICAgICAgICAgICB9CiAgICAgICAgfQogCiAgICAgICAgZG91YmxlIGNvbXBsZXggc2NhbGUgPSB4X25leHRbaWR4XTsKIAogICAgICAgIGZvcihpbnQgaT0wO2k8TjtpKyspCiAgICAgICAgICAgIHhfbmV4dFtpXSAvPSBzY2FsZTsKIAogICAgICAgIGRpZmYgPSAwLjA7CiAKICAgICAgICBmb3IoaW50IGk9MDtpPE47aSsrKQogICAgICAgICAgICBkaWZmICs9IGNhYnMoeF9uZXh0W2ldLXhbaV0pOwogCiAgICAgICAgZm9yKGludCBpPTA7aTxOO2krKykKICAgICAgICAgICAgeFtpXT14X25leHRbaV07CiAKICAgICAgICBsYW1iZGEgPSBsYW1iZGFfaGF0ICsgMS4wL3NjYWxlOwogCiAgICAgICAgcHJpbnRmKCLlj43lvqkgJTJkOiB4ID0gWyglLjRmJSsuNGZpKSwgKCUuNGYlKy40ZmkpLCAoJS40ZiUrLjRmaSksICglLjRmJSsuNGZpKV0sIOWbuuacieWApCA9ICUuNmYlKy42ZmlcbiIsIAogICAgICAgICAgICAgICBpdGVyICsgMSwgCiAgICAgICAgICAgICAgIGNyZWFsKHhbMF0pLCBjaW1hZyh4WzBdKSwKICAgICAgICAgICAgICAgY3JlYWwoeFsxXSksIGNpbWFnKHhbMV0pLAogICAgICAgICAgICAgICBjcmVhbCh4WzJdKSwgY2ltYWcoeFsyXSksCiAgICAgICAgICAgICAgIGNyZWFsKHhbM10pLCBjaW1hZyh4WzNdKSwKICAgICAgICAgICAgICAgY3JlYWwobGFtYmRhKSwgY2ltYWcobGFtYmRhKSk7CiAKICAgICAgICBpdGVyKys7CiAKICAgIH13aGlsZShkaWZmPkVQUyAmJiBpdGVyPE1BWF9JVEVSKTsKIAoKICAgIHByaW50ZigiIOWPjuadn+WbuuacieWApCAgID0gJTEwLjZmICUrLjZmaVxuIiwgY3JlYWwobGFtYmRhKSwgY2ltYWcobGFtYmRhKSk7CiAgICBwcmludGYoIiDlj47mnZ/lm7rmnInjg5njgq/jg4jjg6sgKHgpICA9IFsoJS40ZiUrLjRmaSksICglLjRmJSsuNGZpKSwgKCUuNGYlKy40ZmkpLCAoJS40ZiUrLjRmaSldXG4iLCAKICAgICAgICAgICBjcmVhbCh4WzBdKSwgY2ltYWcoeFswXSksCiAgICAgICAgICAgY3JlYWwoeFsxXSksIGNpbWFnKHhbMV0pLAogICAgICAgICAgIGNyZWFsKHhbMl0pLCBjaW1hZyh4WzJdKSwKICAgICAgICAgICBjcmVhbCh4WzNdKSwgY2ltYWcoeFszXSkpOwogCiAgICBwcmludGYoIlxuIOWbuuacieWvvuOBruaknOiovOODgeOCp+ODg+OCryAoQSp4IC0gbGFtYmRhKngpOlxuIik7CiAgICBmb3IgKGludCBpID0gMDsgaSA8IE47IGkrKykgewogICAgICAgIGRvdWJsZSBjb21wbGV4IGF4X2kgPSAwLjAgKyAwLjAqSTsKICAgICAgICBmb3IgKGludCBqID0gMDsgaiA8IE47IGorKykgewogICAgICAgICAgICBheF9pICs9IEFbaV1bal0gKiB4W2pdOwogICAgICAgIH0KICAgICAgICBkb3VibGUgY29tcGxleCBseF9pID0gbGFtYmRhICogeFtpXTsKICAgICAgICBkb3VibGUgY29tcGxleCByZXNpZHVlID0gYXhfaSAtIGx4X2k7CiAgICAgICAgcHJpbnRmKCIg56ysJWTooYw6IEEqeCA9ICglOS41ZiUrLjVmaSksIGxhbWJkYSp4ID0gKCU5LjVmJSsuNWZpKSwg5q6L5beu57W25a++5YCkID0gJWVcbiIsIAogICAgICAgICAgICAgICBpICsgMSwgCiAgICAgICAgICAgICAgIGNyZWFsKGF4X2kpLCBjaW1hZyhheF9pKSwgCiAgICAgICAgICAgICAgIGNyZWFsKGx4X2kpLCBjaW1hZyhseF9pKSwgCiAgICAgICAgICAgICAgIGNhYnMocmVzaWR1ZSkpOwogICAgfQogICAgcHJpbnRmKCI9PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT1cbiIpOwp9CiAKaW50IG1haW4oKQp7CiAgICBpbnZlcnNlX2l0ZXJhdGlvbigtNS40ICsgMC4wKkkpOwogICAgaW52ZXJzZV9pdGVyYXRpb24oLTIuMSArIDEuNypJKTsKICAgIGludmVyc2VfaXRlcmF0aW9uKC0yLjEgLSAxLjcqSSk7CiAKICAgIHJldHVybiAwOwp9Cg==