// na_rep09_hinagata.c
// ガウスの消去法

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<time.h>

#define N 128										// 行列の最大サイズ（未知数の個数 <= N - 1）
#define RANDOM(n) (rand()%(2*n+1)-n)				// 絶対値 |n| 以下の整数乱数 
#define REALRANDOM(x) (x*((2.0*rand()/RAND_MAX)-1))	// 絶対値 |x| 以下の実数乱数 

void matprint(double (*a)[N], int m, int n);		// 行列を出力する関数
void vecprint(double *v, int n);					// ベクトルを出力する関数

int main()
{
	double a0[N][N], b0[N];	// 連立一次方程式の係数
	double a[N][N];			// 拡大係数行列
	double x[N];			// 未知数ベクトル
	double z[N];			// 検算用ベクトル
	int size;				// 行列のサイズ
	int i, j, k;
	double w, e, f;
	// 変数は適宜定義せよ

///// 乱数の初期化 /////

	srand(time(NULL));	
	
///// 方程式の設定(鶴亀算の場合) /////

	size = 2; 	// 行列のサイズを設定
	a0[0][0] = 1; a0[0][1] = 1; a0[1][0] = 2; a0[1][1] = 4;	b0[0] = 7; b0[1] = 22;

/*
///// 方程式の設定(ランダムな整数係数の場合) /////

	size = 4; 	// 行列のサイズを設定
	for(i = 0; i < size; i++) for(j = 0; j < size; j++)	a0[i][j] = RANDOM(10);
	for(i = 0; i < size; i++) b0[i] = RANDOM(10);

///// 方程式の設定(ランダムな実数係数の場合) /////

	size = 4; 	// 行列のサイズを設定
	for(i = 0; i < size; i++) for(j = 0; j < size; j++) a0[i][j] = REALRANDOM(10);
	for(i = 0; i < size; i++) b0[i] = REALRANDOM(10);
*/

///// 拡大係数行列の設定 /////

	for(i = 0; i < size; i++){
		for(j = 0; j < size; j++) a[i][j] = a0[i][j];
		a[i][size] = b0[i];
	}
	printf("拡大系数行列 :\n"); 
	matprint(a, size, size + 1); printf("\n");

///// 前進消去（ピボット選択無し） /////

	for(i = 0; i < size - 1; i++){
		// k = i + 1,..., size - 1 について
		// k 行目から, i 行目の (a[k][i] / a[i][i])倍を引く

		//////////////////////////////////////////////
		/////     課題その１：ここを完成せよ     /////
		//////////////////////////////////////////////
	}

///// 後退代入 /////

	for(i = size - 1 ; i >= 0; i--){
		// x[i] の計算

		//////////////////////////////////////////////
		/////     課題その２：ここを完成せよ     /////
		//////////////////////////////////////////////
	}

///// 検算: z = Ax /////

	for(i = 0; i < size; i++){
		w = 0.0;
		for(j = 0; j < size; j++) w += a0[i][j] * x[j];
		z[i] = w;
	}

///// 結果出力 /////

	printf("->\n");
	printf("A' :\n"); matprint(a, size, size + 1); 		
	printf("x :\n"); vecprint(x, size); printf("\n");
	printf("検算 Ax :\n"); vecprint(z, size); printf("\n");
	printf("b :\n"); vecprint(b0, size); printf("\n");

	e = 0.0;
	for(i = 0; i < size; i++){
		f = fabs(z[i] - b0[i]);
		if(f > e) e = f;
	}
	printf("最大誤差 = %12.10lf\n\n", e);

	return 0;
}

///// 以下は出力のための関数定義部 /////

void matprint(double (*a)[N], int m, int n)
{
	int i, j;
	for(i = 0; i < m; i++){
		for(j = 0; j < n; j++) printf("%12.7lf",a[i][j]);
		printf("\n");
	}
}

void vecprint(double *v, int n)
{
	int i;
	for(i = 0; i < n; i++) printf("%12.7lf\n", v[i]);
}