#include "Mat_ZZ_p.h"
#include<fstream>
#include<cassert>
Mat_ZZ_p::Mat_ZZ_p(int row, int line) 
{
	if(row <= 0||line <= 0)
	{
		std:: cerr<<"bad Matrix size."<<std::endl;
		exit(1);
	}
	Row = row;
	Col = line;		
	m = new unsigned char[row * line];
	memset(m,0,row * line *sizeof(unsigned char));
}

Mat_ZZ_p::~Mat_ZZ_p(void)
{
	Row = 0;
	Col = 0;
	if(this->m != nullptr)
		delete []this->m;
}

void Mat_ZZ_p::resize(int i, int j)
{
	if(i <= 0||j <= 0)
	{
		std:: cerr<<"bad Matrix resize."<<std::endl;
		exit(1);
	}
	Row = i;
	Col = j;		
	m = new unsigned char[i * j];
	memset(m, 0, i * j *sizeof(unsigned char));
}
int Mat_ZZ_p::getRow()const
{
	return Row;
}
int Mat_ZZ_p::getCol()const
{
	return Col;
}

Mat_ZZ_p::Mat_ZZ_p(const Mat_ZZ_p& Mat)
{
	this->Row = Mat.getRow();
	this->Col = Mat.getCol(); 
	int size = Mat.getRow() * Mat.getCol();
    if(size != 0)
    {
        this->m = new unsigned char[size];

        for(int i = 0;i < size;++i)
            this->m[i] = Mat.m[i];
    }
    else
        m = NULL;
}
Mat_ZZ_p& Mat_ZZ_p::operator=(const Mat_ZZ_p& Mat)
{
	 if(this==&Mat) 
        return *this;
	 this->resize(Mat.getRow(),Mat.getCol());//很重要！

    for(int i = 0;i < (Row * Col);i++)
        m[i] = Mat.m[i];  
    return *this;
}

