逆行列を取得するInverse()の実装等，基本機能の拡充(#26510, #26511)

diff --git a/Karinto/Numeric/Matrix.cs b/Karinto/Numeric/Matrix.cs
index bbc622b..64c4e44 100755 (executable)
--- a/Karinto/Numeric/Matrix.cs
+++ b/Karinto/Numeric/Matrix.cs
@@ -18,7 +18,6 @@ namespace Karinto.Numeric
     public class Matrix : Matrix<double>\r
     {\r
         #region constructors\r
-\r
         public Matrix(int row, int column)\r
             : base(row, column)\r
         {\r
@@ -33,6 +32,7 @@ namespace Karinto.Numeric
             : base(elements)\r
         {\r
         }\r
+\r
         #endregion\r
 \r
         #region operators\r
@@ -43,6 +43,11 @@ namespace Karinto.Numeric
             return (Matrix)((Matrix<double>)left + right);\r
         }\r
 \r
+        public static Matrix operator -(\r
+                                    Matrix left, Matrix<double> right)\r
+        {\r
+            return (Matrix)((Matrix<double>)left - right);\r
+        }\r
 \r
         public static Matrix operator *(Matrix left, Matrix<double> right)\r
         {\r
@@ -59,7 +64,31 @@ namespace Karinto.Numeric
             return (Matrix)((Matrix<double>)left * right);\r
         }\r
 \r
+        public static Matrix operator /(Matrix left, Matrix<double> right)\r
+        {\r
+            return (Matrix)(left * right.Inverse());\r
+        }\r
+\r
+        #endregion\r
+\r
+        #region public methods\r
+        public override Object Clone()\r
+        {\r
+            return new Matrix(e);\r
+        }\r
+\r
+        public new Matrix Transpose()\r
+        {\r
+            return (Matrix)((Matrix<double>)this).Transpose();\r
+        }\r
+\r
+        public new Matrix Inverse()\r
+        {\r
+            return (Matrix)((Matrix<double>)this).Inverse();\r
+        }\r
+\r
         #endregion\r
+\r
         protected override Matrix<double> CreateSkeleton(int row, int column)\r
         {\r
             return new Matrix(row, column);\r
@@ -131,7 +160,6 @@ namespace Karinto.Numeric
             col = elements.GetLength(1);\r
             e = (T[,])elements.Clone();\r
         }\r
-\r
         #endregion\r
 \r
         #region operators\r
@@ -242,6 +270,11 @@ namespace Karinto.Numeric
             }\r
             return m;\r
         }\r
+\r
+        public static Matrix<T> operator /(Matrix<T> left, Matrix<T> right)\r
+        {\r
+            return left * right.Inverse();\r
+        }\r
         #endregion\r
 \r
         #region properties\r
@@ -271,6 +304,7 @@ namespace Karinto.Numeric
         /// <summary>\r
         /// 列ベクトルの配列\r
         /// </summary>\r
