﻿using System;
using System.Collections.Generic;
using MatLib.Numerics;
using System.Text.RegularExpressions;

/// <summary>
/// Contains support for matrix manipulation
/// </summary>
namespace MatLib
{
    /// <summary>
    /// Contains support for matrices and vectors
    /// </summary>
    namespace ArrayTypes
    {


        /// <summary>
        /// Enumerator that stores the different types of matrices to aid in calculation
        /// </summary>
        [Flags]
        internal enum MatType
        {
            /// <summary>
            /// A matrix that has no outstanding properties
            /// </summary>
            NoProperty = 0x00,

            /// <summary>
            /// identity a_ij=0 if i!=j and 1 elsewhere
            /// </summary>
            Identity = 0x01,

            /// <summary>
            /// Zero if all elements are zero
            /// </summary>
            Zero = 0x02,

            /// <summary>
            ///square if _m=_n
            /// </summary>
            Square = 0x04,

            /// <summary>
            ///diagonal if all elements not on the leading diagonal are zero
            /// </summary>
            Diagonal = 0x08,

            /// <summary>
            ///Upper Triangular if all nonzero elements are above or on the leading diagonal
            /// </summary>
            UT = 0x10,

            /// <summary>
            ///Lower Triangular if all nonzero elements are below or on the leading diagonal
            /// </summary>
            LT = 0x20

        }

        /// <summary>
        /// Provides interoperability between matrices and vectors. NOT IMPLEMENTED
        /// </summary>
        /// <typeparam name="T">A type inheriting Real and new()</typeparam>
        public class Interop
        {
            /// <summary>
            /// Extract a column Vector from a Matrix
            /// </summary>
            /// <param name="InP">The matrix to shave from</param>
            /// <param name="OutP">The vector to return</param>
            /// <param name="m">The zero based index of the column to extract, the last index by default</param>
            /// <returns>'InP' without the column at 'm'</returns>
            public static Matrix Slice(Matrix InP, out Vector OutVec, int m = -1)
            {
                //validation
                if ((m < -1) || (m >= InP.Length))
                { throw new IndexOutOfRangeException("Column index is out of range"); }
                //setting variables
                if (m == -1) { m = InP.Length - 1; }
                //formatting the variables
                var OutMat = new Matrix(InP.Height, InP.Length - 1);
                OutVec = new Vector(InP.Height);
                //counter for matrix columns
                int c = 0;

                //for every column
                for (int j = 0; j < InP.Length; j++)
                {
                    //if on the vector's column 
                    if (j == m)
                    {
                        //extract the vector
                        for (int i = 0; i < InP.Height; i++)
                        {
                            OutVec[i] = InP[i, j];
                        }
                    }
                    else
                    {
                        //add the column to the matrix
                        for (int i = 0; i < InP.Height; i++)
                        {
                            OutMat[i,c] = InP[i, j];
                        }
                        c++;
                    }

                }
                return OutMat;
            }


        }

        /// <summary>
        /// A matrix of type T where T is derived from Real
        /// </summary>
        /// <typeparam name="T">A type that is derived from Real (Frac|Floating)</typeparam>
        public class Matrix
        {
            //Matrix specific types, variables and methods
            #region
            /// <summary>
            /// Tests the flag "y" against type to see if it's triggered
            /// </summary>
            /// <param name="y">The flag to be tested for</param>
            /// <returns>A boolean indicating whether "y" is triggered in type</returns>
            private bool TestForFlag(MatType y) { return (type & y) == y; }

            /// <summary>
            /// Stores the rows of the matrix as arrays of whatever type the matrix was called as
            /// </summary>
            private List<Real[]> storage;

            /// <summary>
            /// Stores the number of rows
            /// </summary>
            private int _m;

            /// <summary>
            /// Stores the number of columns
            /// </summary>
            private int _n;

            /// <summary>
            /// Stores the type type of the matrix
            /// </summary>
            private MatType type;

            /// <summary>
            /// Checks the type of a matrix.
            /// </summary>
            private void CheckMatType()
            {
                //Clean the type
                type = MatType.NoProperty;

                //assumed to be of every matrix type (even though that's impossible) and the types are steadily eliminated
                type |= (MatType.Diagonal | MatType.Identity | MatType.LT | MatType.Square | MatType.UT | MatType.Zero);

                //testing for squareness is relatively easy, also tests for identity as well
                if (_m != _n) { type ^= MatType.Square; type ^= MatType.Identity; }

                //Row number, position along a column
                int i = 0;
                //Column number, position along a row
                int j = 0;

                //for every row in the matrix
                foreach (Real[] Row in storage)
                {
                    foreach (Real Element in Row)
                    {
                        //general tests for all elements
                        //zero test
                        if ((TestForFlag(MatType.Zero) && !Element.IsZero())) { type ^= MatType.Zero; }

                        //tests for elements on the main diagonal
                        if (i == j)
                        {
                            //identity test
                            if ((TestForFlag(MatType.Identity)) && (Element.GetValue() != 1))
                            { type ^= MatType.Identity; }
                        }
                        //tests for elements not on the main diagonal
                        else
                        {
                            //diagonal test
                            if ((TestForFlag(MatType.Diagonal)) && !Element.IsZero()) { type ^= MatType.Diagonal; }

                            //identity test
                            if ((TestForFlag(MatType.Identity)) && !Element.IsZero()) { type ^= MatType.Identity; }

                            //Lower Triangular test
                            //tests for elements above the main diagonal (more advanced along a row than along a column)
                            if (TestForFlag(MatType.LT) && j > i && !Element.IsZero()) { type ^= MatType.LT; }

                            //Upper Triangular test
                            //tests for elements below the main diagonal (more advanced along a column than along a row)
                            if (TestForFlag(MatType.UT) && j < i && !Element.IsZero()) { type ^= MatType.UT; }
                        }
                        j++;
                    }
                    j = 0;
                    i++;
                }

                //Ironing out inconsistencies
                /*
                 * If matrix is a zero matrix, it would also pass the diagonal, UT and LT tests
                 * The same goes for an identity matrix which would pass the diagonal, UT and LT tests
                 * Lastly a diagonal matrix would also pass the UT and LT tests
                 * While this isn't neccessarily incorrect, it could present problems later when doing calculations
                 * The following tests make sure only one tag is triggered
                 */

                //using the property that enums can be cast to integers
                if ((int)type == (int)(MatType.Zero | MatType.Diagonal | MatType.UT | MatType.LT))
                { type ^= (MatType.Diagonal | MatType.UT | MatType.LT); }
                if (TestForFlag(MatType.Identity))
                { type ^= (MatType.Diagonal | MatType.UT | MatType.LT); }
                if (TestForFlag(MatType.Diagonal))
                { type ^= (MatType.UT | MatType.LT); }


            }

