矩阵计算库

/*
---- From XDU's mzb
*/
#include <bits/stdc++.h>
using namespace std;
using ll = long long int;
namespace matrix
  {
  	template <typename T>
  	struct matrix : vector<vector<T>>
  	  {
  	  	using value_type = T;
  	  	matrix(ll n = 0)
  	  	  :matrix(n,n)
  	  	    {}
  	  	matrix(ll n,ll m)
  	  	  :vector<vector<value_type>>(n,vector<value_type>(m))
  	  	    {}
  	  	ll row() const
  	  	  {
  	  	  	return this -> size();
			  }
		ll column() const
		  {
		  	return (this -> begin()) -> size();
		  }
  	  	void check(ll a,ll b) const
  	  	  {
  	  	  	if (a < 0 or b < 0 or a >= row() or b >= column())
  	  	  	  throw std::out_of_range("matirx::matrix::operator(ll,ll)::out_of_range");
			  }
		value_type const& operator () (ll a,ll b) const
		  {
		  	check(a,b);
		  	return (*this)[a][b];
		  }
		matrix& one()
		  {
		  	for (ll i = 0;i < row() and i < column();i++)
		  	  (*this)(i,i) = 1;
		  	return *this;
		   } 
		value_type& operator () (ll a,ll b)
		  {
		  	check(a,b);
		  	return (*this)[a][b];
		  }
		};
	template<typename T>
	ostream& operator << (ostream& os,matrix<T> const& rhs)
	  {
	  	for (ll i = 0;i < rhs.row();i++)
	  	  {
	  	  	for (ll k = 0;k < rhs.column();k++)
	  	  	  os << setw(5) << (fabs(rhs(i,k)) <= 1e-8 ? 0 : rhs(i,k)) << " ";
	  	  	os << "\n";
			}
		return os;
	  }
	template<typename T>
	matrix<T> operator * (matrix<T> const& a,matrix<T> const& b)
	  {
	  	if (a.column() != b.row())
	  	  throw std::out_of_range("matrix::operator * ()::wrong");
	  	matrix<T> ret(a.row(),b.column());
	  	for (ll i = 0;i < a.row();i++)
	  	  for (ll j = 0;j < b.column();j++)
	  	    for (ll k = 0;k < a.column();k++)
	  	      {
	  	      	ret(i,j) = ret(i,j) + a(i,k) * b(k,j);
				}
	  	return ret;
	  }
	template<typename T>
	vector<ll> lup_decomposition(matrix<T> & a)
	  {
	  	vector<ll> p;
	  	for (ll i = 0;i < a.row();i++)
	  	  {
	  	  	p.push_back(i);
			}
		for (ll i = 0;i < a.row();i++)
		  {
		  	auto pos = i;
		  	for (ll j = i + 1;j < a.row();j++)
		  	  {
		  	  	if (abs(a(j,i)) > abs(a(pos,i)))
		  	  	  {
		  	  	  	pos = j;
					  }
				}
			if (a(pos,i) == 0)
			  {
			  	return vector<ll>();
			  }
			swap(p[i],p[pos]);  
			for (ll j = 0;j < a.row();j++)
			  {
			  	swap(a(i,j),a(pos,j));
			  }
			for (ll j = i + 1;j < a.row();j++)
			  {
			  	a(j,i) /= a(i,i);
			  	for (ll k = i + 1;k < a.row();k++)
			  	  {
			  	  	a(j,k) -= a(j,i) * a(i,k);
					}
			  }
			}
		return p;
	  }
	template<typename T>
	vector<T> LUP_gen(matrix<T> a,vector<T> raw_y)
	  {
		auto p = lup_decomposition(a);
	  	auto y = raw_y;
	  	for (ll i = 0;i < ll(p.size());i++)
	  	  {
	  	  	y[i] = raw_y[p[i]];
			}
		/*cout << "after P" << "\n"; 
		for (auto it : y)
		  cout << setw(5) << it << " ";cout << "\n";*/
		for (ll i = 0;i < ll(y.size());i++)
		  {
			for (ll j = 0;j < i;j++)
			  {
			  	y[i] -= a(i,j) * y[j];
			   } 
		  }	
	/*	cout << "after L" << "\n";
		for (auto it : y)
		  cout << setw(5) << it << " ";cout << "\n";*/
		for (ll i = ll(y.size()) - 1;i >= 0;i--)
		  {
			for (ll j = i + 1;j < ll(y.size());j++)
			  {
			  	y[i] -= a[i][j] * y[j];
			  } 
			y[i] /= a[i][i];
		  }
		/*cout << "after U" << "\n";
		for (auto it : y)
		  cout << setw(5) << it << " ";cout << "\n";*/
	  	return y;
	  }
	template<typename T>
	void show_lup(matrix<T> a)
	  {
	  	auto p = lup_decomposition(a);
	  	if (p.size())
	  	  {
			cout << "L = \n";
			for (ll i = 0;i < a.row();i++)
	  	  	  {
	  	  	  	for (ll k = 0;k < a.row();k++)
	  	  	  	  {
	  	  	  	  	cout << setw(5) << (i >= k ? (i == k ? 1 : a(i,k)) : 0);
						}
				cout << "\n";
				  }
			cout << "U = \n";
			for (ll i = 0;i < a.row();i++)
	  	  	  {
	  	  	  	for (ll k = 0;k < a.row();k++)
	  	  	  	  {
	  	  	  	  	cout << setw(5) << (k >= i ? a(i,k) : 0);
						}
				cout << "\n";
				  }
			cout << "P = \n";
			for (auto it : p)
			  cout << setw(5) << it << " ";
			cout << "\n";
			
			}
		else
		  {
		  	cout << "det = 0" << "\n";
		  	return;
		  }
	  }
	template<typename T>
	matrix<T> inverse(matrix<T> const& rhs) // 转置 
	  {
	  	matrix<T> ret(rhs.row());
	  	for (ll i = 0;i < rhs.row();i++)
	  	  {
	  	  	vector<double> b(rhs.row());
	  	  	b[i] = 1.0;
	  	  	auto v = LUP_gen(rhs,b);
	  	  	for (ll k = 0;k < ll(v.size());k++)
	  	  	  {
	  	  	  	ret(k,i) = v[k];
				  }
			}
	  	return ret;
	  }
	template<typename T>
	matrix<T> transpose(matrix<T> const& rhs) // 转置 
	  {
	  	matrix<T> ret(rhs.column(),rhs.row());
	  	for (ll i = 0;i < rhs.row();i++)
	  	  for (ll k = 0;k < rhs.column();k++)
	  	    {
	  	    	ret(k,i) = rhs(i,k);
			  }
	  	return ret;
	  }
  };
int main()
  {
  	matrix::matrix<double> a(5,3),y(5,1);
  	a[0] = {1,-1,1};
  	a[1] = {1,1,1};
  	a[2] = {1,2,4};
  	a[3] = {1,3,9};
  	a[4] = {1,5,25};
  	
  	y[0][0] = 2;
  	y[1][0] = 1;
  	y[2][0] = 1;
  	y[3][0] = 0;
  	y[4][0] = 3;
  	
  	auto at = matrix::transpose(a);
  	auto fin = matrix::inverse(at * a) * at;
  	
  	cout << fin * y;
	return 0;
  }