Matrix Class supported Generic in C#

Matrix.cs

using System;
using GenericOperator;

namespace NeuralNetwork.Matrix
{
    public static class Matrix
    {
        public static Matrix<int> Zero(int n)
        {
            if (n <= 0) throw new ArgumentException();
                return new Matrix<int>(n, n);
        }
        
            public static Matrix<int> Identity(int n)
            {
                if (n <= 0) throw new ArgumentException();
                    var ret = new Matrix<int>(n, n);
                        for (int i = 0; i < n; ++i)
                            ret[i, i] = 1;
                                return ret;
            }
        
            public static Matrix<double> ZeroD(int n)
            {
                if (n <= 0) throw new ArgumentException();
                    return new Matrix<double>(n, n);
            }
        
            public static Matrix<double> IdentityD(int n)
            {
                if (n <= 0) throw new ArgumentException();
                    var ret = new Matrix<double>(n, n);
                        for (int i = 0; i < n; ++i)
                            ret[i, i] = 1;
                                return ret;
            }
        
            public static int Trace(Matrix<int> mat)
            {
                if (mat.Row != mat.Col) throw new ArithmeticException();
                    var ret = 0;
                        for (int i = 0; i < mat.Row; ++i)
                            ret += mat[i, i];
                                return ret;
            }
        
            public static double Trace(Matrix<double> mat)
            {
                if (mat.Row != mat.Col) throw new ArithmeticException();
                    var ret = 0.0;
                        for (int i = 0; i < mat.Row; ++i)
                            ret += mat[i, i];
                                return ret;
            }
        
            public static Matrix<double> Inverse(Matrix<double> mat)
            {
                if (mat.Col != mat.Row) throw new ArithmeticException();
                    if (mat.Col == 2) return Inverse2(mat);
                        
                            var N = mat.Row;
                            var ret = new Matrix<double>(N, N);
                            var identify = IdentityD(N);
                            
                            for (int i = 0; i < N; ++i)
                            {
                                if (Math.Abs(mat[i, i]) < double.Epsilon) return null;
                                    var inv = 1.0 / mat[i, i];
                                        
                                        for (int j = 0; j < N; j++)
                                        {
                                            mat[i, j] *= inv;
                                            ret[i, j] *= inv;
                                        }
                                
                                    for (int j = 0; j < N; j++)
                                    {
                                        if (i != j)
                                        {
                                            inv = mat[j, i];
                                            for (int k = 0; k < N; k++)
                                            {
                                                mat[j, k] -= mat[i, k] * inv;
                                                ret[k, j] -= ret[i, k] * inv;
                                            }
                                        }
                                    }
                            }
                
                    return ret;
            }
        
            public static Matrix<double> Inverse2(Matrix<int> mat)
            {
                var det = Determinant2(mat);
                    if (det == 0) return null;
                        return new Matrix<double>(2, 2, mat[1, 1] / det, -mat[0, 1] / det, -mat[1, 0] / det, mat[0, 0] / det);
            }
        
            public static Matrix<double> Inverse2(Matrix<double> mat)
            {
                var det = Determinant2(mat);
                    if (Math.Abs(det) < double.Epsilon) return null;
                        return new Matrix<double>(2, 2, mat[1, 1] / det, -mat[0, 1] / det, -mat[1, 0] / det, mat[0, 0] / det);
            }
        
            public static double Determinate(Matrix<int> mat)
            {
                if (mat.Col != mat.Row) throw new ArithmeticException();
                    if (mat.Col == 2) return Determinant2(mat);
                        if (mat.Col == 3) return Determinant3(mat);
                            
                                Matrix<double> L, U;
                                LUDecomposition(mat, out L, out U);
                                double ret = 1;
                                for (int i = 0; i < U.Row; ++i)
                                    ret *= U[i, i];
                                        return ret;
            }
        
            public static double Determinate(Matrix<double> mat)
            {
                if (mat.Col != mat.Row) throw new ArithmeticException();
                    if (mat.Col == 2) return Determinant2(mat);
                        if (mat.Col == 3) return Determinant3(mat);
                            
                                Matrix<double> L, U;
                                LUDecomposition(mat, out L, out U);
                                double ret = 1;
                                for (int i = 0; i < U.Row; ++i)
                                    ret *= U[i, i];
                                        return ret;
            }
        
            public static int Determinant2(Matrix<int> mat)
            {
                if (mat.Col != 2 || mat.Row != 2) throw new ArithmeticException();
                    return mat[0, 0] * mat[1, 1] - mat[0, 1] * mat[1, 0];
            }
        
            public static double Determinant2(Matrix<double> mat)
            {
                if (mat.Col != 2 || mat.Row != 2) throw new ArithmeticException();
                    return mat[0, 0] * mat[1, 1] - mat[0, 1] * mat[1, 0];
            }
        