            /// <summary>
            /// Recursively parse a string to see whether it represents a nested list
            /// </summary>
            /// <param name="phrase">string to be parsed</param>
            /// <param name="outList">List to contain the output of string if there is any</param>
            /// <param name="pointer">The current read position of the algorithm</param>
            /// <param name="nested">Whether this call is a nested call or not</param>
            /// <returns>Whether 'phrase' can be parsed as a nested list</returns>
            private static bool RecursiveParse(string phrase, out List<object> outList, ref int pointer, bool nested = false)
            {
                //list to output
                outList = new List<object>();
                //whether an element has ended
                bool endElem = true;
                //a real value passed as an output variable
                Real tmpVal = new Fraction();
                //list of elements to be passed as an output variable
                List<object> elementList = new List<object>();
                //an accumulator string
                string tmpStr = "";

                //for the length of the string
                while (pointer < phrase.Length)
                {
                    switch (phrase[pointer])
                    {
                        //in the case of an open bracket, increase the read pointer and recurse
                        case '[':
                            pointer++;
                            //if an element has not ended or if the parse fails on the nested list
                            //clear the output list and return false
                            if (!(endElem && RecursiveParse(phrase, out elementList, ref pointer, true))) { outList.Clear(); return false; }
                            //if it succeeds add the returned value to the output list
                            outList.Add(elementList);

                            break;
                        //close brackets indicate the end of a successful parse but only if an element has ended, increment the pointer and return true
                        case ']':
                            if (!endElem) { outList.Clear(); return false; }
                            pointer++;
                            return true;
                        //semicolons indicate the end of an element, attempt to parse it and add it to the list
                        case ';':
                            //an element has finished
                            endElem = true;
                            //an empty temporary string implies that the item just parsed was a list and the parse can skip this
                            if (tmpStr == "") { pointer++; break; }
                            //try parsing the string as a fraction or a floating
                            
                            if (!(Real.TryParse(tmpStr, out tmpVal)))
                            { outList.Clear(); return false; }
                            //if it succeeds add the returned value to the output list
                            outList.Add(tmpVal);
                            //reset the temp string and value
                            tmpStr = "";
                            tmpVal = new Fraction();
                            //increment the pointer
                            pointer++;
                            break;
                        default:
                            //an element has begun
                            endElem = false;
                            //add the element to the parse string and increment the pointer
                            tmpStr += phrase[pointer];
                            pointer++;
                            break;

                    }
                }
                //if the elements are not finalised, return false
                if (!endElem || nested) { outList.Clear(); return false; }
                return true;

            }
            #endregion

            //Matrix suboperations
            #region
            /// <summary>
            /// Swap two rows
            /// </summary>
            /// <param name="index1">Index of the first row to be swapped</param>
            /// <param name="index2">Index of the second row to be swapped</param>
            private void Swap(int index1, int index2)
            {
                Real[] tmp = storage[index1];
                storage[index1] = storage[index2];
                storage[index2] = tmp;
            }

            /// <summary>
            /// Multiply a row by a scalar
            /// </summary>
            /// <param name="index">Row to be multiplied</param>
            /// <param name="multiplier">Real value to multiply the row by</param>
            private void MultiplyByScalar(int index, Real multiplier)
            {
                //for every column
                for (int i = 0; i < _n; i++)
                {
                    //multiply the value in the given row by the multiplier
                    storage[index][i] = storage[index][i] * multiplier;
                }
            }

            /// <summary>
            /// Subtract a scalar multiple of a row from another
            /// </summary>
            /// <param name="index1">The row to be subtracted from</param>
            /// <param name="index2">The row to be subtracted</param>
            /// <param name="multiplier">The scalar multiple to multiply the row at "index2" by</param>
            private void SubtractRow(int index1, int index2, Real multiplier)
            {
                //perform the subtraction
                for (int i = 0; i < _n; i++)
                { storage[index1][i] = storage[index1][i] - (multiplier * storage[index2][i]); }

            }
            #endregion

            //Constructors and Casts
            #region
            /// <summary>
            /// Construct a matrix using a 2D array of values
            /// </summary>
            /// <param name="Values">Array of values to be stored in the matrix</param>
            public Matrix(Real[,] Values)
            {
                //instantiate storage
                storage = new List<Real[]>();

                //set the height and width
                _m = Values.GetLength(0);
                _n = Values.GetLength(1);

                //if either the length or height of the matrix is zero
                //throw an exception
                if ((_m == 0) || (_n == 0))
                { throw new Exception("Cannot have a matrix with no elements", new ArgumentNullException()); }

                //if either the length or height of the matrix is one
                //throw an exception
                if ((_m == 1) || (_n == 1))
                { throw new Exception("Cannot have a matrix with one row or column, please use a matrix instead", new ArgumentNullException()); }

                //for every row
                for (int i = 0; i < _m; i++)
                {
                    //copy a row
                    Real[] row = new Real[_n];
                    for (int j = 0; j < _n; j++)
                    { row[j] = Values[i, j]; }
                    //add that row
                    storage.Add(row);
                }


                CheckMatType();
            }

