public class Matrix : Matrix<double>\r
{\r
#region constructors\r
-\r
public Matrix(int row, int column)\r
: base(row, column)\r
{\r
: base(elements)\r
{\r
}\r
+\r
#endregion\r
\r
#region operators\r
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
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
col = elements.GetLength(1);\r
e = (T[,])elements.Clone();\r
}\r
-\r
#endregion\r
\r
#region operators\r
}\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
/// <summary>\r
/// 列ベクトルの配列\r
/// </summary>\r
+ /// <remarks>取得した配列,要素に破壊的操作をしてはならない</remarks>\r
public virtual T[][] Rows\r
{\r
get\r
/// <summary>\r
/// 行ベクトルの配列\r
/// </summary>\r
+ /// <remarks>取得した配列,要素に破壊的操作をしてはならない</remarks>\r
public virtual T[][] Columns\r
{\r
get\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
+ 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
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
{\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
#endregion\r
\r
#region public methods\r
- public Object Clone()\r
+ public virtual Object Clone()\r
{\r
return new Matrix<T>(e);\r
}\r
/// <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
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
{\r
return new Matrix<T>(row, column);\r
}\r
+\r
private T[] GetColumn(int columnIndex)\r
{\r
if (cs == null)\r
}\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
/// <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
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
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
{\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
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
{\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
}\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
[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
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