summaryrefslogtreecommitdiff
path: root/ruby_1_8_6/lib/matrix.rb
diff options
context:
space:
mode:
Diffstat (limited to 'ruby_1_8_6/lib/matrix.rb')
-rw-r--r--ruby_1_8_6/lib/matrix.rb1272
1 files changed, 1272 insertions, 0 deletions
diff --git a/ruby_1_8_6/lib/matrix.rb b/ruby_1_8_6/lib/matrix.rb
new file mode 100644
index 0000000000..c62acdf9aa
--- /dev/null
+++ b/ruby_1_8_6/lib/matrix.rb
@@ -0,0 +1,1272 @@
+#!/usr/local/bin/ruby
+#--
+# matrix.rb -
+# $Release Version: 1.0$
+# $Revision: 1.11 $
+# $Date: 1999/10/06 11:01:53 $
+# Original Version from Smalltalk-80 version
+# on July 23, 1985 at 8:37:17 am
+# by Keiju ISHITSUKA
+#++
+#
+# = matrix.rb
+#
+# An implementation of Matrix and Vector classes.
+#
+# Author:: Keiju ISHITSUKA
+# Documentation:: Gavin Sinclair (sourced from <i>Ruby in a Nutshell</i> (Matsumoto, O'Reilly))
+#
+# See classes Matrix and Vector for documentation.
+#
+
+
+require "e2mmap.rb"
+
+module ExceptionForMatrix # :nodoc:
+ extend Exception2MessageMapper
+ def_e2message(TypeError, "wrong argument type %s (expected %s)")
+ def_e2message(ArgumentError, "Wrong # of arguments(%d for %d)")
+
+ def_exception("ErrDimensionMismatch", "\#{self.name} dimension mismatch")
+ def_exception("ErrNotRegular", "Not Regular Matrix")
+ def_exception("ErrOperationNotDefined", "This operation(%s) can\\'t defined")
+end
+
+#
+# The +Matrix+ class represents a mathematical matrix, and provides methods for creating
+# special-case matrices (zero, identity, diagonal, singular, vector), operating on them
+# arithmetically and algebraically, and determining their mathematical properties (trace, rank,
+# inverse, determinant).
+#
+# Note that although matrices should theoretically be rectangular, this is not
+# enforced by the class.
+#
+# Also note that the determinant of integer matrices may be incorrectly calculated unless you
+# also <tt>require 'mathn'</tt>. This may be fixed in the future.
+#
+# == Method Catalogue
+#
+# To create a matrix:
+# * <tt> Matrix[*rows] </tt>
+# * <tt> Matrix.[](*rows) </tt>
+# * <tt> Matrix.rows(rows, copy = true) </tt>
+# * <tt> Matrix.columns(columns) </tt>
+# * <tt> Matrix.diagonal(*values) </tt>
+# * <tt> Matrix.scalar(n, value) </tt>
+# * <tt> Matrix.scalar(n, value) </tt>
+# * <tt> Matrix.identity(n) </tt>
+# * <tt> Matrix.unit(n) </tt>
+# * <tt> Matrix.I(n) </tt>
+# * <tt> Matrix.zero(n) </tt>
+# * <tt> Matrix.row_vector(row) </tt>
+# * <tt> Matrix.column_vector(column) </tt>
+#
+# To access Matrix elements/columns/rows/submatrices/properties:
+# * <tt> [](i, j) </tt>
+# * <tt> #row_size </tt>
+# * <tt> #column_size </tt>
+# * <tt> #row(i) </tt>
+# * <tt> #column(j) </tt>
+# * <tt> #collect </tt>
+# * <tt> #map </tt>
+# * <tt> #minor(*param) </tt>
+#
+# Properties of a matrix:
+# * <tt> #regular? </tt>
+# * <tt> #singular? </tt>
+# * <tt> #square? </tt>
+#
+# Matrix arithmetic:
+# * <tt> *(m) </tt>
+# * <tt> +(m) </tt>
+# * <tt> -(m) </tt>
+# * <tt> #/(m) </tt>
+# * <tt> #inverse </tt>
+# * <tt> #inv </tt>
+# * <tt> ** </tt>
+#
+# Matrix functions:
+# * <tt> #determinant </tt>
+# * <tt> #det </tt>
+# * <tt> #rank </tt>
+# * <tt> #trace </tt>
+# * <tt> #tr </tt>
+# * <tt> #transpose </tt>
+# * <tt> #t </tt>
+#
+# Conversion to other data types:
+# * <tt> #coerce(other) </tt>
+# * <tt> #row_vectors </tt>
+# * <tt> #column_vectors </tt>
+# * <tt> #to_a </tt>
+#
+# String representations:
+# * <tt> #to_s </tt>
+# * <tt> #inspect </tt>
+#
+class Matrix
+ @RCS_ID='-$Id: matrix.rb,v 1.11 1999/10/06 11:01:53 keiju Exp keiju $-'
+
+# extend Exception2MessageMapper
+ include ExceptionForMatrix
+
+ # instance creations
+ private_class_method :new
+
+ #
+ # Creates a matrix where each argument is a row.
+ # Matrix[ [25, 93], [-1, 66] ]
+ # => 25 93
+ # -1 66
+ #
+ def Matrix.[](*rows)
+ new(:init_rows, rows, false)
+ end
+
+ #
+ # Creates a matrix where +rows+ is an array of arrays, each of which is a row
+ # to the matrix. If the optional argument +copy+ is false, use the given
+ # arrays as the internal structure of the matrix without copying.
+ # Matrix.rows([[25, 93], [-1, 66]])
+ # => 25 93
+ # -1 66
+ def Matrix.rows(rows, copy = true)
+ new(:init_rows, rows, copy)
+ end
+
+ #
+ # Creates a matrix using +columns+ as an array of column vectors.
+ # Matrix.columns([[25, 93], [-1, 66]])
+ # => 25 -1
+ # 93 66
+ #
+ #
+ def Matrix.columns(columns)
+ rows = (0 .. columns[0].size - 1).collect {
+ |i|
+ (0 .. columns.size - 1).collect {
+ |j|
+ columns[j][i]
+ }
+ }
+ Matrix.rows(rows, false)
+ end
+
+ #
+ # Creates a matrix where the diagonal elements are composed of +values+.
+ # Matrix.diagonal(9, 5, -3)
+ # => 9 0 0
+ # 0 5 0
+ # 0 0 -3
+ #
+ def Matrix.diagonal(*values)
+ size = values.size
+ rows = (0 .. size - 1).collect {
+ |j|
+ row = Array.new(size).fill(0, 0, size)
+ row[j] = values[j]
+ row
+ }
+ rows(rows, false)
+ end
+
+ #
+ # Creates an +n+ by +n+ diagonal matrix where each diagonal element is
+ # +value+.
+ # Matrix.scalar(2, 5)
+ # => 5 0
+ # 0 5
+ #
+ def Matrix.scalar(n, value)
+ Matrix.diagonal(*Array.new(n).fill(value, 0, n))
+ end
+
+ #
+ # Creates an +n+ by +n+ identity matrix.
+ # Matrix.identity(2)
+ # => 1 0
+ # 0 1
+ #
+ def Matrix.identity(n)
+ Matrix.scalar(n, 1)
+ end
+ class << Matrix
+ alias unit identity
+ alias I identity
+ end
+
+ #
+ # Creates an +n+ by +n+ zero matrix.
+ # Matrix.zero(2)
+ # => 0 0
+ # 0 0
+ #
+ def Matrix.zero(n)
+ Matrix.scalar(n, 0)
+ end
+
+ #
+ # Creates a single-row matrix where the values of that row are as given in
+ # +row+.
+ # Matrix.row_vector([4,5,6])
+ # => 4 5 6
+ #
+ def Matrix.row_vector(row)
+ case row
+ when Vector
+ Matrix.rows([row.to_a], false)
+ when Array
+ Matrix.rows([row.dup], false)
+ else
+ Matrix.rows([[row]], false)
+ end
+ end
+
+ #
+ # Creates a single-column matrix where the values of that column are as given
+ # in +column+.
+ # Matrix.column_vector([4,5,6])
+ # => 4
+ # 5
+ # 6
+ #
+ def Matrix.column_vector(column)
+ case column
+ when Vector
+ Matrix.columns([column.to_a])
+ when Array
+ Matrix.columns([column])
+ else
+ Matrix.columns([[column]])
+ end
+ end
+
+ #
+ # This method is used by the other methods that create matrices, and is of no
+ # use to general users.
+ #
+ def initialize(init_method, *argv)
+ self.send(init_method, *argv)
+ end
+
+ def init_rows(rows, copy)
+ if copy
+ @rows = rows.collect{|row| row.dup}
+ else
+ @rows = rows
+ end
+ self
+ end
+ private :init_rows
+
+ #
+ # Returns element (+i+,+j+) of the matrix. That is: row +i+, column +j+.
+ #
+ def [](i, j)
+ @rows[i][j]
+ end
+
+ #
+ # Returns the number of rows.
+ #
+ def row_size
+ @rows.size
+ end
+
+ #
+ # Returns the number of columns. Note that it is possible to construct a
+ # matrix with uneven columns (e.g. Matrix[ [1,2,3], [4,5] ]), but this is
+ # mathematically unsound. This method uses the first row to determine the
+ # result.
+ #
+ def column_size
+ @rows[0].size
+ end
+
+ #
+ # Returns row vector number +i+ of the matrix as a Vector (starting at 0 like
+ # an array). When a block is given, the elements of that vector are iterated.
+ #
+ def row(i) # :yield: e
+ if block_given?
+ for e in @rows[i]
+ yield e
+ end
+ else
+ Vector.elements(@rows[i])
+ end
+ end
+
+ #
+ # Returns column vector number +j+ of the matrix as a Vector (starting at 0
+ # like an array). When a block is given, the elements of that vector are
+ # iterated.
+ #
+ def column(j) # :yield: e
+ if block_given?
+ 0.upto(row_size - 1) do
+ |i|
+ yield @rows[i][j]
+ end
+ else
+ col = (0 .. row_size - 1).collect {
+ |i|
+ @rows[i][j]
+ }
+ Vector.elements(col, false)
+ end
+ end
+
+ #
+ # Returns a matrix that is the result of iteration of the given block over all
+ # elements of the matrix.
+ # Matrix[ [1,2], [3,4] ].collect { |i| i**2 }
+ # => 1 4
+ # 9 16
+ #
+ def collect # :yield: e
+ rows = @rows.collect{|row| row.collect{|e| yield e}}
+ Matrix.rows(rows, false)
+ end
+ alias map collect
+
+ #
+ # Returns a section of the matrix. The parameters are either:
+ # * start_row, nrows, start_col, ncols; OR
+ # * col_range, row_range
+ #
+ # Matrix.diagonal(9, 5, -3).minor(0..1, 0..2)
+ # => 9 0 0
+ # 0 5 0
+ #
+ def minor(*param)
+ case param.size
+ when 2
+ from_row = param[0].first
+ size_row = param[0].end - from_row
+ size_row += 1 unless param[0].exclude_end?
+ from_col = param[1].first
+ size_col = param[1].end - from_col
+ size_col += 1 unless param[1].exclude_end?
+ when 4
+ from_row = param[0]
+ size_row = param[1]
+ from_col = param[2]
+ size_col = param[3]
+ else
+ Matrix.Raise ArgumentError, param.inspect
+ end
+
+ rows = @rows[from_row, size_row].collect{
+ |row|
+ row[from_col, size_col]
+ }
+ Matrix.rows(rows, false)
+ end
+
+ #--
+ # TESTING -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
+ #++
+
+ #
+ # Returns +true+ if this is a regular matrix.
+ #
+ def regular?
+ square? and rank == column_size
+ end
+
+ #
+ # Returns +true+ is this is a singular (i.e. non-regular) matrix.
+ #
+ def singular?
+ not regular?
+ end
+
+ #
+ # Returns +true+ is this is a square matrix. See note in column_size about this
+ # being unreliable, though.
+ #
+ def square?
+ column_size == row_size
+ end
+
+ #--
+ # OBJECT METHODS -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
+ #++
+
+ #
+ # Returns +true+ if and only if the two matrices contain equal elements.
+ #
+ def ==(other)
+ return false unless Matrix === other
+
+ other.compare_by_row_vectors(@rows)
+ end
+ alias eql? ==
+
+ #
+ # Not really intended for general consumption.
+ #
+ def compare_by_row_vectors(rows)
+ return false unless @rows.size == rows.size
+
+ 0.upto(@rows.size - 1) do
+ |i|
+ return false unless @rows[i] == rows[i]
+ end
+ true
+ end
+
+ #
+ # Returns a clone of the matrix, so that the contents of each do not reference
+ # identical objects.
+ #
+ def clone
+ Matrix.rows(@rows)
+ end
+
+ #
+ # Returns a hash-code for the matrix.
+ #
+ def hash
+ value = 0
+ for row in @rows
+ for e in row
+ value ^= e.hash
+ end
+ end
+ return value
+ end
+
+ #--
+ # ARITHMETIC -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
+ #++
+
+ #
+ # Matrix multiplication.
+ # Matrix[[2,4], [6,8]] * Matrix.identity(2)
+ # => 2 4
+ # 6 8
+ #
+ def *(m) # m is matrix or vector or number
+ case(m)
+ when Numeric
+ rows = @rows.collect {
+ |row|
+ row.collect {
+ |e|
+ e * m
+ }
+ }
+ return Matrix.rows(rows, false)
+ when Vector
+ m = Matrix.column_vector(m)
+ r = self * m
+ return r.column(0)
+ when Matrix
+ Matrix.Raise ErrDimensionMismatch if column_size != m.row_size
+
+ rows = (0 .. row_size - 1).collect {
+ |i|
+ (0 .. m.column_size - 1).collect {
+ |j|
+ vij = 0
+ 0.upto(column_size - 1) do
+ |k|
+ vij += self[i, k] * m[k, j]
+ end
+ vij
+ }
+ }
+ return Matrix.rows(rows, false)
+ else
+ x, y = m.coerce(self)
+ return x * y
+ end
+ end
+
+ #
+ # Matrix addition.
+ # Matrix.scalar(2,5) + Matrix[[1,0], [-4,7]]
+ # => 6 0
+ # -4 12
+ #
+ def +(m)
+ case m
+ when Numeric
+ Matrix.Raise ErrOperationNotDefined, "+"
+ when Vector
+ m = Matrix.column_vector(m)
+ when Matrix
+ else
+ x, y = m.coerce(self)
+ return x + y
+ end
+
+ Matrix.Raise ErrDimensionMismatch unless row_size == m.row_size and column_size == m.column_size
+
+ rows = (0 .. row_size - 1).collect {
+ |i|
+ (0 .. column_size - 1).collect {
+ |j|
+ self[i, j] + m[i, j]
+ }
+ }
+ Matrix.rows(rows, false)
+ end
+
+ #
+ # Matrix subtraction.
+ # Matrix[[1,5], [4,2]] - Matrix[[9,3], [-4,1]]
+ # => -8 2
+ # 8 1
+ #
+ def -(m)
+ case m
+ when Numeric
+ Matrix.Raise ErrOperationNotDefined, "-"
+ when Vector
+ m = Matrix.column_vector(m)
+ when Matrix
+ else
+ x, y = m.coerce(self)
+ return x - y
+ end
+
+ Matrix.Raise ErrDimensionMismatch unless row_size == m.row_size and column_size == m.column_size
+
+ rows = (0 .. row_size - 1).collect {
+ |i|
+ (0 .. column_size - 1).collect {
+ |j|
+ self[i, j] - m[i, j]
+ }
+ }
+ Matrix.rows(rows, false)
+ end
+
+ #
+ # Matrix division (multiplication by the inverse).
+ # Matrix[[7,6], [3,9]] / Matrix[[2,9], [3,1]]
+ # => -7 1
+ # -3 -6
+ #
+ def /(other)
+ case other
+ when Numeric
+ rows = @rows.collect {
+ |row|
+ row.collect {
+ |e|
+ e / other
+ }
+ }
+ return Matrix.rows(rows, false)
+ when Matrix
+ return self * other.inverse
+ else
+ x, y = other.coerce(self)
+ rerurn x / y
+ end
+ end
+
+ #
+ # Returns the inverse of the matrix.
+ # Matrix[[1, 2], [2, 1]].inverse
+ # => -1 1
+ # 0 -1
+ #
+ def inverse
+ Matrix.Raise ErrDimensionMismatch unless square?
+ Matrix.I(row_size).inverse_from(self)
+ end
+ alias inv inverse
+
+ #
+ # Not for public consumption?
+ #
+ def inverse_from(src)
+ size = row_size - 1
+ a = src.to_a
+
+ for k in 0..size
+ if (akk = a[k][k]) == 0
+ i = k
+ begin
+ Matrix.Raise ErrNotRegular if (i += 1) > size
+ end while a[i][k] == 0
+ a[i], a[k] = a[k], a[i]
+ @rows[i], @rows[k] = @rows[k], @rows[i]
+ akk = a[k][k]
+ end
+
+ for i in 0 .. size
+ next if i == k
+ q = a[i][k] / akk
+ a[i][k] = 0
+
+ (k + 1).upto(size) do
+ |j|
+ a[i][j] -= a[k][j] * q
+ end
+ 0.upto(size) do
+ |j|
+ @rows[i][j] -= @rows[k][j] * q
+ end
+ end
+
+ (k + 1).upto(size) do
+ |j|
+ a[k][j] /= akk
+ end
+ 0.upto(size) do
+ |j|
+ @rows[k][j] /= akk
+ end
+ end
+ self
+ end
+ #alias reciprocal inverse
+
+ #
+ # Matrix exponentiation. Defined for integer powers only. Equivalent to
+ # multiplying the matrix by itself N times.
+ # Matrix[[7,6], [3,9]] ** 2
+ # => 67 96
+ # 48 99
+ #
+ def ** (other)
+ if other.kind_of?(Integer)
+ x = self
+ if other <= 0
+ x = self.inverse
+ return Matrix.identity(self.column_size) if other == 0
+ other = -other
+ end
+ z = x
+ n = other - 1
+ while n != 0
+ while (div, mod = n.divmod(2)
+ mod == 0)
+ x = x * x
+ n = div
+ end
+ z *= x
+ n -= 1
+ end
+ z
+ elsif other.kind_of?(Float) || defined?(Rational) && other.kind_of?(Rational)
+ Matrix.Raise ErrOperationNotDefined, "**"
+ else
+ Matrix.Raise ErrOperationNotDefined, "**"
+ end
+ end
+
+ #--
+ # MATRIX FUNCTIONS -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
+ #++
+
+ #
+ # Returns the determinant of the matrix. If the matrix is not square, the
+ # result is 0.
+ # Matrix[[7,6], [3,9]].determinant
+ # => 63
+ #
+ def determinant
+ return 0 unless square?
+
+ size = row_size - 1
+ a = to_a
+
+ det = 1
+ k = 0
+ begin
+ if (akk = a[k][k]) == 0
+ i = k
+ begin
+ return 0 if (i += 1) > size
+ end while a[i][k] == 0
+ a[i], a[k] = a[k], a[i]
+ akk = a[k][k]
+ det *= -1
+ end
+ (k + 1).upto(size) do
+ |i|
+ q = a[i][k] / akk
+ (k + 1).upto(size) do
+ |j|
+ a[i][j] -= a[k][j] * q
+ end
+ end
+ det *= akk
+ end while (k += 1) <= size
+ det
+ end
+ alias det determinant
+
+ #
+ # Returns the rank of the matrix. Beware that using Float values, with their
+ # usual lack of precision, can affect the value returned by this method. Use
+ # Rational values instead if this is important to you.
+ # Matrix[[7,6], [3,9]].rank
+ # => 2
+ #
+ def rank
+ if column_size > row_size
+ a = transpose.to_a
+ a_column_size = row_size
+ a_row_size = column_size
+ else
+ a = to_a
+ a_column_size = column_size
+ a_row_size = row_size
+ end
+ rank = 0
+ k = 0
+ begin
+ if (akk = a[k][k]) == 0
+ i = k
+ exists = true
+ begin
+ if (i += 1) > a_column_size - 1
+ exists = false
+ break
+ end
+ end while a[i][k] == 0
+ if exists
+ a[i], a[k] = a[k], a[i]
+ akk = a[k][k]
+ else
+ i = k
+ exists = true
+ begin
+ if (i += 1) > a_row_size - 1
+ exists = false
+ break
+ end
+ end while a[k][i] == 0
+ if exists
+ k.upto(a_column_size - 1) do
+ |j|
+ a[j][k], a[j][i] = a[j][i], a[j][k]
+ end
+ akk = a[k][k]
+ else
+ next
+ end
+ end
+ end
+ (k + 1).upto(a_row_size - 1) do
+ |i|
+ q = a[i][k] / akk
+ (k + 1).upto(a_column_size - 1) do
+ |j|
+ a[i][j] -= a[k][j] * q
+ end
+ end
+ rank += 1
+ end while (k += 1) <= a_column_size - 1
+ return rank
+ end
+
+ #
+ # Returns the trace (sum of diagonal elements) of the matrix.
+ # Matrix[[7,6], [3,9]].trace
+ # => 16
+ #
+ def trace
+ tr = 0
+ 0.upto(column_size - 1) do
+ |i|
+ tr += @rows[i][i]
+ end
+ tr
+ end
+ alias tr trace
+
+ #
+ # Returns the transpose of the matrix.
+ # Matrix[[1,2], [3,4], [5,6]]
+ # => 1 2
+ # 3 4
+ # 5 6
+ # Matrix[[1,2], [3,4], [5,6]].transpose
+ # => 1 3 5
+ # 2 4 6
+ #
+ def transpose
+ Matrix.columns(@rows)
+ end
+ alias t transpose
+
+ #--
+ # CONVERTING -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
+ #++
+
+ #
+ # FIXME: describe #coerce.
+ #
+ def coerce(other)
+ case other
+ when Numeric
+ return Scalar.new(other), self
+ else
+ raise TypeError, "#{self.class} can't be coerced into #{other.class}"
+ end
+ end
+
+ #
+ # Returns an array of the row vectors of the matrix. See Vector.
+ #
+ def row_vectors
+ rows = (0 .. row_size - 1).collect {
+ |i|
+ row(i)
+ }
+ rows
+ end
+
+ #
+ # Returns an array of the column vectors of the matrix. See Vector.
+ #
+ def column_vectors
+ columns = (0 .. column_size - 1).collect {
+ |i|
+ column(i)
+ }
+ columns
+ end
+
+ #
+ # Returns an array of arrays that describe the rows of the matrix.
+ #
+ def to_a
+ @rows.collect{|row| row.collect{|e| e}}
+ end
+
+ #--
+ # PRINTING -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
+ #++
+
+ #
+ # Overrides Object#to_s
+ #
+ def to_s
+ "Matrix[" + @rows.collect{
+ |row|
+ "[" + row.collect{|e| e.to_s}.join(", ") + "]"
+ }.join(", ")+"]"
+ end
+
+ #
+ # Overrides Object#inspect
+ #
+ def inspect
+ "Matrix"+@rows.inspect
+ end
+
+ # Private CLASS
+
+ class Scalar < Numeric # :nodoc:
+ include ExceptionForMatrix
+
+ def initialize(value)
+ @value = value
+ end
+
+ # ARITHMETIC
+ def +(other)
+ case other
+ when Numeric
+ Scalar.new(@value + other)
+ when Vector, Matrix
+ Scalar.Raise WrongArgType, other.class, "Numeric or Scalar"
+ when Scalar
+ Scalar.new(@value + other.value)
+ else
+ x, y = other.coerce(self)
+ x + y
+ end
+ end
+
+ def -(other)
+ case other
+ when Numeric
+ Scalar.new(@value - other)
+ when Vector, Matrix
+ Scalar.Raise WrongArgType, other.class, "Numeric or Scalar"
+ when Scalar
+ Scalar.new(@value - other.value)
+ else
+ x, y = other.coerce(self)
+ x - y
+ end
+ end
+
+ def *(other)
+ case other
+ when Numeric
+ Scalar.new(@value * other)
+ when Vector, Matrix
+ other.collect{|e| @value * e}
+ else
+ x, y = other.coerce(self)
+ x * y
+ end
+ end
+
+ def / (other)
+ case other
+ when Numeric
+ Scalar.new(@value / other)
+ when Vector
+ Scalar.Raise WrongArgType, other.class, "Numeric or Scalar or Matrix"
+ when Matrix
+ self * _M.inverse
+ else
+ x, y = other.coerce(self)
+ x / y
+ end
+ end
+
+ def ** (other)
+ case other
+ when Numeric
+ Scalar.new(@value ** other)
+ when Vector
+ Scalar.Raise WrongArgType, other.class, "Numeric or Scalar or Matrix"
+ when Matrix
+ other.powered_by(self)
+ else
+ x, y = other.coerce(self)
+ x ** y
+ end
+ end
+ end
+end
+
+
+#
+# The +Vector+ class represents a mathematical vector, which is useful in its own right, and
+# also constitutes a row or column of a Matrix.
+#
+# == Method Catalogue
+#
+# To create a Vector:
+# * <tt> Vector.[](*array) </tt>
+# * <tt> Vector.elements(array, copy = true) </tt>
+#
+# To access elements:
+# * <tt> [](i) </tt>
+#
+# To enumerate the elements:
+# * <tt> #each2(v) </tt>
+# * <tt> #collect2(v) </tt>
+#
+# Vector arithmetic:
+# * <tt> *(x) "is matrix or number" </tt>
+# * <tt> +(v) </tt>
+# * <tt> -(v) </tt>
+#
+# Vector functions:
+# * <tt> #inner_product(v) </tt>
+# * <tt> #collect </tt>
+# * <tt> #map </tt>
+# * <tt> #map2(v) </tt>
+# * <tt> #r </tt>
+# * <tt> #size </tt>
+#
+# Conversion to other data types:
+# * <tt> #covector </tt>
+# * <tt> #to_a </tt>
+# * <tt> #coerce(other) </tt>
+#
+# String representations:
+# * <tt> #to_s </tt>
+# * <tt> #inspect </tt>
+#
+class Vector
+ include ExceptionForMatrix
+
+ #INSTANCE CREATION
+
+ private_class_method :new
+
+ #
+ # Creates a Vector from a list of elements.
+ # Vector[7, 4, ...]
+ #
+ def Vector.[](*array)
+ new(:init_elements, array, copy = false)
+ end
+
+ #
+ # Creates a vector from an Array. The optional second argument specifies
+ # whether the array itself or a copy is used internally.
+ #
+ def Vector.elements(array, copy = true)
+ new(:init_elements, array, copy)
+ end
+
+ #
+ # For internal use.
+ #
+ def initialize(method, array, copy)
+ self.send(method, array, copy)
+ end
+
+ #
+ # For internal use.
+ #
+ def init_elements(array, copy)
+ if copy
+ @elements = array.dup
+ else
+ @elements = array
+ end
+ end
+
+ # ACCESSING
+
+ #
+ # Returns element number +i+ (starting at zero) of the vector.
+ #
+ def [](i)
+ @elements[i]
+ end
+
+ #
+ # Returns the number of elements in the vector.
+ #
+ def size
+ @elements.size
+ end
+
+ #--
+ # ENUMERATIONS -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
+ #++
+
+ #
+ # Iterate over the elements of this vector and +v+ in conjunction.
+ #
+ def each2(v) # :yield: e1, e2
+ Vector.Raise ErrDimensionMismatch if size != v.size
+ 0.upto(size - 1) do
+ |i|
+ yield @elements[i], v[i]
+ end
+ end
+
+ #
+ # Collects (as in Enumerable#collect) over the elements of this vector and +v+
+ # in conjunction.
+ #
+ def collect2(v) # :yield: e1, e2
+ Vector.Raise ErrDimensionMismatch if size != v.size
+ (0 .. size - 1).collect do
+ |i|
+ yield @elements[i], v[i]
+ end
+ end
+
+ #--
+ # COMPARING -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
+ #++
+
+ #
+ # Returns +true+ iff the two vectors have the same elements in the same order.
+ #
+ def ==(other)
+ return false unless Vector === other
+
+ other.compare_by(@elements)
+ end
+ alias eqn? ==
+
+ #
+ # For internal use.
+ #
+ def compare_by(elements)
+ @elements == elements
+ end
+
+ #
+ # Return a copy of the vector.
+ #
+ def clone
+ Vector.elements(@elements)
+ end
+
+ #
+ # Return a hash-code for the vector.
+ #
+ def hash
+ @elements.hash
+ end
+
+ #--
+ # ARITHMETIC -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
+ #++
+
+ #
+ # Multiplies the vector by +x+, where +x+ is a number or another vector.
+ #
+ def *(x)
+ case x
+ when Numeric
+ els = @elements.collect{|e| e * x}
+ Vector.elements(els, false)
+ when Matrix
+ Matrix.column_vector(self) * x
+ else
+ s, x = x.coerce(self)
+ s * x
+ end
+ end
+
+ #
+ # Vector addition.
+ #
+ def +(v)
+ case v
+ when Vector
+ Vector.Raise ErrDimensionMismatch if size != v.size
+ els = collect2(v) {
+ |v1, v2|
+ v1 + v2
+ }
+ Vector.elements(els, false)
+ when Matrix
+ Matrix.column_vector(self) + v
+ else
+ s, x = v.coerce(self)
+ s + x
+ end
+ end
+
+ #
+ # Vector subtraction.
+ #
+ def -(v)
+ case v
+ when Vector
+ Vector.Raise ErrDimensionMismatch if size != v.size
+ els = collect2(v) {
+ |v1, v2|
+ v1 - v2
+ }
+ Vector.elements(els, false)
+ when Matrix
+ Matrix.column_vector(self) - v
+ else
+ s, x = v.coerce(self)
+ s - x
+ end
+ end
+
+ #--
+ # VECTOR FUNCTIONS -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
+ #++
+
+ #
+ # Returns the inner product of this vector with the other.
+ # Vector[4,7].inner_product Vector[10,1] => 47
+ #
+ def inner_product(v)
+ Vector.Raise ErrDimensionMismatch if size != v.size
+
+ p = 0
+ each2(v) {
+ |v1, v2|
+ p += v1 * v2
+ }
+ p
+ end
+
+ #
+ # Like Array#collect.
+ #
+ def collect # :yield: e
+ els = @elements.collect {
+ |v|
+ yield v
+ }
+ Vector.elements(els, false)
+ end
+ alias map collect
+
+ #
+ # Like Vector#collect2, but returns a Vector instead of an Array.
+ #
+ def map2(v) # :yield: e1, e2
+ els = collect2(v) {
+ |v1, v2|
+ yield v1, v2
+ }
+ Vector.elements(els, false)
+ end
+
+ #
+ # Returns the modulus (Pythagorean distance) of the vector.
+ # Vector[5,8,2].r => 9.643650761
+ #
+ def r
+ v = 0
+ for e in @elements
+ v += e*e
+ end
+ return Math.sqrt(v)
+ end
+
+ #--
+ # CONVERTING
+ #++
+
+ #
+ # Creates a single-row matrix from this vector.
+ #
+ def covector
+ Matrix.row_vector(self)
+ end
+
+ #
+ # Returns the elements of the vector in an array.
+ #
+ def to_a
+ @elements.dup
+ end
+
+ #
+ # FIXME: describe Vector#coerce.
+ #
+ def coerce(other)
+ case other
+ when Numeric
+ return Scalar.new(other), self
+ else
+ raise TypeError, "#{self.class} can't be coerced into #{other.class}"
+ end
+ end
+
+ #--
+ # PRINTING -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
+ #++
+
+ #
+ # Overrides Object#to_s
+ #
+ def to_s
+ "Vector[" + @elements.join(", ") + "]"
+ end
+
+ #
+ # Overrides Object#inspect
+ #
+ def inspect
+ str = "Vector"+@elements.inspect
+ end
+end
+
+
+# Documentation comments:
+# - Matrix#coerce and Vector#coerce need to be documented