#include <iostream>
using namespace std;
 
// 行列出力用関数，nは行列サイズ，Mは先頭ポインタ
void printM(int n, double* M) {
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
printf("%lf ", M[i * n + j]);
}
printf("\n");
}
}
// ベクトル出力用関数，nはベクトルサイズ，Vは先頭ポインタ
void printV(int n, double* V) {
for (int j = 0; j < n; j++) {
printf("%5.3f ", V[j]);
}
printf("\n");
}
 
int main() {
// 乱数セット
srand(29);
// 行列サイズ設定
int n = 3;
// A,x,bを格納するメモリの確保
double* A = (double*)malloc(n * n * sizeof(double));
double* x = (double*)malloc(n * sizeof(double));
double* xAns = (double*)malloc(n * sizeof(double));
double* b = (double*)malloc(n * sizeof(double));
 
// ベクトルの初期値を入れる
for (int i = 0; i < n; i++) {
x[i] = 0;
b[i] = 0;
xAns[i] = 5 + i; // 最終的に答えはx=5,6,7...となる
}
 
// Aの行列を生成
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
A[i * n + j] = 1 + (double)(rand() % 100) / 100; // 1~2の乱数で埋める
}
// 対角成分はその行の要素の総和にする(対角成分が大きい方が安定するため)
for (int j = 0; j < n; j++) {
double tmp = 0;
tmp += A[i * n + j];
A[i * n + i] = tmp;
}
}
 
// A,xAnsからbベクトルを逆算
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
b[i] += A[i * n + j] * xAns[j];
}
}
 
/*------ここより下でAx=bを解く------*/
	double* L = (double*)malloc(n * n * sizeof(double));
	double* U = (double*)malloc(n * n * sizeof(double));
	double* L2 = (double*)malloc(n * n * sizeof(double));
	double* U2 = (double*)malloc(n * n * sizeof(double));
	double* Z = (double*)malloc(n * sizeof(double));
	for(int i = 0; i < n; i++){
		for(int j = 0; j < n; j++){
			L[i * n + j] = 0;
			U[i * n + j] = A[i * n + j];
			L2[i * n + j] = 0;
			U2[i * n + j] = 0;
		}
		Z[i] = b[i];
	}

	for(int i = 0; i < n; i++){
		L[i * n + i] = U[i * n + i];
		for(int j = i; j < n; j++){
			U[i * n + j] /= L[i * n + i];
		}
		for(int j = i+1; j < n; j++){
			L[j * n + i] = U[j * n + i];
			for(int k = i; k < n; k++){
			U[j * n + k] -= U[i * n + k] * L[j * n + i];
			}
		}
	}
	for(int i = 0; i < n; i++){
		for(int j = 0; j < n; j++){
			L2[i * n + j] = L[i * n + j];
			U2[i * n + j] = U[i * n + j];
		}
	}
	for(int i = 0; i < n; i++){
		double ii = L2[i * n + i];
		L2[i * n + i] /= ii;
		Z[i] /= ii;
		for(int j = i+1; j < n; j++){
			double ji = L2[j * n + i];
			L2[j * n + i] -= L2[i * n + i] * ji;
			Z[j] -= Z[i] * ji;
		}
	}
	for(int i = 0; i < n; i++){
		x[i] = Z[i];
	}
	for(int i = n-1; i > 0; i--){
		for(int j = i-1; j >= 0; j--){
			double ji = U2[j * n + i];
			U2[j * n + i] -= U2[i * n + i] * ji;
			x[j] -= x[i] * ji;
		}
	}

/*------ここより上でAx=bを解く------*/
 
// 最後に残ったA,b,x,xの答えを出力
printf("A:\n");
printM(n, A);
printf("b:\n");
printV(n, b);
printf("L:\n");
printM(n, L);
printf("U:\n");
printM(n, U);
printf("Z:\n");
printV(n, Z);
printf("x:\n");
printV(n, x);
printf("xAns:\n");
printV(n, xAns);
}

