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
