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. nonsingular) matrix.


#


def regular?


square? and rank == column_size


not singular?


end


#


# Returns +true+ is this is a singular (i.e. nonregular) 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 GaussBareiss algorithm


determinant_bareiss


end


end


alias_method :det, :determinant


#


# Private. Use Matrix#determinant


#


# Returns the determinant of the matrix, using


# Bareiss' multistep integerpreserving 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 integerpreserving 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 »