            public static int Determinant3(Matrix<int> mat)
            {
                if (mat.Col != 3 || mat.Row != 3) throw new ArithmeticException();
                    return mat[0, 0] * mat[1, 1] * mat[2, 2] + mat[0, 1] * mat[1, 2] * mat[2, 0] + mat[0, 2] * mat[1, 0] * mat[2, 1]
                        - mat[0, 2] * mat[1, 1] * mat[2, 0] - mat[1, 2] * mat[2, 1] * mat[0, 0] - mat[2, 2] * mat[0, 1] * mat[1, 0];
            }
        
            public static double Determinant3(Matrix<double> mat)
            {
                if (mat.Col != 3 || mat.Row != 3) throw new ArithmeticException();
                    return mat[0, 0] * mat[1, 1] * mat[2, 2] + mat[0, 1] * mat[1, 2] * mat[2, 0] + mat[0, 2] * mat[1, 0] * mat[2, 1]
                        - mat[0, 2] * mat[1, 1] * mat[2, 0] - mat[1, 2] * mat[2, 1] * mat[0, 0] - mat[2, 2] * mat[0, 1] * mat[1, 0];
            }
        
            public static Matrix<int> Pow(Matrix<int> mat, int n)
            {
                if (mat.Col != mat.Row) throw new ArithmeticException();
                    if (n < 0) throw new ArithmeticException();
                        if (n == 1) return mat;
                            if (n % 2 == 0)
                            {
                                return Pow(mat * mat, n / 2);
                            }
                            else
                            {
                                return Pow(mat, n - 1) * mat;
                            }
            }
        
            public static Matrix<double> Pow(Matrix<double> mat, int n)
            {
                if (mat.Col != mat.Row) throw new ArithmeticException();
                    if (n < 0) throw new ArithmeticException();
                        if (n == 1) return mat;
                            if (n % 2 == 0)
                            {
                                return Pow(mat * mat, n / 2);
                            }
                            else
                            {
                                return Pow(mat, n - 1) * mat;
                            }
            }
        
            public static Matrix<int> Transposition(Matrix<int> mat)
            {
                var ret = new Matrix<int>(mat.Col, mat.Row);
                    for (int i = 0; i < mat.Row; ++i)
                    {
                        for (int j = 0; j < mat.Col; ++j)
                        {
                            ret[j, i] = mat[i, j];
                        }
                    }
                return ret;
            }
        
            public static Matrix<double> Transposition(Matrix<double> mat)
            {
                var ret = new Matrix<double>(mat.Col, mat.Row);
                    for (int i = 0; i < mat.Row; ++i)
                    {
                        for (int j = 0; j < mat.Col; ++j)
                        {
                            ret[j, i] = mat[i, j];
                        }
                    }
                return ret;
            }
        
            public static Matrix<double> ToDouble(Matrix<int> mat)
            {
                var ret = new Matrix<double>(mat.Row, mat.Col);
                    for (int i = 0; i < mat.Row; ++i)
                    {
                        for (int j = 0; j < mat.Col; ++j)
                        {
                            ret[i, j] = (double)mat[i, j];
                        }
                    }
                return ret;
            }
        
            public static void LUDecomposition(Matrix<int> mat, out Matrix<double> L, out Matrix<double> U)
            {
                if (mat.Col != mat.Row) throw new ArithmeticException();
                    
                        var N = mat.Col;
                        
                        L = IdentityD(N);
                        U = ZeroD(N);
                        
                        var sum = 0.0;
                        for (int i = 0; i < N; i++)
                        {
                            for (int j = 0; j < N; j++)
                            {
                                if (i > j)
                                {
                                    sum = 0.0;
                                    for (int k = 0; k < j; k++)
                                    {
                                        sum += L[i, k] * U[k, j];
                                    }
                                    L[i, j] = (mat[i, j] - sum) / U[j, j];
                                }
                                else
                                {
                                    sum = 0.0;
                                    for (int k = 0; k < i; k++)
                                    {
                                        sum += L[i, k] * U[k, j];
                                    }
                                    U[i, j] = mat[i, j] - sum;
                                }
                            }
                        }
            }
        