+        /// <remarks>取得した配列，要素に破壊的操作をしてはならない</remarks>\r
         public virtual T[][] Rows\r
         {\r
             get\r
@@ -287,6 +321,7 @@ namespace Karinto.Numeric
         /// <summary>\r
         /// 行ベクトルの配列\r
         /// </summary>\r
+        /// <remarks>取得した配列，要素に破壊的操作をしてはならない</remarks>\r
         public virtual T[][] Columns\r
         {\r
             get\r
@@ -308,19 +343,16 @@ namespace Karinto.Numeric
             get\r
             {\r
                 CheckSquare(this);\r
-                T one = Operator<T>.One;\r
                 T zero = Operator<T>.Zero;\r
-                int l = row;\r
-                T[][] rows = new T[l][];\r
-                for (int r = 0; r < l; ++r)\r
+                T[][] rows = new T[row][];\r
+                for (int r = 0; r < row; ++r)\r
                 {\r
                     rows[r] = (T[])this[r].Clone();\r
                 }\r
-                int last = l - 1;\r
                 bool sign = false;\r
                 T pivot = Operator<T>.One;\r
                 T prevPivot = pivot;\r
-                for (int i = 0; i < l; ++i)\r
+                for (int i = 0; i < col; ++i)\r
                 {\r
                     prevPivot = pivot;\r
                     pivot = rows[i][i];\r
@@ -329,27 +361,25 @@ namespace Karinto.Numeric
                         int r = 0;\r
                         while (rows[r][i].Equals(zero))\r
                         {\r
-                            if (++r >= l)\r
+                            if (++r >= row)\r
                             {\r
                                 r = 0;\r
                                 break;\r
                             }\r
                         }\r
-                        T[] t = rows[r];\r
-                        rows[r] = rows[i];\r
-                        rows[i] = t;\r
+                        SwapRows(rows, r, i);\r
                         pivot = rows[i][i];\r
                         sign = !sign;\r
                     }\r
-                    for (int j = i + 1; j < l; ++j)\r
+                    for (int m = i + 1; m < row; ++m)\r
                     {\r
-                        T[] rowVec = rows[j];\r
-                        for (int k = i + 1; k < l; ++k)\r
+                        T[] rowVec = rows[m];\r
+                        for (int n = i + 1; n < col; ++n)\r
                         {\r
-                            T x = Operator<T>.Mul(pivot, rowVec[k]);\r
-                            T y = Operator<T>.Mul(rowVec[i], rows[i][k]);\r
+                            T x = Operator<T>.Mul(pivot, rowVec[n]);\r
+                            T y = Operator<T>.Mul(rowVec[i], rows[i][n]);\r
                             T z = Operator<T>.Sub(x, y);\r
-                            rowVec[k] = Operator<T>.Div(z, prevPivot);\r
+                            rowVec[n] = Operator<T>.Div(z, prevPivot);\r
                         }\r
                     }\r
                 }\r
@@ -364,7 +394,50 @@ namespace Karinto.Numeric
         {\r
             get\r
             {\r
-                throw new NotImplementedException();\r
+                T zero = Operator<T>.Zero;\r
+                T[][] rows = new T[row][];\r
+                for (int r = 0; r < row; ++r)\r
+                {\r
+                    rows[r] = (T[])this[r].Clone();\r
+                }\r
+                T pivot = Operator<T>.One;\r
+                T prevPivotR = pivot;\r
+                int pivotIndex = 0;\r
+                for (int i = 0; i < col; ++i)\r
+                {\r
+                    int r = pivotIndex;\r
+                    while (rows[r][i].Equals(zero))\r
+                    {\r
+                        if (++r >= row)\r
+                        {\r
+                            break;\r
+                        }\r
+                    }\r
+                    if (r == row)\r
+                    {\r
+                        continue;\r
+                    }\r
+                    SwapRows(rows, r, pivotIndex);\r
+                    pivot = rows[pivotIndex][i];\r
+                    for (int m = pivotIndex + 1; m < row; ++m)\r
+                    {\r
+                        T[] rowVec = rows[m];\r
+                        for (int n = i + 1; n < col; ++n)\r
+                        {\r
+                            T x = Operator<T>.Mul(pivot, rowVec[n]);\r
+                            T y = Operator<T>.Mul(rowVec[i],\r
+                                                        rows[pivotIndex][n]);\r
+                            T z = Operator<T>.Sub(x, y);\r
+                            rowVec[n] = Operator<T>.Mul(z, prevPivotR);\r
+                        }\r
+                    }\r
+                    if (++pivotIndex == row)\r
+                    {\r
+                        return row;\r
+                    }\r
+                    prevPivotR = Operator<T>.Div(Operator<T>.One, pivot);\r
+                }\r
+                return pivotIndex;\r
             }\r
         }\r
 \r
@@ -410,7 +483,7 @@ namespace Karinto.Numeric
         #endregion\r
 \r
         #region public methods\r
-        public Object Clone()\r
+        public virtual Object Clone()\r
         {\r
             return new Matrix<T>(e);\r
         }\r
@@ -436,7 +509,91 @@ namespace Karinto.Numeric
         /// <returns>逆行列</returns>\r
         public virtual Matrix<T> Inverse()\r
         {\r
-            throw new NotImplementedException();\r
+            CheckSquare(this);\r
+            T one = Operator<T>.One;\r
+            T zero = Operator<T>.Zero;\r
+            T[][] inv = new T[row][];\r
+            T[][] rows = new T[row][];\r
+            if (!new T().Equals(zero))\r
+            {\r
+                for (int r = 0; r < row; ++r)\r
+                {\r
+                    T[] vector = new T[col];\r
+                    for (int c = 0; c < col; ++c)\r
+                    {\r
+                        vector[c] = zero;\r
+                    }\r
+                    vector[r] = one;\r
+                    inv[r] = vector;\r
+                    rows[r] = (T[])this[r];\r
+                }\r
+            }\r
+            else\r
+            {\r
+                for (int r = 0; r < row; ++r)\r
+                {\r
+                    T[] vector = new T[col];\r
+                    vector[r] = one;\r
+                    inv[r] = vector;\r
+                    rows[r] = (T[])this[r];\r
+                }\r
+            }\r
+\r
+            for (int i = 0; i < row; ++i)\r
+            {\r
+                int pivotIndex = i;\r
+                T max = Operator<T>.Abs(rows[i][i]);\r
+                for (int r = i + 1; r < row; ++r)\r
+                {\r
+                    T mag = Operator<T>.Abs(rows[r][i]);\r
+                    if (mag.CompareTo(max) > 0)\r
+                    {\r
+                        pivotIndex = r;\r
+                        max = mag;\r
+                    }\r
+                }\r
+                if (max.Equals(zero))\r
+                {\r
+                    throw new InvalidOperationException("not reqular");\r
+                }\r
+                if (i != pivotIndex)\r
+                {\r
+                    SwapRows(rows, i, pivotIndex);\r
+                    SwapRows(inv, i, pivotIndex);\r
+                }\r
+                T pivotR = Operator<T>.Div(one, rows[i][i]);\r
+                for (int r = 0; r < row; ++r)\r
+                {\r
+                    if (r == i)\r
+                    {\r
+                        continue;\r
+                    }\r
+                    T q = Operator<T>.Mul(rows[r][i], pivotR);\r
+                    rows[r][i] = zero;\r
+                    for (int c = i + 1; c < col; ++c)\r
+                    {\r
+                        T x = Operator<T>.Mul(rows[i][c], q);\r
+                        T y = Operator<T>.Sub(rows[r][c], x);\r
+                        rows[r][c] = y;\r
+                    }\r
+                    for (int c = 0; c < col; ++c)\r
+                    {\r
+                        T x = Operator<T>.Mul(inv[i][c], q);\r
+                        T y = Operator<T>.Sub(inv[r][c], x);\r
+                        inv[r][c] = y;\r
+                    }\r
+                }\r
+                rows[i][i] = one;\r
+                ScaleRow(rows[i], i + 1, col, pivotR);\r
+                ScaleRow(inv[i], 0, col, pivotR);\r
+\r
+            }\r
+            Matrix<T> m = CreateSkeleton(row, col);\r
+            for (int r = 0; r < row; ++r)\r
+            {\r
+                m.SetRow(r, inv[r]);\r
+            }\r
+            return m;\r
         }\r
 \r
         /// <summary>\r
@@ -514,6 +671,34 @@ namespace Karinto.Numeric
             return e.GetHashCode() ^ row ^ col;\r
         }\r
 \r
+        public override string ToString()\r
+        {\r
+            int length = Math.Max(Math.Max(e[0, 0].ToString().Length,\r
+                                        e[row - 1, col - 1].ToString().Length),\r
+                                    Math.Max(e[row - 1, 0].ToString().Length,\r
+                                        e[0, col - 1].ToString().Length));\r
+            int capacity = row * col * (length + 2) + (row * 3) + 64;\r
+            StringBuilder sb = new StringBuilder(capacity);\r
+            sb.Append(GetType().Name);\r
+            sb.Append("{");\r
+            for (int r = 0; r < row; ++r)\r
+            {\r
+                if (r > 0)\r
+                {\r
+                    sb.Append("}, ");\r
+                }\r
+                sb.Append("{");\r
+                sb.Append(e[r, 0].ToString());\r
+                for (int c = 1; c < col; ++c)\r
+                {\r
+                    sb.Append(", ");\r
+                    sb.Append(e[r, c].ToString());\r
+                }\r
+            }\r
+            sb.Append("}}");\r
+            return sb.ToString();\r
+        }\r
+\r
         #endregion\r
 \r
         #region private or protected methods\r
@@ -521,6 +706,7 @@ namespace Karinto.Numeric
         {\r
             return new Matrix<T>(row, column);\r
         }\r
+\r
         private T[] GetColumn(int columnIndex)\r
         {\r
             if (cs == null)\r
@@ -585,6 +771,21 @@ namespace Karinto.Numeric
                 }\r
             }\r
         }\r