            /// <summary>
            /// Construct the "n" by "n" identity matrix
            /// </summary>
            /// <param name="n">The size of the identity matrix</param>
            public Matrix(int n)
            {
                //Instance of T to call the factory method
                Real instance = new Fraction();
                List<Real[]> OutP = new List<Real[]>();
                //add the value
                for (int i = 0; i < n; i++)
                {
                    Real[] tmp = new Real[n];
                    for (int j = 0; j < n; j++)
                    {
                        tmp[j] = (i == j) ? instance.FromInt(1) : instance.FromInt(0);
                    }
                    OutP.Add(tmp);
                }

                //copy the elements over
                storage = OutP;
                //assign properties values
                _n = n; _m = n;

                //setting the type
                type |= MatType.Identity | MatType.Square;
            }

            /// <summary>
            /// Constructor for an "m" by "n" zero matrix
            /// </summary>
            /// <param name="m">The height of the matrix</param>
            /// <param name="n">The width of the matrix</param>
            public Matrix(int m, int n)
            {
                _n = n;
                _m = m;
                storage = new List<Real[]>();
                if ((_m == 1) || (_n == 1))
                {
                    throw new Exception(
                        String.Format("Matrix is of dimensions {0} by {1}, use {2} instead.", _m, _n, (_m == _n) ? "MatLib.ArrayTypes.Vector" : "Real"),
                        new ArgumentException());
                }

                //instance of T to call the factory method
                Real instance = new Fraction();
                //rows
                for (int i = 0; i < m; i++)
                {
                    //row list
                    Real[] temp = new Real[n];
                    //columns
                    for (int j = 0; j < n; j++)
                    {
                        temp[j] = instance.FromInt(0);
                    }
                    //add to storage
                    storage.Add(temp);
                }
                //types
                type = MatType.Zero;
                if (_n == _m) { type |= MatType.Square; }
            }

            /// <summary>
            /// Test whether a string is parsable as a matrix and return that matrix if it is
            /// </summary>
            /// <param name="phrase">string to be parsed</param>
            /// <param name="OutP">Matrix to be returned, returns 1*1 indentity matrix if string cannot be parsed</param>
            /// <returns></returns>
            public static bool TryParse(string phrase, out Matrix OutP)
            {
                //OutP is initially a 1*1 matrix (the error value)
                OutP = new Matrix(1);

                //if a recursive parse fails, return false
                var x = 0;
                if (!RecursiveParse(phrase, out var parseList, ref x)) { return false; }

                var outList = (parseList[0] as List<object>);

                //outer list is redundant, inner list should be list of lists
                if ((outList is null) || ((outList[0] as List<object>) is null)) { return false; }

                //no empty or single lists or sublists
                if (outList.Count < 2 || (outList[0] as List<object>).Count < 2) { return false; }

                //storage matrix
                Real[,] outVals = new Real[outList.Count, (outList[0] as List<object>).Count];

                //for every object in the first list
                for (int i = 0; i < outVals.GetLength(0); i++)
                {
                    //element as a list
                    var iRow = outList[i] as List<object>;

                    //if it can't be cast as a list of objects or doesn't match the list length, return false
                    if ((iRow is null) || (iRow.Count != outVals.GetLength(1))) { return false; }

                    //add the items
                    for (int j = 0; j < outVals.GetLength(1); j++)
                    {
                        //item as a T
                        var tmpVal = iRow[j] as Real;
                        //if the cast fails return false
                        if (tmpVal is null) { return false; }
                        //add the item
                        outVals[i, j] = tmpVal;
                    }
                }
                //generate the matrix
                OutP = new Matrix(outVals);
                return true;

            }
            #endregion

            //Accessors
            #region
            /// <summary>
            /// Get the element at row "i", column "j"
            /// </summary>
            /// <param name="i">The zero-based row number of the element to be found</param>
            /// <param name="j">The zero-based column number of the element to be found</param>
            /// <returns>The element at row "i", column "j" of this matrix</returns>            
            public Real this[int i, int j]
            {
                get
                {
                    if (i < 0 || i >= _m) { throw new Exception("Row number is out of range", new ArgumentOutOfRangeException()); }
                    if (j < 0 || j >= _n) { throw new Exception("Column number is out of range", new ArgumentOutOfRangeException()); }
                    else { return storage[i][j]; }
                }
                internal set
                {
                    if (i < 0 || i >= _m) { throw new Exception("Row number is out of range", new ArgumentOutOfRangeException()); }
                    if (j < 0 || j >= _n) { throw new Exception("Column number is out of range", new ArgumentOutOfRangeException()); }
                    else { storage[i][j] = value; }
                }
            }

            /// <summary>
            /// Get the height of the Matrix
            /// </summary>
            public int Height
            { get { return _m; } }

            /// <summary>
            /// Get the length of the Matrix
            /// </summary>
            public int Length
            { get { return _n; } }
            #endregion

