<ItemGroup>\r
<Reference Include="System" />\r
<Reference Include="System.Data" />\r
+ <Reference Include="System.Drawing" />\r
<Reference Include="System.Xml" />\r
</ItemGroup>\r
<ItemGroup>\r
<Compile Include="HiokiHicorderData.cs" />\r
<Compile Include="HiokiHicorderDataReader.cs" />\r
<Compile Include="Numeric\Complex.cs" />\r
+ <Compile Include="Numeric\Matrix.cs" />\r
+ <Compile Include="Numeric\Matrix2x2.cs" />\r
<Compile Include="Numeric\Operator.cs" />\r
<Compile Include="Numeric\Rational.cs" />\r
<Compile Include="RegexSet.cs" />\r
--- /dev/null
+/*\r
+ * Karinto Library Project\r
+ *\r
+ * This software is distributed under a zlib-style license.\r
+ * See license.txt for more information.\r
+ */\r
+\r
+using System;\r
+using System.Collections.Generic;\r
+using System.Text;\r
+\r
+namespace Karinto.Numeric\r
+{\r
+ #region Matrix of Double\r
+ /// <summary>\r
+ /// 倍精度浮動小数点数の行列\r
+ /// </summary>\r
+ public class Matrix : Matrix<double>\r
+ {\r
+ #region constructors\r
+\r
+ public Matrix(int row, int column)\r
+ : base(row, column)\r
+ {\r
+ }\r
+\r
+ public Matrix(int row, int column, params double[] elements)\r
+ : base(row, column, elements)\r
+ {\r
+ }\r
+\r
+ public Matrix(double[,] elements)\r
+ : base(elements)\r
+ {\r
+ }\r
+ #endregion\r
+\r
+ #region operators\r
+\r
+ public static Matrix operator +(\r
+ Matrix left, Matrix<double> right)\r
+ {\r
+ return (Matrix)((Matrix<double>)left + right);\r
+ }\r
+\r
+\r
+ public static Matrix operator *(Matrix left, Matrix<double> right)\r
+ {\r
+ return (Matrix)((Matrix<double>)left * right);\r
+ }\r
+\r
+ public static Matrix operator *(double left, Matrix right)\r
+ {\r
+ return (Matrix)(left * (Matrix<double>)right);\r
+ }\r
+\r
+ public static Matrix operator *(Matrix left, double right)\r
+ {\r
+ return (Matrix)((Matrix<double>)left * right);\r
+ }\r
+\r
+ #endregion\r
+ protected override Matrix<double> CreateSkeleton(int row, int column)\r
+ {\r
+ return new Matrix(row, column);\r
+ }\r
+\r
+ public static new Matrix I(int size)\r
+ {\r
+ Matrix m = new Matrix(size, size);\r
+ for (int i = 0; i < size; ++i)\r
+ {\r
+ m.e[i, i] = 1.0;\r
+ }\r
+ return m;\r
+ }\r
+ }\r
+ #endregion\r
+\r
+\r
+\r
+\r
+\r
+ /// <summary>\r
+ /// 行列\r
+ /// </summary>\r
+ /// <param name="T">要素の型</param>\r
+ public class Matrix<T> : ICloneable where T : IComparable, new()\r
+ {\r
+ #region fields\r
+ protected int row;\r
+ protected int col;\r
+ protected T[,] e;\r
+ protected T[][] rs;\r
+ protected T[][] cs;\r
+ #endregion\r
+\r
+ #region constructors\r
+ public Matrix(int row, int column)\r
+ {\r
+ this.row = row;\r
+ col = column;\r
+ e = new T[row, col];\r
+ if (!e[0, 0].Equals(Operator<T>.Zero))\r
+ {\r
+ FillWithZero(this);\r
+ }\r
+ }\r
+\r
+ public Matrix(int row, int column, params T[] elements)\r
+ : this(row, column)\r
+ {\r
+ if (elements.Length != (row * column))\r
+ {\r
+ throw new ArgumentException(\r
+ "The Length of 'elements' is inconsistent with 'row' and 'column'");\r
+ }\r
+ int i = 0;\r
+ for (int r = 0; r < row; ++r)\r
+ {\r
+ for (int c = 0; c < col; ++c)\r
+ {\r
+ e[r, c] = elements[i++];\r
+ }\r
+ }\r
+ }\r
+\r
+ public Matrix(T[,] elements)\r
+ {\r
+ row = elements.GetLength(0);\r
+ col = elements.GetLength(1);\r
+ e = (T[,])elements.Clone();\r
+ }\r
+\r
+ #endregion\r
+\r
+ #region operators\r
+ public T[] this[int rowIndex]\r
+ {\r
+ get\r
+ {\r
+ if (rs == null)\r
+ {\r
+ rs = new T[row][];\r
+ }\r
+ if (rs[rowIndex] == null)\r
+ {\r
+ T[] vector = new T[col];\r
+ for (int c = 0; c < col; ++c)\r
+ {\r
+ vector[c] = e[rowIndex, c];\r
+ }\r
+ rs[rowIndex] = vector;\r
+ }\r
+ return rs[rowIndex];\r
+ }\r
+ }\r
+\r
+ public T this[int rowIndex, int columnIndex]\r
+ {\r
+ get\r
+ {\r
+ return e[rowIndex, columnIndex];\r
+ }\r
+ set\r
+ {\r
+ e[rowIndex, columnIndex] = value;\r
+ if (rs != null)\r
+ {\r
+ if (rs[rowIndex] != null)\r
+ {\r
+ rs[rowIndex][columnIndex] = value;\r
+ }\r
+ }\r
+ if (cs != null)\r
+ {\r
+ if (cs[columnIndex] != null)\r
+ {\r
+ cs[columnIndex][rowIndex] = value;\r
+ }\r
+ };\r
+ }\r
+ }\r
+\r
+ public static Matrix<T> operator +(Matrix<T> left, Matrix<T> right)\r
+ {\r
+ return Operate(left, right, Operator<T>.Add);\r
+ }\r
+\r
+ public static Matrix<T> operator -(Matrix<T> left, Matrix<T> right)\r
+ {\r
+ return Operate(left, right, Operator<T>.Sub);\r
+ }\r
+\r
+ public static Matrix<T> operator *(Matrix<T> left, Matrix<T> right)\r
+ {\r
+ int l = left.col;\r
+ if (l != right.row)\r
+ {\r
+ throw new Exception("size mismatch");\r
+ }\r
+ T zero = Operator<T>.Zero;\r
+ Matrix<T> m = left.CreateSkeleton(left.row, right.col);\r
+ for (int r = 0; r < left.row; ++r)\r
+ {\r
+ for (int c = 0; c < right.col; ++c)\r
+ {\r
+ T element = zero;\r
+ for (int i = 0; i < l; ++i)\r
+ {\r
+ T dot = Operator<T>.Mul(left.e[r, i], right.e[i, c]);\r
+ element = Operator<T>.Add(element, dot);\r
+ }\r
+ m.e[r, c] = element;\r
+ }\r
+ }\r
+ return m;\r
+ }\r
+\r
+ public static Matrix<T> operator *(Matrix<T> left, T right)\r
+ {\r
+ Matrix<T> m = left.CreateSkeleton(left.row, left.col);\r
+ for (int r = 0; r < m.row; ++r)\r
+ {\r
+ for (int c = 0; c < m.col; ++c)\r
+ {\r
+ m.e[r, c] = Operator<T>.Mul(left[r, c], right);\r
+ }\r
+ }\r
+ return m;\r
+ }\r
+\r
+ public static Matrix<T> operator *(T left, Matrix<T> right)\r
+ {\r
+ Matrix<T> m = right.CreateSkeleton(right.row, right.col);\r
+ for (int r = 0; r < m.row; ++r)\r
+ {\r
+ for (int c = 0; c < m.col; ++c)\r
+ {\r
+ m.e[r, c] = Operator<T>.Mul(left, right[r, c]);\r
+ }\r
+ }\r
+ return m;\r
+ }\r
+ #endregion\r
+\r
+ #region properties\r
+\r
+ /// <summary>\r
+ /// 行数\r
+ /// </summary>\r
+ public int RowSize\r
+ {\r
+ get\r
+ {\r
+ return row;\r
+ }\r
+ }\r
+\r
+ /// <summary>\r
+ /// 列数\r
+ /// </summary>\r
+ public int ColumnSize\r
+ {\r
+ get\r
+ {\r
+ return col;\r
+ }\r
+ }\r
+\r
+ /// <summary>\r
+ /// 列ベクトルの配列\r
+ /// </summary>\r
+ public virtual T[][] Rows\r
+ {\r
+ get\r
+ {\r
+ for (int r = 0; r < row; ++r)\r
+ {\r
+ T[] vector = this[r];\r
+ rs[r] = vector;\r
+ }\r
+ return rs;\r
+ }\r
+ }\r
+\r
+ /// <summary>\r
+ /// 行ベクトルの配列\r
+ /// </summary>\r
+ public virtual T[][] Columns\r
+ {\r
+ get\r
+ {\r
+ for (int c = 0; c < col; ++c)\r
+ {\r
+ T[] vector = GetColumn(c);\r
+ cs[c] = vector;\r
+ }\r
+ return cs;\r
+ }\r
+ }\r
+\r
+ /// <summary>\r
+ /// 行列式\r
+ /// </summary>\r
+ public virtual T Determinant\r
+ {\r
+ 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
+ {\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
+ {\r
+ prevPivot = pivot;\r
+ pivot = rows[i][i];\r
+ if (pivot.Equals(zero))\r
+ {\r
+ int r = 0;\r
+ while (rows[r][i].Equals(zero))\r
+ {\r
+ if (++r >= l)\r
+ {\r
+ r = 0;\r
+ break;\r
+ }\r
+ }\r
+ T[] t = rows[r];\r
+ rows[r] = rows[i];\r
+ rows[i] = t;\r
+ pivot = rows[i][i];\r
+ sign = !sign;\r
+ }\r
+ for (int j = i + 1; j < l; ++j)\r
+ {\r
+ T[] rowVec = rows[j];\r
+ for (int k = i + 1; k < l; ++k)\r
+ {\r
+ T x = Operator<T>.Mul(pivot, rowVec[k]);\r
+ T y = Operator<T>.Mul(rowVec[i], rows[i][k]);\r
+ T z = Operator<T>.Sub(x, y);\r
+ rowVec[k] = Operator<T>.Div(z, prevPivot);\r
+ }\r
+ }\r
+ }\r
+ return sign ? Operator<T>.Sub(zero, pivot) : pivot;\r
+ }\r
+ }\r
+\r
+ /// <summary>\r
+ /// 階数\r
+ /// </summary>\r
+ public virtual int Rank\r
+ {\r
+ get\r
+ {\r
+ throw new NotImplementedException();\r
+ }\r
+ }\r
+\r
+ /// <summary>\r
+ /// 正則行列か否か\r
+ /// </summary>\r
+ public virtual bool Regular\r
+ {\r
+ get\r
+ {\r
+ return !Determinant.Equals(Operator<T>.Zero);\r
+ }\r
+ }\r
+\r
+ /// <summary>\r
+ /// 正方行列か否か\r
+ /// </summary>\r
+ public bool Square\r
+ {\r
+ get\r
+ {\r
+ return row == col;\r
+ }\r
+ }\r
+\r
+ /// <summary>\r
+ /// 跡（トレース）\r
+ /// </summary>\r
+ public virtual T Trace\r
+ {\r
+ get\r
+ {\r
+ CheckSquare(this);\r
+ T t = e[0, 0];\r
+ for (int i = 1; i < row; ++i)\r
+ {\r
+ t = Operator<T>.Add(t, e[i, i]);\r
+ }\r
+ return t;\r
+ }\r
+ }\r
+\r
+ #endregion\r
+\r
+ #region public methods\r
+ public Object Clone()\r
+ {\r
+ return new Matrix<T>(e);\r
+ }\r
+\r
+ public static Matrix<T> I(int size)\r
+ {\r
+ T one = Operator<T>.One;\r
+ Matrix<T> m = new Matrix<T>(size, size);\r
+ if (m.e[0, 0].Equals(Operator<T>.Zero))\r
+ {\r
+ FillWithZero(m);\r
+ }\r
+ for (int i = 0; i < size; ++i)\r
+ {\r
+ m[i, i] = one;\r
+ }\r
+ return m;\r
+ }\r
+\r
+ /// <summary>\r
+ /// 逆行列を生成する\r
+ /// </summary>\r
+ /// <returns>逆行列</returns>\r
+ public virtual Matrix<T> Inverse()\r
+ {\r
+ throw new NotImplementedException();\r
+ }\r
+\r
+ /// <summary>\r
+ /// 転置行列を生成する\r
+ /// </summary>\r
+ /// <returns>転置行列</returns>\r
+ public virtual Matrix<T> Transpose()\r
+ {\r
+ Matrix<T> m = CreateSkeleton(col, row);\r
+ for (int r = 0; r < m.row; ++r)\r
+ {\r
+ for (int c = 0; c < m.col; ++c)\r
+ {\r
+ m.e[r, c] = e[c, r];\r
+ }\r
+ }\r
+ return m;\r
+ }\r
+\r
+ public virtual void SetRow(int index, T[] vector)\r
+ {\r
+ if (vector.Length != col)\r
+ {\r
+ throw new ArgumentException("size mismatch");\r
+ }\r
+ if (rs != null)\r
+ {\r
+ rs[index] = null;\r
+ }\r
+ for (int c = 0; c < col; ++c)\r
+ {\r
+ e[index, c] = vector[c];\r
+ }\r
+ }\r
+\r
+ public virtual void SetColumn(int index, T[] vector)\r
+ {\r
+ if (vector.Length != row)\r
+ {\r
+ throw new ArgumentException("size mismatch");\r
+ }\r
+ if (cs != null)\r
+ {\r
+ cs[index] = null;\r
+ }\r
+ for (int r = 0; r < row; ++r)\r
+ {\r
+ e[r, index] = vector[r];\r
+ }\r
+ }\r
+\r
+ public override bool Equals(object other)\r
+ {\r
+ Matrix<T> m = other as Matrix<T>;\r
+ if (m == null)\r
+ {\r
+ throw new ArgumentException();\r
+ }\r
+ CheckSize(this, m);\r
+ for (int r = 0; r < m.row; ++r)\r
+ {\r
+ for (int c = 0; c < m.col; ++c)\r
+ {\r
+ if (!e[r, c].Equals(m.e[r, c]))\r
+ {\r
+ return false;\r
+ }\r
+ }\r
+ }\r
+ return true;\r
+ }\r
+\r
+ public override int GetHashCode()\r
+ {\r
+ return e.GetHashCode() ^ row ^ col;\r
+ }\r
+\r
+ #endregion\r
+\r
+ #region private or protected methods\r
+ protected virtual Matrix<T> CreateSkeleton(int row, int column)\r
+ {\r
+ return new Matrix<T>(row, column);\r
+ }\r
+ private T[] GetColumn(int columnIndex)\r
+ {\r
+ if (cs == null)\r
+ {\r
+ cs = new T[col][];\r
+ }\r
+ if (cs[columnIndex] == null)\r
+ {\r
+ T[] vector = new T[row];\r
+ for (int r = 0; r < row; ++r)\r
+ {\r
+ vector[r] = e[r, columnIndex];\r
+ }\r
+ cs[columnIndex] = vector;\r
+ }\r
+ return cs[columnIndex];\r
+ }\r
+\r
+ private static void CheckSize(Matrix<T> a, Matrix<T> b)\r
+ {\r
+ if (a.row != b.row)\r
+ {\r
+ throw new ArgumentException("row size mismatch");\r
+ }\r
+ if (a.col != b.col)\r
+ {\r
+ throw new ArgumentException("colmun size mismatch");\r
+ }\r
+ }\r
+\r
+ private static void CheckSquare(Matrix<T> m)\r
+ {\r
+ if (m.row != m.col)\r
+ {\r
+ throw new InvalidOperationException("non-square matrix");\r
+ }\r
+ }\r
+\r
+ private static Matrix<T> Operate(Matrix<T> left, Matrix<T> right,\r
+ Operator<T>.BinaryOperator op)\r
+ {\r
+ CheckSize(left, right);\r
+ Matrix<T> m = left.CreateSkeleton(left.row, left.col);\r
+ for (int r = 0; r < m.row; ++r)\r
+ {\r
+ for (int c = 0; c < m.col; ++c)\r
+ {\r
+ m.e[r, c] = op(left[r, c], right[r, c]);\r
+ }\r
+ }\r
+ return m;\r
+ }\r
+\r
+ private static void FillWithZero(Matrix<T> m)\r
+ {\r
+ T zero = Operator<T>.Zero;\r
+ for (int r = 0; r < m.row; ++r)\r
+ {\r
+ for (int c = 0; c < m.col; ++c)\r
+ {\r
+ m.e[r, c] = zero;\r
+ }\r
+ }\r
+ }\r
+ #endregion\r
+ }\r
+}\r
--- /dev/null
+/*\r
+ * Karinto Library Project\r
+ *\r
+ * This software is distributed under a zlib-style license.\r
+ * See license.txt for more information.\r
+ */\r
+\r
+using System;\r
+using System.Collections.Generic;\r
+using System.Text;\r
+\r
+namespace Karinto.Numeric\r
+{\r
+ /// <summary>\r
+ /// 2×2の倍精度浮動小数点行列\r
+ /// </summary>\r
+ public class Matrix2x2 : Matrix\r
+ {\r
+ #region constructors\r
+ public Matrix2x2()\r
+ : base(2, 2)\r
+ {\r
+ }\r
+\r
+ public Matrix2x2(double e11, double e12, double e21, double e22)\r
+ : base(2, 2)\r
+ {\r
+ e[0, 0] = e11;\r
+ e[0, 1] = e12;\r
+ e[1, 0] = e21;\r
+ e[1, 1] = e22;\r
+ }\r
+\r
+ public Matrix2x2(double[,] elements)\r
+ : base(elements)\r
+ {\r
+ if (row != 2 || col != 2)\r
+ {\r
+ throw new ArgumentException("size mismatch");\r
+ }\r
+ }\r
+ #endregion\r
+\r
+ #region operators\r
+ public static Matrix2x2 operator +(Matrix2x2 left, Matrix2x2 right)\r
+ {\r
+ Matrix2x2 m = new Matrix2x2();\r
+ m.e[0, 0] = left.e[0, 0] + right.e[0, 0];\r
+ m.e[0, 1] = left.e[0, 1] + right.e[0, 1];\r
+ m.e[1, 0] = left.e[1, 0] + right.e[1, 0];\r
+ m.e[1, 1] = left.e[1, 1] + right.e[1, 1];\r
+ return m;\r
+ }\r
+\r
+ public static Matrix2x2 operator -(Matrix2x2 left, Matrix2x2 right)\r
+ {\r
+ Matrix2x2 m = new Matrix2x2();\r
+ m.e[0, 0] = left.e[0, 0] - right.e[0, 0];\r
+ m.e[0, 1] = left.e[0, 1] - right.e[0, 1];\r
+ m.e[1, 0] = left.e[1, 0] - right.e[1, 0];\r
+ m.e[1, 1] = left.e[1, 1] - right.e[1, 1];\r
+ return m;\r
+ }\r
+\r
+ public static Matrix2x2 operator *(Matrix2x2 left, Matrix2x2 right)\r
+ {\r
+ Matrix2x2 m = new Matrix2x2();\r
+ double l00r00 = left.e[0, 0] * right.e[0, 0];\r
+ double l00r01 = left.e[0, 0] * right.e[0, 1];\r
+ double l10r00 = left.e[1, 0] * right.e[0, 0];\r
+ double l10r01 = left.e[1, 0] * right.e[0, 1];\r
+ m.e[0, 0] = l00r00 + left.e[0, 1] * right.e[1, 0];\r
+ m.e[0, 1] = l00r01 + left.e[0, 1] * right.e[1, 1];\r
+ m.e[1, 0] = l10r00 + left.e[1, 1] * right.e[1, 0];\r
+ m.e[1, 1] = l10r01 + left.e[1, 1] * right.e[1, 1];\r
+ return m;\r
+ }\r
+ #endregion\r
+\r
+ public static new Matrix2x2 I(int size)\r
+ {\r
+ if (size != 2)\r
+ {\r
+ throw new ArgumentException("size must be 2");\r
+ }\r
+ return Matrix2x2.I();\r
+ }\r
+\r
+ public static Matrix2x2 I()\r
+ {\r
+ return new Matrix2x2(1.0, 0.0, 0.0, 1.0);\r
+ }\r
+ }\r
+}\r
<ItemGroup>\r
<Compile Include="CsvWriterTest.cs" />\r
<Compile Include="NumericTest\ComplexTest.cs" />\r
+ <Compile Include="NumericTest\Matrix2x2Test.cs" />\r
<Compile Include="NumericTest\OperatorTest.cs" />\r
+ <Compile Include="NumericTest\MatrixTest.cs" />\r
<Compile Include="PointListTest.cs" />\r
<Compile Include="PointTest.cs" />\r
<Compile Include="RangeTest.cs" />\r
--- /dev/null
+/*\r
+ * Karinto Library Project\r
+ *\r
+ * This software is distributed under a zlib-style license.\r
+ * See license.txt for more information.\r
+ */\r
+\r
+using System;\r
+using System.Collections.Generic;\r
+using System.Text;\r
+using System.Diagnostics;\r
+using Karinto.Numeric;\r
+using NUnit.Framework;\r
+using M = Karinto.Numeric.Matrix2x2;\r
+\r
+namespace KarintoTest.NumericTest\r
+{\r
+ [TestFixture]\r
+ public class Matrix2x2Test\r
+ {\r
+ private static readonly double[] d;\r
+ static Matrix2x2Test()\r
+ {\r
+ d = new double[]{\r
+ 1.0,\r
+ -2.0,\r
+ 0.5,\r
+ -0.25,\r
+ 1.5,\r
+ -3.5,\r
+ 4.0,\r
+ -6.125,\r
+ 16.0,\r
+ -32.0\r
+ };\r
+ }\r
+\r
+ [Test]\r
+ public void Construct()\r
+ {\r
+ M m = new M();\r
+ Assert.AreEqual(2, m.RowSize);\r
+ Assert.AreEqual(2, m.ColumnSize);\r
+ Assert.AreEqual(0.0, m[0, 0]);\r
+\r
+ double[,] elements = new double[,] {\r
+ { d[0], d[1]},\r
+ { d[2], d[3]}\r
+ };\r
+ m = new M(elements);\r
+ Assert.AreEqual(2, m.RowSize);\r
+ Assert.AreEqual(2, m.ColumnSize);\r
+ Assert.AreEqual(d[0], m[0, 0]);\r
+ Assert.AreEqual(d[1], m[0, 1]);\r
+ Assert.AreEqual(d[2], m[1, 0]);\r
+ Assert.AreEqual(d[3], m[1, 1]);\r
+\r
+ m = new M(d[0], d[1], d[2], d[3]);\r
+ Assert.AreEqual(2, m.RowSize);\r
+ Assert.AreEqual(2, m.ColumnSize);\r
+ Assert.AreEqual(d[0], m[0, 0]);\r
+ Assert.AreEqual(d[1], m[0, 1]);\r
+ Assert.AreEqual(d[2], m[1, 0]);\r
+ Assert.AreEqual(d[3], m[1, 1]);\r
+\r
+ }\r
+\r
+ [Test]\r
+ public void Add()\r
+ {\r
+ const int n = 100000;\r
+ DateTime start;\r
+ start = DateTime.Now;\r
+ M a = new M();\r
+ M x = new M(1.0, -1.0, 0.5, -0.5);\r
+ for (int i = 0; i < n; ++i)\r
+ {\r
+ x += a;\r
+ }\r
+ TimeSpan dt1 = DateTime.Now - start;\r
+\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
+ for (int i = 0; i < n; ++i)\r
+ {\r
+ y += b;\r
+ }\r
+ TimeSpan dt2 = DateTime.Now - start;\r
+ \r
+ Assert.AreEqual(y, x);\r
+ Assert.LessOrEqual(dt1, dt2);\r
+ }\r
+\r
+ [Test]\r
+ public void Mul()\r
+ {\r
+ const int n = 100000;\r
+ DateTime start;\r
+ start = DateTime.Now;\r
+ M a = M.I();\r
+ M x = new M(1.0, -1.0, 0.5, -0.5);\r
+ for (int i = 0; i < n; ++i)\r
+ {\r
+ x *= a;\r
+ }\r
+ TimeSpan dt1 = DateTime.Now - start;\r
+\r
+ start = DateTime.Now;\r
+ Matrix b = Matrix.I(2);\r
+ Matrix y = new Matrix(2, 2, 1.0, -1.0, 0.5, -0.5);\r
+ for (int i = 0; i < n; ++i)\r
+ {\r
+ y *= b;\r
+ }\r
+ TimeSpan dt2 = DateTime.Now - start;\r
+\r
+ Assert.AreEqual(y, x);\r
+ Assert.LessOrEqual(dt1, dt2);\r
+ }\r
+ }\r
+}\r
--- /dev/null
+/*\r
+ * Karinto Library Project\r
+ *\r
+ * This software is distributed under a zlib-style license.\r
+ * See license.txt for more information.\r
+ */\r
+\r
+using System;\r
+using System.Collections.Generic;\r
+using System.Text;\r
+using System.Diagnostics;\r
+using Karinto.Numeric;\r
+using NUnit.Framework;\r
+using M = Karinto.Numeric.Matrix;\r
+\r
+namespace KarintoTest.NumericTest\r
+{\r
+ [TestFixture]\r
+ public class MatrixTest\r
+ {\r
+ private const double err = 3e-16;\r
+ private static readonly double[] d;\r
+ private static readonly M m1x1;\r
+ private static readonly M m1x3;\r
+ private static readonly M m3x1;\r
+ private static readonly M m2x2;\r
+ private static readonly M m3x3;\r
+ private static readonly M m4x4;\r
+ private static readonly M m10x10;\r
+\r
+ static MatrixTest()\r
+ {\r
+ d = new double[]{\r
+ 1.0,\r
+ -2.0,\r
+ 0.5,\r
+ -0.25,\r
+ 1.5,\r
+ -3.5,\r
+ 4.0,\r
+ -6.125,\r
+ 16.0,\r
+ -32.0\r
+ };\r
+\r
+ m1x1 = new M(1, 1, d[0]);\r
+ m1x3 = new M(1, 3, d[0], d[1], d[2]);\r
+ m3x1 = new M(3, 1, d[0], d[1], d[2]);\r
+ m2x2 = new M(2, 2, d[0], d[1], d[2], d[3]);\r
+ m3x3 = new M(3, 3,\r
+ d[0], d[1], d[2],\r
+ d[3], d[4], d[5],\r
+ d[6], d[7], d[8]);\r
+ }\r
+\r
+ [Test]\r
+ public void Construct()\r
+ {\r
+ M m = new M(1, 3);\r
+ Assert.AreEqual(1, m.RowSize);\r
+ Assert.AreEqual(3, m.ColumnSize);\r
+ Assert.AreEqual(0.0, m[0, 0]);\r
+\r
+ double[,] elements = new double[,] {\r
+ { d[0], d[1]},\r
+ { d[2], d[3]}\r
+ };\r
+ m = new M(elements);\r
+ Assert.AreEqual(2, m.RowSize);\r
+ Assert.AreEqual(2, m.ColumnSize);\r
+ Assert.AreEqual(d[0], m[0, 0]);\r
+ Assert.AreEqual(d[1], m[0, 1]);\r
+ Assert.AreEqual(d[2], m[1, 0]);\r
+ Assert.AreEqual(d[3], m[1, 1]);\r
+\r
+ double[] elements2 = new double[] { d[0], d[1], d[2], d[3] };\r
+ m = new M(2, 2, elements2);\r
+ Assert.AreEqual(2, m.RowSize);\r
+ Assert.AreEqual(2, m.ColumnSize);\r
+ Assert.AreEqual(d[0], m[0, 0]);\r
+ Assert.AreEqual(d[1], m[0, 1]);\r
+ Assert.AreEqual(d[2], m[1, 0]);\r
+ Assert.AreEqual(d[3], m[1, 1]);\r
+\r
+ }\r
+\r
+ [Test]\r
+ public void ChangeElement()\r
+ {\r
+ M m = new M(3, 3);\r
+ Assert.AreEqual(0.0, m[1, 2]);\r
+ Assert.AreEqual(0.0, m[1][2]);\r
+ Assert.AreEqual(0.0, m.Columns[2][1]);\r
+ m[1, 2] = 1.0;\r
+ Assert.AreEqual(1.0, m[1, 2]);\r
+ Assert.AreEqual(1.0, m[1][2]);\r
+ Assert.AreEqual(1.0, m.Columns[2][1]);\r
+ m[1, 2] = -1.0;\r
+ Assert.AreEqual(-1.0, m[1, 2]);\r
+ Assert.AreEqual(-1.0, m[1][2]);\r
+ Assert.AreEqual(-1.0, m.Columns[2][1]);\r
+ }\r
+\r
+ [Test]\r
+ public void Add()\r
+ {\r
+ M m1x1a = new M(1, 1, d[0] + d[0]);\r
+ Assert.AreEqual(m1x1a, m1x1 + m1x1);\r
+ }\r
+\r
+ [Test]\r
+ public void Sub()\r
+ {\r
+ Assert.AreEqual(m1x1, m1x1 + m1x1 - m1x1);\r
+ Assert.AreEqual(m1x3, m1x3 + m1x3 - m1x3);\r
+ Assert.AreEqual(m3x1, m3x1 + m3x1 - m3x1);\r
+ Assert.AreEqual(m2x2, m2x2 + m2x2 - m2x2);\r
+ Assert.AreEqual(m3x3, m3x3 + m3x3 - m3x3);\r
+ //Assert.AreEqual(m4x4, m4x4 + m4x4 - m4x4);\r
+ //Assert.AreEqual(m10x10, m10x10 + m10x10 - m10x10);\r
+ }\r
+\r
+ [Test]\r
+ public void Mul()\r
+ {\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
+ 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
+ Assert.AreEqual(3, (m3x1 * m1x3).RowSize);\r
+ Assert.AreEqual(3, (m3x1 * m1x3).ColumnSize);\r
+\r
+ Assert.AreEqual(m1x1 + m1x1, m1x1 * 2.0);\r
+ Assert.AreEqual(m1x3 + m1x3, m1x3 * 2.0);\r
+ Assert.AreEqual(m3x1 + m3x1, m3x1 * 2.0);\r
+ Assert.AreEqual(m2x2 + m2x2, m2x2 * 2.0);\r
+ Assert.AreEqual(m1x1, 1.0 * m1x1);\r
+ Assert.AreEqual(m1x3, 1.0 * m1x3);\r
+ Assert.AreEqual(m3x1, 1.0 * m3x1);\r
+ Assert.AreEqual(m2x2, 1.0 * m2x2);\r
+ }\r
+\r
+ [Test]\r
+ public void Determinant()\r
+ {\r
+ Assert.AreEqual(d[0], m1x1.Determinant);\r
+ Assert.AreEqual(d[0] * d[3] - d[1] * d[2], m2x2.Determinant);\r
+ Assert.AreEqual(d[0] * d[4] * d[8] - d[0] * d[5] * d[7]\r
+ + d[1] * d[5] * d[6] - d[1] * d[3] * d[8]\r
+ + d[2] * d[3] * d[7] - d[2] * d[4] * d[6]\r
+ , m3x3.Determinant);\r
+ }\r
+\r
+ [Test]\r
+ public void Trace()\r
+ {\r
+ Assert.AreEqual(d[0], m1x1.Trace);\r
+ Assert.AreEqual(d[0] + d[3], m2x2.Trace);\r
+ Assert.AreEqual(d[0] + d[4] + d[8], m3x3.Trace);\r
+ }\r
+\r
+ [Test]\r
+ public void Inverse()\r
+ {\r
+ }\r
+\r
+ [Test]\r
+ public void Transpose()\r
+ {\r
+ Assert.AreEqual(m1x3, m3x1.Transpose());\r
+ Assert.AreEqual(m3x1, m1x3.Transpose());\r
+ }\r
+ }\r
+}\r
OSDN Git Service
