#ifndef LINEAR_2025_04_13
#define LINEAR_2025_04_13

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

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

class IntVector;
class Vector;
class Matrix;

// 整数ベクトルクラス
class IntVector {
private:
	class IntVectorPtr {
	public:
		int *Array;
		int Len;
		int RefCount;
	};
	IntVectorPtr *Ptr;
	void Free();
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;
	int Size() const;
	void Clear();
	IntVector Copy() const;
};

// ベクトルクラス
class Vector {
	friend Matrix;
private:
	class VectorPtr {
	public:
		double *Array;
		int Len;
		int RefCount;
	};
	VectorPtr *Ptr;
	void Free();
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) const;
	Vector Shift(int k) const;
	double Norm1() const;
	double Norm2() const;
	double NormInfty() const;
	int Size() const;
	void Show(int n) const;
	void Show() const;
	void Clear();
	Vector Copy() const;
};

// 行列クラス
class Matrix {
	friend Vector;
private:
	class SubIndexClass {
	public:
		int ColLen;
		double *Array;
		mutable int RowIndex;
		double& operator[](int i);
		const double& operator[](int i) const;
	};
	SubIndexClass SubIndex;
	class WorkingMemory {// ワーキングメモリ管理クラス
	public:
		double *Array;
		int size;
		void Keep(int n);
		WorkingMemory();
		~WorkingMemory();
	};
	static WorkingMemory WM;
	class MatrixPtr {
	public:
		double *Array;
		int RowLen;
		int ColLen;
		int RefCount;
	};
	MatrixPtr *Ptr;
	void Free();
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& b) const;
	int RowSize() const;
	int ColSize() const;
	void Show(int n) const;
	void Show() const;
	void Clear();
	Matrix Copy() const;
};

Matrix::WorkingMemory Matrix::WM;

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

// メモリ解放
void IntVector::Free()
{
	if(Ptr != NULL){
		Ptr->RefCount--;
		if(Ptr->RefCount <= 0){
			delete[] Ptr->Array;
			delete Ptr;
		}
	}
}
// デフォルトコンストラクタ
IntVector::IntVector()
{
	Ptr = NULL;
}
// コンストラクタ
IntVector::IntVector(int n)
{
	if(n<=0){
		Error("Illegal size");
	}
	Ptr = new IntVectorPtr;
	Ptr->RefCount = 1;
	Ptr->Len = n;
	Ptr->Array = new int[n+1];
}
// デストラクタ
IntVector::~IntVector()
{
	Free();
}
// コピーコンストラクタ
IntVector::IntVector(const IntVector& obj)
{
	Ptr = obj.Ptr;
	if(Ptr != NULL){
		Ptr->RefCount++;
	}
}
// 代入演算子のオーバーロード
IntVector IntVector::operator =(const IntVector& obj)
{
	if(Ptr != obj.Ptr){
		Free();
		Ptr = obj.Ptr;
		if(Ptr != NULL){
			Ptr->RefCount++;
		}
	}
	return *this;
}
// []演算子のオーバーロード
int& IntVector::operator[](int i)
{
	if(Ptr==NULL || Ptr->Len<=0){
		Error("Null vector");
	}
	if(i<0 || i>Ptr->Len){
		Error("Out of range");
	}
	return Ptr->Array[i];
}
// []演算子のオーバーロード（constバージョン）
const int& IntVector::operator[](int i) const
{
	if(Ptr==NULL || Ptr->Len<=0){
		Error("Null vector");
	}
	if(i<0 || i>Ptr->Len){
		Error("Out of range");
	}
	return Ptr->Array[i];
}
// 整数ベクトルのサイズ
int IntVector::Size() const
{
	if(Ptr==NULL || Ptr->Len<=0){
		return 0;
	}
	return Ptr->Len;
}
// ゼロクリア
void IntVector::Clear()
{
	if(Ptr==NULL || Ptr->Len<=0){
		return;
	}
	int n = Ptr->Len;
	for(int i=0; i<=n; i++){
		Ptr->Array[i] = 0;
	}
}
// 複製
IntVector IntVector::Copy() const
{
	if(Ptr == NULL){
		return *this;
	}
	int n = Ptr->Len;
	IntVector v(n);
	for(int i=0; i<=n; i++){
		v.Ptr->Array[i] = Ptr->Array[i];
	}
	return v;
}

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

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

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

