#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}
};
 
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 = [(%.3f%+.3fi), (%.3f%+.3fi), (%.3f%+.3fi), (%.3f%+.3fi)], 固有値 = %.4f%+.4fi\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), (%.5f%+.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+CiNkZWZpbmUgTiA0CiNkZWZpbmUgRVBTIDFlLTgKI2RlZmluZSBNQVhfSVRFUiAxMDAKIApkb3VibGUgY29tcGxleCBBW05dW05dID0gewogICAgezEuMCwgMi4wLCAzLjAsIDUuMH0sCiAgICB7My4wLCAtNS4wLCAxLjAsIDQuMH0sCiAgICB7NS4wLCA5LjAsIDIuMCwgLTYuMH0sCiAgICB7MS4wLCA3LjAsIDQuMCwgMS4wfQp9OwoKdm9pZCBsdV9kZWNvbXBvc2l0aW9uKGRvdWJsZSBjb21wbGV4IEJbTl1bTl0sCiAgICAgICAgICAgICAgICAgICAgICBkb3VibGUgY29tcGxleCBMW05dW05dLAogICAgICAgICAgICAgICAgICAgICAgZG91YmxlIGNvbXBsZXggVVtOXVtOXSkKewogICAgZm9yKGludCBpPTA7aTxOO2krKyl7CiAgICAgICAgZm9yKGludCBqPTA7ajxOO2orKyl7CiAgICAgICAgICAgIExbaV1bal09KGk9PWopPzEuMCswLjAqSTowLjArMC4wKkk7CiAgICAgICAgICAgIFVbaV1bal09MC4wKzAuMCpJOwogICAgICAgIH0KICAgIH0KIAogICAgZm9yKGludCBpPTA7aTxOO2krKyl7CiAKICAgICAgICBmb3IoaW50IGo9aTtqPE47aisrKXsKICAgICAgICAgICAgZG91YmxlIGNvbXBsZXggc3VtID0gMC4wICsgMC4wKkk7CiAKICAgICAgICAgICAgZm9yKGludCBrPTA7azxpO2srKykKICAgICAgICAgICAgICAgIHN1bSArPSBMW2ldW2tdKlVba11bal07CiAKICAgICAgICAgICAgVVtpXVtqXSA9IEJbaV1bal0gLSBzdW07CiAgICAgICAgfQogCiAgICAgICAgZm9yKGludCBqPWkrMTtqPE47aisrKXsKICAgICAgICAgICAgZG91YmxlIGNvbXBsZXggc3VtID0gMC4wICsgMC4wKkk7CiAKICAgICAgICAgICAgZm9yKGludCBrPTA7azxpO2srKykKICAgICAgICAgICAgICAgIHN1bSArPSBMW2pdW2tdKlVba11baV07CiAKICAgICAgICAgICAgTFtqXVtpXSA9IChCW2pdW2ldLXN1bSkvVVtpXVtpXTsKICAgICAgICB9CiAgICB9Cn0KIAovKiDliY3pgLLku6PlhaUgKi8Kdm9pZCBmb3J3YXJkX3N1YnN0aXR1dGlvbihkb3VibGUgY29tcGxleCBMW05dW05dLAogICAgICAgICAgICAgICAgICAgICAgICAgIGRvdWJsZSBjb21wbGV4IHhbTl0sCiAgICAgICAgICAgICAgICAgICAgICAgICAgZG91YmxlIGNvbXBsZXggeVtOXSkKewogICAgZm9yKGludCBpPTA7aTxOO2krKyl7CiAKICAgICAgICBkb3VibGUgY29tcGxleCBzdW0gPSAwLjAgKyAwLjAqSTsKIAogICAgICAgIGZvcihpbnQgaj0wO2o8aTtqKyspCiAgICAgICAgICAgIHN1bSArPSBMW2ldW2pdKnlbal07CiAKICAgICAgICB5W2ldID0geFtpXSAtIHN1bTsKICAgIH0KfQogCi8qIOW+jOmAgOS7o+WFpSAqLwp2b2lkIGJhY2t3YXJkX3N1YnN0aXR1dGlvbihkb3VibGUgY29tcGxleCBVW05dW05dLAogICAgICAgICAgICAgICAgICAgICAgICAgICBkb3VibGUgY29tcGxleCB5W05dLAogICAgICAgICAgICAgICAgICAgICAgICAgICBkb3VibGUgY29tcGxleCB4W05dKQp7CiAgICBmb3IoaW50IGk9Ti0xO2k+PTA7aS0tKXsKIAogICAgICAgIGRvdWJsZSBjb21wbGV4IHN1bSA9IDAuMCArIDAuMCpJOwogCiAgICAgICAgZm9yKGludCBqPWkrMTtqPE47aisrKQogICAgICAgICAgICBzdW0gKz0gVVtpXVtqXSp4W2pdOwogCiAgICAgICAgeFtpXSA9ICh5W2ldLXN1bSkvVVtpXVtpXTsKICAgIH0KfQogCnZvaWQgaW52ZXJzZV9pdGVyYXRpb24oZG91YmxlIGNvbXBsZXggbGFtYmRhX2hhdCkKewogICAgZG91YmxlIGNvbXBsZXggQltOXVtOXTsKICAgIGRvdWJsZSBjb21wbGV4IExbTl1bTl07CiAgICBkb3VibGUgY29tcGxleCBVW05dW05dOwogCiAgICBkb3VibGUgY29tcGxleCB4W05dID0KICAgIHsKICAgICAgICAxLjArMC4wKkksCiAgICAgICAgMS4wKzAuMCpJLAogICAgICAgIDEuMCswLjAqSSwKICAgICAgICAxLjArMC4wKkkKICAgIH07CiAKICAgIGRvdWJsZSBjb21wbGV4IHlbTl07CiAgICBkb3VibGUgY29tcGxleCB4X25leHRbTl07CiAKICAgIGZvcihpbnQgaT0wO2k8TjtpKyspewogICAgICAgIGZvcihpbnQgaj0wO2o8TjtqKyspewogCiAgICAgICAgICAgIEJbaV1bal0gPSBBW2ldW2pdOwogCiAgICAgICAgICAgIGlmKGk9PWopCiAgICAgICAgICAgICAgICBCW2ldW2pdIC09IGxhbWJkYV9oYXQ7CiAgICAgICAgfQogICAgfQogCiAgICBsdV9kZWNvbXBvc2l0aW9uKEIsTCxVKTsKIAogICAgcHJpbnRmKCJcbiAgICAgICAgICAgICDov5HkvLzlm7rmnInlgKQgbGFtYmRhX2hhdCA9ICUuM2YgJSsuM2ZpKVxuIiwKICAgICAgICAgICBjcmVhbChsYW1iZGFfaGF0KSwgY2ltYWcobGFtYmRhX2hhdCkpOwoKIAogICAgaW50IGl0ZXI9MDsKICAgIGRvdWJsZSBkaWZmOwogICAgZG91YmxlIGNvbXBsZXggbGFtYmRhID0gMC4wICsgMC4wKkk7CiAKICAgIGRvewogCiAgICAgICAgZm9yd2FyZF9zdWJzdGl0dXRpb24oTCx4LHkpOwogICAgICAgIGJhY2t3YXJkX3N1YnN0aXR1dGlvbihVLHkseF9uZXh0KTsKIAogICAgICAgIGRvdWJsZSBtYXhfbWFnID0gMC4wOwogICAgICAgIGludCBpZHggPSAwOwogCiAgICAgICAgZm9yKGludCBpPTA7aTxOO2krKyl7CiAgICAgICAgICAgIGlmKGNhYnMoeF9uZXh0W2ldKSA+IG1heF9tYWcpewogICAgICAgICAgICAgICAgbWF4X21hZyA9IGNhYnMoeF9uZXh0W2ldKTsKICAgICAgICAgICAgICAgIGlkeCA9IGk7CiAgICAgICAgICAgIH0KICAgICAgICB9CiAKICAgICAgICBkb3VibGUgY29tcGxleCBzY2FsZSA9IHhfbmV4dFtpZHhdOwogCiAgICAgICAgZm9yKGludCBpPTA7aTxOO2krKykKICAgICAgICAgICAgeF9uZXh0W2ldIC89IHNjYWxlOwogCiAgICAgICAgZGlmZiA9IDAuMDsKIAogICAgICAgIGZvcihpbnQgaT0wO2k8TjtpKyspCiAgICAgICAgICAgIGRpZmYgKz0gY2Ficyh4X25leHRbaV0teFtpXSk7CiAKICAgICAgICBmb3IoaW50IGk9MDtpPE47aSsrKQogICAgICAgICAgICB4W2ldPXhfbmV4dFtpXTsKIAogICAgICAgIGxhbWJkYSA9IGxhbWJkYV9oYXQgKyAxLjAvc2NhbGU7CiAKICAgICAgICBwcmludGYoIuWPjeW+qSAlMmQ6IHggPSBbKCUuM2YlKy4zZmkpLCAoJS4zZiUrLjNmaSksICglLjNmJSsuM2ZpKSwgKCUuM2YlKy4zZmkpXSwg5Zu65pyJ5YCkID0gJS40ZiUrLjRmaVxuIiwgCiAgICAgICAgICAgICAgIGl0ZXIgKyAxLCAKICAgICAgICAgICAgICAgY3JlYWwoeFswXSksIGNpbWFnKHhbMF0pLAogICAgICAgICAgICAgICBjcmVhbCh4WzFdKSwgY2ltYWcoeFsxXSksCiAgICAgICAgICAgICAgIGNyZWFsKHhbMl0pLCBjaW1hZyh4WzJdKSwKICAgICAgICAgICAgICAgY3JlYWwoeFszXSksIGNpbWFnKHhbM10pLAogICAgICAgICAgICAgICBjcmVhbChsYW1iZGEpLCBjaW1hZyhsYW1iZGEpKTsKIAogICAgICAgIGl0ZXIrKzsKIAogICAgfXdoaWxlKGRpZmY+RVBTICYmIGl0ZXI8TUFYX0lURVIpOwogCgogICAgcHJpbnRmKCIg5Y+O5p2f5Zu65pyJ5YCkICAgPSAlMTAuNmYgJSsuNmZpXG4iLCBjcmVhbChsYW1iZGEpLCBjaW1hZyhsYW1iZGEpKTsKICAgIHByaW50ZigiIOWPjuadn+WbuuacieODmeOCr+ODiOODqyAoeCkgID0gWyglLjRmJSsuNGZpKSwgKCUuNGYlKy40ZmkpLCAoJS40ZiUrLjRmaSksICglLjVmJSsuNGZpKV1cbiIsIAogICAgICAgICAgIGNyZWFsKHhbMF0pLCBjaW1hZyh4WzBdKSwKICAgICAgICAgICBjcmVhbCh4WzFdKSwgY2ltYWcoeFsxXSksCiAgICAgICAgICAgY3JlYWwoeFsyXSksIGNpbWFnKHhbMl0pLAogICAgICAgICAgIGNyZWFsKHhbM10pLCBjaW1hZyh4WzNdKSk7CiAKICAgIHByaW50ZigiXG4g5Zu65pyJ5a++44Gu5qSc6Ki844OB44Kn44OD44KvIChBKnggLSBsYW1iZGEqeCk6XG4iKTsKICAgIGZvciAoaW50IGkgPSAwOyBpIDwgTjsgaSsrKSB7CiAgICAgICAgZG91YmxlIGNvbXBsZXggYXhfaSA9IDAuMCArIDAuMCpJOwogICAgICAgIGZvciAoaW50IGogPSAwOyBqIDwgTjsgaisrKSB7CiAgICAgICAgICAgIGF4X2kgKz0gQVtpXVtqXSAqIHhbal07CiAgICAgICAgfQogICAgICAgIGRvdWJsZSBjb21wbGV4IGx4X2kgPSBsYW1iZGEgKiB4W2ldOwogICAgICAgIGRvdWJsZSBjb21wbGV4IHJlc2lkdWUgPSBheF9pIC0gbHhfaTsKICAgICAgICBwcmludGYoIiDnrKwlZOihjDogQSp4ID0gKCU5LjVmJSsuNWZpKSwgbGFtYmRhKnggPSAoJTkuNWYlKy41ZmkpLCDmrovlt67ntbblr77lgKQgPSAlZVxuIiwgCiAgICAgICAgICAgICAgIGkgKyAxLCAKICAgICAgICAgICAgICAgY3JlYWwoYXhfaSksIGNpbWFnKGF4X2kpLCAKICAgICAgICAgICAgICAgY3JlYWwobHhfaSksIGNpbWFnKGx4X2kpLCAKICAgICAgICAgICAgICAgY2FicyhyZXNpZHVlKSk7CiAgICB9CiAgICBwcmludGYoIj09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PVxuIik7Cn0KIAppbnQgbWFpbigpCnsKICAgIGludmVyc2VfaXRlcmF0aW9uKC01LjQgKyAwLjAqSSk7CiAgICBpbnZlcnNlX2l0ZXJhdGlvbigtMi4xICsgMS43KkkpOwogICAgaW52ZXJzZV9pdGVyYXRpb24oLTIuMSAtIDEuNypJKTsKIAogICAgcmV0dXJuIDA7Cn0K