This is an alternative site for discovering Elm packages. You may be looking for the official Elm package site instead.

# Gauss

The Gauss and Gauss-Jordan Elimination Algorithms, on type `Matrix` from the `eeue56/elm-flat-matrix` package, over `Rational`s from the `imeckler/ratio` package.

gauss : Int -> Matrix Rational -> Matrix Rational

`gauss 3` on matrix

``````[ [ 3, 4,  0, 2 ]
, [ 1, 0, -3, 1 ]
, [ 0, 2, -1, 5 ]
]
``````

gives matrix

``````[ [ 1, 4/3,   0,  2/3 ]
, [ 0,   1, 9/4, -1/4 ]
, [ 0,   0,   1,   -1 ]
]
``````
gaussJordan : Int -> Matrix Rational -> Matrix Rational

`gaussJordan 3` on matrix

``````[ [ 3, 4,  0, 2 ]
, [ 1, 0, -3, 1 ]
, [ 0, 2, -1, 5 ]
]
``````

gives matrix

``````[ [ 1, 0, 0, -2 ]
, [ 0, 1, 0,  2 ]
, [ 0, 0, 1, -1 ]
]
``````
``````module Gauss exposing (gauss, gaussJordan)

{-| The Gauss and Gauss-Jordan Elimination Algorithms, on type `Matrix`
from the [`eeue56/elm-flat-matrix` package](http://package.elm-lang.org/packages/eeue56/elm-flat-matrix/latest),
over `Rational`s from the [`imeckler/ratio` package](http://package.elm-lang.org/packages/imeckler/ratio/latest).

@docs gauss, gaussJordan
-}

import Matrix exposing (Matrix)
import Ratio exposing (Rational)
import Array
import List.Extra as List

{-| `gauss 3` on matrix

[ [ 3, 4,  0, 2 ]
, [ 1, 0, -3, 1 ]
, [ 0, 2, -1, 5 ]
]

gives matrix

[ [ 1, 4/3,   0,  2/3 ]
, [ 0,   1, 9/4, -1/4 ]
, [ 0,   0,   1,   -1 ]
]
-}
gauss : Int -> Matrix Rational -> Matrix Rational
gauss n =
let
swapIfNecessary i j m =
if i == j then
m
else
List.foldl (\c -> Matrix.set c j (safeGet c i m) >> Matrix.set c i (safeGet c j m)) m (List.range i n)

forward i m =
case Matrix.getColumn i m of
Nothing ->
Debug.crash "IMPOSSIBLE!"

Just column ->
case List.find (\( _, v ) -> Ratio.numerator v /= 0) (List.drop i (Array.toIndexedList column)) of
Nothing ->
m

Just ( j, v ) ->
let
m_ =
swapIfNecessary i j m
|> Matrix.set i i one
|> flip (List.foldl (\c -> Matrix.update c i (flip Ratio.divide v))) (List.range (i + 1) n)
in
List.foldl
(\r ->
let
f =
Ratio.negate (safeGet i r m_)
in
Matrix.set i r zero
>> flip (List.foldl (\c -> Matrix.update c r (Ratio.add (Ratio.multiply f (safeGet c i m_))))) (List.range (i + 1) n)
)
m_
(List.range (i + 1) (n - 1))
in
flip (List.foldl forward) (List.range 0 (n - 1))

{-| `gaussJordan 3` on matrix

[ [ 3, 4,  0, 2 ]
, [ 1, 0, -3, 1 ]
, [ 0, 2, -1, 5 ]
]

gives matrix

[ [ 1, 0, 0, -2 ]
, [ 0, 1, 0,  2 ]
, [ 0, 0, 1, -1 ]
]
-}
gaussJordan : Int -> Matrix Rational -> Matrix Rational
gaussJordan n =
let
backward i m =
if Ratio.numerator (safeGet i i m) == 0 then
m
else
List.foldl
(\r ->
let
k =
Ratio.negate (safeGet i r m)
in
Matrix.set i r zero
>> flip (List.foldl (\c -> Matrix.update c r (Ratio.add (Ratio.multiply k (safeGet c i m))))) (List.range (i + 1) n)
)
m
(List.range 0 (i - 1))
in
gauss n >> flip (List.foldl backward) (List.reverse (List.range 1 (n - 1)))

safeGet i j m =
case Matrix.get i j m of
Nothing ->
Debug.crash "IMPOSSIBLE!"

Just x ->
x

zero =
Ratio.fromInt 0

one =
Ratio.fromInt 1
```
```