Nice.
Note though that this implementation doesn't seem to be as robust for
trickier matrices.
matrix_2
2 0 _1 0 0
1 0 0 _1 0
3 0 0 _2 _1
0 1 0 0 _2
0 1 _1 0 0
gauss_jordan matrix_2
1 0 0 0 _1
0 1 0 0 _2
0 0 1 0 _2
0 0 0 1 _1
0 0 0 0 0
RREF matrix_2
1 0 0 0 0
0 1 0 0 0
0 0 1 0 0
0 0 0 1 0
0 0 0 0 1
matrix_3
1 2 3 4 3 1
2 4 6 2 6 2
3 6 18 9 9 _6
4 8 12 10 12 4
5 10 24 11 15 _4
gauss_jordan matrix_3
1 2 0 0 3 0
0 0 1 0 0 0
0 0 0 1 0 0
0 0 0 0 0 1
0 0 0 0 0 0
gauss_jordan x: matrix_3
1 2 0 0 3 4
0 0 1 0 0 _1
0 0 0 1 0 0
0 0 0 0 0 0
0 0 0 0 0 0
RREF matrix_3
|stack error: REF
| 0,. REF(}.ks),}."1}.y
RREF x: matrix_3
|stack error: REF
| 0,. REF(}.ks),}."1}.y
On Sun, Apr 22, 2012 at 5:43 AM, Aai <agroeneveld...@gmail.com> wrote:
> In search of a functional version of subject not using indices I found a
> Haskell program implemented by David Amos (for those interested see [1];
> I also added my slightly modified version).
>
> The idea is to improve my attempt of this approach.
>
> BTW: the J contribution at rosetta code (see [2]) uses the gauss_jordan
> methode of 'math/misc/linear'
>
> An almost one to one translation of the Haskell program results in:
>
> REF=: 3 :0
> if. 0=#y do. i.0 0 return. end.
> if. 0~:k=.{.ks=.{.y do.
> r1=. k%~}.ks
> (1,r1),0,. REF r1 (}.@]-(*{.))"1 }.y
> else. if. 0=#z=. (#~0~:{."1) }.y do.
> 0,. REF (}. ks),}."1 }.y
> else. REF (}.y),~ks+{.z
> end.
> end.
> )
>
> reduceStep=: 3 :0
> select. {.ks=.{.y
> case. 1 do. ks, 0,. (}.-(}.ks)*{.)"1 }.y
> case. 0 do. ({."1 y),"_1 reduceStep }."1 y
> case. do. y
> end.
> )
>
> RREF =: ([:{.&> (reduceStep@}.&.>^:(1<#&>)^:a:@ <@reduceStep))&.|. @ REF
>
>
> Tests:
> =======
>
> rcTest=: 1 2 _1 _4,2 3 _1 _11,:_2 0 _3 22
> test =: 0 3 _6 6 4 _5, 3 _7 8 _5 8 9,: 3 _9 12 _9 6 15
>
> RREF rcTest
> 1 0 0 _8
> 0 1 0 1
> 0 0 1 _2
>
> RREF test
> 1 0 _2 3 0 _24
> 0 1 _2 2 0 _7
> 0 0 0 0 1 4
>
>
>
>
>
>
> [1] http://tinyurl.com/7tfthg9
>
> -- D. Amos
> -- HaskelForMath Linear.algebra
>
> import Data.List
> import Control.Monad
> import Control.Arrow
>
> rowEchelonForm [] = []
> rowEchelonForm (hs@(x:xs):rs) =
> if x /= 0
> then let r' = map (/x) xs
> in (1:r') : map (0:) (rowEchelonForm [zipWith (-) ys (map (y*) r')
> | (y:ys) <- rs])
> else case filter ((/= 0).head) rs of
> [] -> map (0:) (rowEchelonForm $ xs : map tail rs)
> r:_ -> rowEchelonForm ((zipWith (+) hs r) : rs)
> --rowEchelonForm zs@([]:_) = zs
>
> reducedRowEchelonForm :: (Eq a, Fractional a) => [[a]] -> [[a]]
> reducedRowEchelonForm m = reverse $ reduce $ reverse $ rowEchelonForm m
> where
> reduce = unfoldr (\m -> if null m then Nothing else Just ((head &&&
> tail) $ reduceStep m) )
> reduceStep ((1:xs):rs) = (1:xs) : [ 0: zipWith (-) ys (map (y*) xs) |
> y:ys <- rs]
> reduceStep rs@((0:_):_) = zipWith (:) (map head rs) (reduceStep $ map
> tail rs)
> reduceStep rs = rs
>
>
> [2] http://rosettacode.org/wiki/Reduced_row_echelon_form#J
>
>
>
> --
> Met vriendelijke groet,
> @@i = Arie Groeneveld
>
> ----------------------------------------------------------------------
> For information about J forums see http://www.jsoftware.com/forums.htm
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
