Daniel-Diaz / matrix

A Haskell native implementation of matrices and their operations.
BSD 3-Clause "New" or "Revised" License
35 stars 31 forks source link

rref bug #42

Open rzil opened 8 years ago

rzil commented 8 years ago

The following commands result in unexpected behaviour with matrix-0.3.5.1

let mat = fromList 9 10 [(-1) % 1,1 % 1,0 % 1,0 % 1,0 % 1,0 % 1,0 % 1,0 % 1,1 % 1,1 % 1,1 % 1,(-1) % 1,0 % 1,0 % 1,0 % 1,0 % 1,1 % 1,0 % 1,1 % 1,1 % 1,1 % 1,1 % 1,(-1) % 1,0 % 1,0 % 1,0 % 1,0 % 1,1 % 1,0 % 1,1 % 1,0 % 1,0 % 1,1 % 1,(-1) % 1,0 % 1,0 % 1,0 % 1,0 % 1,0 % 1,1 % 1,0 % 1,0 % 1,0 % 1,1 % 1,(-1) % 1,1 % 1,0 % 1,0 % 1,0 % 1,1 % 1,0 % 1,0 % 1,0 % 1,0 % 1,0 % 1,(-1) % 1,1 % 1,0 % 1,0 % 1,1 % 1,0 % 1,0 % 1,0 % 1,0 % 1,0 % 1,0 % 1,(-1) % 1,0 % 1,1 % 1,1 % 1,0 % 1,0 % 1,0 % 1,0 % 1,0 % 1,0 % 1,0 % 1,(-1) % 1,0 % 1,1 % 1,0 % 1,0 % 1,0 % 1,0 % 1,0 % 1,0 % 1,0 % 1,0 % 1,(-1) % 1,1 % 1]

rref mat

Right ┌ *** Exception: getElem: Trying to get the (6,1) element from a 5x6 matrix.
CallStack (from HasCallStack):
  error, called at ./Data/Matrix.hs:434:16 in matrix-0.3.5.1-1pdOyADkVl9H4vOxaM7loN:Data.Matrix

NB: the correct row reduced form is given below (computed with sage)

[ 1  1  0  0 -1  0  0  0  0  1]
[ 0  2  0  0 -1  0  0  0  0  3]
[ 0  0  1  0 -1  0  0  0  0  5]
[ 0  0  0  1 -1  0  0  0  0  4]
[ 0  0  0  0  0  1  0  0  0  3]
[ 0  0  0  0  0  0  1  0  0  4]
[ 0  0  0  0  0  0  0  1  0  5]
[ 0  0  0  0  0  0  0  0  1  5]
[ 0  0  0  0  0  0  0  0  0  6]
wchresta commented 6 years ago
octave:8> output_precision(1)
octave:9> A=[-1,1,0,0,0,0,0,0,1,1;1,-1,0,0,0,0,1,0,1,1;1,1,-1,0,0,0,0,1,0,1;0,0,1,-1,0,0,0,0,0,1;0,0,0,1,-1,1,0,0,0,1;0,0,0,0,0,-1,1,0,0,1;0,0,0,0,0,0,-1,0,1,1;0,0,0,0,0,0,0,-1,0,1;0,0,0,0,0,0,0,0,-1,1]
A =

  -1   1   0   0   0   0   0   0   1   1
   1  -1   0   0   0   0   1   0   1   1
   1   1  -1   0   0   0   0   1   0   1
   0   0   1  -1   0   0   0   0   0   1
   0   0   0   1  -1   1   0   0   0   1
   0   0   0   0   0  -1   1   0   0   1
   0   0   0   0   0   0  -1   0   1   1
   0   0   0   0   0   0   0  -1   0   1
   0   0   0   0   0   0   0   0  -1   1

octave:10> rref(A)
ans =

   1.0   0.0   0.0   0.0  -0.5   0.0   0.0   0.0   0.0   0.0
   0.0   1.0   0.0   0.0  -0.5   0.0   0.0   0.0   0.0   0.0
   0.0   0.0   1.0   0.0  -1.0   0.0   0.0   0.0   0.0   0.0
   0.0   0.0   0.0   1.0  -1.0   0.0   0.0   0.0   0.0   0.0
   0.0   0.0   0.0   0.0   0.0   1.0   0.0  -0.0   0.0   0.0
   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0
   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0  -0.0   0.0
   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0
   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0

Seems octave doesn't agree with sage. Reason is, what you provided is the echelon form; not the reduced echelon form. One can clearly see that in your matrix, the last column's elements of rows 1-9 can be easily eliminated by the last row multiplying by ?/6.