#ifndef LINEAR_2023_05_01
#define LINEAR_2023_05_01

#include <iostream>
#include <iomanip>
#include <math.h>

// エラー処理
void Error(const char *s)
{
	std::cerr << s << std::endl;
	std::exit(EXIT_FAILURE);
}

class Matrix;

// 整数ベクトルクラス
class IntVector {
private:
	int *Array;
	int Len;
public:
	IntVector();
	IntVector(int n);
	~IntVector();
	IntVector(const IntVector& obj);
	IntVector operator =(const IntVector& obj);
	int& operator[](int i);
	const int& operator[](int i) const;
};

// ベクトルクラス
class Vector {
	friend Matrix;
private:
	double *Array;
	int Len;
	void New(int n);
public:
	Vector();
	Vector(int n);
	~Vector();
	Vector(const Vector& obj);
	Vector operator =(const Vector& obj);
	double& operator[](int i);
	const double& operator[](int i) const;
	Vector operator +(const Vector& obj) const;
	Vector operator +() const;
	Vector operator -(const Vector& obj) const;
	Vector operator -() const;
	Vector operator +=(const Vector& obj);
	Vector operator -=(const Vector& obj);
	double operator *(const Vector& obj) const;
	Vector operator *(double a) const;
	Vector operator *=(double a);
	friend Vector operator *(double a, const Vector& obj);
	Vector operator /(double a) const;
	Vector operator /=(double a);
	Vector operator *(const Matrix& obj);
	Vector Shift(int k) const;
	double Norm1() const;
	double Norm2() const;
	double NormInfty() const;
	int Size() const;
	void Show(int n) const;
};

// 行列クラス
class Matrix {
	friend Vector;
private:
	class SubIndexClass {
	public:
		mutable int RowIndex;
		int ColLen;
		double *Array;
		double& operator[](int i);
		const double& operator[](int i) const;
	};
	SubIndexClass SubIndex;
	class WorkingMemory {// ワーキングメモリ管理クラス
	public:
		double *Array;
		int size;
		void Keep(int n);
		Vector v;
		WorkingMemory();
		~WorkingMemory();
	};
	static WorkingMemory WM;
	double *Array;
	int RowLen;
	int ColLen;
	void New(int m, int n);
public:
	Matrix();
	Matrix(int m, int n);
	~Matrix();
	Matrix(const Matrix& obj);
	Matrix operator =(const Matrix& obj);
	SubIndexClass& operator[](int i);
	const SubIndexClass& operator[](int i) const;
	Matrix operator +(const Matrix& obj) const;
	Matrix operator +() const;
	Matrix operator -(const Matrix& obj) const;
	Matrix operator -() const;
	Matrix operator +=(const Matrix& obj);
	Matrix operator -=(const Matrix& obj);
	Matrix operator *(double a) const;
	Matrix operator *=(double a);
	friend Matrix operator *(double a, const Matrix& obj);
	Matrix operator *(const Matrix& obj) const;
	Vector operator *(const Vector& obj) const;
	Matrix operator /(double a) const;
	Matrix operator /=(double a);
	Matrix Transpose() const;
	Vector RowVector(int k) const;
	Vector ColVector(int k) const;
	void SetRowVector(int k, const Vector& v);
	void SetColVector(int k, const Vector& v);
	Vector Gauss(const Vector& x);
	int RowSize() const;
	int ColSize() const;
	void Show(int n) const;
};

Matrix::WorkingMemory Matrix::WM;

// 整数ベクトルクラスのメンバ関数

