Feature #2772 » matrix_bareiss.diff
| lib/matrix.rb | ||
|---|---|---|
|
#++
|
||
|
#
|
||
|
# Returns the determinant of the matrix. If the matrix is not square, the
|
||
|
# result is 0. 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.
|
||
|
# Returns the determinant of the matrix.
|
||
|
# This method's algorithm is Bareiss algorithm, which has the nice properties
|
||
|
# of being exact for integers and minimizing the rounding error for floats.
|
||
|
#
|
||
|
# Matrix[[7,6], [3,9]].determinant
|
||
|
# => 45.0
|
||
|
# => 45
|
||
|
#
|
||
|
def determinant
|
||
|
return 0 unless square?
|
||
|
Matrix.Raise ErrDimensionMismatch unless square?
|
||
|
size = row_size
|
||
|
last = row_size - 1
|
||
|
a = to_a
|
||
|
det = 1
|
||
|
last_akk = 1
|
||
|
size.times do |k|
|
||
|
if (akk = a[k][k]) == 0
|
||
|
i = (k+1 ... size).find {|ii|
|
||
| ... | ... | |
|
return 0 if i.nil?
|
||
|
a[i], a[k] = a[k], a[i]
|
||
|
akk = a[k][k]
|
||
|
det *= -1
|
||
|
det = -det
|
||
|
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 |ii|
|
||
|
last.downto(k) do |j|
|
||
|
a[ii][j] = (akk * a[ii][j] - a[ii][k] * a[k][j]) / last_akk
|
||
|
end
|
||
|
end
|
||
|
det *= akk
|
||
|
det = det / last_akk * akk
|
||
|
last_akk = akk
|
||
|
end
|
||
|
det
|
||
|
end
|
||