void Mat_ZZ_p::put(int i,int j,unsigned char a)
{
	if(i < 0 || Row <= i|| j < 0 || Col <= j)
	{
		std::cerr<<"Matrix index out of range."<<std::endl;
		exit(1);
	}
	m[i * Col + j] = a % 3;
}
unsigned char& Mat_ZZ_p::elem(int i,int j)//引用返回的目的是返回值可以做左值
{
	if(i < 0 || Row <= i|| j < 0 || Col <= j)
	{
		std::cerr<<"Matrix index out of range."<<std::endl;
		exit(1);
	}
	return *(m + i * Col + j);
}
unsigned char Mat_ZZ_p::elem(int i,int j)const
{
	if(i < 0 || Row <= i|| j < 0 || Col <= j)
	{
		std::cerr<<"Matrix index out of range."<<std::endl;
		exit(1);
	}
	return *(m + i * Col + j);
}
Mat_ZZ_p& Mat_ZZ_p::operator+=(const Mat_ZZ_p &A)
{
	if(!A.m) return *this;
	if(A.getRow() != this->Row || A.getCol() != this->Col)
	{
		std::cerr << "error matrix add equal due to dim different！"<< std::endl;
		std::abort();
	}
	int size = A.getRow() * A.getCol();
	for (int i = 0;i < size;i++)
		this->m[i] = (this->m[i] + A.m[i]) % THREE;
	return *this;
}
Mat_ZZ_p& Mat_ZZ_p::operator-=(const Mat_ZZ_p &A)
{
	if(!A.m) return *this;
	if(A.getRow() != Row || A.getCol() != Col)
	{
		std::cerr << "error matrix sub equal due to dim different！"<< std::endl;
		std::abort();
	}
	int size = A.getRow() * A.getCol();
	for (int i = 0;i < size;i++)
		this->m[i] = (THREE + this->m[i] - A.m[i]) % THREE;
	return *this;
}
Mat_ZZ_p operator+(const Mat_ZZ_p &A,const Mat_ZZ_p &B)
{
	Mat_ZZ_p tmp = A;
    tmp += B;
    return tmp;
}
Mat_ZZ_p operator-(const Mat_ZZ_p &A,const Mat_ZZ_p &B)
{
	Mat_ZZ_p tmp = A;
    tmp -= B;
    return tmp;
	//if(A.getRow() != B.getRow() || A.getCol() != B.getCol())
	//{
	//	std::cerr << "error matrix sub due to dim different！！"<< std::endl;
	//	std::abort();
	//}
	//int size = A.getRow() * A.getCol();
	//Mat_ZZ_p tmp(A.getRow(), A.getCol());
	//for (int i = 0;i < size;i++)
	//	tmp.m[i] = (3 + A.m[i] - B.m[i]) % 3;
	//return tmp;
}
Mat_ZZ_p operator*(int a,const Mat_ZZ_p& Mat)
{
	Mat_ZZ_p Res(Mat.getRow(),Mat.getCol());
	if(a == 0)
		return Res;
	else if(a == 1)
		return Mat;
	int size = Mat.getRow() * Mat.getCol();;
	for(int i = 0;i < size;i++)
		Res.m[i] = (a * Mat.m[i]) % 3;
	return Res;
}
Mat_ZZ_p operator*(const Mat_ZZ_p& Mat,int a)
{
	Mat_ZZ_p Res(Mat.getRow(),Mat.getCol());
	if(a == 0)
		return Res;
	else if(a == 1)
		return Mat;
	int size = Mat.getRow() * Mat.getCol();;
	for(int i = 0;i < size;i++)
		Res.m[i] = (a * Mat.m[i]) % 3;
	return Res;
}
Vec_ZZ_p operator*(const Mat_ZZ_p& Mat,const Vec_ZZ_p& Vec)
{
	if(Mat.getCol() != Vec.getlength())
	{
		std::cerr<<"bad multiply Matrix with Vector."<<std::endl;
		exit(1);
	}
	Vec_ZZ_p Res(Mat.getRow());
	for(int i = 0;i < Mat.getRow();i++)
	{
		Res[i] = 0;
		for(int j = 0;j < Mat.getCol();j++)
			Res[i] = ( Res[i] + Mat.m[i * Mat.getCol() + j] * Vec[j]) % THREE;
	}
	return Res;
}
void Mat_ZZ_p::clear()
{
	int size = Row * Col;
	memset(m,0,(size * sizeof(unsigned char)));
}

std::ostream & operator<<(std::ostream &os,const Mat_ZZ_p& Mat)
{
	int i,j;
	os << "[";
	for(i = 0;i < Mat.getRow();i++)
	{
		os << "[";
		for(j = 0;j < Mat.getCol() - 1;j++)
			os << (int)(Mat.elem(i, j)) << " ";
		os << (int)( Mat.elem(i, Mat.getCol() - 1));
		os << "]";
	}
	os << "]";
	return os;
}
bool Mat_ZZ_p::SaveToTxt(const char *filename)
{
	std::ofstream file(filename,std::ios::app);
    if(!file)
        return false;
    for(int i = 0;i < Row * Col;i++)
        file << (int)(*(m + i));
    file.close();
    return true;
}
bool Mat_ZZ_p::LoadFromTxt(const char *filename)
{
	std::ifstream file(filename, std::ios::in);
	assert(file);
	if(!file)
		return false;
	for(int i = 0;i < Row * Col;i++)
	{
		char ch;
		file >> ch;
		m[i] = ch - '0';
	}
	file.close();
    return true;
}