// デフォルトコンストラクタ
IntVector::IntVector()
{
	Len = 0;
	Array = NULL;
}
// コンストラクタ
IntVector::IntVector(int n)
{
	if(n<=0){
		Error("Illegal size");
	}
	Len = n;
	Array = new int[n+1];
	for(int i=0; i<=Len; i++){
		Array[i] = 0;
	}
}
// デストラクタ
IntVector::~IntVector()
{
	delete[] Array;
}
// コピーコンストラクタ
IntVector::IntVector(const IntVector& obj)
{
	Len = obj.Len;
	if(Len<=0){
		Len = 0;
		Array = NULL;
	} else {
		Array = new int[Len+1];
		for(int i=0; i<=Len; i++){
			Array[i] = obj.Array[i];
		}
	}
}
// 代入演算子のオーバーロード
IntVector IntVector::operator =(const IntVector& obj)
{
	if(obj.Len<=0){
		delete[] Array;
		Len = 0;
		Array = NULL;
		return *this;
	}
	if(obj.Len!=Len){
		delete[] Array;
		Len = obj.Len;
		Array = new int[Len+1];
	}
	for(int i=0; i<=Len; i++){
		Array[i] = obj.Array[i];
	}
	return *this;
}
// []演算子のオーバーロード
int& IntVector::operator[](int i)
{
	if(Len<=0){
		Error("Null vector");
	}
	if(i<0 || i>Len){
		Error("Out of range");
	}
	return Array[i];
}
// []演算子のオーバーロード（constバージョン）
const int& IntVector::operator[](int i) const
{
	if(Len<=0){
		Error("Null vector");
	}
	if(i<0 || i>Len){
		Error("Out of range");
	}
	return Array[i];
}

// ベクトルクラスのメンバ関数