            //Matrix operations
            #region
            /// <summary>
            /// Add two matrices together
            /// </summary>
            /// <param name="other">The matrix to be added</param>
            /// <returns>The sum of this matrix and "other"</returns>
            private Matrix Add(Matrix other)
            {
                //size mismatch
                if (!((_m == other._m) && (_n == other._n)))
                { throw new Exception(String.Format("Matrix size mismatch, first matrix is of size {0}x{1} while second matrix is of size {2}x{3}", _m, _n, other._m, other._n), new ArgumentException()); }

                //efficiency checks
                if (TestForFlag(MatType.Zero))
                { return other; }
                if (other.TestForFlag(MatType.Zero))
                { return this; }

                //create an output matrix of the required size
                Matrix OutP = new Matrix(_m, _n);

                //for every row
                for (int i = 0; i < _m; i++)
                {
                    //for every column
                    for (int j = 0; j < _n; j++)
                    { OutP[i, j] = this[i, j] + other[i, j]; }
                }

                //check the matrix type
                OutP.CheckMatType();

                //return the output matrix
                return OutP;
            }

            /// <summary>
            /// Find the difference of two matrices
            /// </summary>
            /// <param name="other">The matrix to be subtracted</param>
            /// <returns>The difference between this matrtix and "other"</returns>
            private Matrix Subtract(Matrix other)
            {
                //size mismatch
                if (!((_m == other._m) && (_n == other._n)))
                { throw new Exception(String.Format("Matrix size mismatch, first matrix is of size {0}x{1} while second matrix is of size {2}x{3}", _m, _n, other._m, other._n), new ArgumentException()); }

                //efficiency checks
                if (TestForFlag(MatType.Zero))
                { return other.Negate(); }
                if (other.TestForFlag(MatType.Zero))
                { return this; }

                //create an output matrix of the required size
                Matrix OutP = new Matrix(_m, _n);

                //for every row
                for (int i = 0; i < _m; i++)
                {
                    //for every column
                    for (int j = 0; j < _n; j++)
                    { OutP[i, j] = this[i, j] - other[i, j]; }
                }

                //check the matrix type
                OutP.CheckMatType();

                //return the output matrix
                return OutP;
            }

            /// <summary>
            /// Multiply two matrices together
            /// </summary>
            /// <param name="other">The matrix to be right multiplied by</param>
            /// <returns>This matrix right multiplied by "other"</returns>
            private Matrix Multiply(Matrix other)
            {
                //size mismatch
                if (_n != other._m)
                { throw new Exception(String.Format("Matrix size mismatch, first matrix is of size {0}x{1} while second matrix is of size {2}x{3}, {1} != {2}", _m, _n, other._m, other._n), new ArgumentException()); }

                //efficiency checks
                //identity
                if (TestForFlag(MatType.Identity))
                { return other; }
                if (other.TestForFlag(MatType.Identity))
                { return this; }

                //zero
                if (TestForFlag(MatType.Zero) | other.TestForFlag(MatType.Zero))
                { return new Matrix(_m, other._n); }

                //create an output matrix of the required size
                Matrix OutP = new Matrix(_m, other._n);

                //diagonal matrices can be multiplied efficiently
                //left multiplication by diagonal
                if (TestForFlag(MatType.Diagonal))
                {
                    //for every row
                    for (int i = 0; i < OutP._m; i++)
                    {
                        //for every column
                        for (int j = 0; j < OutP._n; j++)
                        { OutP[i, j] = this[i, i] * other[i, j]; }
                    }
                    //check and return
                    OutP.CheckMatType();
                    return OutP;
                }
                //right multiplication by diagonal
                if (other.TestForFlag(MatType.Diagonal))
                {
                    //for every row
                    for (int i = 0; i < OutP._m; i++)
                    {
                        //for every column
                        for (int j = 0; j < OutP._n; j++)
                        { OutP[i, j] = this[i, j] * other[j, j]; }
                    }
                    //check and returm
                    OutP.CheckMatType();
                    return OutP;
                }

                //General matrix multiplication
                //for every row in the output
                for (int i = 0; i < _m; i++)
                {
                    //and for every column
                    for (int j = 0; j < other._n; j++)
                    {
                        //first product
                        OutP[i, j] = this[i, 0] * other[0, j];
                        //sum products over a row and column
                        for (int m = 1; m < _n; m++)
                        { OutP[i, j] = OutP[i, j] + this[i, m] * other[m, j]; }
                    }
                }
                //check and return
                OutP.CheckMatType();
                return OutP;

            }

            /// <summary>
            /// Multiply a matrix and a scalar together
            /// </summary>
            /// <param name="Multiplier">Scalar to multiply by</param>
            /// <returns>This matrix but with every element multiplied by 'Multiplier'</returns>
            private Matrix ScalarMultiply(Real Multiplier)
            {
                //create an output matrix of the required size
                Matrix OutP = new Matrix(_m, _n);

                for (int i = 0; i < _m; i++)
                {
                    for (int j = 0; j < _n; j++)
                    {
                        OutP[i, j] = Multiplier * this[i, j];
                    }
                }
                OutP.CheckMatType();
                return OutP;
            }

            /// <summary>
            /// Concatenate two matrices
            /// </summary>
            /// <param name="other">The matrix to be right concatenated</param>
            /// <returns>This matrix right concatenated with "other"</returns>
            public Matrix Concat(Matrix other)
            {
                if (_m != other._m) { throw new Exception(String.Format("Height mismatch, first matrix is of height {0} while second is of {1}", _m, other._m), new ArgumentException()); }
                Matrix OutP = new Matrix(_m, _n + other._n);
                //for every row
                for (int i = 0; i < _m; i++)
                {
                    //add the values of the first matrix
                    for (int j = 0; j < _n; j++)
                    { OutP[i, j] = this[i, j]; }
                    //and then the second
                    for (int j = 0; j < other._n; j++)
                    { OutP[i, j + _n] = other[i, j]; }
                }
                //return the combined matrix
                OutP.CheckMatType();
                return OutP;
            }