+\r
+        private static void SwapRows(T[][] rows, int index1, int index2)\r
+        {\r
+            T[] t = rows[index1];\r
+            rows[index1] = rows[index2];\r
+            rows[index2] = t;\r
+        }\r
+\r
+        private static void ScaleRow(T[] row, int start, int end, T scaler)\r
+        {\r
+            for (int c = start; c < end; ++c)\r
+            {\r
+                row[c] = Operator<T>.Mul(row[c], scaler);\r
+            }\r
+        }\r
         #endregion\r
     }\r
 }\r

diff --git a/Karinto/Numeric/Matrix2x2.cs b/Karinto/Numeric/Matrix2x2.cs
index e2605b7..71dafc4 100755 (executable)
--- a/Karinto/Numeric/Matrix2x2.cs
+++ b/Karinto/Numeric/Matrix2x2.cs
@@ -14,7 +14,7 @@ namespace Karinto.Numeric
     /// <summary>\r
     /// 2×2の倍精度浮動小数点行列\r
     /// </summary>\r
-    public class Matrix2x2 : Matrix\r
+    public class Matrix2x2 : Matrix, IComparable\r
     {\r
         #region constructors\r
         public Matrix2x2()\r
@@ -39,6 +39,32 @@ namespace Karinto.Numeric
                 throw new ArgumentException("size mismatch");\r
             }\r
         }\r