            public static void LUDecomposition(Matrix<double> mat, out Matrix<double> L, out Matrix<double> U)
            {
                if (mat.Col != mat.Row) throw new ArithmeticException();
                    
                        var N = mat.Row;
                        
                        L = IdentityD(N);
                        U = ZeroD(N);
                        
                        var sum = 0.0;
                        for (int i = 0; i < N; i++)
                        {
                            for (int j = 0; j < N; j++)
                            {
                                if (i > j)
                                {
                                    sum = 0.0;
                                        for (int k = 0; k < j; k++)
                                        {
                                            sum += L[i, k] * U[k, j];
                                        }
                                    L[i, j] = (mat[i, j] - sum) / U[j, j];
                                }
                                else
                                {
                                    sum = 0.0;
                                        for (int k = 0; k < i; k++)
                                        {
                                            sum += L[i, k] * U[k, j];
                                        }
                                    U[i, j] = mat[i, j] - sum;
                                }
                            }
                        }
            }
        
            public static Matrix<double> LForwardsubs(Matrix<double> L, Matrix<double> b)
            {
                if (L.Row != b.Row && b.Col != 1) throw new ArithmeticException();
                    var N = L.Row;
                        
                        var ret = new Matrix<double>(N, 1);
                        for (int i = 0; i < N; ++i)
                            ret[i, 0] = b[i, 0];
                                
                                for (int i = 0; i < N; i++)
                                {
                                    ret[i, 0] /= L[i, i];
                                    for (int j = i + 1; j < N; j++)
                                    {
                                        ret[j, 0] -= ret[i, 0] * L[j, i];
                                    }
                                }
                
                    return ret;
            }
        
            public static Matrix<double> UBackwardsubs(Matrix<double> U, Matrix<double> y)
            {
                if (U.Row != y.Row && y.Col != 1) throw new ArithmeticException();
                    var N = U.Row;
                        
                        var ret = new Matrix<double>(N, 1);
                for (int i = 0; i < N; i++)
                {
                    ret[i, 0] = y[i, 0];
                }
                for (int i = N - 1; i >= 0; i--)
                {
                    ret[i, 0] /= U[i, i];
                    for (int j = i - 1; j >= 0; j--)
                    {
                        ret[j, 0] -= ret[i, 0] * U[j, i];
                    }
                }
                
                    return ret;
            }
        
            public static Matrix<double> SolveSimultaneousEquations(Matrix<double> mat, Matrix<double> b)
            {
                var N = mat.Row;
                    Matrix<double> L, U;
                    L = IdentityD(N);
                    U = ZeroD(N);
                    
                    var y = LForwardsubs(L, b);
                    return UBackwardsubs(U, y);
            }
    }
    
        /// <summary>
        /// ジェネリック型 T の行列クラス
        /// </summary>
        /// <typeparam name="T">対象となる型。</typeparam>
        public class Matrix<T> where T : IEquatable<T>
        {
            private T[,] _Matrix;
                public int Row
                {
                    private set;
                        get;
                }
            public int Col
            {
                private set;
                    get;
            }
            
                public Matrix(T[,] mat)
                {
                    if (mat == null) throw new ArgumentNullException();
                        
                            _Matrix = mat;
                            Row = mat.GetLength(0);
                            Col = mat.GetLength(1);
                }
            
                /// <summary>
                /// Matrix [row * col]
                /// </summary>
                /// <param name="row">Row</param>
                /// <param name="col">Col</param>
                public Matrix(int row, int col)
                {
                    if (row <= 0 || col <= 0) throw new ArgumentException();
                        
                            _Matrix = new T[row, col];
                            Row = row;
                            Col = col;
                }
            
                /// <summary>
                /// Matrix [row * col]
                /// </summary>
                /// <param name="row">Row</param>
                /// <param name="col">Col</param>
                /// <param name="args">Matrix elements</param>
                public Matrix(int row, int col, params T[] args)
                {
                    if (row <= 0 || col <= 0) throw new ArgumentException();
                        if (args.Length != row * col) throw new ArgumentException();
                            
                                _Matrix = new T[row, col];
                                Row = row;
                                Col = col;
                                
                                for (int i = 0; i < row; ++i)
                                {
                                    for (int j = 0; j < col; ++j)
                                    {
                                        _Matrix[i, j] = args[i * col + j];
                                    }
                                }
                }

            public T this[int i, int j]
            {
                set
                {
                    if (i < 0 || i >= Row) throw new IndexOutOfRangeException();
                        if (j < 0 || j >= Col) throw new IndexOutOfRangeException();
                            _Matrix[i, j] = value;
                }
                get
                {
                    if (i < 0 || i >= Row) throw new IndexOutOfRangeException();
                        if (j < 0 || j >= Col) throw new IndexOutOfRangeException();
                            return _Matrix[i, j];
                }
            }
            
                public Matrix<T> GetRow(int r)
                {
                    if (r < 0 || r >= Row) throw new IndexOutOfRangeException();
                        var ret = new Matrix<T>(1, Col);
                            
                            for (int i = 0; i < Col; ++i)
                            {
                                ret[0, i] = this[0, i];
                            }
                    return ret;
                }
            
