// INPUT : Just the matrix. (No dimension) // For a sample, use: /* 1+i 2+2i i+3 4-2i 5+5i 6+6i 7+7i 8+8i 9+9i */ /* This program calculates the row echelon form of a matrix, and additionally the reduced row echelon form and the rank of the matrix. The elements can be integers, fractions or complex numbers (represented by a+bi form, where a and b can themselves be integers or fractions). Along with the results, you_get the complete set of elementary row operations you need to apply to get the reduced form yourself, and the intermediate results...*/ // It will be very kind of you if you can // leave your suggestions about code // practices, readability and/or efficiency // of the code. Doubts are welcome! // If your matrix has more than 2 columns, // switch to landscape mode for the output. // It may get messy in portrait mode. // Created by Kinshuk Vasisht, 2019 #include // atan2, sqrt #include // get, forward_as_tuple #include // numeric_limits #include // string #include // vector #include // stringstream #include // setw #include // ostream, istream #include // exception, etc. /// Namespace for storing classes and /// functions for Mathematical Purposes. namespace math { /// Returns the floor of a number - constexpr mode. constexpr double floor(double num) { return static_cast((static_cast(static_cast(num))==num)?static_cast(num):static_cast(num)-((num<0.0)?1:0)); } /// Returns the ceil of a number - constexpr mode. constexpr double ceil(double num) { return static_cast((static_cast(static_cast(num))==num)?static_cast(num):static_cast(num)+((num>0.0)?1:0)); } /// Returns the number rounded off to the nearest integer - constexpr mode. constexpr double round(double num) { if((math::ceil(num)-num)<(math::floor(num)-num)) return math::ceil(num); else return math::floor(num); } /// Class to represent fractions. class Fraction { long long num, den; /// Simplifies the fraction, removing any irregularities. constexpr void simplify() noexcept { if(den==0) { den=1; } if(num<0&&den<0) { num=-num; den=-den; } else if(den<0) { num=-num; den=-den; } long long gcd = hcf(labs(num),labs(den)); num/=gcd; den/=gcd; } /// Returns HCF of two numbers. constexpr long long hcf(long long n, long long d) const { return ((d==0)?n:hcf(d,n%d)); } /// Returns LCM of two numbers. constexpr long long lcm(long long n, long long d) const { long long gcd = hcf(n,d); return (n/gcd)*(d/gcd); } public: /// Default Constructor : Assigns fraction with a default value, equal to 0/1. constexpr Fraction() : num(0), den(1) {} /// Parametrized Constructor : Assigns fraction with numerator and denominator values. constexpr Fraction(long long n, long long d) : num(n), den(d) { this->simplify(); } /// Parametrized Constructor : Assigns fraction the fractional representation of a double. constexpr Fraction(long double dec) : num(0), den(1) { this->assign(dec); } /// Copy Constructor : Assigns fraction values from another Fraction constexpr Fraction(const Fraction& f) : num(f.num), den(f.den) {} /// Returns the fraction converted to a double. explicit constexpr operator long double() const noexcept { return num/static_cast(den); } /// Returns the fraction converted to a double. explicit constexpr operator double() const noexcept { return num/static_cast(den); } /// Returns the fraction converted to a double. explicit constexpr operator float() const noexcept { return num/static_cast(den); } /// Assigns the fraction from the value of a double. constexpr Fraction& assign(const long double& d) noexcept { double integral = math::floor(d), frac = d-integral; const long long precision = 10000000ULL; den = precision; num = 0; num = math::round(frac*precision); this->simplify(); num += (den*integral); return *this; } /// Returns the fraction, negated. constexpr Fraction operator-() const noexcept { return {-num,den}; } /// Returns the sum of two fractions. constexpr Fraction operator+(const Fraction& f2) const noexcept { return {(num*f2.den)+(f2.num*den),den*f2.den}; } /// Returns the difference of two fractions. constexpr Fraction operator-(const Fraction& f2) const noexcept { return {(num*f2.den)-(f2.num*den),den*f2.den}; } /// Returns the product of two fractions. constexpr Fraction operator*(const Fraction& f2) const noexcept { return {num*f2.num,den*f2.den}; } /// Returns the quotient of two fractions. constexpr Fraction operator/(const Fraction& f2) const noexcept { return {num*f2.den,den*f2.num}; } /// Returns the sum of a fraction and an integer. constexpr Fraction operator+(const long long& n2) const noexcept { return {num+(n2*den),den}; } /// Returns the difference of a fraction and an integer. constexpr Fraction operator-(const long long& n2) const noexcept { return {num-(n2*den),den}; } /// Returns the product of a fraction and an integer. constexpr Fraction operator*(const long long& n2) const noexcept { return {num*n2,den}; } /// Returns the quotient of a fraction and an integer. constexpr Fraction operator/(const long long& n2) const noexcept { return {num,den*n2}; } /// Returns the Fraction raised to the power of an integer. constexpr Fraction operator^(const long long& power) const noexcept { Fraction res{1,1}; if(power==0) return res; else { for(long long i=0;i (const Fraction& F1, const Fraction& F2) { return F1.num*F2.den>F1.den*F2.num; } /// Compares two Fraction Objects, checking if the Fraction object is greater than or equal to the other. friend constexpr bool operator>=(const Fraction& F1, const Fraction& F2) { return F1.num*F2.den>=F2.num*F1.den; } /// Compares a Fraction Object and an integer for equality. (integer on the right). friend constexpr bool operator==(const Fraction& F1, long long F2) { return F1.num==F2&&F1.den==1LL; } /// Compares a Fraction Object and an integer for inequality. (integer on the right). friend constexpr bool operator!=(const Fraction& F1, long long F2) { return F1.num!=F2||F1.den!=1LL; } /// Compares a Fraction Object and an integer, checking if the Object is smaller than the integer. (integer on the right). friend constexpr bool operator< (const Fraction& F1, long long F2) { return F1.den*F2>F1.num; } /// Compares a Fraction Object and an integer, checking if the Object is smaller than or equal to the integer. (integer on the right). friend constexpr bool operator<=(const Fraction& F1, long long F2) { return F1.den*F2>=F1.num; } /// Compares a Fraction Object and an integer, checking if the Object is greater than the integer. (integer on the right). friend constexpr bool operator> (const Fraction& F1, long long F2) { return F1.den*F2=(const Fraction& F1, long long F2) { return F1.den*F2<=F1.num; } /// Compares a Fraction Object and an integer for equality. (integer on the left). friend constexpr bool operator==(long long F1, const Fraction& F2) { return F1==F2.num&&F2.den==1LL; } /// Compares a Fraction Object and an integer for inequality. (integer on the left). friend constexpr bool operator!=(long long F1, const Fraction& F2) { return F1!=F2.num||F2.den!=1LL; } /// Compares a Fraction Object and an integer, checking if the integer is smaller than the Object. (integer on the left). friend constexpr bool operator< (long long F1, const Fraction& F2) { return F1*F2.den (long long F1, const Fraction& F2) { return F1*F2.den>F2.num; } /// Compares a Fraction Object and an integer, checking if the integer is greater than or equal to the Object. (integer on the left). friend constexpr bool operator>=(long long F1, const Fraction& F2) { return F1*F2.den>=F2.num; } /// Reads a fraction from an input stream. friend std::istream& operator>>(std::istream& is, Fraction& f) { char delim; is>>f.num; if(is.peek()=='/') { is>>delim>>f.den; } else { f.den=1; } f.simplify(); return is; } /// Writes a fraction to an output stream. friend std::ostream& operator<<(std::ostream& os, const Fraction& f) { std::stringstream ss; ss< : Checks if the Complex Object is greater than a scalar. constexpr bool operator> (const Fraction& val) const noexcept { return re(*this)>val; } /// operator >= : Checks if the Complex Object is greater than or equal to a scalar. constexpr bool operator>=(const Fraction& val) const noexcept { return (*this>val)||(*this==val); } /// operator == : Checks if the Complex Object and the scalar are equal. constexpr bool operator==(const long long& val) const noexcept { return val==real&&imag==0; } /// operator != : Checks if the Complex Object and the scalar are unequal. constexpr bool operator!=(const long long& val) const noexcept { return !(*this==val); } /// operator < : Checks if the Complex Object is less than a scalar. constexpr bool operator< (const long long& val) const noexcept { return re(*this) : Checks if the Complex Object is greater than a scalar. constexpr bool operator> (const long long& val) const noexcept { return re(*this)>val; } /// operator >= : Checks if the Complex Object is greater than or equal to a scalar. constexpr bool operator>=(const long long& val) const noexcept { return (*this>val)||(*this==val); } /// operator + : Returns sum of two Complex Objects. constexpr Complex operator+(const Complex& C2) const noexcept { return { real+C2.real, imag+C2.imag }; } /// operator - : Returns difference of two Complex Objects. constexpr Complex operator-(const Complex& C2) const noexcept { return { real-C2.real, imag-C2.imag }; } /// operator * : Returns product of two Complex Objects. constexpr Complex operator*(const Complex& C2) const noexcept { return { (real*C2.real)-(imag*C2.imag), (imag*C2.real)+(C2.imag*real) }; } /// operator / : Returns quotient of two Complex Objects. constexpr Complex operator/(const Complex& C2) const noexcept { return *this * inv(C2); } /// operator = : Assigns another Complex Object to Complex Object. constexpr Complex& operator= (const Complex& C2) noexcept { if(this!=&C2) { real=C2.real; imag=C2.imag; } return *this; } /// operator += : Assigns Complex object with the sum of two Complex objects. constexpr Complex& operator+=(const Complex& C2) noexcept { return *this=*this+C2; } /// operator -= : Assigns Complex object with the difference of two Complex objects. constexpr Complex& operator-=(const Complex& C2) noexcept { return *this=*this-C2; } /// operator *= : Assigns Complex object with the product of two Complex objects. constexpr Complex& operator*=(const Complex& C2) noexcept { return *this=*this*C2; } /// operator /= : Assigns Complex object with the quotient of two Complex objects. constexpr Complex& operator/=(const Complex& C2) noexcept { return *this=*this/C2; } /// operator == : Checks if the Complex Objects are equal. constexpr bool operator==(const Complex& C2) const noexcept { return (real==C2.real)&&(imag==C2.imag); } /// operator != : Checks if the Complex Objects are unequal. constexpr bool operator!=(const Complex& C2) const noexcept { return !(*this==C2); } /// operator < : Checks if the Complex Object is less than another Complex object. constexpr bool operator< (const Complex& C2) const noexcept { return abs(*this) : Checks if the Complex Object is greater than another Complex object. constexpr bool operator> (const Complex& C2) const noexcept { return abs(*this)>abs(C2); } /// operator >= : Checks if the Complex Object is greater than or equal to another Complex object. constexpr bool operator>=(const Complex& C2) const noexcept { return (*this>C2)||(*this==C2); } /// Returns the Real value of the Complex Object. friend constexpr Fraction re(const Complex& c) { return c.real; } /// Returns the Imaginary value of the Complex Object. friend constexpr Fraction im(const Complex& c) { return c.imag; } /// Returns the Principal Argument of the Complex Object. friend constexpr Fraction arg (const Complex& c) { return atan2(static_cast(c.imag), static_cast(c.real)); } /// Returns the Modulus of the Complex Object, or |z|. friend constexpr Fraction abs (const Complex& c) { return sqrt(static_cast((c.real*c.real)+(c.imag*c.imag))); } /// Returns the Normal of the Complex Object, or |z|^2. friend constexpr Fraction norm(const Complex& c) { return (c.real*c.real)+(c.imag*c.imag); } /// Returns the Conjugate of the Complex Object. friend constexpr Complex conj(const Complex& c) { return { c.real, -c.imag }; } /// Returns the Inverse of a Complex Object, or 1/z. friend constexpr Complex inv (const Complex& c) { return conj(c)/norm(c); } /// operator+ : Adds constant to Complex Object. friend constexpr Complex operator+ (const Fraction& rs, const Complex& c) { return { c.real+rs, c.imag }; } /// operator- : Subtracts constant from Complex Object. friend constexpr Complex operator- (const Fraction& rs, const Complex& c) { return { rs-c.real,-c.imag }; } /// operator* : Multiplies Complex Object by a constant. friend constexpr Complex operator* (const Fraction& rs, const Complex& c) { return { c.real*rs, c.imag*rs }; } /// operator/ : Divides Complex Object by a constant. friend constexpr Complex operator/ (const Fraction& rs, const Complex& c) { return { rs * inv(c) }; } /// operator == : Checks if the Complex Object and the scalar are equal. friend constexpr bool operator==(const Fraction& rs, const Complex& c) { return c==rs; } /// operator != : Checks if the Complex Object and the scalar are unequal. friend constexpr bool operator!=(const Fraction& rs, const Complex& c) { return !(c==rs); } /// operator < : Checks if the Complex Object is less than a scalar. friend constexpr bool operator< (const Fraction& rs, const Complex& c) { return c>rs; } /// operator <= : Checks if the Complex Object is less than or equal to a scalar. friend constexpr bool operator<=(const Fraction& rs, const Complex& c) { return c>rs||c==rs; } /// operator > : Checks if the Complex Object is greater than a scalar. friend constexpr bool operator> (const Fraction& rs, const Complex& c) { return c= : Checks if the Complex Object is less than or equal to a scalar. friend constexpr bool operator>=(const Fraction& rs, const Complex& c) { return crs; } /// operator <= : Checks if the Complex Object is less than or equal to a scalar. friend constexpr bool operator<=(long long rs, const Complex& c) { return c>rs||c==rs; } /// operator > : Checks if the Complex Object is greater than a scalar. friend constexpr bool operator> (long long rs, const Complex& c) { return c= : Checks if the Complex Object is less than or equal to a scalar. friend constexpr bool operator>=(long long rs, const Complex& c) { return c>(std::istream&, Complex&); friend std::ostream& operator<<(std::ostream&, const Complex&); }; /// operator >> : Reads Complex Object from specified stream. std::istream& operator>>(std::istream& is, Complex& c) { // discard trailing whitespace, similar to read operations of normal integers and all. while(is.peek()==' '||is.peek()=='\t'||is.peek()=='\n') is.get(); // approach : peek, follow the path discovered, save work so far, and continue till the end. char c1 = is.peek(),c2; c.real = 0.0; c.imag = 0.0; if(c1=='+'||c1=='-') // number in form +... or -... { c2 = c1; c1 = is.get(); c1 = is.peek(); if(c1=='i'||c1=='j') // number in form +i... or -j... { c.imag = ((c2=='-')?-1.0:1.0); is.get(); c1=is.peek(); if(c1=='+'||c1=='-') is>>c.real; // number in form +i-4 or -j+3 else {} // number in form +i or -j } else // number in form +3... or -4... { is>>c.real; if(c2=='-') { c.real=-c.real; } c1=is.peek(); if(c1=='+'||c1=='-') // number in form +3-... or -4+... { c2 = c1; is.get(); c1=is.peek(); if(c1=='i'||c1=='j') { c.imag = ((c2=='-')?-1.0:1.0); is.get(); } // number in form +3-i or -4+j else { is>>c.imag>>c1; if(c2=='-') c.imag=-c.imag; } // number in form +3+4i or -4+6j } else if(c1=='i'||c1=='j') // number in form +3i... or -4j... { std::swap(c.real,c.imag); is.get(); c1=is.peek(); if(c1=='+'||c1=='-') { is>>c.real; } // number in form +3i+7 or -4j-5 else {} // number in form +3i or -4j } else { } // number in form +3 or -4 } } else if(c1=='i'||c1=='j') // number in form i... or j... { c.imag = 1.0; is.get(); c1=is.peek(); if(c1=='+'||c1=='-') { is>>c.real; } // number in form i+3 or j-7 else {} // number in form i or j } else // number in form 3... or 4... { is>>c.real; c1=is.peek(); if(c1=='+'||c1=='-') // number in form 3+... or 4-... { c2 = c1; is.get(); c1=is.peek(); if(c1=='i'||c1=='j') { c.imag = ((c2=='-')?-1.0:1.0); is.get(); } // number in form 3+i or 4-j else { is>>c.imag>>c1; if(c2=='-') c.imag=-c.imag; } // number in form 3+7i or 4-5j } else if(c1=='i'||c1=='j') // number in form 3i... or 4j... { std::swap(c.real,c.imag); is.get(); c1=is.peek(); if(c1=='+'||c1=='-') { is>>c.real; } // number in form 3i+7 or 4j-5 else {} // number in form 3i or 4j } else {} // number in form 3 or 4 } return is; } /// operator << : Prints Complex Object to specified stream. std::ostream& operator<<(std::ostream& os, const Complex& c) { std::stringstream ss; ss.setf(os.flags()); ss.precision(os.precision()); os.unsetf(std::ios::showpos); ss.unsetf(std::ios::showpos); if(c.real!=0||c.imag==0) ss<0?(c.real==0?"":"+"):"-"); if(abs(c.imag)==1) ss<<"i"; else ss<<(abs(c.imag))<<"i"; } os< class Matrix; template std::ostream& operator<<(std::ostream&, const Matrix&); std::ostream& operator<<(std::ostream&, const Matrix&); /// Generic Class to represent Matrices and perform operations on them. template class Matrix { size_t r,c,_r,_c; M** mat; static M* parse(std::istream& is, std::string& s, size_t& c) { bool inc = c == 0; size_t i = 0; s+=" "; std::stringstream ss(s); ss.setf(is.flags()); M dummy; if(inc) { while(ss>>dummy) { c++; } } ss.clear(); ss.str(s); M* res = new M[c]; while(i>res[i++]) { } while(i(0); } return res; } std::pair print_args() const noexcept; public: struct row_iterator { using iterator_category = std::random_access_iterator_tag; using value_type = M*; using difference_type = std::ptrdiff_t; using pointer = value_type*; using reference = value_type; row_iterator() : p(nullptr) {} row_iterator(pointer x) : p(x) {} row_iterator(const row_iterator& rit) : p(rit.p) { } row_iterator(row_iterator&& rit) : p(rit.p) { rit.p = nullptr; } value_type operator*() { return *p; } row_iterator operator+ (difference_type n) {return p+n;} row_iterator operator- (difference_type n) {return p-n;} row_iterator& operator+=(difference_type n) {*this=*this+n; return *this;} row_iterator& operator-=(difference_type n) {*this=*this-n; return *this;} difference_type operator-(const row_iterator& q) {return static_cast((p-q.p));} friend row_iterator operator+(difference_type n, const row_iterator& c) {return c.p+n;} value_type operator[](difference_type n) { return *(p+n); } operator bool() const noexcept { return p==nullptr; } pointer operator->() { return p; } row_iterator& operator++() {p++; return *this;} row_iterator operator++(int) {row_iterator tmp(*this); ++*this; return tmp;} row_iterator& operator--() {p--; return *this;} row_iterator operator--(int) {row_iterator tmp(*this); --*this; return tmp;} row_iterator& operator=(const row_iterator& ri) { if(&ri!=this){this->p=ri.p; } return *this; } row_iterator& operator=(row_iterator&& ri) { if(&ri!=this){this->p=ri.p; ri.p = nullptr; } return *this; } bool operator==(const row_iterator& rhs) const {return p==rhs.p;} bool operator!=(const row_iterator& rhs) const {return p!=rhs.p;} bool operator< (const row_iterator& rhs) const {return p< rhs.p;} bool operator<=(const row_iterator& rhs) const {return p<=rhs.p;} bool operator> (const row_iterator& rhs) const {return p> rhs.p;} bool operator>=(const row_iterator& rhs) const {return p>=rhs.p;} protected: pointer p; }; using const_row_iterator = const row_iterator; using reverse_row_iterator = std::reverse_iterator; using const_reverse_row_iterator = const reverse_row_iterator; using iterator = row_iterator; using reverse_iterator = reverse_row_iterator; using const_iterator = const iterator; using const_reverse_iterator = const reverse_iterator; iterator begin() { return Matrix::iterator(&mat[0]); } iterator end() { return Matrix::iterator(&mat[r]); } const_iterator cbegin() const { return Matrix::const_iterator(&mat[0]); } const_iterator cend() const { return Matrix::const_iterator(&mat[r]); } reverse_iterator rbegin() { return Matrix::reverse_iterator(this->end()); } reverse_iterator rend() { return Matrix::reverse_iterator(this->begin()); } const_reverse_iterator crbegin() const { return Matrix::const_reverse_iterator(this->end()); } const_reverse_iterator crend() const { return Matrix::const_reverse_iterator(this->begin()); } using scalar_type = M; enum class Operation : bool { Row, Column }; struct flags { static const std::ios_base::fmtflags noformat = static_cast(0x100000), nobrackets = static_cast(0x200000); }; Matrix(unsigned short row=1, unsigned short col=1, M&& val=M(0)); Matrix(std::initializer_list>&&); Matrix(const Matrix&); Matrix(Matrix&&); ~Matrix(); static Matrix identity(size_t); static std::string format(size_t, std::vector, Operation); static std::ostream& nobrackets(std::ostream& os) { os.setf(Matrix::flags::nobrackets); os.unsetf(Matrix::flags::noformat); return os; } static std::ostream& noformat(std::ostream& os) { os.setf(Matrix::flags::noformat); os.unsetf(Matrix::flags::nobrackets); return os; } static std::ostream& normal(std::ostream& os) { os.unsetf(Matrix::flags::noformat); os.unsetf(Matrix::flags::nobrackets);return os; } void clear() noexcept; bool empty() { return r==0&&c==0; } Matrix& transform(size_t, M*, Operation = Operation::Row); /// Returns the Number of Rows of the Matrix object. size_t row_size() const noexcept { return r; } /// Returns the Number of Columns of the Matrix object. size_t col_size() const noexcept { return c; } /// Returns the size of the Matrix Object as a pair. std::pair size() const noexcept { return std::make_pair(r,c); } /// Returns the number of elements present in the matrix. size_t count() { return r*c; } /// Returns the number of allocated rows for the Matrix object. size_t row_capacity() const noexcept { return _r; } /// Returns the number of allocated columns for the Matrix object. size_t col_capacity() const noexcept { return _c; } /// Returns the capacity of the Matrix Object as a pair. std::pair capacity() const noexcept { return std::make_pair(_r,_c); } bool operator==(const Matrix&) const noexcept; bool operator!=(const Matrix&) const noexcept; M& operator()(size_t, size_t) const; Matrix& operator= (Matrix&&); Matrix& operator= (const Matrix&); template friend size_t rank (const Matrix& m); template friend Matrix ref (const Matrix& m); template friend Matrix rref (const Matrix& m); template friend std::vector>> rref_stepwise (const Matrix& m); template friend void multi_print(std::ostream& os, std::string opsym, Args... args); template friend Matrix operator*(const MT&, const Matrix&); template friend std::istream& operator>>(std::istream&, Matrix&); friend std::ostream& operator<< <>(std::ostream& os, const Matrix& m); friend std::ostream& operator<<(std::ostream& os, const Matrix& m); }; /// Returns the Row Echelon Form of the Matrix object. template Matrix ref(const Matrix& m) { Matrix res(m); M* arr = new M[m.r]; size_t x = 0; for(size_t i=0;i(0)) { for(size_t j=i+1;j(0)) { for(size_t k=0;k(0); } arr[j]=static_cast(1); res.transform(i,arr); break; } } if(res.mat[i][x]!=static_cast(0)) { if(res.mat[i][x]!=static_cast(1)) { for(size_t j=0;j(0); arr[i] = static_cast(1)/res.mat[i][x]; res.transform(i,arr); } for(size_t j=i+1;j(0)) { for(size_t k=0;k(0); arr[j] = static_cast(1); arr[i] = static_cast(-1)*res.mat[j][x]; res.transform(j,arr); } x++; break; } } } delete[] arr; return res; } /// Returns the Reduced Row Echelon Form of the Matrix object. template Matrix rref(const Matrix& m) { Matrix res(m); M* arr = new M[m.r]; size_t x = 0; for(size_t i=0;i(0)) { for(size_t j=i+1;j(0)) { for(size_t k=0;k(0); } arr[j]=static_cast(1); res.transform(i,arr); break; } } if(res.mat[i][x]!=static_cast(0)) { if(res.mat[i][x]!=static_cast(1)) { for(size_t j=0;j(0); arr[i] = static_cast(1)/res.mat[i][x]; res.transform(i,arr); } for(size_t j=0;j(0)) { for(size_t k=0;k(0); arr[j] = static_cast(1); arr[i] = static_cast(-1)*res.mat[j][x]; res.transform(j,arr); } x++; break; } } } delete[] arr; return res; } /// Returns a vector of steps required to convert a Matrix object into its reduced row echelon form. template std::vector>> rref_stepwise(const Matrix& m) { std::vector>> rstep; std::vector vec; Matrix res(m); M* arr = new M[m.r]; size_t x = 0; for(size_t i=0;i(0)) { //std::cout<<"Value is 0. Swap May be performed...\n"; for(size_t j=i+1;j(0)) { for(size_t k=0;k(0); } arr[j]=static_cast(1); //std::cout<<"Swapping row "<> "<::format(i,vec,Matrix::Operation::Row)<<"\n"; res.transform(i,arr); break; } //std::cout<<"Swap check complete.\n"; } if(res.mat[i][x]!=static_cast(0)) { if(res.mat[i][x]!=static_cast(1)) { //std::cout<<"Factorizing row "<(0); arr[i] = static_cast(1)/res.mat[i][x]; vec.assign(arr,arr+m.r); rstep.emplace_back(i,vec); res.transform(i,arr); //std::cout<<" >> "<::format(i,vec,Matrix::Operation::Row)<<"\n"; //std::cout<<"Factorization complete.\n"; } //std::cout<(0)) { for(size_t k=0;k(0); arr[j] = static_cast(1); arr[i] = static_cast(-1)*res.mat[j][x]; vec.assign(arr,arr+m.r); rstep.emplace_back(j,vec); res.transform(j,arr); //std::cout<<" >> "<::format(i,vec,Matrix::Operation::Row)<<"\n"; } //std::cout<<"Normalization of column complete.\n\n"; std::cout< size_t rank(const Matrix& m) { size_t x = 0, rank_ct = 0; for(const auto& r : rref(m)) { for(size_t i=x;i=m.col_size()) break; } return rank_ct; } #if __cplusplus>=201703 /// Prints multiple Matrix Objects in the same rows (one side by side with others), along with additional operators between them. template void multi_print(std::ostream& os, std::string opsym, Args... args) { int cct = 0; opsym = " "+opsym; auto tsz = opsym.size(); std::string out; std::stringstream ss; auto tmp = std::get<0>(std::forward_as_tuple(args...)); Matrix max_m(1,1); auto find_max = [&max_m](const auto& a){ if(a.r>max_m.r) max_m = a; }; (find_max(args),...); size_t i=0, mid = (max_m.r)/2; for(;i=mid&&row_idx<(mid+arg.r)) ss<> : Reads the Matrix from an istream object, without the need of entering rows or columns. template std::istream& operator>>(std::istream& is, Matrix& m) { for(size_t i=0;i arr; while(getline(is,str)&&!str.empty()) { arr.push_back(str); } m.r = arr.size(); m.mat = new M*[m.r]; for(size_t i=0;i::parse(is,arr[i],m.c); // m.reserve(50,50); while(getline(is,str)&&!str.empty()) { m.push_back(math::Matrix::parse(is,str,m.c)); } return is; return is; } /// operator<< : Prints the Matrix to an ostream object, with perfect alignment and appropriate brackets. template std::ostream& operator<<(std::ostream& os, const Matrix& m) { if((os.flags()&Matrix::flags::noformat)==0x100000) { for(size_t i=0;i::flags::nobrackets)==0x200000) { std::stringstream ss; size_t mwidth = 0; bool neg = false; for(size_t i=0;imwidth) mwidth = ss.str().size(); ss.str(""); } } if(neg) { os.setf(std::ios::showpos); } for(size_t i=0;imwidth) mwidth = ss.str().size(); ss.str(""); } } if(neg) { os.setf(std::ios::showpos); } os<& m) { if((os.flags()&Matrix::flags::noformat)==0x100000) { for(size_t i=0;i::flags::nobrackets)==0x200000) { std::stringstream ss; ss.setf(os.flags()); size_t mwidth = 0; for(size_t i=0;imwidth) mwidth = ss.str().size(); ss.str(""); } } for(size_t i=0;imwidth) mwidth = ss.str().size(); ss.str(""); } } os< std::pair math::Matrix::print_args() const noexcept { std::stringstream ss; size_t mwidth = 0; bool neg = false; for(size_t i=0;imwidth) mwidth = ss.str().size(); ss.str(""); } } return std::make_pair(mwidth,neg); } /// Matrix Constructor 0 - Creates a Matrix of size 'r'x'c', and assign all its elements with 'val' template math::Matrix::Matrix(unsigned short row, unsigned short col, M&& val) : r(row),c(col), _r(r), _c(c) { mat = new M*[r]; for(size_t i=0;i math::Matrix::Matrix(std::initializer_list>&& ls) : r(ls.size()), c(ls.begin()->size()), _r(r), _c(c) { mat = new M*[r]; for(auto it=ls.begin(); it!=ls.end(); it++) { mat[std::distance(ls.begin(),it)] = new M[c]; for(auto iit = it->begin(); iit!=it->end(); iit++) mat[std::distance(ls.begin(),it)][std::distance(it->begin(),iit)] = *(iit); } } /// Matrix Constructor 2 - Copy Constructor for deep/logical copying. template math::Matrix::Matrix(const Matrix& M2) : r(M2.r), c(M2.c) { mat = new M*[r]; for(size_t i=0;i math::Matrix::Matrix(math::Matrix&& M2) : r(M2.r), c(M2.c), _r(r), _c(c), mat(M2.mat) { M2.r = 0ULL; M2.c = 0ULL; M2._r = 0ULL; M2._c=0ULL; M2.mat = nullptr; } /// Matrix Destructor - For Deallocation of reserved memory. template math::Matrix::~Matrix() { for(size_t i=0;i math::Matrix math::Matrix::identity(size_t order) { math::Matrix eye(order,order); for(size_t i=0;i std::string math::Matrix::format(size_t idx, std::vector opr, math::Matrix::Operation type) { std::stringstream ss; char c = (type==math::Matrix::Operation::Row?'R':'C'); ss<0?"":"- "); if(abs(opr[k])!=1) { ss<0?" +":" -"); if(abs(opr[k])!=1) { ss<<" "< std::string math::Matrix::format(size_t idx, std::vector opr, math::Matrix::Operation type) { std::stringstream ss; char c = (type==math::Matrix::Operation::Row?'R':'C'); ss<0?"":"- "); if(abs(re(opr[k]))!=1) { ss<0?" +":" -"); if(abs(re(opr[k])!=1)) { ss<<" "< math::Matrix& math::Matrix::transform(size_t idx, M* opr, math::Matrix::Operation type) { size_t ones = 0, zeros = 0, id = 0; for(size_t i=0;i<(type==math::Matrix::Operation::Row?r:c);i++) { if(opr[i]==static_cast(0)) zeros++; else if(opr[i]==static_cast(1)) { ones++; id = i; } } if(ones==1&&zeros==(type==math::Matrix::Operation::Row?r-1:c-1)&&id!=idx) { switch(type) { case math::Matrix::Operation::Row: { M tmp; for(size_t i=0;i::Operation::Column: { M tmp; for(size_t i=0;i::Operation::Row: { M* tmp = new M[c]; for(size_t i=0;i(0); for(size_t i=0;i::Operation::Column: { M* tmp = new M[c]; for(size_t i=0;i(0); for(size_t i=0;i bool math::Matrix::operator==(const math::Matrix& M2) const noexcept { if(M2.r==r&&M2.c==c) { for(size_t i=0;i bool math::Matrix::operator!=(const math::Matrix& M2) const noexcept { return !(*this==M2); } /// operator() : Returns a reference to the element at i,j in the Matrix object. template inline M& math::Matrix::operator()(size_t i, size_t j) const { if(i<0||i>r-1||j<0||j>c-1) { throw std::out_of_range("MatrixError : operator(size_t,size_t) : Index for access is out of bounds!"); } return mat[i][j]; } /// operator= : Allocates new memory and performs deep/logical assignment. template math::Matrix& math::Matrix::operator= (const math::Matrix& m) { if(this!=&m) { for(size_t i=0;i math::Matrix& math::Matrix::operator= (math::Matrix&& m) { for(size_t i=0;ir = m.r; m.r = 0ULL; this->c = m.c; m.c = 0ULL; this->_r=_r; m._r=0ULL; this->_c=_c; m._c=0ULL; this->mat = m.mat; m.mat = nullptr; } return *this; } int main() { try { math::Matrix M1,M2; std::cin>>M1; M2=M1; std::cout<<"Your Matrix : \n\n"<>> steps = rref_stepwise(M1); for(std::pair>& step : steps) { std::cout<<">> "<::format(step.first,step.second,math::Matrix::Operation::Row)<<"\n"; std::cout<> Result in "<> Matrix, its row echelon and its reduced row echelon : \n\n"; multi_print(std::cout,",,",M2,ref(M2),M1); std::cout<<">> Rank of the Matrix = "<

Matrix to Row Echelon