+\r
+        static Matrix2x2()\r
+        {\r
+            Operator<Matrix2x2>.Zero = new Matrix2x2();\r
+            Operator<Matrix2x2>.One = Matrix2x2.I();\r
+            Operator<Matrix2x2>.Add = delegate(Matrix2x2 l, Matrix2x2 r)\r
+            {\r
+                return l + r;\r
+            };\r
+            Operator<Matrix2x2>.Sub = delegate(Matrix2x2 l, Matrix2x2 r)\r
+            {\r
+                return l - r;\r
+            };\r
+            Operator<Matrix2x2>.Mul = delegate(Matrix2x2 l, Matrix2x2 r)\r
+            {\r
+                return l * r;\r
+            };\r
+            Operator<Matrix2x2>.Div = delegate(Matrix2x2 l, Matrix2x2 r)\r
+            {\r
+                return l / r;\r
+            };\r
+            Operator<Matrix2x2>.Abs = delegate(Matrix2x2 m)\r
+            {\r
+                throw new InvalidOperationException();\r
+            };\r
+        }\r
         #endregion\r
 \r
         #region operators\r
@@ -75,8 +101,86 @@ namespace Karinto.Numeric
             m.e[1, 1] = l10r01 + left.e[1, 1] * right.e[1, 1];\r
             return m;\r
         }\r
+\r
+        public static Matrix2x2 operator /(Matrix2x2 left, Matrix2x2 right)\r
+        {\r
+            double d = 1.0 / (right.e[0, 0] * right.e[1, 1] -\r
+                                                right.e[0, 1] * right.e[1, 0]);\r
+            Matrix2x2 m = new Matrix2x2();\r
+            double r00 = right.e[1, 1] * d;\r
+            double r01 = right.e[0, 1] * -d;\r
+            double r10 = right.e[1, 0] * -d;\r
+            double r11 = right.e[0, 0] * d;\r
+            m.e[0, 0] = left.e[0, 0] * r00 + left.e[0, 1] * r10;\r
+            m.e[0, 1] = left.e[0, 0] * r01 + left.e[0, 1] * r11;\r
+            m.e[1, 0] = left.e[1, 0] * r00 + left.e[1, 1] * r10;\r
+            m.e[1, 1] = left.e[1, 0] * r01 + left.e[1, 1] * r11;\r
+            return m;\r
+        }\r
         #endregion\r
 \r
