[Haskell-cafe] assistance with using haskell to calculate the determinant

[Daniel Fischer]

Am Sonntag, 27. März 2005 07:04 schrieb Carter Tazio Schonwald:
> I've written a simple program to calculate determinants of matrices, and
> while it executes to completion
> it fails to yield a correct answer. (and yes, I know I'm using lists to
> represent matrices)
>
> I would greatly appreciate any assistance in correcting this algorithm
>
> --Carter Schonwald
>
>
> import Prelude
 ^^^^^^^^^^^^^^^
is automatically done

>
>
> first [] = []
> first (a:b) = a
>
> nth::Int->[a]->a
> nth _ [] = error "list is too small"
> nth 1 (a:b) = a
> nth n (a:b)  = nth (n-1) b
>
> takeAfter::Int->[a]->[a]
> takeAfter _ [] = error "list too small"
> takeAfter 0 a= a
> takeAfter  1  (a:b) = b
> takeAfter  n (a:b)  = takeAfter  ( n-1)  b
^^^^^^^^^^
this is just 'drop'
>
>
>
> type Matrix = [[Rational]]
>
>
> pad a  = [x++x| x<- a]
>
>
> time2 []  _ = []
> time2  _ [] = []
> time2 (a:b) (c:d) = (a  * c):(time2 b d)
^^^^^^^^^
zipWith (*)
>
>
> tupleProd (a,b) = a * b
^^^^^^^^^^^
unused
>
>
> altSign a = [b* (-1^num) |  b<-a, num <- [2..]]
^^^^^^^^^
first real problem: this does absolutely not what you want,
take 4 (altSign [1,2])  == [-1,-1,-1,-1].
What you want is altSign [1,2] == [1,-2], which you get by

altSign = zipWith (*) [(-1)^n | n <- [0 .. ]]

or, with less calculation:

altSign = zipWith (*) (cycle [1,-1])


>
> index::Int->Int->Matrix->Rational
> index a b c  = nth b (nth a c) --- ath row, bth column
>
> slice::(Int,Int)->(Int,Int)->Matrix->Matrix
> slice (a,b) (c,d) list =   let rowSliceFront = takeAfter  (a-1) list
> ---  a and c are rows, b and c are columns
>                 in let rowSlice = take  (c-a+1) list
                                                                ^^^^^^
should be rowSliceFront
>                   in let columnSliceFront = map (takeAfter (b-1)) rowSlice
>                     in map (take  (d-b+1)) columnSliceFront
>
this would be more readable with layout:

slice (i,j) (k,l) mat
     = let rowSliceFront = drop (i-1) mat
             rowSlice = take (k-i+1) rowSliceFront
             columnSliceFront = map (drop (j-1)) rowSlice
         in map (take (l-j+1)) columnSliceFront

>
> determinant::Matrix->Rational
> determinant a = det (pad a) (length a)
>
>   -- only called from determinant
> det::Matrix->Int->Rational
> det a 1 = index 1 1 a
> det a size =let
>         coeffs = altSign (first (take 1 a))
>         newsize = size - 1
>         leastToMax = newsize - 1
>         slices = [ slice (2,  i) (size, i + leastToMax) a | i <-
> [2..(size+1)]  ]
>           in
>           sum (time2  coeffs (map (\l-> det l newsize) (map (pad)
> slices))  )
>
here is the next problem, you get the wrong signs here, for the 'slices' for 
the 3x3-matrix

1 2 3
6 5 4
2 4 2

will be 
5 4
4 2 -correct,

4 6
2 2 -which ought to be

6 4
2 2

and, correct again
6 5
2 4.

So some of the subdeterminants are summed with the wrong sign.

I fixed that in the attched code, clinging close to your version.
However I think it's not very efficient.

Hope to have been of assistance,

Daniel
>
>
>
> _______________________________________________
> Haskell-Cafe mailing list
> Haskell-Cafe at haskell.org
> http://www.haskell.org/mailman/listinfo/haskell-cafe
-------------- next part --------------
module MyDet where

type Matrix = [[Rational]]

first :: [[a]] -> [a]
first [] = []
first (xs:_) = xs

pad a = [x ++ x | x <- a]

sig n = case odd n of
           True -> -1
	   False -> 1

ind :: Int -> Int -> Matrix -> Rational
ind i k mat = (mat !! (i-1)) !! (k-1)

slice :: (Int,Int) -> (Int,Int) -> Matrix -> Matrix
slice (i,j) (k,l) mat
    = let rowSliceFront = drop (i-1) mat
          rowSlice      = take (k-i+1) rowSliceFront
	  columnSliceFront = map (drop (j-1)) rowSlice
      in map (take (l-j+1)) columnSliceFront

determinant :: Matrix -> Rational
determinant mat = det (pad mat) (length mat)

det :: Matrix -> Int -> Rational
det mat 1 = ind 1 1 mat
det mat size
    = let coeffs = zipWith (*) (cycle [1,sig newsize]) $ first mat
          newsize = size - 1
	  leastToMax = newsize - 1
	  slices = [slice (2,i) (size,i+leastToMax) mat | i <- [2 .. (size+1)]]
      in sum (zipWith (*) coeffs (map (\l -> det l newsize) (map pad slices)))