// ワーキングメモリー処理
Matrix::WorkingMemory::WorkingMemory()
{
	size = 0;
	Array = NULL;
}
Matrix::WorkingMemory::~WorkingMemory()
{
	delete[] Array;
}
void Matrix::WorkingMemory::Keep(int n)
{
	if(n>size){
		if(size>0){
			delete[] Array;
		}
		Array = new double[n];
		size = n;
	}
}
// メモリ解放
void Matrix::Free()
{
	if(Ptr != NULL){
		Ptr->RefCount--;
		if(Ptr->RefCount <= 0){
			delete[] Ptr->Array;
			delete Ptr;
		}
	}
}
// デフォルトコンストラクタ
Matrix::Matrix()
{
	Ptr = NULL;
}
// コンストラクタ
Matrix::Matrix(int m, int n)
{
	if(m<=0 || n<=0){
		Error("Illegal size");
	}
	Ptr = new MatrixPtr;
	Ptr->RefCount = 1;
	Ptr->RowLen = m;
	Ptr->ColLen = n;
	int g = (m+1)*(n+1);
	Ptr->Array = new double[g];
	SubIndex.ColLen = n;
	SubIndex.Array = Ptr->Array;
}
// デストラクタ
Matrix::~Matrix()
{
	Free();
}
// コピーコンストラクタ
Matrix::Matrix(const Matrix& obj)
{
	Ptr = obj.Ptr;
	if(Ptr != NULL){
		Ptr->RefCount++;
		SubIndex.ColLen = Ptr->ColLen;
		SubIndex.Array = Ptr->Array;
	}
}
// 代入演算子のオーバーロード
Matrix Matrix::operator =(const Matrix& obj)
{
	if(Ptr != obj.Ptr){
		Free();
		Ptr = obj.Ptr;
		if(Ptr != NULL){
			Ptr->RefCount++;
			SubIndex.ColLen = Ptr->ColLen;
			SubIndex.Array = Ptr->Array;
		}
	}
	return *this;
}
// []演算子のオーバーロード（2つ目の添字）
double& Matrix::SubIndexClass::operator[](int i)
{
	if(i<0 || i>ColLen){
		Error("Out of range");
	}
	return Array[RowIndex*(ColLen+1)+i];
}
// []演算子のオーバーロード（2つ目の添字、constバージョン）
const double& Matrix::SubIndexClass::operator[](int i) const
{
	if(i<0 || i>ColLen){
		Error("Out of range");
	}
	return Array[RowIndex*(ColLen+1)+i];
}
// []演算子のオーバーロード
Matrix::SubIndexClass& Matrix::operator[](int i)
{
	if(Ptr==NULL || Ptr->RowLen<=0 || Ptr->ColLen<=0){
		Error("Null matrix");
	}
	if(i<0 || i>Ptr->RowLen){
		Error("Out of range");
	}
	SubIndex.RowIndex = i;
	return SubIndex;
}
// []演算子のオーバーロード（constバージョン）
const Matrix::SubIndexClass& Matrix::operator[](int i) const
{
	if(Ptr==NULL || Ptr->RowLen<=0 || Ptr->ColLen<=0){
		Error("Null matrix");
	}
	if(i<0 || i>Ptr->RowLen){
		Error("Out of range");
	}
	SubIndex.RowIndex = i;
	return SubIndex;
}
// +演算子のオーバーロード
Matrix Matrix::operator +(const Matrix& obj) const
{
	if(Ptr==NULL || obj.Ptr==NULL || Ptr->RowLen<=0 || Ptr->ColLen<=0 || obj.Ptr->RowLen<=0 || obj.Ptr->ColLen<=0){
		Error("Null matrix");
	}
	int row = Ptr->RowLen;
	int col = Ptr->ColLen;
	if(obj.Ptr->RowLen != row || obj.Ptr->ColLen != col){
		Error("Size doesn't much");
	}
	Matrix A(row, col);
	int g = (row+1)*(col+1);
	for(int i=0; i<g; i++){
		A.Ptr->Array[i] = Ptr->Array[i] + obj.Ptr->Array[i];
	}
	return A;
}
// +演算子のオーバーロード（単項演算子）
Matrix Matrix::operator +() const
{
	if(Ptr==NULL || Ptr->RowLen<=0 || Ptr->ColLen<=0){
		Error("Null matrix");
	}
	return *this;
}
// -演算子のオーバーロード
Matrix Matrix::operator -(const Matrix& obj) const
{
	if(Ptr==NULL || obj.Ptr==NULL || Ptr->RowLen<=0 || Ptr->ColLen<=0 || obj.Ptr->RowLen<=0 || obj.Ptr->ColLen<=0){
		Error("Null matrix");
	}
	int row = Ptr->RowLen;
	int col = Ptr->ColLen;
	if(obj.Ptr->RowLen != row || obj.Ptr->ColLen != col){
		Error("Size doesn't much");
	}
	Matrix A(row, col);
	int g = (row+1)*(col+1);
	for(int i=0; i<g; i++){
		A.Ptr->Array[i] = Ptr->Array[i] - obj.Ptr->Array[i];
	}
	return A;
}
// -演算子のオーバーロード（単項演算子）
Matrix Matrix::operator -() const
{
	if(Ptr==NULL || Ptr->RowLen<=0 || Ptr->ColLen<=0){
		Error("Null matrix");
	}
	int row = Ptr->RowLen;
	int col = Ptr->ColLen;
	Matrix A(row, col);
	int g = (row+1)*(col+1);
	for(int i=0; i<g; i++){
		A.Ptr->Array[i] = -Ptr->Array[i];
	}
	return A;
}
// +=演算子のオーバーロード
Matrix Matrix::operator +=(const Matrix& obj)
{
	if(Ptr==NULL || obj.Ptr==NULL || Ptr->RowLen<=0 || Ptr->ColLen<=0 || obj.Ptr->RowLen<=0 || obj.Ptr->ColLen<=0){
		Error("Null matrix");
	}
	int row = Ptr->RowLen;
	int col = Ptr->ColLen;
	if(obj.Ptr->RowLen != row || obj.Ptr->ColLen != col){
		Error("Size doesn't much");
	}
	int g = (row+1)*(col+1);
	for(int i=0; i<g; i++){
		Ptr->Array[i] += obj.Ptr->Array[i];
	}
	return *this;
}
// -=演算子のオーバーロード
Matrix Matrix::operator -=(const Matrix& obj)
{
	if(Ptr==NULL || obj.Ptr==NULL || Ptr->RowLen<=0 || Ptr->ColLen<=0 || obj.Ptr->RowLen<=0 || obj.Ptr->ColLen<=0){
		Error("Null matrix");
	}
	int row = Ptr->RowLen;
	int col = Ptr->ColLen;
	if(obj.Ptr->RowLen != row || obj.Ptr->ColLen != col){
		Error("Size doesn't much");
	}
	int g = (row+1)*(col+1);
	for(int i=0; i<g; i++){
		Ptr->Array[i] -= obj.Ptr->Array[i];
	}
	return *this;
}
// *演算子のオーバーロード（スカラー倍）
Matrix Matrix::operator *(double a) const
{
	if(Ptr==NULL || Ptr->RowLen<=0 || Ptr->ColLen<=0){
		Error("Null matrix");
	}
	int row = Ptr->RowLen;
	int col = Ptr->ColLen;
	Matrix A(row, col);
	int g = (row+1)*(col+1);
	for(int i=0; i<g; i++){
		A.Ptr->Array[i] = Ptr->Array[i] * a;
	}
	return A;
}
// *=演算子のオーバーロード（スカラー倍）
Matrix Matrix::operator *=(double a)
{
	if(Ptr==NULL || Ptr->RowLen<=0 || Ptr->ColLen<=0){
		Error("Null matrix");
	}
	int row = Ptr->RowLen;
	int col = Ptr->ColLen;
	int g = (row+1)*(col+1);
	for(int i=0; i<g; i++){
		Ptr->Array[i] *= a;
	}
	return *this;
}
// *演算子のオーバーロード（左からのスカラー倍）
Matrix operator *(double a, const Matrix& obj)
{
	if(obj.Ptr==NULL || obj.Ptr->RowLen<=0 || obj.Ptr->ColLen<=0){
		Error("Null matrix");
	}
	int row = obj.Ptr->RowLen;
	int col = obj.Ptr->ColLen;
	Matrix A(row, col);
	int g = (row+1)*(col+1);
	for(int i=0; i<g; i++){
		A.Ptr->Array[i] = a * obj.Ptr->Array[i];
	}
	return A;
}
// *演算子のオーバーロード（行列との積）
Matrix Matrix::operator *(const Matrix& obj) const
{
	if(Ptr==NULL || obj.Ptr==NULL || Ptr->RowLen<=0 || Ptr->ColLen<=0 || obj.Ptr->RowLen<=0 || obj.Ptr->ColLen<=0){
		Error("Null matrix");
	}
	int row1 = Ptr->RowLen;
	int col1 = Ptr->ColLen;
	int col2 = obj.Ptr->ColLen;
	if(col1 != obj.Ptr->RowLen){
		Error("Size doesn't much");
	}
	Matrix A(row1, col2);
	int c1 = col1 + 1;
	int c2 = col2 + 1;
	for(int i=0; i<=row1; i++){
		for(int j=0; j<=col2; j++){
			double a = 0.0;
			if(i>0 && j>0){
				for(int k=1; k<=col1; k++){
					a += Ptr->Array[i*c1+k] * obj.Ptr->Array[k*c2+j];
				}
			}
			A.Ptr->Array[i*c2+j] = a;
		}
	}
	return A;
}
// *演算子のオーバーロード（ベクトルとの積）
Vector Matrix::operator *(const Vector& obj) const
{
	if(obj.Ptr==NULL || obj.Ptr->Len<=0){
		Error("Null vector");
	}
	if(Ptr==NULL || Ptr->RowLen<=0 || Ptr->ColLen<=0){
		Error("Null matrix");
	}
	int row = Ptr->RowLen;
	int col = Ptr->ColLen;
	if(col != obj.Ptr->Len){
		Error("Size doesn't much");
	}
	Vector v(row);
	v.Ptr->Array[0] = 0.0;
	int c = col + 1;
	for(int i=1; i<=row; i++){
		double a = 0.0;
		for(int j=1; j<=col; j++){
			a += Ptr->Array[i*c+j] * obj.Ptr->Array[j];
		}
		v.Ptr->Array[i] = a;
	}
	return v;
}
// /演算子のオーバーロード（スカラー除算）
Matrix Matrix::operator /(double a) const
{
	if(Ptr==NULL || Ptr->RowLen<=0 || Ptr->ColLen<=0){
		Error("Null matrix");
	}
	if(a==0.0){
		Error("Divide by zero");
	}
	int row = Ptr->RowLen;
	int col = Ptr->ColLen;
	Matrix A(row, col);
	int g = (row+1)*(col+1);
	for(int i=0; i<g; i++){
		A.Ptr->Array[i] = Ptr->Array[i] / a;
	}
	return A;
}
// /=演算子のオーバーロード（スカラー除算）
Matrix Matrix::operator /=(double a)
{
	if(Ptr==NULL || Ptr->RowLen<=0 || Ptr->ColLen<=0){
		Error("Null matrix");
	}
	if(a==0.0){
		Error("Divide by zero");
	}
	int row = Ptr->RowLen;
	int col = Ptr->ColLen;
	int g = (row+1)*(col+1);
	for(int i=0; i<g; i++){
		Ptr->Array[i] /= a;
	}
	return *this;
}
// 転置行列
Matrix Matrix::Transpose() const
{
	if(Ptr==NULL || Ptr->RowLen<=0 || Ptr->ColLen<=0){
		Error("Null matrix");
	}
	int row = Ptr->RowLen;
	int col = Ptr->ColLen;
	Matrix A(col, row);
	int c1 = col + 1;
	int c2 = row + 1;
	for(int i=0; i<=row; i++){
		for(int j=0; j<=col; j++){
			A.Ptr->Array[j*c2+i] = Ptr->Array[i*c1+j];
		}
	}
	return A;
}
// 横ベクトルの取り出し
Vector Matrix::RowVector(int k) const
{
	if(Ptr==NULL || Ptr->RowLen<=0 || Ptr->ColLen<=0){
		Error("Null matrix");
	}
	int row = Ptr->RowLen;
	int col = Ptr->ColLen;
	if(k<0 || k>row){
		Error("Out of range");
	}
	Vector v(col);
	int c = col + 1;
	for(int i=0; i<=col; i++){
		v.Ptr->Array[i] = Ptr->Array[k*c+i];
	}
	return v;
}
// 縦ベクトルの取り出し
Vector Matrix::ColVector(int k) const
{
	if(Ptr==NULL || Ptr->RowLen<=0 || Ptr->ColLen<=0){
		Error("Null matrix");
	}
	int row = Ptr->RowLen;
	int col = Ptr->ColLen;
	if(k<0 || k>col){
		Error("Out of range");
	}
	Vector v(row);
	int c = col + 1;
	for(int i=0; i<=row; i++){
		v.Ptr->Array[i] = Ptr->Array[i*c+k];
	}
	return v;
}
// 横ベクトルのセット
void Matrix::SetRowVector(int k, const Vector& v)
{
	if(v.Ptr==NULL || v.Ptr->Len<=0){
		Error("Null vector");
	}
	if(Ptr==NULL || Ptr->RowLen<=0 || Ptr->ColLen<=0){
		Error("Null matrix");
	}
	int row = Ptr->RowLen;
	int col = Ptr->ColLen;
	if(k<0 || k>row){
		Error("Out of range");
	}
	if(v.Ptr->Len != col){
		Error("Size doesn't mutch");
	}
	int c = col + 1;
	for(int i=0; i<=col; i++){
		Ptr->Array[k*c+i] = v.Ptr->Array[i];
	}
}
// 縦ベクトルのセット
void Matrix::SetColVector(int k, const Vector& v)
{
	if(v.Ptr==NULL || v.Ptr->Len<=0){
		Error("Null vector");
	}
	if(Ptr==NULL || Ptr->RowLen<=0 || Ptr->ColLen<=0){
		Error("Null matrix");
	}
	int row = Ptr->RowLen;
	int col = Ptr->ColLen;
	if(k<0 || k>col){
		Error("Out of range");
	}
	if(v.Ptr->Len != row){
		Error("Size doesn't mutch");
	}
	int c = col + 1;
	for(int i=0; i<=row; i++){
		Ptr->Array[i*c+k] = v.Ptr->Array[i];
	}
}
// ガウスの消去法による求解
Vector Matrix::Gauss(const Vector& b) const
{
	if(b.Ptr==NULL || b.Ptr->Len<=0){
		Error("Null vector");
	}
	if(Ptr==NULL || Ptr->RowLen<=0 || Ptr->ColLen<=0){
		Error("Null matrix");
	}
	int row = Ptr->RowLen;
	int col = Ptr->ColLen;
	if(row != col){
		Error("Not square matrix");
	}
	if(row != b.Ptr->Len){
		Error("Size doesn't match");
	}
	int g = (row+1)*(col+1);
	WM.Keep(g);
	double *A = WM.Array;
	for(int i=0; i<g; i++){
		A[i] = Ptr->Array[i];
	}
	Vector v = b.Copy();
	double *x = v.Ptr->Array;
	int c = col + 1;
	// リスケーリング
	for(int i=1; i<=row; i++){
		double ma = 0.0;
		for(int j=1; j<=col; j++){
			double a = fabs(A[i*c+j]);
			if(ma<a){ ma = a;}
		}
		for(int j=1; j<=col; j++){
			A[i*c+j] /= ma;
		}
		x[i] /= ma;
	}
	// 前進消去
	for(int k=1; k<=row-1; k++){
		//ピボット選択
		int p = k;
		double pv = fabs(A[k*c+k]);
		for(int i=k+1; i<=row; i++){
			double a = fabs(A[i*c+k]);
			if(a>pv){p = i; pv = a;}
		}
		if(p>k){
			for(int i=k; i<=col; 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<=row; i++){
			if(A[k*c+k]==0.0){
				Error("Matrix is singular");
			}
			double a = A[i*c+k] / A[k*c+k];
			for(int j=k+1; j<=col; j++){
				A[i*c+j] -= a*A[k*c+j];
			}
			x[i] -= a*x[k];
		}
	}
	//後退代入
	for(int k=row; k>=1; k--){
		for(int i=k+1; i<=col; i++){
			x[k] -= A[k*c+i] * x[i];
		}
		x[k] /= A[k*c+k];
	}
	return v;
}
// 行サイズ
int Matrix::RowSize() const
{
	if(Ptr==NULL || Ptr->RowLen<=0){
		return 0;
	}
	return Ptr->RowLen;
}
// 列サイズ
int Matrix::ColSize() const
{
	if(Ptr==NULL || Ptr->ColLen<=0){
		return 0;
	}
	return Ptr->ColLen;
}
// 表示
void Matrix::Show(int n) const
{
	if(Ptr==NULL || Ptr->RowLen<=0 || Ptr->ColLen<=0 || n<=0){
		return;
	}
	int row = Ptr->RowLen;
	int col = Ptr->ColLen;
	const int c = col + 1;
	for(int i=1; i<=row; i++){
		for(int j=1; j<=col; j++){
			std::cout << std::left << std::setw(n+8);
			if(Ptr->Array[i*c+j]==0.0){
				std::cout << std::fixed << std::setprecision(n-1) << Ptr->Array[i*c+j];
				continue;
			}
			int k;
			double b = 1.0;
			double f = fabs(Ptr->Array[i*c+j])*100000.0;
			for(k=0; k<=n+5; k++){
				if(b>f){
					break;
				}
				b *= 10.0;
			}
			if(k==0 || k>n+5){
				std::cout << std::scientific << std::setprecision(n-1) << Ptr->Array[i*c+j];
			} else {
				std::cout << std::fixed << std::setprecision(n+5-k) << Ptr->Array[i*c+j];
			}
		}
		std::cout << std::endl;
	}
}
// 表示（引数無し）
void Matrix::Show() const
{
	Show(16);
}
// ゼロクリア
void Matrix::Clear()
{
	if(Ptr==NULL || Ptr->RowLen<=0 || Ptr->ColLen<=0){
		return;
	}
	int g = (Ptr->RowLen+1)*(Ptr->ColLen+1);
	for(int i=0; i<g; i++){
		Ptr->Array[i] = 0.0;
	}
}
// 複製
Matrix Matrix::Copy() const
{
	if(Ptr == NULL){
		return *this;
	}
	int row = Ptr->RowLen;
	int col = Ptr->ColLen;
	Matrix A(row, col);
	int g = (row+1)*(col+1);
	for(int i=0; i<g; i++){
		A.Ptr->Array[i] = Ptr->Array[i];
	}
	return A;
}

#endif