+        #region properties\r
+        public override double Determinant\r
+        {\r
+            get\r
+            {\r
+                return e[0, 0] * e[1, 1] - e[0, 1] * e[1, 0];\r
+            }\r
+        }\r
+\r
+        public override int Rank\r
+        {\r
+            get\r
+            {\r
+                if (Determinant != 0)\r
+                {\r
+                    return 2;\r
+                }\r
+                return this.Equals(Operator<Matrix2x2>.Zero) ? 0 : 1;\r
+            }\r
+        }\r
+\r
+        public override double Trace\r
+        {\r
+            get\r
+            {\r
+                return e[0, 0] + e[1, 1];\r
+            }\r
+        }\r
+        #endregion\r
+\r
+        #region public methods\r
+        public int CompareTo(Object other)\r
+        {\r
+            Matrix2x2 m = other as Matrix2x2;\r
+            if (m == null)\r
+            {\r
+                throw new ArgumentException();\r
+            }\r
+            return Determinant.CompareTo(m.Determinant);\r
+        }\r
+\r
+        public new Matrix2x2 Inverse()\r
+        {\r
+            double d = 1.0 / (e[0, 0] * e[1, 1] - e[0, 1] * e[1, 0]);\r
+            Matrix2x2 m = new Matrix2x2();\r
+            m.e[0, 0] = e[1, 1] * d;\r
+            m.e[0, 1] = e[0, 1] * -d;\r
+            m.e[1, 0] = e[1, 0] * -d;\r
+            m.e[1, 1] = e[0, 0] * d;\r
+            return m;\r
+        }\r
+\r
+        public new Matrix2x2 Transpose()\r
+        {\r
+            Matrix2x2 m = new Matrix2x2();\r
+            m.e[0, 0] = e[0, 0];\r
+            m.e[0, 1] = e[1, 0];\r
+            m.e[1, 0] = e[0, 1];\r
+            m.e[1, 1] = e[1, 1];\r
+            return m;\r
+        }\r
+\r
         public static new Matrix2x2 I(int size)\r
         {\r
             if (size != 2)\r
@@ -90,5 +194,14 @@ namespace Karinto.Numeric
         {\r
             return new Matrix2x2(1.0, 0.0, 0.0, 1.0);\r
         }\r
+        \r
+        #endregion\r
+\r
+        #region private or protected methods\r
+        protected override Matrix<double> CreateSkeleton(int row, int column)\r
+        {\r
+            return new Matrix2x2();\r
+        }\r
+        #endregion\r
     }\r
 }\r

diff --git a/KarintoTest/NumericTest/Matrix2x2Test.cs b/KarintoTest/NumericTest/Matrix2x2Test.cs
index 94e4a10..c34fc1e 100755 (executable)
--- a/KarintoTest/NumericTest/Matrix2x2Test.cs
+++ b/KarintoTest/NumericTest/Matrix2x2Test.cs
@@ -33,6 +33,8 @@ namespace KarintoTest.NumericTest
                 16.0,\r
                 -32.0\r
             };\r
+            M dummy = new Matrix2x2();\r
+            Matrix dummy2 = new Matrix(2, 2);\r
         }\r
 \r
         [Test]\r
@@ -79,6 +81,7 @@ namespace KarintoTest.NumericTest
             }\r
             TimeSpan dt1 = DateTime.Now - start;\r
 \r
+            Matrix dummy2 = Matrix.I(2);\r
             start = DateTime.Now;\r
             Matrix b = new Matrix(2, 2);\r
             Matrix y = new Matrix(2, 2, 1.0, -1.0, 0.5, -0.5);\r