I2luY2x1ZGUgPGlvc3RyZWFtPgp1c2luZyBuYW1lc3BhY2Ugc3RkOwogCi8vIOihjOWIl+WHuuWKm+eUqOmWouaVsO+8jG7jga/ooYzliJfjgrXjgqTjgrrvvIxN44Gv5YWI6aCt44Od44Kk44Oz44K/CnZvaWQgcHJpbnRNKGludCBuLCBkb3VibGUqIE0pIHsKZm9yIChpbnQgaSA9IDA7IGkgPCBuOyBpKyspIHsKZm9yIChpbnQgaiA9IDA7IGogPCBuOyBqKyspIHsKcHJpbnRmKCIlbGYgIiwgTVtpICogbiArIGpdKTsKfQpwcmludGYoIlxuIik7Cn0KfQovLyDjg5njgq/jg4jjg6vlh7rlipvnlKjplqLmlbDvvIxu44Gv44OZ44Kv44OI44Or44K144Kk44K677yMVuOBr+WFiOmgreODneOCpOODs+OCvwp2b2lkIHByaW50VihpbnQgbiwgZG91YmxlKiBWKSB7CmZvciAoaW50IGogPSAwOyBqIDwgbjsgaisrKSB7CnByaW50ZigiJTUuM2YgIiwgVltqXSk7Cn0KcHJpbnRmKCJcbiIpOwp9CiAKaW50IG1haW4oKSB7Ci8vIOS5seaVsOOCu+ODg+ODiApzcmFuZCgyOSk7Ci8vIOihjOWIl+OCteOCpOOCuuioreWumgppbnQgbiA9IDM7Ci8vIEEseCxi44KS5qC857SN44GZ44KL44Oh44Oi44Oq44Gu56K65L+dCmRvdWJsZSogQSA9IChkb3VibGUqKW1hbGxvYyhuICogbiAqIHNpemVvZihkb3VibGUpKTsKZG91YmxlKiB4ID0gKGRvdWJsZSopbWFsbG9jKG4gKiBzaXplb2YoZG91YmxlKSk7CmRvdWJsZSogeEFucyA9IChkb3VibGUqKW1hbGxvYyhuICogc2l6ZW9mKGRvdWJsZSkpOwpkb3VibGUqIGIgPSAoZG91YmxlKiltYWxsb2MobiAqIHNpemVvZihkb3VibGUpKTsKIAovLyDjg5njgq/jg4jjg6vjga7liJ3mnJ/lgKTjgpLlhaXjgozjgosKZm9yIChpbnQgaSA9IDA7IGkgPCBuOyBpKyspIHsKeFtpXSA9IDA7CmJbaV0gPSAwOwp4QW5zW2ldID0gNSArIGk7IC8vIOacgOe1gueahOOBq+etlOOBiOOBr3g9NSw2LDcuLi7jgajjgarjgosKfQogCi8vIEHjga7ooYzliJfjgpLnlJ/miJAKZm9yIChpbnQgaSA9IDA7IGkgPCBuOyBpKyspIHsKZm9yIChpbnQgaiA9IDA7IGogPCBuOyBqKyspIHsKQVtpICogbiArIGpdID0gMSArIChkb3VibGUpKHJhbmQoKSAlIDEwMCkgLyAxMDA7IC8vIDF+MuOBruS5seaVsOOBp+Wfi+OCgeOCiwp9Ci8vIOWvvuinkuaIkOWIhuOBr+OBneOBruihjOOBruimgee0oOOBrue3j+WSjOOBq+OBmeOCiyjlr77op5LmiJDliIbjgYzlpKfjgY3jgYTmlrnjgYzlronlrprjgZnjgovjgZ/jgoEpCmZvciAoaW50IGogPSAwOyBqIDwgbjsgaisrKSB7CmRvdWJsZSB0bXAgPSAwOwp0bXAgKz0gQVtpICogbiArIGpdOwpBW2kgKiBuICsgaV0gPSB0bXA7Cn0KfQogCi8vIEEseEFuc+OBi+OCiWLjg5njgq/jg4jjg6vjgpLpgIbnrpcKZm9yIChpbnQgaSA9IDA7IGkgPCBuOyBpKyspIHsKZm9yIChpbnQgaiA9IDA7IGogPCBuOyBqKyspIHsKYltpXSArPSBBW2kgKiBuICsgal0gKiB4QW5zW2pdOwp9Cn0KIAovKi0tLS0tLeOBk+OBk+OCiOOCiuS4i+OBp0F4PWLjgpLop6PjgY8tLS0tLS0qLwoJZG91YmxlKiBMID0gKGRvdWJsZSopbWFsbG9jKG4gKiBuICogc2l6ZW9mKGRvdWJsZSkpOwoJZG91YmxlKiBVID0gKGRvdWJsZSopbWFsbG9jKG4gKiBuICogc2l6ZW9mKGRvdWJsZSkpOwoJZG91YmxlKiBMMiA9IChkb3VibGUqKW1hbGxvYyhuICogbiAqIHNpemVvZihkb3VibGUpKTsKCWRvdWJsZSogVTIgPSAoZG91YmxlKiltYWxsb2MobiAqIG4gKiBzaXplb2YoZG91YmxlKSk7Cglkb3VibGUqIFogPSAoZG91YmxlKiltYWxsb2MobiAqIHNpemVvZihkb3VibGUpKTsKCWZvcihpbnQgaSA9IDA7IGkgPCBuOyBpKyspewoJCWZvcihpbnQgaiA9IDA7IGogPCBuOyBqKyspewoJCQlMW2kgKiBuICsgal0gPSAwOwoJCQlVW2kgKiBuICsgal0gPSBBW2kgKiBuICsgal07CgkJCUwyW2kgKiBuICsgal0gPSAwOwoJCQlVMltpICogbiArIGpdID0gMDsKCQl9CgkJWltpXSA9IGJbaV07Cgl9CgoJZm9yKGludCBpID0gMDsgaSA8IG47IGkrKyl7CgkJTFtpICogbiArIGldID0gVVtpICogbiArIGldOwoJCWZvcihpbnQgaiA9IGk7IGogPCBuOyBqKyspewoJCQlVW2kgKiBuICsgal0gLz0gTFtpICogbiArIGldOwoJCX0KCQlmb3IoaW50IGogPSBpKzE7IGogPCBuOyBqKyspewoJCQlMW2ogKiBuICsgaV0gPSBVW2ogKiBuICsgaV07CgkJCWZvcihpbnQgayA9IGk7IGsgPCBuOyBrKyspewoJCQlVW2ogKiBuICsga10gLT0gVVtpICogbiArIGtdICogTFtqICogbiArIGldOwoJCQl9CgkJfQoJfQoJZm9yKGludCBpID0gMDsgaSA8IG47IGkrKyl7CgkJZm9yKGludCBqID0gMDsgaiA8IG47IGorKyl7CgkJCUwyW2kgKiBuICsgal0gPSBMW2kgKiBuICsgal07CgkJCVUyW2kgKiBuICsgal0gPSBVW2kgKiBuICsgal07CgkJfQoJfQoJZm9yKGludCBpID0gMDsgaSA8IG47IGkrKyl7CgkJZG91YmxlIGlpID0gTDJbaSAqIG4gKyBpXTsKCQlMMltpICogbiArIGldIC89IGlpOwoJCVpbaV0gLz0gaWk7CgkJZm9yKGludCBqID0gaSsxOyBqIDwgbjsgaisrKXsKCQkJZG91YmxlIGppID0gTDJbaiAqIG4gKyBpXTsKCQkJTDJbaiAqIG4gKyBpXSAtPSBMMltpICogbiArIGldICogamk7CgkJCVpbal0gLT0gWltpXSAqIGppOwoJCX0KCX0KCWZvcihpbnQgaSA9IDA7IGkgPCBuOyBpKyspewoJCXhbaV0gPSBaW2ldOwoJfQoJZm9yKGludCBpID0gbi0xOyBpID4gMDsgaS0tKXsKCQlmb3IoaW50IGogPSBpLTE7IGogPj0gMDsgai0tKXsKCQkJZG91YmxlIGppID0gVTJbaiAqIG4gKyBpXTsKCQkJVTJbaiAqIG4gKyBpXSAtPSBVMltpICogbiArIGldICogamk7CgkJCXhbal0gLT0geFtpXSAqIGppOwoJCX0KCX0KCi8qLS0tLS0t44GT44GT44KI44KK5LiK44GnQXg9YuOCkuino+OBjy0tLS0tLSovCiAKLy8g5pyA5b6M44Gr5q6L44Gj44GfQSxiLHgseOOBruetlOOBiOOCkuWHuuWKmwpwcmludGYoIkE6XG4iKTsKcHJpbnRNKG4sIEEpOwpwcmludGYoImI6XG4iKTsKcHJpbnRWKG4sIGIpOwpwcmludGYoIkw6XG4iKTsKcHJpbnRNKG4sIEwpOwpwcmludGYoIlU6XG4iKTsKcHJpbnRNKG4sIFUpOwpwcmludGYoIlo6XG4iKTsKcHJpbnRWKG4sIFopOwpwcmludGYoIng6XG4iKTsKcHJpbnRWKG4sIHgpOwpwcmludGYoInhBbnM6XG4iKTsKcHJpbnRWKG4sIHhBbnMpOwp9Cg==