            /// <summary>
            /// Find whether this matrix is equal to another
            /// </summary>
            /// <param name="other">The matrix to compare this one against</param>
            /// <returns>Whether this matrix is equal to "other"</returns>
            private bool CompareTo(Matrix other)
            {
                //must have same dimensions
                if (_m != other._m || _n != other._n) { return false; }

                //must have same elements
                for (int i = 0; i < _m; i++)
                {
                    for (int j = 0; j < _n; j++)
                    {
                        if (this[i, j] != other[i, j]) { return false; }
                    }
                }
                return true;
            }

            /// <summary>
            /// Returns a value indicating whether this instance is equal to a specific object
            /// </summary>
            /// <param name="obj">The object to test equality against</param>
            /// <returns>A boolean value indicating whether this instance is equal to "obj"</returns>
            public override bool Equals(object obj)
            {
                Matrix temp = obj as Matrix;
                if (temp is null) { return false; }
                else { return this.CompareTo(temp); }
            }

            /// <summary>
            /// Returns the hash code for this instance
            /// </summary>
            /// <returns>The hash code for this instance</returns>
            public override int GetHashCode()
            {
                return storage.GetHashCode();
            }
            #endregion

            //Operator overloads
            #region
            //Unary
            public static Matrix operator -(Matrix a)
            { return a.Negate(); }

            //Binary
            public static Matrix operator *(Matrix a, Matrix b)
            { return a.Multiply(b); }
            public static Matrix operator *(Real a, Matrix b)
            { return b.ScalarMultiply(a); }
            public static Matrix operator +(Matrix a, Matrix b)
            { return a.Add(b); }
            public static Matrix operator -(Matrix a, Matrix b)
            { return a.Subtract(b); }
            public static bool operator ==(Matrix a, Matrix b)
            { return a.CompareTo(b); }
            public static bool operator !=(Matrix a, Matrix b)
            { return !a.CompareTo(b); }


            #endregion

            //Matrix methods returning matrices
            #region
            /// <summary>
            /// The Negative of this matrix
            /// </summary>
            /// <returns>The negative of this matrix</returns>
            private Matrix Negate()
            {
                //create an output matrix of the required size
                Matrix OutP = new Matrix(_m, _n);

                //for every row
                for (int i = 0; i < _m; i++)
                {
                    //for every column
                    for (int j = 0; j < _n; j++)
                    { OutP[i, j] = -this[i, j]; }
                }
                //casting back and forth to ensure value is copied by value and not by reference
                OutP.type = (MatType)(int)this.type;
                //return the output matrix
                return OutP;
            }

            /// <summary>
            /// Finds the Row-Echelon form of a matrix
            /// </summary>
            /// <param name="DetMult">The determinant multiplier</param>
            /// <returns>This matrix in Row-Echelon form</returns>
            public Matrix REF(out Real DetMult)
            {
                //an instance of T to issue factory methods
                Real instance = new Fraction();

                //output matrix
                Matrix OutP = (Matrix)this.Clone();

                //initially the multiplier is one
                DetMult = instance.FromInt(1);
                //the current column being processed
                int CurCol = 0;
                //for every row
                for (int m = 0; m < _m; m++)
                {
                    //the largest is initially zero and assumed to be in the mth columns
                    Real Largest = instance.FromInt(0);
                    int i_max = m;
                    //While the largest value is zero, keep proceeding through the matrix to find the column with a non zero element
                    while (Largest.IsZero() && (CurCol < _n))
                    {
                        //find the largest value in the current column
                        for (int n = m; n < _m; n++)
                        {
                            //if the value is larger than the current largest value
                            if (OutP[n, CurCol].Absolute() > Largest.Absolute())
                            {
                                //change the index and the stored value
                                i_max = n;
                                Largest = OutP[n, CurCol];
                            }
                        }
                        CurCol++;
                    }
                    //if the largest value is still zero at this point, then the matrix has been eliminated as far as possible so end the function
                    if (Largest.IsZero()) { return OutP; }
                    //if the largest value is not at the current row level
                    if (i_max != m)
                    {
                        //else, move the largest value to the row being processed
                        OutP.Swap(m, i_max);
                        //multiply the determinant by -1;
                        DetMult *= instance.FromInt(-1);
                    }
                    //Make the current row a unit row
                    OutP.MultiplyByScalar(m, Largest.Reciprocal());
                    DetMult *= Largest.Reciprocal();

                    //sanity check
                    if (!((OutP[m, m] == instance.FromInt(1) || (OutP[m, m] == instance.FromInt(0))))) { throw new Exception("Leading value of a row is not zero or a one"); }

                    //for every row after the current one
                    for (int i = (m + 1); i < _m; i++)
                    {
                        //if the first value is already zero, skip the rest of this loop
                        if (OutP[i, m].IsZero()) { continue; }
                        //subtract the current row to make the leading value nought
                        OutP.SubtractRow(i, m, OutP[i, m]);
                    }
                }
                //check the matrix type
                OutP.CheckMatType();
                return OutP;
            }

            /// <summary>
            /// Finds the Reduced Row-Echelon form of the matrix
            /// </summary>
            /// <param name="DetMult">The determinant multiplier</param>
            /// <returns>This matrix in Reduced Row-Echelon form</returns>
            public Matrix RREF(out Real DetMult)
            {
                //run a GaussJordanElimination to get it into row-echelon form
                Matrix OutP = REF(out DetMult);

                //for every row in reverse order
                for (int i = _m - 1; i >= 0; i--)
                {
                    //the first non-zero element in every row is one, so we only need to store the position
                    //initially I've set it to -1 since we can easily check if a non-zero value was found
                    int c = -1;

                    //find the first non-zero value in the row, values before the row number will be zero due to it being in row-echelon form
                    for (int j = i; j < _n; j++)
                    {
                        //find the first non-zero number and end the loop
                        if (OutP[i, j].Absolute() > 0) { c = j; break; }
                    }
                    //if the array index is still -1, then a non-zero element hasn't been found and the row is a zero row. Skip the row
                    if (c == -1) { break; }

                    //for every row above the current one, turn all values above the first non-zero value into zero
                    for (int j = 0; j < i; j++) { OutP.SubtractRow(j, i, OutP[j, c]); }
                }
                //check the matrix type
                OutP.CheckMatType();
                return OutP;
            }