diff --git a/KarintoTest/NumericTest/MatrixTest.cs b/KarintoTest/NumericTest/MatrixTest.cs
index 3070173..4a9b97f 100755 (executable)
--- a/KarintoTest/NumericTest/MatrixTest.cs
+++ b/KarintoTest/NumericTest/MatrixTest.cs
@@ -31,10 +31,10 @@ namespace KarintoTest.NumericTest
         static MatrixTest()\r
         {\r
             d = new double[]{\r
-                1.0,\r
-                -2.0,\r
+                2.0,\r
+                -1.0,\r
                 0.5,\r
-                -0.25,\r
+                -1.25,\r
                 1.5,\r
                 -3.5,\r
                 4.0,\r
@@ -125,7 +125,10 @@ namespace KarintoTest.NumericTest
         {\r
             M m1x1Sq = new M(1, 1, d[0] * d[0]);\r
             Assert.AreEqual(m1x1Sq, m1x1 * m1x1);\r
-            M m2x2Sq = new M(2, 2, 0.0, -1.5, 0.375, -0.9375);\r
+            M m2x2Sq = new M(2, 2, d[0] * d[0] + d[1] * d[2],\r
+                                    d[0] * d[1] + d[1] * d[3],\r
+                                    d[2] * d[0] + d[3] * d[2],\r
+                                    d[2] * d[1] + d[3] * d[3]);\r
             Assert.AreEqual(m2x2Sq, m2x2 * m2x2);\r
             M m1x3m3x1 = new M(1, 1, d[0] * d[0] + d[1] * d[1] + d[2] * d[2]);\r
             Assert.AreEqual(m1x3m3x1, m1x3 * m3x1);\r
@@ -154,6 +157,23 @@ namespace KarintoTest.NumericTest
         }\r
 \r
         [Test]\r
+        public void Rank()\r
+        {\r
+            Assert.AreEqual(0, new M(1, 1).Rank);\r
+            Assert.AreEqual(1, m1x1.Rank);\r
+            Assert.AreEqual(1, m1x3.Rank);\r
+            Assert.AreEqual(1, m3x1.Rank);\r
+            Assert.AreEqual(2, m2x2.Rank);\r
+            Assert.AreEqual(3, m3x3.Rank);\r
+            M m3x4 = new M(3, 4, \r
+                            1.0, 2.0, 3.0, 4.0,\r
+                            0.5, 1.0, 1.5, 2.0,\r
+                            0.0, 0.0, 0.0, 0.0\r
+                            );\r
+            Assert.AreEqual(1, m3x4.Rank);\r
+        }\r
+\r
+        [Test]\r
         public void Trace()\r
         {\r
             Assert.AreEqual(d[0], m1x1.Trace);\r
@@ -164,6 +184,9 @@ namespace KarintoTest.NumericTest
         [Test]\r
         public void Inverse()\r
         {\r
+            AreEqual(M.I(1), m1x1.Inverse() * m1x1, 3e-16);\r
+            AreEqual(M.I(2), m2x2.Inverse() * m2x2, 3e-16);\r
+            AreEqual(M.I(3), m3x3.Inverse() * m3x3, 3e-10);\r
         }\r
 \r
         [Test]\r
@@ -172,5 +195,28 @@ namespace KarintoTest.NumericTest
             Assert.AreEqual(m1x3, m3x1.Transpose());\r
             Assert.AreEqual(m3x1, m1x3.Transpose());\r
         }\r
+\r
+        [Test]\r
+        public new void ToString()\r
+        {\r
+            string m2x2s = String.Format(\r
+                                        "Matrix{{{{{0}, {1}}}, {{{2}, {3}}}}}",\r
+                                        d[0], d[1], d[2], d[3]);\r
+            Assert.AreEqual(m2x2s, m2x2.ToString());\r
+        }\r
+\r
+        private static void AreEqual(M expected, M actual, double delta)\r
+        {\r
+            Assert.AreEqual(expected.RowSize, actual.RowSize);\r
+            Assert.AreEqual(expected.ColumnSize, actual.ColumnSize);\r
+            for (int r = 0; r < expected.ColumnSize; ++r)\r
+            {\r
+                for (int c = 0; c < expected.ColumnSize; ++c)\r
+                { \r
+                    Assert.AreEqual(expected[r, c], actual[r, c], \r
+                        delta, "Matrix[{0}, {1}]\n{2}\n{3}", r, c, expected, actual);\r
+                }\r
+            }\r
+        }\r
     }\r
 }\r