matrix_bareiss.diff
| b/lib/matrix.rb | ||
|---|---|---|
| 661 | 661 |
#++ |
| 662 | 662 | |
| 663 | 663 |
# |
| 664 |
# Returns the determinant of the matrix. If the matrix is not square, the |
|
| 665 |
# result is 0. This method's algorithm is Gaussian elimination method |
|
| 666 |
# and using Numeric#quo(). Beware that using Float values, with their |
|
| 667 |
# usual lack of precision, can affect the value returned by this method. Use |
|
| 668 |
# Rational values or Matrix#det_e instead if this is important to you. |
|
| 664 |
# Returns the determinant of the matrix. |
|
| 665 |
# This method's algorithm is Bareiss algorithm, which has the nice properties |
|
| 666 |
# of being exact for integers and minimizing the rounding error for floats. |
|
| 669 | 667 |
# |
| 670 | 668 |
# Matrix[[7,6], [3,9]].determinant |
| 671 |
# => 45.0
|
|
| 669 |
# => 45 |
|
| 672 | 670 |
# |
| 673 | 671 |
def determinant |
| 674 |
return 0 unless square?
|
|
| 672 |
Matrix.Raise ErrDimensionMismatch unless square?
|
|
| 675 | 673 | |
| 676 | 674 |
size = row_size |
| 675 |
last = row_size - 1 |
|
| 677 | 676 |
a = to_a |
| 678 | 677 | |
| 679 | 678 |
det = 1 |
| 679 |
last_akk = 1 |
|
| 680 | 680 |
size.times do |k| |
| 681 | 681 |
if (akk = a[k][k]) == 0 |
| 682 | 682 |
i = (k+1 ... size).find {|ii|
|
| ... | ... | |
| 685 | 685 |
return 0 if i.nil? |
| 686 | 686 |
a[i], a[k] = a[k], a[i] |
| 687 | 687 |
akk = a[k][k] |
| 688 |
det *= -1
|
|
| 688 |
det = -det
|
|
| 689 | 689 |
end |
| 690 | 690 | |
| 691 |
(k + 1 ... size).each do |ii| |
|
| 692 |
q = a[ii][k].quo(akk) |
|
| 693 |
(k + 1 ... size).each do |j| |
|
| 694 |
a[ii][j] -= a[k][j] * q |
|
| 691 |
(k+1).upto(last) do |ii| |
|
| 692 |
last.downto(k) do |j| |
|
| 693 |
a[ii][j] = (akk * a[ii][j] - a[ii][k] * a[k][j]) / last_akk |
|
| 695 | 694 |
end |
| 696 | 695 |
end |
| 697 |
det *= akk |
|
| 696 |
det = det / last_akk * akk |
|
| 697 |
last_akk = akk |
|
| 698 | 698 |
end |
| 699 | 699 |
det |
| 700 | 700 |
end |