// デフォルトコンストラクタ
Vector::Vector()
{
	Len = 0;
	Array = NULL;
}
// コンストラクタ
Vector::Vector(int n)
{
	if(n<=0){
		Error("Illegal size");
	}
	Len = n;
	Array = new double[n+1];
	for(int i=0; i<=Len; i++){
		Array[i] = 0;
	}
}
// デストラクタ
Vector::~Vector()
{
	delete[] Array;
}
// コピーコンストラクタ
Vector::Vector(const Vector& obj)
{
	Len = obj.Len;
	if(Len<=0){
		Len = 0;
		Array = NULL;
	} else {
		Array = new double[Len+1];
		for(int i=0; i<=Len; i++){
			Array[i] = obj.Array[i];
		}
	}
}
// 代入演算子のオーバーロード
Vector Vector::operator =(const Vector& obj)
{
	if(obj.Len<=0){
		delete[] Array;
		Len = 0;
		Array = NULL;
		return *this;
	}
	if(obj.Len!=Len){
		delete[] Array;
		Len = obj.Len;
		Array = new double[Len+1];
	}
	for(int i=0; i<=Len; i++){
		Array[i] = obj.Array[i];
	}
	return *this;
}
// 内部用初期化関数
void Vector::New(int n)
{
	Len = n;
	Array = new double[n+1];
}
// []演算子のオーバーロード
double& Vector::operator[](int i)
{
	if(Len<=0){
		Error("Null vector");
	}
	if(i<0 || i>Len){
		Error("Out of range");
	}
	return Array[i];
}
// []演算子のオーバーロード（constバージョン）
const double& Vector::operator[](int i) const
{
	if(Len<=0){
		Error("Null vector");
	}
	if(i<0 || i>Len){
		Error("Out of range");
	}
	return Array[i];
}
// +演算子のオーバーロード
Vector Vector::operator +(const Vector& obj) const
{
	if(Len<=0 || obj.Len<=0){
		Error("Null vector");
	}
	if(obj.Len != Len){
		Error("Size doesn't much");
	}
	Vector v;
	v.New(Len);
	for(int i=0; i<=Len; i++){
		v.Array[i] = Array[i] + obj.Array[i];
	}
	return v;
}
// +演算子のオーバーロード（単項演算子）
Vector Vector::operator +() const
{
	if(Len<=0){
		Error("Null vector");
	}
	return *this;
}
// -演算子のオーバーロード
Vector Vector::operator -(const Vector& obj) const
{
	if(Len<=0 || obj.Len<=0){
		Error("Null vector");
	}
	if(obj.Len != Len){
		Error("Size doesn't much");
	}
	Vector v;
	v.New(Len);
	for(int i=0; i<=Len; i++){
		v.Array[i] = Array[i] - obj.Array[i];
	}
	return v;
}
// -演算子のオーバーロード（単項演算子）
Vector Vector::operator -() const
{
	if(Len<=0){
		Error("Null vector");
	}
	Vector v;
	v.New(Len);
	for(int i=0; i<=Len; i++){
		v.Array[i] = -Array[i];
	}
	return v;
}
// +=演算子のオーバーロード
Vector Vector::operator +=(const Vector& obj)
{
	if(Len<=0 || obj.Len<=0){
		Error("Null vector");
	}
	if(obj.Len != Len){
		Error("Size doesn't much");
	}
	for(int i=0; i<=Len; i++){
		Array[i] += obj.Array[i];
	}
	return *this;
}
// -=演算子のオーバーロード
Vector Vector::operator -=(const Vector& obj)
{
	if(Len<=0 || obj.Len<=0){
		Error("Null vector");
	}
	if(obj.Len != Len){
		Error("Size doesn't much");
	}
	for(int i=0; i<=Len; i++){
		Array[i] -= obj.Array[i];
	}
	return *this;
}
// *演算子のオーバーロード（内積）
double Vector::operator *(const Vector& obj) const
{
	if(Len<=0 || obj.Len<=0){
		Error("Null vector");
	}
	if(obj.Len != Len){
		Error("Size doesn't much");
	}
	double a = 0;
	for(int i=1; i<=Len; i++){
		a += Array[i]*obj.Array[i];
	}
	return a;
}
// *演算子のオーバーロード（スカラー倍）
Vector Vector::operator *(double a) const
{
	if(Len<=0){
		Error("Null vector");
	}
	Vector v;
	v.New(Len);
	for(int i=0; i<=Len; i++){
		v.Array[i] = Array[i] * a;
	}
	return v;
}
// *=演算子のオーバーロード（スカラー倍）
Vector Vector::operator *=(double a)
{
	if(Len<=0){
		Error("Null vector");
	}
	for(int i=0; i<=Len; i++){
		Array[i] *= a;
	}
	return *this;
}
// *演算子のオーバーロード（左からのスカラー倍）
Vector operator *(double a, const Vector& obj)
{
	if(obj.Len<=0){
		Error("Null vector");
	}
	Vector v;
	v.New(obj.Len);
	for(int i=0; i<=obj.Len; i++){
		v.Array[i] = a * obj.Array[i];
	}
	return v;
}
// /演算子のオーバーロード（スカラー除算）
Vector Vector::operator /(double a) const
{
	if(Len<=0){
		Error("Null vector");
	}
	if(a==0){
		Error("Divide by zero");
	}
	Vector v;
	v.New(Len);
	for(int i=0; i<=Len; i++){
		v.Array[i] = Array[i] / a;
	}
	return v;
}
// /=演算子のオーバーロード（スカラー除算）
Vector Vector::operator /=(double a)
{
	if(Len<=0){
		Error("Null vector");
	}
	if(a==0){
		Error("Divide by zero");
	}
	for(int i=0; i<=Len; i++){
		Array[i] /= a;
	}
	return *this;
}
// *演算子のオーバーロード（行列との積）
Vector Vector::operator *(const Matrix& obj)
{
	if(Len<=0){
		Error("Null vector");
	}
	if(obj.RowLen<=0 || obj.ColLen<=0){
		Error("Null matrix");
	}
	if(Len!=obj.RowLen){
		Error("Size doesn't much");
	}
	Vector v;
	v.New(obj.ColLen);
	v.Array[0] = 0;
	const int c = obj.ColLen + 1;
	for(int i=1; i<=obj.ColLen; i++){
		double a = 0;
		for(int j=1; j<=Len; j++){
			a += Array[j] * obj.Array[j*c+i];
		}
		v.Array[i] = a;
	}
	return v;
}
// シフト
Vector Vector::Shift(int k) const
{
	if(Len<=0){
		Error("Null vector");
	}
	Vector v;
	v.New(Len);
	for(int i=0, j=-k; i<=Len; i++,j++){
		if(j>=0 && j<=Len){
			v.Array[i] = Array[j];
		} else {
			v.Array[i] = 0;
		}
	}
	return v;
}
// 1ノルム
double Vector::Norm1() const
{
	if(Len<=0){
		Error("Null vector");
	}
	double r = 0;
	for(int i=1; i<=Len; i++){
		r += fabs(Array[i]);
	}
	return r;
}
// 2ノルム
double Vector::Norm2() const
{
	if(Len<=0){
		Error("Null vector");
	}
	double r = 0;
	for(int i=1; i<=Len; i++){
		r += Array[i]*Array[i];
	}
	return sqrt(r);
}
// 無限大ノルム
double Vector::NormInfty() const
{
	if(Len<=0){
		Error("Null vector");
	}
	double r = 0;
	for(int i=1; i<=Len; i++){
		double a = fabs(Array[i]);
		if(a>r){
			r = a;
		}
	}
	return r;
}
// ベクトルのサイズ
int Vector::Size() const
{
	return Len;
}
// 表示
void Vector::Show(int n) const
{
	if(Len<=0 && n<=0){
		return;
	}
	for(int i=1; i<=Len; i++){
		if(Array[i]==0){
			std::cout << std::fixed << std::setprecision(n-1) << Array[i] << std::endl;
			continue;
		}
		int k;
		double b = 1;
		double f = fabs(Array[i])*100000;
		for(k=0; k<=n+5; k++){
			if(b>f){
				break;
			}
			b *= 10;
		}
		if(k==0 || k>n+5){
			std::cout << std::scientific << std::setprecision(n-1) << Array[i];
		} else {
			std::cout << std::fixed << std::setprecision(n+5-k) << Array[i];
		}
		std::cout << std::endl;
	}
}