            /// <summary>
            /// Returns the inverse of a matrix, returns a zero matrix if a matrix is singular
            /// </summary>
            /// <returns>The inverse of this matrix, the corresponding zero matrix if this matrix does not have an inverse</returns>
            public Matrix Inverse()
            {
                //if it isn't square, can't have an inverse
                if (!TestForFlag(MatType.Square))
                {
                    //Return zero matrix
                    return new Matrix(_m, _n);
                }

                //concatenate this matrix and an identity matrix of the same size and perfrom a Gauss Jordan elimination
                var temp = this.Concat(new Matrix(_m)).RREF(out var DetMult);

                //the left half should end up as an identity matrix, furthermore the Gaussian eliminiation method migrates zero rows to the bottom
                //this allows for an easy check to be made on whether the matrix is singular or not
                if (temp[_m - 1, _m - 1].IsZero())
                {
                    return new Matrix(_m, _n);
                }

                //split the matrix and return the inverse
                temp.Split(_m, out var Inverse);
                //check the matrix type
                Inverse.CheckMatType();
                return Inverse;
            }

            /// <summary>
            /// Create a shallow (reference) copy of this matrix
            /// </summary>
            /// <returns>A shallow copy of this matrix</returns>
            public Matrix Clone()
            {
                var OutP = new Matrix(_m, _n);
                for (int i = 0; i < _m; i++)
                {
                    for (int j = 0; j < _n; j++)
                    {
                        OutP[i, j] = this[i, j].Clone();
                    }
                }
                //copy over the matrix type
                OutP.type = (MatType)(int)this.type;
                return OutP;
            }

            /// <summary>
            /// Splits a matrix into two matrices with the first matrix having width "SplitLength"
            /// </summary>
            /// <param name="SplitLength">The number of columns in the first matrix</param>
            /// <param name="Second">The matrix formed by taking the columns from "SplitLength" to n from this matrix</param>
            /// <returns>The matrix formed by taking columns 0 to "SplitLength-1" from this Matrix </returns>
            public Matrix Split(int SplitLength, out Matrix Second)
            {

                if (SplitLength >= _n) { throw new Exception(String.Format("SplitLength, {0}, is greater than or equal to matrix length, {1}", SplitLength, _n), new ArgumentOutOfRangeException()); }
                if ((SplitLength == _n - 1) || (SplitLength == 1)) { throw new Exception(String.Format("Splitting the matrix with SplitLength {0} would result in a vector. Try using GetColumn", SplitLength), new ArgumentOutOfRangeException()); }
                //storage lists
                Matrix OutP1 = new Matrix(_m, SplitLength); Second = new Matrix(_m, _n - SplitLength);
                //first matrix
                for (int i = 0; i < _m; i++)
                {
                    for (int j = 0; j < SplitLength; j++)
                    { OutP1[i, j] = this[i, j].Clone(); }
                }
                //second matrix
                for (int i = 0; i < _m; i++)
                {
                    for (int j = SplitLength; j < _n; j++)
                    { Second[i, j - SplitLength] = this[i, j].Clone(); }
                }
                //check types
                OutP1.CheckMatType();
                Second.CheckMatType();
                return OutP1;
            }

            /// <summary>
            /// Find the transpose of a matrix. The transpose is in Row-Major form
            /// </summary>
            /// <returns>The transpose of this matrix in Row-Major form</returns>
            public Matrix Transpose()
            {
                //Output matrix
                Matrix OutP = new Matrix(_n, _m);
                //for each column
                for (int j = 0; j < _n; j++)
                {
                    //progress down a column
                    for (int i = 0; i < _m; i++)
                    {
                        OutP[j, i] = this[i, j].Clone();
                    }
                }
                //check type
                OutP.CheckMatType();
                return OutP;
            }

            /// <summary>
            /// Insert a vector into matrix
            /// </summary>
            /// <param name="other">The vector to be inserted</param>
            /// <param name="index">The zero based index specifying where to insert the vector after. A value of -1 indicates the last index for the relevant dimension</param>
            /// <param name="RowMode">Whether to insert the vector as a row or a column</param>
            /// <returns></returns>
            public Matrix Insert(Vector other, int index = -1, bool RowMode = true)
            {
                throw new NotImplementedException();
            }
            #endregion

            //Matrix methods returning non-matrices
            #region
            /// <summary>
            /// Calculate the determinant of a matrix
            /// </summary>
            /// <returns>The determinant of this matrix</returns>
            public Real Determinant()
            {
                //non-square matrices don't have determinants
                if (!TestForFlag(MatType.Square)) { throw new Exception("Matrix is not square", new ArgumentException()); }

                //perform a Gaussian elimination
                Matrix tmp = this.REF(out Real DetMult);

                //The leading diagonal should all be either one or zero
                //If there is a zero, then the product of the diagonals is zero and so the determinant is zero
                for (int i = 0; i < _m; i++)
                { if (tmp[i, i].IsZero()) { return new Fraction().FromInt(0); } }

                //otherwise return the reciprocal of the determinant multiplier which is the determinant
                return (DetMult.Reciprocal());

            }

            /// <summary>
            /// Get the matrix as a string which satisfies the conversion criteria
            /// </summary>
            /// <returns>This matrix as a reconvertible string</returns>
            public string ToParseable()
            {
                //Initially the matrix will begin with an open bracket
                string OutP = "[";
                for (int i = 0; i < _m; i++)
                {
                    //for every row add an open bracket
                    OutP += "[";
                    for (int j = 0; j < _n; j++)
                    {
                        //add every element and a semicolon
                        OutP += this[i, j].ToString();
                        OutP += ";";
                    }
                    //and finish with a close bracket and a semicolon
                    OutP += "];";
                }
                //finish with a close bracket
                OutP += "]";

                return OutP;
            }