int det(const Mat_ZZ_p& mat)//矩阵随机化的时候都产生 0 1 2 的数
{
	if(mat.Row != mat.Col)
	{
		std::cerr<< "nonsquare matrix !"<<std::endl;
		abort();
	}
	if(mat.Row == 1)
		return mat.elem(0,0);
	else if(mat.Row == 2)
		return (6 + mat.elem(0,0) * mat.elem(1,1) - mat.elem(0,1)* mat.elem(1,0)) % 3;
	else
	{
		Mat_ZZ_p mat1 = mat;
		int n = mat.Col;
		int i;
		int det_calc = 1;
		for(int k = 0;k < n;k++)
		{
			//保证第k行第k个元素不为0
			for(i = k;i < n && mat1.elem(i,k) == 0;i++);
			if(i == n)//第k列如果全0,则行列式为0
				return 0;
			if(i != k)
			{
				//j代表列数,交换第i,k行
				for(int j = k; j < n;j++)
					std::swap(mat1.elem(i,j), mat1.elem(k,j));
			}

			//使第k行第k个元素为1
			if(mat1.elem(k,k) == 2)
			{
				for(int s = k;s < n;s++ )
					mat1.elem(k,s) = (TWO * mat1.elem(k,s)) % THREE;
			}
			//下三角变为0
			for(i = k + 1;i < n;i++)
			{
				if(mat1.elem(i,k) == 0)continue;
				else if(mat1.elem(i,k) == 1)
				{
					for(int t = k;t < n;t++)
						mat1.elem(i,t) = (THREE + mat1.elem(i,t) - mat1.elem(k,t)) % THREE;
				}
				else if(mat1.elem(i,k) == 2)
				{
					for(int t = k;t < n;t++)
						mat1.elem(i,t) = (SIX + mat1.elem(i,t) - TWO * mat1.elem(k,t)) % THREE;
				}					
			}
		}
		for(int k = 0;k < n;k++)
			det_calc  = (det_calc * mat1.elem(k,k) ) % 3;
		return det_calc;
	}
}

Mat_ZZ_p inv(Mat_ZZ_p& A)
{
	if(A.Row != A.Col)
	{
		std:: cerr<<"bad Matrix inverse."<<std::endl;
		exit(1);
	}
	int rows = A.Row;
    Mat_ZZ_p C(rows,rows);
    // 构造一个与A行相同的单位阵B = |A : E|
    Mat_ZZ_p B(rows, rows * 2);
    for(int s = 0;s < rows;++s)
        for(int t = 0;t < rows;++t)
			B.put(s,t,A.elem(s,t));
	for(int s = 0;s < rows;++s)
		B.put(s, s + rows, ONE);
	int i;
	for(int k = 0;k < rows;k++)
	{
		//保证第k行第k个元素不为0
		for(i = k;i < rows && B.elem(i,k) == 0;i++);
		if(i == rows)//第k列如果全0,则行列式为0
		{
			std:: cerr<<"Matrix can not inverse."<<std::endl;
			exit(1);
		}
		if(i != k)
		{
			//j代表列数,交换第i,k行
			for(int j = k; j < 2 * rows;j++)
				std::swap(B.elem(i,j), B.elem(k,j));
		}

		//使第k行第k个元素为1
		if(B.elem(k,k) == 2)
		{
			for(int s = k;s < 2 * rows;s++ )
				B.elem(k,s) = (TWO * B.elem(k,s)) % THREE;
		}

		//上三角变为0
		for(i = 0;i < k;i++)
		{
			if(B.elem(i,k) == 0)continue;
			else if(B.elem(i,k) == ONE)
			{
				for(int t = k;t < 2 * rows;t++)
					B.elem(i,t) = (THREE + B.elem(i,t) - B.elem(k,t)) % THREE;
			}
			else if(B.elem(i,k) == TWO)
			{
				for(int t = k;t < 2 * rows;t++)
					B.elem(i,t) = (SIX + B.elem(i,t) - TWO * B.elem(k,t)) % THREE;
			}					
		}

		//下三角变为0
		for(i = k + 1;i < rows;i++)
		{
			if(B.elem(i,k) == 0)continue;
			else if(B.elem(i,k) == ONE)
			{
				for(int t = k;t < 2 * rows;t++)
					B.elem(i,t) = (THREE + B.elem(i,t) - B.elem(k,t)) % THREE;
			}
			else if(B.elem(i,k) == TWO)
			{
				for(int t = k;t < 2 * rows;t++)
					B.elem(i,t) = (SIX + B.elem(i,t) - TWO * B.elem(k,t)) % THREE;
			}					
		}
	}

	for(int s = 0;s < rows;++s)
        for(int t = 0;t < rows;++t)
			C.put(s, t, B.elem(s, t + rows));
	return C;
}
	