#include <stdio.h>
#include <stdlib.h>
 
// 行列出力用関数，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を解く------*/
for (int k = 0; k < n; k++) {
	double a1 = A[k * n + k];
    for (int j = k; j < n; j++) {
        A[k * n + j] /= a1;
    }
    b[k] /= a1;
    for (int i = k + 1; i < n; i++) {
    	double a2 = A[i * n + k];
        for (int j = k; j < n; j++) {
            A[i * n + j] -= a2 * A[k * n + j];
        }
        b[i] -= b[k] * a2;
    }
}
for (int k = n - 1; k > 0; k--) {
    for (int j = k - 1; j >= 0; j--) {
    	double a3 = A[j * n + k];
        for (int i = k; i < n; i++) {
            A[j * n + i] -= a3 * A[k * n + i];
        }
        b[j] -= b[k] * a3;
    }
}
 
/*------ここより上でAx=bを解く------*/
 
// 最後に残ったA,b,x,xの答えを出力
printf("A:\n");
printM(n, A);
printf("b:\n");
printV(n, b);
printf("A:\n");
printV(n, x);
printf("xAns:\n");
printV(n, xAns);
}
 
 

I2luY2x1ZGUgPHN0ZGlvLmg+CiNpbmNsdWRlIDxzdGRsaWIuaD4KIAovLyDooYzliJflh7rlipvnlKjplqLmlbDvvIxu44Gv6KGM5YiX44K144Kk44K677yMTeOBr+WFiOmgreODneOCpOODs+OCvwp2b2lkIHByaW50TShpbnQgbiwgZG91YmxlKiBNKSB7CmZvciAoaW50IGkgPSAwOyBpIDwgbjsgaSsrKSB7CmZvciAoaW50IGogPSAwOyBqIDwgbjsgaisrKSB7CnByaW50ZigiJWxmICIsIE1baSAqIG4gKyBqXSk7Cn0KcHJpbnRmKCJcbiIpOwp9Cn0KLy8g44OZ44Kv44OI44Or5Ye65Yqb55So6Zai5pWw77yMbuOBr+ODmeOCr+ODiOODq+OCteOCpOOCuu+8jFbjga/lhYjpoK3jg53jgqTjg7Pjgr8Kdm9pZCBwcmludFYoaW50IG4sIGRvdWJsZSogVikgewpmb3IgKGludCBqID0gMDsgaiA8IG47IGorKykgewpwcmludGYoIiU1LjNmICIsIFZbal0pOwp9CnByaW50ZigiXG4iKTsKfQogCmludCBtYWluKCkgewovLyDkubHmlbDjgrvjg4Pjg4gKc3JhbmQoMjkpOwovLyDooYzliJfjgrXjgqTjgrroqK3lrpoKaW50IG4gPSAzOwovLyBBLHgsYuOCkuagvOe0jeOBmeOCi+ODoeODouODquOBrueiuuS/nQpkb3VibGUqIEEgPSAoZG91YmxlKiltYWxsb2MobiAqIG4gKiBzaXplb2YoZG91YmxlKSk7CmRvdWJsZSogeCA9IChkb3VibGUqKW1hbGxvYyhuICogc2l6ZW9mKGRvdWJsZSkpOwpkb3VibGUqIHhBbnMgPSAoZG91YmxlKiltYWxsb2MobiAqIHNpemVvZihkb3VibGUpKTsKZG91YmxlKiBiID0gKGRvdWJsZSopbWFsbG9jKG4gKiBzaXplb2YoZG91YmxlKSk7CiAKLy8g44OZ44Kv44OI44Or44Gu5Yid5pyf5YCk44KS5YWl44KM44KLCmZvciAoaW50IGkgPSAwOyBpIDwgbjsgaSsrKSB7CnhbaV0gPSAwOwpiW2ldID0gMDsKeEFuc1tpXSA9IDUgKyBpOyAvLyDmnIDntYLnmoTjgavnrZTjgYjjga94PTUsNiw3Li4u44Go44Gq44KLCn0KIAovLyBB44Gu6KGM5YiX44KS55Sf5oiQCmZvciAoaW50IGkgPSAwOyBpIDwgbjsgaSsrKSB7CmZvciAoaW50IGogPSAwOyBqIDwgbjsgaisrKSB7CkFbaSAqIG4gKyBqXSA9IDEgKyAoZG91YmxlKShyYW5kKCkgJSAxMDApIC8gMTAwOyAvLyAxfjLjga7kubHmlbDjgafln4vjgoHjgosKfQovLyDlr77op5LmiJDliIbjga/jgZ3jga7ooYzjga7opoHntKDjga7nt4/lkozjgavjgZnjgoso5a++6KeS5oiQ5YiG44GM5aSn44GN44GE5pa544GM5a6J5a6a44GZ44KL44Gf44KBKQpmb3IgKGludCBqID0gMDsgaiA8IG47IGorKykgewpkb3VibGUgdG1wID0gMDsKdG1wICs9IEFbaSAqIG4gKyBqXTsKQVtpICogbiArIGldID0gdG1wOwp9Cn0KIAovLyBBLHhBbnPjgYvjgoli44OZ44Kv44OI44Or44KS6YCG566XCmZvciAoaW50IGkgPSAwOyBpIDwgbjsgaSsrKSB7CmZvciAoaW50IGogPSAwOyBqIDwgbjsgaisrKSB7CmJbaV0gKz0gQVtpICogbiArIGpdICogeEFuc1tqXTsKfQp9CiAKLyotLS0tLS3jgZPjgZPjgojjgorkuIvjgadBeD1i44KS6Kej44GPLS0tLS0tKi8KZm9yIChpbnQgayA9IDA7IGsgPCBuOyBrKyspIHsKCWRvdWJsZSBhMSA9IEFbayAqIG4gKyBrXTsKICAgIGZvciAoaW50IGogPSBrOyBqIDwgbjsgaisrKSB7CiAgICAgICAgQVtrICogbiArIGpdIC89IGExOwogICAgfQogICAgYltrXSAvPSBhMTsKICAgIGZvciAoaW50IGkgPSBrICsgMTsgaSA8IG47IGkrKykgewogICAgCWRvdWJsZSBhMiA9IEFbaSAqIG4gKyBrXTsKICAgICAgICBmb3IgKGludCBqID0gazsgaiA8IG47IGorKykgewogICAgICAgICAgICBBW2kgKiBuICsgal0gLT0gYTIgKiBBW2sgKiBuICsgal07CiAgICAgICAgfQogICAgICAgIGJbaV0gLT0gYltrXSAqIGEyOwogICAgfQp9CmZvciAoaW50IGsgPSBuIC0gMTsgayA+IDA7IGstLSkgewogICAgZm9yIChpbnQgaiA9IGsgLSAxOyBqID49IDA7IGotLSkgewogICAgCWRvdWJsZSBhMyA9IEFbaiAqIG4gKyBrXTsKICAgICAgICBmb3IgKGludCBpID0gazsgaSA8IG47IGkrKykgewogICAgICAgICAgICBBW2ogKiBuICsgaV0gLT0gYTMgKiBBW2sgKiBuICsgaV07CiAgICAgICAgfQogICAgICAgIGJbal0gLT0gYltrXSAqIGEzOwogICAgfQp9CgovKi0tLS0tLeOBk+OBk+OCiOOCiuS4iuOBp0F4PWLjgpLop6PjgY8tLS0tLS0qLwogCi8vIOacgOW+jOOBq+aui+OBo+OBn0EsYix4LHjjga7nrZTjgYjjgpLlh7rlipsKcHJpbnRmKCJBOlxuIik7CnByaW50TShuLCBBKTsKcHJpbnRmKCJiOlxuIik7CnByaW50VihuLCBiKTsKcHJpbnRmKCJBOlxuIik7CnByaW50VihuLCB4KTsKcHJpbnRmKCJ4QW5zOlxuIik7CnByaW50VihuLCB4QW5zKTsKfQoK