Feature #2772 » newdet.patch
| lib/matrix.rb | ||
|---|---|---|
|
end
|
||
|
#
|
||
|
# Returns +true+ if this is a regular matrix.
|
||
|
# Returns +true+ if this is a regular (i.e. non-singular) matrix.
|
||
|
#
|
||
|
def regular?
|
||
|
square? and rank == column_size
|
||
|
not singular?
|
||
|
end
|
||
|
#
|
||
|
# Returns +true+ is this is a singular (i.e. non-regular) matrix.
|
||
|
# Returns +true+ is this is a singular matrix.
|
||
|
#
|
||
|
def singular?
|
||
|
not regular?
|
||
|
determinant == 0
|
||
|
end
|
||
|
#
|
||
| ... | ... | |
|
#
|
||
|
# Returns the determinant of the matrix.
|
||
|
# This method's algorithm is Gaussian elimination method
|
||
|
# and using Numeric#quo(). Beware that using Float values, with their
|
||
|
# usual lack of precision, can affect the value returned by this method. Use
|
||
|
# Rational values or Matrix#det_e instead if this is important to you.
|
||
|
#
|
||
|
# Beware that using Float values can yield erroneous results
|
||
|
# because of their lack of precision.
|
||
|
# Consider using exact types like Rational or BigDecimal instead.
|
||
|
#
|
||
|
# Matrix[[7,6], [3,9]].determinant
|
||
|
# => 45.0
|
||
|
# => 45
|
||
|
#
|
||
|
def determinant
|
||
|
Matrix.Raise ErrDimensionMismatch unless square?
|
||
|
m = @rows
|
||
|
case row_size
|
||
|
# Up to 4x4, give result using Laplacian expansion by minors.
|
||
|
# This will typically be faster, as well as giving good results
|
||
|
# in case of Floats
|
||
|
when 0
|
||
|
+1
|
||
|
when 1
|
||
|
+ m[0][0]
|
||
|
when 2
|
||
|
+ m[0][0] * m[1][1] - m[0][1] * m[1][0]
|
||
|
when 3
|
||
|
m0 = m[0]; m1 = m[1]; m2 = m[2]
|
||
|
+ m0[0] * m1[1] * m2[2] - m0[0] * m1[2] * m2[1] \
|
||
|
- m0[1] * m1[0] * m2[2] + m0[1] * m1[2] * m2[0] \
|
||
|
+ m0[2] * m1[0] * m2[1] - m0[2] * m1[1] * m2[0]
|
||
|
when 4
|
||
|
m0 = m[0]; m1 = m[1]; m2 = m[2]; m3 = m[3]
|
||
|
+ m0[0] * m1[1] * m2[2] * m3[3] - m0[0] * m1[1] * m2[3] * m3[2] \
|
||
|
- m0[0] * m1[2] * m2[1] * m3[3] + m0[0] * m1[2] * m2[3] * m3[1] \
|
||
|
+ m0[0] * m1[3] * m2[1] * m3[2] - m0[0] * m1[3] * m2[2] * m3[1] \
|
||
|
- m0[1] * m1[0] * m2[2] * m3[3] + m0[1] * m1[0] * m2[3] * m3[2] \
|
||
|
+ m0[1] * m1[2] * m2[0] * m3[3] - m0[1] * m1[2] * m2[3] * m3[0] \
|
||
|
- m0[1] * m1[3] * m2[0] * m3[2] + m0[1] * m1[3] * m2[2] * m3[0] \
|
||
|
+ m0[2] * m1[0] * m2[1] * m3[3] - m0[2] * m1[0] * m2[3] * m3[1] \
|
||
|
- m0[2] * m1[1] * m2[0] * m3[3] + m0[2] * m1[1] * m2[3] * m3[0] \
|
||
|
+ m0[2] * m1[3] * m2[0] * m3[1] - m0[2] * m1[3] * m2[1] * m3[0] \
|
||
|
- m0[3] * m1[0] * m2[1] * m3[2] + m0[3] * m1[0] * m2[2] * m3[1] \
|
||
|
+ m0[3] * m1[1] * m2[0] * m3[2] - m0[3] * m1[1] * m2[2] * m3[0] \
|
||
|
- m0[3] * m1[2] * m2[0] * m3[1] + m0[3] * m1[2] * m2[1] * m3[0]
|
||
|
else
|
||
|
# For bigger matrices, use an efficient and general algorithm.
|
||
|
# Currently, we use the Gauss-Bareiss algorithm
|
||
|
determinant_bareiss
|
||
|
end
|
||
|
end
|
||
|
alias_method :det, :determinant
|
||
|
#
|
||
|
# Private. Use Matrix#determinant
|
||
|
#
|
||
|
# Returns the determinant of the matrix, using
|
||
|
# Bareiss' multistep integer-preserving gaussian elimination.
|
||
|
# It has the same computational cost order O(n^3) as standard Gaussian elimination.
|
||
|
# Intermediate results are fraction free and of lower complexity.
|
||
|
# A matrix of Integers will have thus intermediate results that are also Integers,
|
||
|
# with smaller bignums (if any), while a matrix of Float will usually have more
|
||
|
# precise intermediate results.
|
||
|
#
|
||
|
def determinant_bareiss
|
||
|
size = row_size
|
||
|
last = size - 1
|
||
|
a = to_a
|
||
|
det = 1
|
||
|
no_pivot = Proc.new{ return 0 }
|
||
|
sign = +1
|
||
|
pivot = 1
|
||
|
size.times do |k|
|
||
|
if (akk = a[k][k]) == 0
|
||
|
i = (k+1 ... size).find {|ii|
|
||
|
a[ii][k] != 0
|
||
|
previous_pivot = pivot
|
||
|
if (pivot = a[k][k]) == 0
|
||
|
switch = (k+1 ... size).find(no_pivot) {|row|
|
||
|
a[row][k] != 0
|
||
|
}
|
||
|
return 0 if i.nil?
|
||
|
a[i], a[k] = a[k], a[i]
|
||
|
akk = a[k][k]
|
||
|
det *= -1
|
||
|
a[switch], a[k] = a[k], a[switch]
|
||
|
pivot = a[k][k]
|
||
|
sign = -sign
|
||
|
end
|
||
|
(k + 1 ... size).each do |ii|
|
||
|
q = a[ii][k].quo(akk)
|
||
|
(k + 1 ... size).each do |j|
|
||
|
a[ii][j] -= a[k][j] * q
|
||
|
(k+1).upto(last) do |i|
|
||
|
ai = a[i]
|
||
|
(k+1).upto(last) do |j|
|
||
|
ai[j] = (pivot * ai[j] - ai[k] * a[k][j]) / previous_pivot
|
||
|
end
|
||
|
end
|
||
|
det *= akk
|
||
|
end
|
||
|
det
|
||
|
sign * pivot
|
||
|
end
|
||
|
alias det determinant
|
||
|
private :determinant_bareiss
|
||
|
#
|
||
|
# Returns the determinant of the matrix.
|
||
|
# This method's algorithm is Gaussian elimination method.
|
||
|
# This method uses Euclidean algorithm. If all elements are integer,
|
||
|
# really exact value. But, if an element is a float, can't return
|
||
|
# exact value.
|
||
|
#
|
||
|
# Matrix[[7,6], [3,9]].determinant
|
||
|
# => 63
|
||
|
# deprecated; use Matrix#determinant
|
||
|
#
|
||
|
def determinant_e
|
||
|
Matrix.Raise ErrDimensionMismatch unless square?
|
||
|
size = row_size
|
||
|
a = to_a
|
||
|
det = 1
|
||
|
size.times do |k|
|
||
|
if a[k][k].zero?
|
||
|
i = (k+1 ... size).find {|ii|
|
||
|
a[ii][k] != 0
|
||
|
}
|
||
|
return 0 if i.nil?
|
||
|
a[i], a[k] = a[k], a[i]
|
||
|
det *= -1
|
||
|
end
|
||
|
(k + 1 ... size).each do |ii|
|
||
|
q = a[ii][k].quo(a[k][k])
|
||
|
(k ... size).each do |j|
|
||
|
a[ii][j] -= a[k][j] * q
|
||
|
end
|
||
|
unless a[ii][k].zero?
|
||
|
a[ii], a[k] = a[k], a[ii]
|
||
|
det *= -1
|
||
|
redo
|
||
|
end
|
||
|
end
|
||
|
det *= a[k][k]
|
||
|
end
|
||
|
det
|
||
|
warn "#{caller(1)[0]}: warning: Matrix#determinant_e is deprecated; use #determinant"
|
||
|
rank
|
||
|
end
|
||
|
alias det_e determinant_e
|
||
|
#
|
||
|
# Returns the rank of the matrix. Beware that using Float values,
|
||
|
# probably return faild value. Use Rational values or Matrix#rank_e
|
||
|
# for getting exact result.
|
||
|
# Returns the rank of the matrix.
|
||
|
# Beware that using Float values can yield erroneous results
|
||
|
# because of their lack of precision.
|
||
|
# Consider using exact types like Rational or BigDecimal instead.
|
||
|
#
|
||
|
# 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
|
||
|
# We currently use Bareiss' multistep integer-preserving gaussian elimination
|
||
|
# (see comments on determinant)
|
||
|
a = to_a
|
||
|
last_column = column_size - 1
|
||
|
last_row = row_size - 1
|
||
|
rank = 0
|
||
|
a_column_size.times do |k|
|
||
|
if (akk = a[k][k]) == 0
|
||
|
i = (k+1 ... a_row_size).find {|ii|
|
||
|
a[ii][k] != 0
|
||
|
}
|
||
|
if i
|
||
|
a[i], a[k] = a[k], a[i]
|
||
|
akk = a[k][k]
|
||
|
else
|
||
|
i = (k+1 ... a_column_size).find {|ii|
|
||
|
a[k][ii] != 0
|
||
|
}
|
||
|
next if i.nil?
|
||
|
(k ... a_column_size).each do |j|
|
||
|
a[j][k], a[j][i] = a[j][i], a[j][k]
|
||
|
end
|
||
|
akk = a[k][k]
|
||
|
end
|
||
|
pivot_row = 0
|
||
|
previous_pivot = 1
|
||
|
0.upto(last_column) do |k|
|
||
|
switch_row = (pivot_row .. last_row).find {|row|
|
||
|
a[row][k] != 0
|
||
|
}
|
||
|
if switch_row
|
||
|
a[switch_row], a[pivot_row] = a[pivot_row], a[switch_row] unless pivot_row == switch_row
|
||
|
pivot = a[pivot_row][k]
|
||
|
(pivot_row+1).upto(last_row) do |i|
|
||
|
ai = a[i]
|
||
|
(k+1).upto(last_column) do |j|
|
||
|
ai[j] = (pivot * ai[j] - ai[k] * a[pivot_row][j]) / previous_pivot
|
||
|
end
|
||
|
end
|
||
|
pivot_row += 1
|
||
|
previous_pivot = pivot
|
||
|
end
|
||
|
(k + 1 ... a_row_size).each do |ii|
|
||
|
q = a[ii][k].quo(akk)
|
||
|
(k + 1... a_column_size).each do |j|
|
||
|
a[ii][j] -= a[k][j] * q
|
||
|
end
|
||
|
end
|
||
|
rank += 1
|
||
|
end
|
||
|
return rank
|
||
|
pivot_row
|
||
|
end
|
||
|
#
|
||
|
# Returns the rank of the matrix. This method uses Euclidean
|
||
|
# algorithm. If all elements are integer, really exact value. But, if
|
||
|
# an element is a float, can't return exact value.
|
||
|
#
|
||
|
# Matrix[[7,6], [3,9]].rank
|
||
|
# => 2
|
||
|
# deprecated; use Matrix#rank
|
||
|
#
|
||
|
def rank_e
|
||
|
a = to_a
|
||
|
a_column_size = column_size
|
||
|
a_row_size = row_size
|
||
|
pi = 0
|
||
|
a_column_size.times do |j|
|
||
|
if i = (pi ... a_row_size).find{|i0| !a[i0][j].zero?}
|
||
|
if i != pi
|
||
|
a[pi], a[i] = a[i], a[pi]
|
||
|
end
|
||
|
(pi + 1 ... a_row_size).each do |k|
|
||
|
q = a[k][j].quo(a[pi][j])
|
||
|
(pi ... a_column_size).each do |j0|
|
||
|
a[k][j0] -= q * a[pi][j0]
|
||
|
end
|
||
|
if k > pi && !a[k][j].zero?
|
||
|
a[k], a[pi] = a[pi], a[k]
|
||
|
redo
|
||
|
end
|
||
|
end
|
||
|
pi += 1
|
||
|
end
|
||
|
end
|
||
|
pi
|
||
|
warn "#{caller(1)[0]}: warning: Matrix#rank_e is deprecated; use #rank"
|
||
|
rank
|
||
|
end
|
||
- « Previous
- 1
- 2
- 3
- Next »