#include <stdio.h>
#include <math.h>
#define N 4 // 行列の次元を設定
 
void householderHB(double A[N][N], double H[N][N]) {
    int i, j, k;
    double norm;   
    double Cop[N][N];
 
    // A マトリックス表示
        for (i = 0; i < N; i++) {
           for (j = 0; j < N; j++) {
                printf("A1[%d][%d]=%lf ", i, j, A[i][j]);
            }
         putchar('\n');        
        }
 
 
for (k = 0; k < N-1; k++) {
 
    // H マトリックスを単位行列に初期化
        for (i = 0; i < N; i++) {
           for (j = 0; j < N; j++) {
             H[i][j] = (i == j) ? 1.0 : 0.0;
                printf("H[%d][%d]=%lf ", i, j, H[i][j]);
            }
         putchar('\n');        
        }
 
    // x ベクトルと v ベクトルの計算
        double x[N]={0.0}, v[N]={0.0}, u[N]={0.0}, alpha, beta;
        norm = 0.0;
 
        for (i = k+1; i < N; i++) {
            x[i] = A[i][k];
            norm += x[i]*x[i];
        }    
        printf("k=%d ", k);
        norm = sqrt(norm);
        printf("norm=%lf ", norm);
 
    // alpha の計算
        alpha = (x[k+1] > 0) ? -norm : norm;
        printf("alpha=%lf ", alpha);
 
    // v ベクトルの計算
        for (i = 0; i < N ; i++) {
            if (i==k+1)
            v[i] = x[i]-alpha;
            else if (i<=k)
            v[i] = v[i];
            else 
            v[i] = x[i];
            x[i] = 0.0;
            printf("v[%d]=%lf ", i, v[i]);
        }
         putchar('\n');
 
    // u ベクトルの計算 (正規化された v)
        beta = 1.0/sqrt(alpha*(-v[k+1]));
          printf("beta=%lf ", beta);
 
        for (j = 0; j <= k; j++) {
                u[j] += beta * v[j];
              printf("u[%d]=%lf ", j, u[j]);
            v[j]=0.0;}
        for (j = k+1; j < N; j++) {
                u[j] += beta * v[j];
              printf("u[%d]=%lf ", j, u[j]);
            v[j]=0.0;}
            putchar('\n');        
 
    // H の更新
        for (i = 0; i < N; i++) {
            for (int j = 0; j < N; j++) {
                H[i][j]-=u[i]*u[j];
            }
        }
 
    // H の表示
        for (i = 0; i < N; i++) {
            for (j = 0; j < N; j++) {
                printf("H1[%d][%d]=%lf ", i, j, H[i][j]);
            }
            putchar('\n');
        }
 
 
    // A の更新
        for (i = 0; i < N; i++) {
            for (j = 0; j < N; j++) {
                for (int l = 0; l < N; l++) {
                    Cop[i][j] += H[i][l]*A[l][j];
                }
            }
        }
 
    // A の初期化
        for (i = k; i < N; i++) {
            for (j = k; j < N; j++) {
                A[i][j] = Cop[i][j];
                printf("A[%d][%d]=%lf ", i, j, A[i][j]);
                Cop[i][j] = 0.0;                
            }
            putchar('\n');
        }
 
    // H の初期化
        for (i = 0; i < N; i++) {
            for (j = 0; j < N; j++) {
                printf("H[%d][%d]=%lf ", i, j, H[i][j]);
                H[i][j] = 0.0;
            }
            putchar('\n');
        }
    }
}
 
void gram_schmidt_qr(double A[N][N], double Q[N][N], double R[N][N]) {
     double H[N][N];
 
    householderHB(A,H);
 
    int i, j, k;
    double temp;
 
    for (j = 0; j < N; j++) {
        for (i = 0; i < N; i++) {
            Q[i][j] = A[i][j];
        }
        for (k = 0; k < j; k++) {
            temp = 0.0;
            for (i = 0; i < N; i++) {
                temp += Q[i][j] * Q[i][k];
            }
            for (i = 0; i < N; i++) {
                Q[i][j] -= temp * Q[i][k];
            }
        }
        temp = 0.0;
        for (i = 0; i < N; i++) {
            temp += Q[i][j] * Q[i][j];
        }
        temp = sqrt(temp);
        for (i = 0; i < N; i++) {
            Q[i][j] /= temp;
        }
        for (i = 0; i <= j; i++) {
            temp = 0.0;
            for (k = 0; k < N; k++) {
                temp += A[k][j] * Q[k][i];
            }
            R[i][j] = temp;
        }
        for (i = j + 1; i < N; i++) {
            R[i][j] = 0.0;
        }
    }
}
 