            /// <summary>
            /// Get the matrix as a multiline string
            /// </summary>
            /// <returns>This matrix as a string that spans several lines</returns>
            public override string ToString()
            {
                //stores the values as strings
                string[,] TempArr = new string[_m, _n];
                int LongestLength = 0;
                string OutP = "";

                for (int i = 0; i < _m; i++)
                {
                    for (int j = 0; j < _n; j++)
                    {
                        TempArr[i, j] = this[i, j].ToString();
                        //check for whether the string is longer than the current longest accepted length
                        LongestLength = (TempArr[i, j].Length > LongestLength) ? TempArr[i, j].Length : LongestLength;
                    }
                }

                for (int i = 0; i < _m; i++)
                {
                    for (int j = 0; j < _n - 1; j++)
                    {
                        //pad and add to the end of the string
                        OutP += (TempArr[i, j] + ',').PadRight(LongestLength + 2);
                    }
                    OutP += TempArr[i, _n - 1];
                    OutP += System.Environment.NewLine;
                }
                return OutP;
            }

            /// <summary>
            /// Get a string representing the matrix type of a matrix
            /// </summary>
            /// <returns>The matrix type of this matrix</returns>
            public string GetMatType()
            { return type.ToString(); }
            #endregion
        }

        /// <summary>
        /// A vector of type T where T is derived from Real
        /// </summary>
        /// <typeparam name="T">A type that is derived from Real (Frac|Floating)</typeparam>
        public class Vector
        {
            //Vector specific variables and types
            #region
            /// <summary>
            /// storage for the values of the vector
            /// </summary>
            private Real[] storage;

            /// <summary>
            /// The dimension of the vector
            /// </summary>
            private int _n;
            #endregion

            //Constructors
            #region
            /// <summary>
            /// Construct a vector from an array of Real values
            /// </summary>
            /// <param name="Values">The values to form the matrix from</param>
            public Vector(Real[] Values)
            {
                storage = Values;
                _n = storage.Length;
            }

            /// <summary>
            /// Construct a zero vector of dimension "n"
            /// </summary>
            /// <param name="n">The dimension of the new vector</param>
            public Vector(int n)
            {
                _n = n;
                Real instance = new Fraction();
                storage = new Real[_n];
                for (int i = 0; i < _n; i++)
                { storage[i] = instance.FromInt(0); }
            }

            /// <summary>
            /// Try to parse a string into a vector
            /// </summary>
            /// <param name="phrase">String to be parsed</param>
            /// <param name="OutP">Stores the output value if there is one</param>
            /// <returns>A boolean indicating whether "Phrase" is the valid </returns>
            public static bool TryParse(string phrase, out Vector OutP)
            {
                string tmpstr = "";
                bool elemEnd = false;
                Real tmpVal = new Fraction();
                var items = new List<Real>();
                OutP = new Vector(1);

                for (int i = 0; i < phrase.Length; i++)
                {
                    switch (phrase[i])
                    {
                        //first character is an open bracket
                        case '[':
                            if (i != 0) { return false; }
                            break;
                        //close bracket is the last character, all items should have ended by then
                        case ']':
                            if ((!elemEnd) || (i != phrase.Length - 1)) { return false; }
                            break;
                        //try and parse the element that has just ended
                        case ';':
                            elemEnd = true;
                            if (!Real.TryParse(tmpstr, out tmpVal)) { return false; }
                            items.Add(tmpVal);
                            tmpstr = "";
                            tmpVal = new Fraction();
                            break;
                        //add it to the accumulator
                        default:
                            elemEnd = false;
                            tmpstr += phrase[i];
                            break;
                    }
                }
                OutP = new Vector(items.ToArray());
                return true;

            }
            #endregion

            //Accessors
            #region
            /// <summary>
            /// Get the ith value of the vector
            /// </summary>
            /// <param name="i">The index of the value to be retrieved</param>
            /// <returns>The value at index i</returns>
            public Real this[int i]
            {
                get
                {
                    if (i < 0 || i >= _n) { throw new Exception("Index out of range", new ArgumentOutOfRangeException()); }
                    return storage[i];
                }
                internal set
                {
                    if (i < 0 || i >= _n) { throw new Exception("Index out of range", new ArgumentOutOfRangeException()); }
                    storage[i] = value;
                }
            }

            /// <summary>
            /// Get the dimension of the vector
            /// </summary>
            /// <returns>An intger representing the dimension of the vector</returns>
            public int Height
            { get { return _n; } }
            #endregion

            //Vector Operations
            #region
            /// <summary>
            /// Vector addition
            /// </summary>
            /// <param name="other">Vector to be to added</param>
            /// <returns>The sum of this vector and "other"</returns>
            private Vector Add(Vector other)
            {
                if (_n != other._n)
                { throw new Exception(String.Format("Dimension mismatch, first vector is of dimension {0}, second is of dimension {1}", _n, other._n), new ArgumentException()); }
                //Output vector
                Vector OutP = new Vector(_n);
                //perform addition
                for (int i = 0; i < _n; i++)
                { OutP[i] = this[i] + other[i]; }
                return OutP;
            }

            /// <summary>
            /// Vector subtraction
            /// </summary>
            /// <param name="other">Vector to be subtracted</param>
            /// <returns>This vector minus "other"</returns>
            private Vector Subtract(Vector other)
            {
                if (_n != other._n)
                { throw new Exception(String.Format("Dimension mismatch, first vector is of dimension {0}, second is of dimension {1}", _n, other._n), new ArgumentException()); }
                //Output vector
                Vector OutP = new Vector(_n);
                //perform addition
                for (int i = 0; i < _n; i++)
                { OutP[i] = this[i] - other[i]; }
                return OutP;
            }