// 行列クラスのメンバ関数

// インデックス処理
double& Matrix::SubIndexClass::operator[](int i)
{
	if(i<0 || i>ColLen){
		Error("Out of range");
	}
	return Array[RowIndex*(ColLen+1)+i];
}
const double& Matrix::SubIndexClass::operator[](int i) const
{
	if(i<0 || i>ColLen){
		Error("Out of range");
	}
	return Array[RowIndex*(ColLen+1)+i];
}

// ワーキングメモリー処理
Matrix::WorkingMemory::WorkingMemory()
{
	size = 0;
	Array = NULL;
}
Matrix::WorkingMemory::~WorkingMemory()
{
	delete[] Array;
}
void Matrix::WorkingMemory::Keep(int n)
{
	if(n>size){
		delete[] Array;
		Array = new double[n];
		size = n;
	}
}

// デフォルトコンストラクタ
Matrix::Matrix()
{
	RowLen = 0;
	ColLen = 0;
	Array = NULL;
}
// コンストラクタ
Matrix::Matrix(int m, int n)
{
	if(m<=0 || n<=0){
		Error("Illegal size");
	}
	RowLen = m;
	ColLen = n;
	const int g = (m+1) * (n+1);
	Array = new double[g];
	for(int i=0; i<g; i++){
		Array[i] = 0;
	}
	SubIndex.ColLen = ColLen;
	SubIndex.Array = Array;
}
// デストラクタ
Matrix::~Matrix()
{
	delete[] Array;
}
// コピーコンストラクタ
Matrix::Matrix(const Matrix& obj)
{
	RowLen = obj.RowLen;
	ColLen = obj.ColLen;
	if(RowLen<=0 || ColLen<=0){
		RowLen = 0;
		ColLen = 0;
		Array = NULL;
	} else {
		const int g = (RowLen+1)*(ColLen+1);
		Array = new double[g];
		for(int i=0; i<g; i++){
			Array[i] = obj.Array[i];
		}
		SubIndex.ColLen = ColLen;
		SubIndex.Array = Array;
	}
}
// 代入演算子のオーバーロード
Matrix Matrix::operator =(const Matrix& obj)
{
	if(obj.RowLen<=0 || obj.ColLen<=0){
		delete[] Array;
		RowLen = 0;
		ColLen = 0;
		Array = NULL;
		return *this;
	}
	ColLen = obj.ColLen;
	RowLen = obj.RowLen;
	const int g = (obj.RowLen+1)*(obj.ColLen+1);
	if((RowLen+1)*(ColLen+1)!=g){
		delete[] Array;
		Array = new double[g];
	}
	for(int i=0; i<g; i++){
		Array[i] = obj.Array[i];
	}
	SubIndex.ColLen = ColLen;
	SubIndex.Array = Array;
	return *this;
}
// 内部用初期化関数
void Matrix::New(int m, int n)
{
	RowLen = m;
	ColLen = n;
	Array = new double[(m+1) * (n+1)];
	SubIndex.ColLen = ColLen;
	SubIndex.Array = Array;
}
// []演算子のオーバーロード
Matrix::SubIndexClass& Matrix::operator[](int i)
{
	if(RowLen<=0 || ColLen<=0){
		Error("Null matrix");
	}
	if(i<0 || i>ColLen){
		Error("Out of range");
	}
	SubIndex.RowIndex = i;
	return SubIndex;
}
// []演算子のオーバーロード（constバージョン）
const Matrix::SubIndexClass& Matrix::operator[](int i) const
{
	if(RowLen<=0 || ColLen<=0){
		Error("Null matrix");
	}
	if(i<0 || i>ColLen){
		Error("Out of range");
	}
	SubIndex.RowIndex = i;
	return SubIndex;
}
// +演算子のオーバーロード
Matrix Matrix::operator +(const Matrix& obj) const
{
	if(RowLen<=0 || ColLen<=0 || obj.RowLen<=0 || obj.ColLen<=0){
		Error("Null matrix");
	}
	if(obj.RowLen != RowLen || obj.ColLen != ColLen){
		Error("Size doesn't much");
	}
	Matrix m;
	m.New(RowLen, ColLen);
	const int g = (RowLen+1)*(ColLen+1);
	for(int i=0; i<g; i++){
		m.Array[i] = Array[i] + obj.Array[i];
	}
	return m;
}
// +演算子のオーバーロード（単項演算子）
Matrix Matrix::operator +() const
{
	if(RowLen<=0 || ColLen<=0){
		Error("Null matrix");
	}
	return *this;
}
// -演算子のオーバーロード
Matrix Matrix::operator -(const Matrix& obj) const
{
	if(RowLen<=0 || ColLen<=0 || obj.RowLen<=0 || obj.ColLen<=0){
		Error("Null matrix");
	}
	if(obj.RowLen != RowLen || obj.ColLen != ColLen){
		Error("Size doesn't much");
	}
	Matrix m;
	m.New(RowLen, ColLen);
	const int g = (RowLen+1)*(ColLen+1);
	for(int i=0; i<g; i++){
		m.Array[i] = Array[i] - obj.Array[i];
	}
	return m;
}
// -演算子のオーバーロード（単項演算子）
Matrix Matrix::operator -() const
{
	if(RowLen<=0 || ColLen<=0){
		Error("Null matrix");
	}
	Matrix m;
	m.New(RowLen, ColLen);
	const int g = (RowLen+1)*(ColLen+1);
	for(int i=0; i<g; i++){
		m.Array[i] = -Array[i];
	}
	return m;
}
// +=演算子のオーバーロード
Matrix Matrix::operator +=(const Matrix& obj)
{
	if(RowLen<=0 || ColLen<=0 || obj.RowLen<=0 || obj.ColLen<=0){
		Error("Null matrix");
	}
	if(obj.RowLen != RowLen || obj.ColLen != ColLen){
		Error("Size doesn't much");
	}
	const int g = (RowLen+1)*(ColLen+1);
	for(int i=0; i<g; i++){
		Array[i] += obj.Array[i];
	}
	return *this;
}
// -=演算子のオーバーロード
Matrix Matrix::operator -=(const Matrix& obj)
{
	if(RowLen<=0 || ColLen<=0 || obj.RowLen<=0 || obj.ColLen<=0){
		Error("Null matrix");
	}
	if(obj.RowLen != RowLen || obj.ColLen != ColLen){
		Error("Size doesn't much");
	}
	const int g = (RowLen+1)*(ColLen+1);
	for(int i=0; i<g; i++){
		Array[i] -= obj.Array[i];
	}
	return *this;
}
// *演算子のオーバーロード（スカラー倍）
Matrix Matrix::operator *(double a) const
{
	if(RowLen<=0 || ColLen<=0){
		Error("Null matrix");
	}
	Matrix m;
	m.New(RowLen, ColLen);
	const int g = (RowLen+1)*(ColLen+1);
	for(int i=0; i<g; i++){
		m.Array[i] = Array[i] * a;
	}
	return m;
}
// *=演算子のオーバーロード（スカラー倍）
Matrix Matrix::operator *=(double a)
{
	if(RowLen<=0 || ColLen<=0){
		Error("Null matrix");
	}
	const int g = (RowLen+1)*(ColLen+1);
	for(int i=0; i<g; i++){
		Array[i] *= a;
	}
	return *this;
}
// *演算子のオーバーロード（左からのスカラー倍）
Matrix operator *(double a, const Matrix& obj)
{
	if(obj.RowLen<=0 || obj.ColLen<=0){
		Error("Null matrix");
	}
	Matrix m;
	m.New(obj.RowLen, obj.ColLen);
	const int g = (obj.RowLen+1)*(obj.ColLen+1);
	for(int i=0; i<g; i++){
		m.Array[i] = a * obj.Array[i];
	}
	return m;
}
// *演算子のオーバーロード（行列との積）
Matrix Matrix::operator *(const Matrix& obj) const
{
	if(RowLen<=0 || ColLen<=0 || obj.RowLen<=0 || obj.ColLen<=0){
		Error("Null matrix");
	}
	if(ColLen!=obj.RowLen){
		Error("Size doesn't much");
	}
	Matrix A;
	A.New(RowLen, obj.ColLen);
	const int c1 = ColLen + 1;
	const int c2 = obj.ColLen + 1;
	for(int i=0; i<=RowLen; i++){
		for(int j=0; j<=obj.ColLen; j++){
			double a = 0;
			if(i>0 && j>0){
				for(int k=1; k<=ColLen; k++){
					a += Array[i*c1+k] * obj.Array[k*c2+j];
				}
			}
			A.Array[i*c2+j] = a;
		}
	}
	return A;
}
// *演算子のオーバーロード（ベクトルとの積）
Vector Matrix::operator *(const Vector& obj) const
{
	if(obj.Len<=0){
		Error("Null vector");
	}
	if(RowLen<=0 || ColLen<=0){
		Error("Null matrix");
	}
	if(ColLen!=obj.Len){
		Error("Size doesn't much");
	}
	Vector v;
	v.New(RowLen);
	v.Array[0] = 0;
	const int c = ColLen + 1;
	for(int i=1; i<=RowLen; i++){
		double a = 0;
		for(int j=1; j<=ColLen; j++){
			a += Array[i*c+j] * obj.Array[j];
		}
		v.Array[i] = a;
	}
	return v;
}
// /演算子のオーバーロード（スカラー除算）
Matrix Matrix::operator /(double a) const
{
	if(RowLen<=0 || ColLen<=0){
		Error("Null matrix");
	}
	if(a==0){
		Error("Divide by zero");
	}
	Matrix m;
	m.New(RowLen, ColLen);
	const int g = (RowLen+1)*(ColLen+1);
	for(int i=0; i<g; i++){
		m.Array[i] = Array[i] / a;
	}
	return m;
}
// /=演算子のオーバーロード（スカラー除算）
Matrix Matrix::operator /=(double a)
{
	if(RowLen<=0 || ColLen<=0){
		Error("Null matrix");
	}
	if(a==0){
		Error("Divide by zero");
	}
	const int g = (RowLen+1)*(ColLen+1);
	for(int i=0; i<g; i++){
		Array[i] /= a;
	}
	return *this;
}
// 転置行列
Matrix Matrix::Transpose() const
{
	if(RowLen<=0 || ColLen<=0){
		Error("Null matrix");
	}
	Matrix A;
	A.New(ColLen, RowLen);
	const int c1 = ColLen + 1;
	const int c2 = RowLen + 1;
	for(int i=0; i<=RowLen; i++){
		for(int j=0; j<=ColLen; j++){
			A.Array[j*c2+i] = Array[i*c1+j];
		}
	}
	return A;
}
// 横ベクトルの取り出し
Vector Matrix::RowVector(int k) const
{
	if(RowLen<=0 || ColLen<=0){
		Error("Null matrix");
	}
	if(k<0 || k>RowLen){
		Error("Out of range");
	}
	Vector v;
	v.New(ColLen);
	const int c = ColLen+1;
	for(int i=0; i<=ColLen; i++){
		v.Array[i] = Array[k*c+i];
	}
	return v;
}
// 縦ベクトルの取り出し
Vector Matrix::ColVector(int k) const
{
	if(RowLen<=0 || ColLen<=0){
		Error("Null matrix");
	}
	if(k<0 || k>ColLen){
		Error("Out of range");
	}
	Vector v;
	v.New(RowLen);
	const int c = ColLen+1;
	for(int i=0; i<=RowLen; i++){
		v.Array[i] = Array[i*c+k];
	}
	return v;
}
// 横ベクトルのセット
void Matrix::SetRowVector(int k, const Vector& v)
{
	if(v.Len<=0){
		Error("Null vector");
	}
	if(RowLen<=0 || ColLen<=0){
		Error("Null matrix");
	}
	if(k<0 || k>RowLen){
		Error("Out of range");
	}
	if(v.Len!=ColLen){
		Error("Size doesn't mutch");
	}
	const int c = ColLen+1;
	for(int i=0; i<=ColLen; i++){
		Array[k*c+i] = v.Array[i];
	}
}
// 縦ベクトルのセット
void Matrix::SetColVector(int k, const Vector& v)
{
	if(v.Len<=0){
		Error("Null vector");
	}
	if(RowLen<=0 || ColLen<=0){
		Error("Null matrix");
	}
	if(k<0 || k>ColLen){
		Error("Out of range");
	}
	if(v.Len!=RowLen){
		Error("Size doesn't mutch");
	}
	const int c = ColLen+1;
	for(int i=0; i<=RowLen; i++){
		Array[i*c+k] = v.Array[i];
	}
}
// ガウスの消去法による求解
Vector Matrix::Gauss(const Vector& b)
{
	if(b.Len<=0){
		Error("Null vector");
	}
	if(RowLen<=0 || ColLen<=0){
		Error("Null matrix");
	}
	if(RowLen != ColLen){
		Error("Not square matrix");
	}
	if(RowLen != b.Len){
		Error("Size doesn't match");
	}

	int g = (RowLen+1)*(ColLen+1);
	WM.Keep(g);
	for(int i=0; i<g; i++){
		WM.Array[i] = Array[i];
	}
	WM.v = b;
	double *A = WM.Array;
	double *x = WM.v.Array;
	
	const int c = ColLen+1;
	// リスケーリング
	for(int i=1; i<=RowLen; i++){
		double ma = 0;
		for(int j=1; j<=ColLen; j++){
			double a = fabs(A[i*c+j]);
			if(ma<a){ ma = a;}
		}
		for(int j=1; j<=ColLen; j++){
			A[i*c+j] /= ma;
		}
		x[i] /= ma;
	}
	// 前進消去
	for(int k=1; k<=RowLen-1; k++){
		//ピボット選択
		int p = k;
		double pv = fabs(A[k*c+k]);
		for(int i=k+1; i<=RowLen; i++){
			double a = fabs(A[i*c+k]);
			if(a>pv){p = i; pv = a;}
		}
		if(p>k){
			for(int i=k; i<=ColLen; i++){
				double a = A[k*c+i]; A[k*c+i] = A[p*c+i]; A[p*c+i] = a;
			}
			double a = x[k]; x[k] = x[p]; x[p] = a;
		}
		for(int i=k+1; i<=RowLen; i++){
			if(A[k*c+k]==0){
				Error("Matrix is singular");
			}
			double a = A[i*c+k] / A[k*c+k];
			for(int j=k+1; j<=ColLen; j++){
				A[i*c+j] -= a*A[k*c+j];
			}
			x[i] -= a*x[k];
		}
	}
	//後退代入
	for(int k=RowLen; k>=1; k--){
		for(int i=k+1; i<=ColLen; i++){
			x[k] -= A[k*c+i] * x[i];
		}
		x[k] /= A[k*c+k];
	}
	return WM.v;
}
// 行サイズ
int Matrix::RowSize() const
{
	return RowLen;
}
// 列サイズ
int Matrix::ColSize() const
{
	return ColLen;
}
// 表示
void Matrix::Show(int n) const
{
	if(RowLen>0 && ColLen>0 && n>0){
		const int c = ColLen+1;
		for(int i=1; i<=RowLen; i++){
			for(int j=1; j<=ColLen; j++){
				std::cout << std::left << std::setw(n+8);
				if(Array[i*c+j]==0){
					std::cout << std::fixed << std::setprecision(n-1) << Array[i*c+j];
					continue;
				}
				int k;
				double b = 1;
				double f = fabs(Array[i*c+j])*100000;
				for(k=0; k<=n+5; k++){
					if(b>f){
						break;
					}
					b *= 10;
				}
				if(k==0 || k>n+5){
					std::cout << std::scientific << std::setprecision(n-1) << Array[i*c+j];
				} else {
					std::cout << std::fixed << std::setprecision(n+5-k) << Array[i*c+j];
				}
			}
			std::cout << std::endl;
		}
	}
}

#endif