                public Matrix<T> GetCol(int c)
                {
                    if (c < 0 || c >= Col) throw new IndexOutOfRangeException();
                        var ret = new Matrix<T>(Row, 1);
                            
                            for (int i = 0; i < Row; ++i)
                            {
                                ret[i, 0] = this[i, 0];
                            }
                    return ret;
                }
            
                public Matrix<T> GetLowerTriangularMatrix()
                {
                    if (Row != Col) throw new ArithmeticException();
                        var ret = new Matrix<T>(Row, Col);
                            
                            for (int i = 0; i < Row; ++i)
                            {
                                for (int j = 0; j <= i; ++j)
                                {
                                    ret[i, j] = this[i, j];
                                }
                            }
                    
                        return ret;
                }
            
                
                public Matrix<T> GetUpperTriangularMatrix()
                {
                    if (Row != Col) throw new ArithmeticException();
                        var ret = new Matrix<T>(Row, Col);
                            
                            for (int i = 0; i < Row; ++i)
                            {
                                for (int j = i; j < Col; ++j)
                                {
                                    ret[i, j] = this[i, j];
                                }
                            }
                    
                        return ret;
                }
            
                public static Matrix<T> operator +(Matrix<T> a, Matrix<T> b)
                {
                    //if (a == null || b == null) throw new ArgumentNullException();
                    if (a.Row != b.Row || a.Col != b.Col) throw new ArithmeticException();
                        var row = a.Row;
                            var col = a.Col;
                            var ret = new Matrix<T>(row, col);
                            
                            for (int i = 0; i < row; ++i)
                            {
                                for (int j = 0; j < col; ++j)
                                {
                                    ret[i, j] = ArithmeticOperator<T>.Add(a[i, j], b[i, j]);
                                }
                            }
                    
                        return ret;
                }
            
                public static Matrix<T> operator -(Matrix<T> a, Matrix<T> b)
                {
                    if (a.Row != b.Row || a.Col != b.Col) throw new ArithmeticException();
                        var row = a.Row;
                            var col = a.Col;
                            var ret = new Matrix<T>(row, col);
                            
                            for (int i = 0; i < row; ++i)
                            {
                                for (int j = 0; j < col; ++j)
                                {
                                    ret[i, j] = ArithmeticOperator<T>.Subtract(a[i, j], b[i, j]);
                                }
                            }
                    return ret;
                }
            
                public static Matrix<T> operator *(T x, Matrix<T> mat)
                {
                    var ret = new Matrix<T>(mat.Row, mat.Col);
                    for (int i = 0; i < mat.Row; ++i)
                    {
                        for (int j = 0; j < mat.Col; ++j)
                        {
                            ret[i, j] = ArithmeticOperator<T>.Multiply(x, mat[i, j]);
                        }
                    }

                    return ret;
                }
            
                public static Matrix<T> operator *(Matrix<T> mat, T x)
                {
                    var ret = new Matrix<T>(mat.Row, mat.Col);
                        for (int i = 0; i < mat.Row; ++i)
                        {
                            for (int j = 0; j < mat.Col; ++j)
                            {
                                ret[i, j] = ArithmeticOperator<T>.Multiply(x, mat[i, j]);
                            }
                        }
                    
                        return ret;
                }
            
                public static Matrix<T> operator *(Matrix<T> a, Matrix<T> b)
                {
                    if (a.Col != b.Row) throw new ArithmeticException();
                        var row = a.Row;
                            var col = b.Col;
                            var m = a.Col;
                            var ret = new Matrix<T>(row, col);
                            
                            for (int i = 0; i < row; ++i)
                            {
                                for (int j = 0; j < col; ++j)
                                {
                                    for (int k = 0; k < m; ++k)
                                    {
                                        ret[i, j] = ArithmeticOperator<T>.Add(ArithmeticOperator<T>.Multiply(a[i, k], b[k, j]), ret[i, j]);
                                    }
                                }
                            }
                    
                        return ret;
                }
            
                public static bool operator ==(Matrix<T> a, Matrix<T> b)
                {
                    if (a.Row != b.Row || a.Col != b.Col) return false;
                        
                            for (int i = 0; i < a.Row; ++i)
                            {
                                for (int j = 0; j < a.Col; ++j)
                                {
                                    if (!a[i, j].Equals(b[i, j])) return false;
                                }
                            }
                    
                        return true;
                }
            
                public static bool operator !=(Matrix<T> a, Matrix<T> b)
                {
                    return !(a == b);
                }
            
                public override bool Equals(object obj)
                {
                    var mat = obj as Matrix<T>;
                    return this == mat;
                }
            
                public override int GetHashCode()
                {
                    return base.GetHashCode();
                }
        }
}