            /// <summary>
            /// Calculate the dot product of this vector and another
            /// </summary>
            /// <param name="other">The vector to be multiplied by</param>
            /// <returns>The dot product of this vector and "other" as a Real</returns>
            private Real Dot(Vector other)
            {
                if (_n != other._n)
                { throw new Exception(String.Format("Dimension mismatch, first vector is of dimension {0}, second is of dimension {1}", _n, other._n), new ArgumentException()); }

                //Output variable
                Real OutP = new Fraction().FromInt(0);
                //Sum of products
                for (int i = 0; i < _n; i++)
                { OutP += this[i] * other[i]; }
                //return value
                return OutP;
            }

            /// <summary>
            /// Calcualate the cross product of this vector and another
            /// </summary>
            /// <param name="other">The vector to be right multiplied by</param>
            /// <returns>The cross product of this vector and "other"</returns>
            private Vector Cross(Vector other)
            {
                if (_n != 3)
                { throw new Exception(String.Format("Cross product is not defined for vectors of dimension != 3", _n, other._n), new ArgumentException()); }
                if (_n != other._n)
                { throw new Exception(String.Format("Dimension mismatch, first vector is of dimension {0}, second is of dimension {1}", _n, other._n), new ArgumentException()); }
                //output vector
                Vector OutP = new Vector(3);
                //calculate vector product
                OutP[0] = this[1] * other[2] - this[2] * other[1];
                OutP[1] = this[2] * other[0] - this[0] * other[2];
                OutP[2] = this[0] * other[1] - this[1] * other[0];
                //return the value calculated
                return OutP;
            }

            /// <summary>
            /// Calculate the vector produced when this vector is right multiplied by a matrix
            /// </summary>
            /// <param name="RightMatrix">The matrix to right multiply this vector by</param>
            /// <returns>This vector right multiplied by RightMatrix</returns>
            private Vector RightProduct(Matrix RightMatrix)
            {
                if (_n != RightMatrix.Height)
                { throw new Exception(String.Format("Vector is of dimension {0} and matrix is of input dimension {1}", _n, RightMatrix.Height), new ArgumentException()); }

                Vector OutP = new Vector(RightMatrix.Length);

                for (int i = 0; i < OutP._n; i++)
                {
                    for (int j = 0; j < _n; j++)
                    {
                        OutP[i] = OutP[i] + this[j] * RightMatrix[j, i];
                    }
                }
                return OutP;
            }

            /// <summary>
            /// Calculate the vector produced when this vector is left multiplied by a matrix
            /// </summary>
            /// <param name="LeftMatrix">The matrix to left multiply this vector by</param>
            /// <returns>This vector left multiplied by LeftMatrix</returns>
            private Vector LeftProduct(Matrix LeftMatrix)
            {
                if (_n != LeftMatrix.Length)
                { throw new Exception(String.Format("Vector is of dimension {0} and matrix is of input dimension {1}", _n, LeftMatrix.Length), new ArgumentException()); }

                Vector OutP = new Vector(LeftMatrix.Height);

                for (int i = 0; i < OutP._n; i++)
                {
                    for (int j = 0; j < _n; j++)
                    {
                        OutP[i] = OutP[i] + this[j] * LeftMatrix[i, j];
                    }
                }
                return OutP;
            }
            
            
            /// <summary>
            /// Returns a value indicating whether this instance is equal to a specific object
            /// </summary>
            /// <param name="obj">The object to test equality against</param>
            /// <returns>A boolean value indicating whether this instance is equal to "obj"</returns>
            public override bool Equals(object obj)
            {
                var tmp = obj as Vector;
                if (tmp is null || tmp._n != _n) { return false; }
                for (int i = 0; i < _n; i++)
                {
                    if (tmp[i] != this[i]) { return false; }
                }
                return true;
            }
                        
            /// <summary>
            /// Returns the hash code for this instance
            /// </summary>
            /// <returns>The hash code for this instance</returns>
            public override int GetHashCode()
            {
                return storage.GetHashCode();
            }
            #endregion

            //Operator overloads
            #region
            public static Vector operator -(Vector a)
            { return a.Negate(); }
            public static Vector operator *(Vector a, Vector b)
            { return a.Cross(b); }
            public static Vector operator *(Vector a, Matrix B)
            { return a.RightProduct(B); }
            public static Vector operator *(Matrix A, Vector b)
            { return b.LeftProduct(A); }
            public static Real operator %(Vector a, Vector b)
            { return a.Dot(b); }
            public static Vector operator +(Vector a, Vector b)
            { return a.Add(b); }
            public static Vector operator -(Vector a, Vector b)
            { return a.Subtract(b); }
            public static bool operator ==(Vector a, Vector b)
            { return a.Equals(b); }
            public static bool operator !=(Vector a, Vector b)
            { return !a.Equals(b); }
            #endregion

            //Vector methods returning vectors
            #region
            private Vector Negate()
            {
                Vector OutP = new Vector(_n);
                for (int i = 0; i < _n; i++)
                { OutP[i] = -this[i]; }
                return OutP;
            }
            #endregion

            //Vector methods returning non-vectors
            #region
            /// <summary>
            /// Find the String representation of the vector
            /// </summary>
            /// <returns>Returns a string representation of the vector</returns>
            public override string ToString()
            {
                string OutP = "[";
                for (int i = 0; i < _n; i++)
                { OutP += this[i].ToString(); OutP += ";"; }
                OutP += "]";
                return OutP;
            }
            #endregion
        }

    }
}