Thanks Ric for uncovering this problem by testing to the limit.
I tried the Haskell function reducedRowEchelonForm with matrix_3 and it
produced the same answer as gauss_jordan x: m.
The problem is that I didn't implement the final pattern of
'rowEchelonForm'. If I omit this pattern in Haskell then it breaks too.
I'll have to puzzle on that one.
On 22-04-12 01:46, Ric Sherlock wrote:
> 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
--
Met vriendelijke groet,
@@i = Arie Groeneveld
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