int main() {
    double A[N][N] = {{16.0, -1.0, 1.0, 2.0},
                      {2.0, 12.0, 1.0, -1.0},
                      {1.0, 3.0, -24.0, 2.0},
                      {4.0, -2.0, 1.0, 20.0}};
     double R[N][N];
     double Q[N][N];
 
gram_schmidt_qr(A,Q,R);
 
 
    // ハウスホルダー変換後の行列 A を表示
    printf("Householder変換後の行列 A:\n");
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            printf("%lf ", A[i][j]);
        }
        printf("\n");
    }
 
    // QR分解後の行列 Q を表示
    printf("QR分解後の行列 Q:\n");
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            printf("Q[%d][%d]=%lf ",i,j, Q[i][j]);
        }
        printf("\n");
    }
     // QR分解後の行列 R を表示
    printf("QR分解後の行列 R:\n");
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            printf("R[%d][%d]=%lf ",i,j, R[i][j]);
        }
        printf("\n");
    }
    return 0;
}

I2luY2x1ZGUgPHN0ZGlvLmg+CiNpbmNsdWRlIDxtYXRoLmg+CiNkZWZpbmUgTiA0IC8vIOihjOWIl+OBruasoeWFg+OCkuioreWumgoKdm9pZCBob3VzZWhvbGRlckhCKGRvdWJsZSBBW05dW05dLCBkb3VibGUgSFtOXVtOXSkgewogICAgaW50IGksIGosIGs7CiAgICBkb3VibGUgbm9ybTsgICAKICAgIGRvdWJsZSBDb3BbTl1bTl07CgogICAgLy8gQSDjg57jg4jjg6rjg4Pjgq/jgrnooajnpLoKICAgICAgICBmb3IgKGkgPSAwOyBpIDwgTjsgaSsrKSB7CiAgICAgICAgICAgZm9yIChqID0gMDsgaiA8IE47IGorKykgewogICAgICAgICAgICAgICAgcHJpbnRmKCJBMVslZF1bJWRdPSVsZiAiLCBpLCBqLCBBW2ldW2pdKTsKICAgICAgICAgICAgfQogICAgICAgICBwdXRjaGFyKCdcbicpOyAgICAgICAgCiAgICAgICAgfQoKCmZvciAoayA9IDA7IGsgPCBOLTE7IGsrKykgewoKICAgIC8vIEgg44Oe44OI44Oq44OD44Kv44K544KS5Y2Y5L2N6KGM5YiX44Gr5Yid5pyf5YyWCiAgICAgICAgZm9yIChpID0gMDsgaSA8IE47IGkrKykgewogICAgICAgICAgIGZvciAoaiA9IDA7IGogPCBOOyBqKyspIHsKICAgICAgICAgICAgIEhbaV1bal0gPSAoaSA9PSBqKSA/IDEuMCA6IDAuMDsKICAgICAgICAgICAgICAgIHByaW50ZigiSFslZF1bJWRdPSVsZiAiLCBpLCBqLCBIW2ldW2pdKTsKICAgICAgICAgICAgfQogICAgICAgICBwdXRjaGFyKCdcbicpOyAgICAgICAgCiAgICAgICAgfQoKICAgIC8vIHgg44OZ44Kv44OI44Or44GoIHYg44OZ44Kv44OI44Or44Gu6KiI566XCiAgICAgICAgZG91YmxlIHhbTl09ezAuMH0sIHZbTl09ezAuMH0sIHVbTl09ezAuMH0sIGFscGhhLCBiZXRhOwogICAgICAgIG5vcm0gPSAwLjA7CgogICAgICAgIGZvciAoaSA9IGsrMTsgaSA8IE47IGkrKykgewogICAgICAgICAgICB4W2ldID0gQVtpXVtrXTsKICAgICAgICAgICAgbm9ybSArPSB4W2ldKnhbaV07CiAgICAgICAgfSAgICAKICAgICAgICBwcmludGYoIms9JWQgIiwgayk7CiAgICAgICAgbm9ybSA9IHNxcnQobm9ybSk7CiAgICAgICAgcHJpbnRmKCJub3JtPSVsZiAiLCBub3JtKTsKCiAgICAvLyBhbHBoYSDjga7oqIjnrpcKICAgICAgICBhbHBoYSA9ICh4W2srMV0gPiAwKSA/IC1ub3JtIDogbm9ybTsKICAgICAgICBwcmludGYoImFscGhhPSVsZiAiLCBhbHBoYSk7CgogICAgLy8gdiDjg5njgq/jg4jjg6vjga7oqIjnrpcKICAgICAgICBmb3IgKGkgPSAwOyBpIDwgTiA7IGkrKykgewogICAgICAgICAgICBpZiAoaT09aysxKQogICAgICAgICAgICB2W2ldID0geFtpXS1hbHBoYTsKICAgICAgICAgICAgZWxzZSBpZiAoaTw9aykKICAgICAgICAgICAgdltpXSA9IHZbaV07CiAgICAgICAgICAgIGVsc2UgCiAgICAgICAgICAgIHZbaV0gPSB4W2ldOwogICAgICAgICAgICB4W2ldID0gMC4wOwogICAgICAgICAgICBwcmludGYoInZbJWRdPSVsZiAiLCBpLCB2W2ldKTsKICAgICAgICB9CiAgICAgICAgIHB1dGNoYXIoJ1xuJyk7CiAgICAKICAgIC8vIHUg44OZ44Kv44OI44Or44Gu6KiI566XICjmraPopo/ljJbjgZXjgozjgZ8gdikKICAgICAgICBiZXRhID0gMS4wL3NxcnQoYWxwaGEqKC12W2srMV0pKTsKICAgICAgICAgIHByaW50ZigiYmV0YT0lbGYgIiwgYmV0YSk7CgogICAgICAgIGZvciAoaiA9IDA7IGogPD0gazsgaisrKSB7CiAgICAgICAgICAgICAgICB1W2pdICs9IGJldGEgKiB2W2pdOwogICAgICAgICAgICAgIHByaW50ZigidVslZF09JWxmICIsIGosIHVbal0pOwogICAgICAgICAgICB2W2pdPTAuMDt9CiAgICAgICAgZm9yIChqID0gaysxOyBqIDwgTjsgaisrKSB7CiAgICAgICAgICAgICAgICB1W2pdICs9IGJldGEgKiB2W2pdOwogICAgICAgICAgICAgIHByaW50ZigidVslZF09JWxmICIsIGosIHVbal0pOwogICAgICAgICAgICB2W2pdPTAuMDt9CiAgICAgICAgICAgIHB1dGNoYXIoJ1xuJyk7ICAgICAgICAKCiAgICAvLyBIIOOBruabtOaWsAogICAgICAgIGZvciAoaSA9IDA7IGkgPCBOOyBpKyspIHsKICAgICAgICAgICAgZm9yIChpbnQgaiA9IDA7IGogPCBOOyBqKyspIHsKICAgICAgICAgICAgICAgIEhbaV1bal0tPXVbaV0qdVtqXTsKICAgICAgICAgICAgfQogICAgICAgIH0KCiAgICAvLyBIIOOBruihqOekugogICAgICAgIGZvciAoaSA9IDA7IGkgPCBOOyBpKyspIHsKICAgICAgICAgICAgZm9yIChqID0gMDsgaiA8IE47IGorKykgewogICAgICAgICAgICAgICAgcHJpbnRmKCJIMVslZF1bJWRdPSVsZiAiLCBpLCBqLCBIW2ldW2pdKTsKICAgICAgICAgICAgfQogICAgICAgICAgICBwdXRjaGFyKCdcbicpOwogICAgICAgIH0KCgogICAgLy8gQSDjga7mm7TmlrAKICAgICAgICBmb3IgKGkgPSAwOyBpIDwgTjsgaSsrKSB7CiAgICAgICAgICAgIGZvciAoaiA9IDA7IGogPCBOOyBqKyspIHsKICAgICAgICAgICAgICAgIGZvciAoaW50IGwgPSAwOyBsIDwgTjsgbCsrKSB7CiAgICAgICAgICAgICAgICAgICAgQ29wW2ldW2pdICs9IEhbaV1bbF0qQVtsXVtqXTsKICAgICAgICAgICAgICAgIH0KICAgICAgICAgICAgfQogICAgICAgIH0KCiAgICAvLyBBIOOBruWIneacn+WMlgogICAgICAgIGZvciAoaSA9IGs7IGkgPCBOOyBpKyspIHsKICAgICAgICAgICAgZm9yIChqID0gazsgaiA8IE47IGorKykgewogICAgICAgICAgICAgICAgQVtpXVtqXSA9IENvcFtpXVtqXTsKICAgICAgICAgICAgICAgIHByaW50ZigiQVslZF1bJWRdPSVsZiAiLCBpLCBqLCBBW2ldW2pdKTsKICAgICAgICAgICAgICAgIENvcFtpXVtqXSA9IDAuMDsgICAgICAgICAgICAgICAgCiAgICAgICAgICAgIH0KICAgICAgICAgICAgcHV0Y2hhcignXG4nKTsKICAgICAgICB9CgogICAgLy8gSCDjga7liJ3mnJ/ljJYKICAgICAgICBmb3IgKGkgPSAwOyBpIDwgTjsgaSsrKSB7CiAgICAgICAgICAgIGZvciAoaiA9IDA7IGogPCBOOyBqKyspIHsKICAgICAgICAgICAgICAgIHByaW50ZigiSFslZF1bJWRdPSVsZiAiLCBpLCBqLCBIW2ldW2pdKTsKICAgICAgICAgICAgICAgIEhbaV1bal0gPSAwLjA7CiAgICAgICAgICAgIH0KICAgICAgICAgICAgcHV0Y2hhcignXG4nKTsKICAgICAgICB9CiAgICB9Cn0KCnZvaWQgZ3JhbV9zY2htaWR0X3FyKGRvdWJsZSBBW05dW05dLCBkb3VibGUgUVtOXVtOXSwgZG91YmxlIFJbTl1bTl0pIHsKICAgICBkb3VibGUgSFtOXVtOXTsKICAgIAogICAgaG91c2Vob2xkZXJIQihBLEgpOwoKICAgIGludCBpLCBqLCBrOwogICAgZG91YmxlIHRlbXA7CgogICAgZm9yIChqID0gMDsgaiA8IE47IGorKykgewogICAgICAgIGZvciAoaSA9IDA7IGkgPCBOOyBpKyspIHsKICAgICAgICAgICAgUVtpXVtqXSA9IEFbaV1bal07CiAgICAgICAgfQogICAgICAgIGZvciAoayA9IDA7IGsgPCBqOyBrKyspIHsKICAgICAgICAgICAgdGVtcCA9IDAuMDsKICAgICAgICAgICAgZm9yIChpID0gMDsgaSA8IE47IGkrKykgewogICAgICAgICAgICAgICAgdGVtcCArPSBRW2ldW2pdICogUVtpXVtrXTsKICAgICAgICAgICAgfQogICAgICAgICAgICBmb3IgKGkgPSAwOyBpIDwgTjsgaSsrKSB7CiAgICAgICAgICAgICAgICBRW2ldW2pdIC09IHRlbXAgKiBRW2ldW2tdOwogICAgICAgICAgICB9CiAgICAgICAgfQogICAgICAgIHRlbXAgPSAwLjA7CiAgICAgICAgZm9yIChpID0gMDsgaSA8IE47IGkrKykgewogICAgICAgICAgICB0ZW1wICs9IFFbaV1bal0gKiBRW2ldW2pdOwogICAgICAgIH0KICAgICAgICB0ZW1wID0gc3FydCh0ZW1wKTsKICAgICAgICBmb3IgKGkgPSAwOyBpIDwgTjsgaSsrKSB7CiAgICAgICAgICAgIFFbaV1bal0gLz0gdGVtcDsKICAgICAgICB9CiAgICAgICAgZm9yIChpID0gMDsgaSA8PSBqOyBpKyspIHsKICAgICAgICAgICAgdGVtcCA9IDAuMDsKICAgICAgICAgICAgZm9yIChrID0gMDsgayA8IE47IGsrKykgewogICAgICAgICAgICAgICAgdGVtcCArPSBBW2tdW2pdICogUVtrXVtpXTsKICAgICAgICAgICAgfQogICAgICAgICAgICBSW2ldW2pdID0gdGVtcDsKICAgICAgICB9CiAgICAgICAgZm9yIChpID0gaiArIDE7IGkgPCBOOyBpKyspIHsKICAgICAgICAgICAgUltpXVtqXSA9IDAuMDsKICAgICAgICB9CiAgICB9Cn0KCmludCBtYWluKCkgewogICAgZG91YmxlIEFbTl1bTl0gPSB7ezE2LjAsIC0xLjAsIDEuMCwgMi4wfSwKICAgICAgICAgICAgICAgICAgICAgIHsyLjAsIDEyLjAsIDEuMCwgLTEuMH0sCiAgICAgICAgICAgICAgICAgICAgICB7MS4wLCAzLjAsIC0yNC4wLCAyLjB9LAogICAgICAgICAgICAgICAgICAgICAgezQuMCwgLTIuMCwgMS4wLCAyMC4wfX07CiAgICAgZG91YmxlIFJbTl1bTl07CiAgICAgZG91YmxlIFFbTl1bTl07CgpncmFtX3NjaG1pZHRfcXIoQSxRLFIpOwoKCiAgICAvLyDjg4/jgqbjgrnjg5vjg6vjg4Djg7zlpInmj5vlvozjga7ooYzliJcgQSDjgpLooajnpLoKICAgIHByaW50ZigiSG91c2Vob2xkZXLlpInmj5vlvozjga7ooYzliJcgQTpcbiIpOwogICAgZm9yIChpbnQgaSA9IDA7IGkgPCBOOyBpKyspIHsKICAgICAgICBmb3IgKGludCBqID0gMDsgaiA8IE47IGorKykgewogICAgICAgICAgICBwcmludGYoIiVsZiAiLCBBW2ldW2pdKTsKICAgICAgICB9CiAgICAgICAgcHJpbnRmKCJcbiIpOwogICAgfQogICAgCiAgICAvLyBRUuWIhuino+W+jOOBruihjOWIlyBRIOOCkuihqOekugogICAgcHJpbnRmKCJRUuWIhuino+W+jOOBruihjOWIlyBROlxuIik7CiAgICBmb3IgKGludCBpID0gMDsgaSA8IE47IGkrKykgewogICAgICAgIGZvciAoaW50IGogPSAwOyBqIDwgTjsgaisrKSB7CiAgICAgICAgICAgIHByaW50ZigiUVslZF1bJWRdPSVsZiAiLGksaiwgUVtpXVtqXSk7CiAgICAgICAgfQogICAgICAgIHByaW50ZigiXG4iKTsKICAgIH0KICAgICAvLyBRUuWIhuino+W+jOOBruihjOWIlyBSIOOCkuihqOekugogICAgcHJpbnRmKCJRUuWIhuino+W+jOOBruihjOWIlyBSOlxuIik7CiAgICBmb3IgKGludCBpID0gMDsgaSA8IE47IGkrKykgewogICAgICAgIGZvciAoaW50IGogPSAwOyBqIDwgTjsgaisrKSB7CiAgICAgICAgICAgIHByaW50ZigiUlslZF1bJWRdPSVsZiAiLGksaiwgUltpXVtqXSk7CiAgICAgICAgfQogICAgICAgIHByaW50ZigiXG4iKTsKICAgIH0KICAgIHJldHVybiAwOwp9