{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE CPP #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE Trustworthy #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TypeOperators #-}

---------------------------------------------------------------------------

-- |

-- Copyright   :  (C) 2012-2015 Edward Kmett

-- License     :  BSD-style (see the file LICENSE)

--

-- Maintainer  :  Edward Kmett <ekmett@gmail.com>

-- Stability   :  experimental

-- Portability :  non-portable

--

-- Simple matrix operation for low-dimensional primitives.

---------------------------------------------------------------------------

module Linear.Matrix
  ( (!*!), (!+!), (!-!), (!*), (*!), (!!*), (*!!), (!!/)
  , column
  , adjoint
  , M22, M23, M24, M32, M33, M34, M42, M43, M44
  , m33_to_m44, m43_to_m44
  , det22, det33, det44, inv22, inv33, inv44
  , identity
  , Trace(..)
  , translation
  , transpose
  , fromQuaternion
  , mkTransformation
  , mkTransformationMat
  , _m22, _m23, _m24
  , _m32, _m33, _m34
  , _m42, _m43, _m44
  , lu
  , luFinite
  , forwardSub
  , forwardSubFinite
  , backwardSub
  , backwardSubFinite
  , luSolve
  , luSolveFinite
  , luInv
  , luInvFinite
  , luDet
  , luDetFinite
  ) where

import Control.Lens hiding (index)
import Control.Lens.Internal.Context
import Data.Distributive
import Data.Foldable as Foldable
import Data.Functor.Rep
import GHC.TypeLits
import Linear.Quaternion
import Linear.V
import Linear.V2
import Linear.V3
import Linear.V4
import Linear.Vector
import Linear.Conjugate
import Linear.Trace

-- $setup

-- >>> import Control.Lens hiding (index)

-- >>> import Data.Complex (Complex (..))

-- >>> import Linear.V2

-- >>> import Linear.V3

-- >>> import Linear.V

-- >>> import qualified Data.IntMap as IntMap

-- >>> import Debug.SimpleReflect.Vars


-- | This is a generalization of 'Control.Lens.inside' to work over any corepresentable 'Functor'.

--

-- @

-- 'column' :: 'Representable' f => 'Lens' s t a b -> 'Lens' (f s) (f t) (f a) (f b)

-- @

--

-- In practice it is used to access a column of a matrix.

--

-- >>> V2 (V3 1 2 3) (V3 4 5 6) ^._x

-- V3 1 2 3

--

-- >>> V2 (V3 1 2 3) (V3 4 5 6) ^.column _x

-- V2 1 4

column :: Representable f => LensLike (Context a b) s t a b -> Lens (f s) (f t) (f a) (f b)
column :: forall (f :: * -> *) a b s t.
Representable f =>
LensLike (Context a b) s t a b -> Lens (f s) (f t) (f a) (f b)
column LensLike (Context a b) s t a b
l f a -> f (f b)
f f s
es = f b -> f t
o forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> f a -> f (f b)
f f a
i where
   go :: s -> Context a b t
go = LensLike (Context a b) s t a b
l (forall a b t. (b -> t) -> a -> Context a b t
Context forall a. a -> a
id)
   i :: f a
i = forall (f :: * -> *) a. Representable f => (Rep f -> a) -> f a
tabulate forall a b. (a -> b) -> a -> b
$ \ Rep f
e -> forall (w :: * -> * -> * -> *) a c t.
IndexedComonadStore w =>
w a c t -> a
ipos forall a b. (a -> b) -> a -> b
$ s -> Context a b t
go (forall (f :: * -> *) a. Representable f => f a -> Rep f -> a
index f s
es Rep f
e)
   o :: f b -> f t
o f b
eb = forall (f :: * -> *) a. Representable f => (Rep f -> a) -> f a
tabulate forall a b. (a -> b) -> a -> b
$ \ Rep f
e -> forall (w :: * -> * -> * -> *) c a t.
IndexedComonadStore w =>
c -> w a c t -> t
ipeek (forall (f :: * -> *) a. Representable f => f a -> Rep f -> a
index f b
eb Rep f
e) (s -> Context a b t
go (forall (f :: * -> *) a. Representable f => f a -> Rep f -> a
index f s
es Rep f
e))

infixl 7 !*!
-- | Matrix product. This can compute any combination of sparse and dense multiplication.

--

-- >>> V2 (V3 1 2 3) (V3 4 5 6) !*! V3 (V2 1 2) (V2 3 4) (V2 4 5)

-- V2 (V2 19 25) (V2 43 58)

--

-- >>> V2 (IntMap.fromList [(1,2)]) (IntMap.fromList [(2,3)]) !*! IntMap.fromList [(1,V3 0 0 1), (2, V3 0 0 5)]

-- V2 (V3 0 0 2) (V3 0 0 15)

(!*!) :: (Functor m, Foldable t, Additive t, Additive n, Num a) => m (t a) -> t (n a) -> m (n a)
m (t a)
f !*! :: forall (m :: * -> *) (t :: * -> *) (n :: * -> *) a.
(Functor m, Foldable t, Additive t, Additive n, Num a) =>
m (t a) -> t (n a) -> m (n a)
!*! t (n a)
g = forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (\ t a
f' -> forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
Foldable.foldl' forall (f :: * -> *) a. (Additive f, Num a) => f a -> f a -> f a
(^+^) forall (f :: * -> *) a. (Additive f, Num a) => f a
zero forall a b. (a -> b) -> a -> b
$ forall (f :: * -> *) a b c.
Additive f =>
(a -> b -> c) -> f a -> f b -> f c
liftI2 forall (f :: * -> *) a. (Functor f, Num a) => a -> f a -> f a
(*^) t a
f' t (n a)
g) m (t a)
f

infixl 6 !+!
-- | Entry-wise matrix addition.

--

-- >>> V2 (V3 1 2 3) (V3 4 5 6) !+! V2 (V3 7 8 9) (V3 1 2 3)

-- V2 (V3 8 10 12) (V3 5 7 9)

(!+!) :: (Additive m, Additive n, Num a) => m (n a) -> m (n a) -> m (n a)
m (n a)
as !+! :: forall (m :: * -> *) (n :: * -> *) a.
(Additive m, Additive n, Num a) =>
m (n a) -> m (n a) -> m (n a)
!+! m (n a)
bs = forall (f :: * -> *) a.
Additive f =>
(a -> a -> a) -> f a -> f a -> f a
liftU2 forall (f :: * -> *) a. (Additive f, Num a) => f a -> f a -> f a
(^+^) m (n a)
as m (n a)
bs

infixl 6 !-!
-- | Entry-wise matrix subtraction.

--

-- >>> V2 (V3 1 2 3) (V3 4 5 6) !-! V2 (V3 7 8 9) (V3 1 2 3)

-- V2 (V3 (-6) (-6) (-6)) (V3 3 3 3)

(!-!) :: (Additive m, Additive n, Num a) => m (n a) -> m (n a) -> m (n a)
m (n a)
as !-! :: forall (m :: * -> *) (n :: * -> *) a.
(Additive m, Additive n, Num a) =>
m (n a) -> m (n a) -> m (n a)
!-! m (n a)
bs = forall (f :: * -> *) a.
Additive f =>
(a -> a -> a) -> f a -> f a -> f a
liftU2 forall (f :: * -> *) a. (Additive f, Num a) => f a -> f a -> f a
(^-^) m (n a)
as m (n a)
bs

infixl 7 !*
-- | Matrix * column vector

--

-- >>> V2 (V3 1 2 3) (V3 4 5 6) !* V3 7 8 9

-- V2 50 122

(!*) :: (Functor m, Foldable r, Additive r, Num a) => m (r a) -> r a -> m a
m (r a)
m !* :: forall (m :: * -> *) (r :: * -> *) a.
(Functor m, Foldable r, Additive r, Num a) =>
m (r a) -> r a -> m a
!* r a
v = forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (\r a
r -> forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
Foldable.sum forall a b. (a -> b) -> a -> b
$ forall (f :: * -> *) a b c.
Additive f =>
(a -> b -> c) -> f a -> f b -> f c
liftI2 forall a. Num a => a -> a -> a
(*) r a
r r a
v) m (r a)
m

infixl 7 *!
-- | Row vector * matrix

--

-- >>> V2 1 2 *! V2 (V3 3 4 5) (V3 6 7 8)

-- V3 15 18 21


-- (*!) :: (Metric r, Additive n, Num a) => r a -> r (n a) -> n a

-- f *! g = dot f <$> distribute g


(*!) :: (Num a, Foldable t, Additive f, Additive t) => t a -> t (f a) -> f a
t a
f *! :: forall a (t :: * -> *) (f :: * -> *).
(Num a, Foldable t, Additive f, Additive t) =>
t a -> t (f a) -> f a
*! t (f a)
g = forall (f :: * -> *) (v :: * -> *) a.
(Foldable f, Additive v, Num a) =>
f (v a) -> v a
sumV forall a b. (a -> b) -> a -> b
$ forall (f :: * -> *) a b c.
Additive f =>
(a -> b -> c) -> f a -> f b -> f c
liftI2 forall (f :: * -> *) a. (Functor f, Num a) => a -> f a -> f a
(*^) t a
f t (f a)
g

infixl 7 *!!
-- | Scalar-matrix product

--

-- >>> 5 *!! V2 (V2 1 2) (V2 3 4)

-- V2 (V2 5 10) (V2 15 20)

(*!!) :: (Functor m, Functor r, Num a) => a -> m (r a) -> m (r a)
a
s *!! :: forall (m :: * -> *) (r :: * -> *) a.
(Functor m, Functor r, Num a) =>
a -> m (r a) -> m (r a)
*!! m (r a)
m = forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (a
s forall (f :: * -> *) a. (Functor f, Num a) => a -> f a -> f a
*^) m (r a)
m
{-# INLINE (*!!) #-}

infixl 7 !!*
-- | Matrix-scalar product

--

-- >>> V2 (V2 1 2) (V2 3 4) !!* 5

-- V2 (V2 5 10) (V2 15 20)

(!!*) :: (Functor m, Functor r, Num a) => m (r a) -> a -> m (r a)
!!* :: forall (m :: * -> *) (r :: * -> *) a.
(Functor m, Functor r, Num a) =>
m (r a) -> a -> m (r a)
(!!*) = forall a b c. (a -> b -> c) -> b -> a -> c
flip forall (m :: * -> *) (r :: * -> *) a.
(Functor m, Functor r, Num a) =>
a -> m (r a) -> m (r a)
(*!!)
{-# INLINE (!!*) #-}

infixl 7 !!/
-- | Matrix-scalar division

(!!/) :: (Functor m, Functor r, Fractional a) => m (r a) -> a -> m (r a)
m (r a)
m !!/ :: forall (m :: * -> *) (r :: * -> *) a.
(Functor m, Functor r, Fractional a) =>
m (r a) -> a -> m (r a)
!!/ a
s = forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (forall (f :: * -> *) a.
(Functor f, Fractional a) =>
f a -> a -> f a
^/ a
s) m (r a)
m
{-# INLINE (!!/) #-}

-- | Hermitian conjugate or conjugate transpose

--

-- >>> adjoint (V2 (V2 (1 :+ 2) (3 :+ 4)) (V2 (5 :+ 6) (7 :+ 8)))

-- V2 (V2 (1.0 :+ (-2.0)) (5.0 :+ (-6.0))) (V2 (3.0 :+ (-4.0)) (7.0 :+ (-8.0)))

adjoint :: (Functor m, Distributive n, Conjugate a) => m (n a) -> n (m a)
adjoint :: forall (m :: * -> *) (n :: * -> *) a.
(Functor m, Distributive n, Conjugate a) =>
m (n a) -> n (m a)
adjoint = forall (g :: * -> *) (f :: * -> *) a b.
(Distributive g, Functor f) =>
(a -> g b) -> f a -> g (f b)
collect (forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap forall a. Conjugate a => a -> a
conjugate)
{-# INLINE adjoint #-}

-- * Matrices

--

-- Matrices use a row-major representation.


-- | A 2x2 matrix with row-major representation

type M22 a = V2 (V2 a)
-- | A 2x3 matrix with row-major representation

type M23 a = V2 (V3 a)
-- | A 2x4 matrix with row-major representation

type M24 a = V2 (V4 a)
-- | A 3x2 matrix with row-major representation

type M32 a = V3 (V2 a)
-- | A 3x3 matrix with row-major representation

type M33 a = V3 (V3 a)
-- | A 3x4 matrix with row-major representation

type M34 a = V3 (V4 a)
-- | A 4x2 matrix with row-major representation

type M42 a = V4 (V2 a)
-- | A 4x3 matrix with row-major representation

type M43 a = V4 (V3 a)
-- | A 4x4 matrix with row-major representation

type M44 a = V4 (V4 a)

-- | Build a rotation matrix from a unit 'Quaternion'.

fromQuaternion :: Num a => Quaternion a -> M33 a
fromQuaternion :: forall a. Num a => Quaternion a -> M33 a
fromQuaternion (Quaternion a
w (V3 a
x a
y a
z)) =
  forall a. a -> a -> a -> V3 a
V3 (forall a. a -> a -> a -> V3 a
V3 (a
1forall a. Num a => a -> a -> a
-a
2forall a. Num a => a -> a -> a
*(a
y2forall a. Num a => a -> a -> a
+a
z2)) (a
2forall a. Num a => a -> a -> a
*(a
xyforall a. Num a => a -> a -> a
-a
zw)) (a
2forall a. Num a => a -> a -> a
*(a
xzforall a. Num a => a -> a -> a
+a
yw)))
     (forall a. a -> a -> a -> V3 a
V3 (a
2forall a. Num a => a -> a -> a
*(a
xyforall a. Num a => a -> a -> a
+a
zw)) (a
1forall a. Num a => a -> a -> a
-a
2forall a. Num a => a -> a -> a
*(a
x2forall a. Num a => a -> a -> a
+a
z2)) (a
2forall a. Num a => a -> a -> a
*(a
yzforall a. Num a => a -> a -> a
-a
xw)))
     (forall a. a -> a -> a -> V3 a
V3 (a
2forall a. Num a => a -> a -> a
*(a
xzforall a. Num a => a -> a -> a
-a
yw)) (a
2forall a. Num a => a -> a -> a
*(a
yzforall a. Num a => a -> a -> a
+a
xw)) (a
1forall a. Num a => a -> a -> a
-a
2forall a. Num a => a -> a -> a
*(a
x2forall a. Num a => a -> a -> a
+a
y2)))
  where x2 :: a
x2 = a
xforall a. Num a => a -> a -> a
*a
x
        y2 :: a
y2 = a
yforall a. Num a => a -> a -> a
*a
y
        z2 :: a
z2 = a
zforall a. Num a => a -> a -> a
*a
z
        xy :: a
xy = a
xforall a. Num a => a -> a -> a
*a
y
        xz :: a
xz = a
xforall a. Num a => a -> a -> a
*a
z
        xw :: a
xw = a
xforall a. Num a => a -> a -> a
*a
w
        yz :: a
yz = a
yforall a. Num a => a -> a -> a
*a
z
        yw :: a
yw = a
yforall a. Num a => a -> a -> a
*a
w
        zw :: a
zw = a
zforall a. Num a => a -> a -> a
*a
w
{-# INLINE fromQuaternion #-}

-- | Build a transformation matrix from a rotation matrix and a

-- translation vector.

mkTransformationMat :: Num a => M33 a -> V3 a -> M44 a
mkTransformationMat :: forall a. Num a => M33 a -> V3 a -> M44 a
mkTransformationMat (V3 V3 a
r1 V3 a
r2 V3 a
r3) (V3 a
tx a
ty a
tz) =
  forall a. a -> a -> a -> a -> V4 a
V4 (forall {a}. V3 a -> a -> V4 a
snoc3 V3 a
r1 a
tx) (forall {a}. V3 a -> a -> V4 a
snoc3 V3 a
r2 a
ty) (forall {a}. V3 a -> a -> V4 a
snoc3 V3 a
r3 a
tz) (forall a. a -> a -> a -> a -> V4 a
V4 a
0 a
0 a
0 a
1)
  where snoc3 :: V3 a -> a -> V4 a
snoc3 (V3 a
x a
y a
z) = forall a. a -> a -> a -> a -> V4 a
V4 a
x a
y a
z
{-# INLINE mkTransformationMat #-}

-- |Build a transformation matrix from a rotation expressed as a

-- 'Quaternion' and a translation vector.

mkTransformation :: Num a => Quaternion a -> V3 a -> M44 a
mkTransformation :: forall a. Num a => Quaternion a -> V3 a -> M44 a
mkTransformation = forall a. Num a => M33 a -> V3 a -> M44 a
mkTransformationMat forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. Num a => Quaternion a -> M33 a
fromQuaternion
{-# INLINE mkTransformation #-}

-- | Convert from a 4x3 matrix to a 4x4 matrix, extending it with the @[ 0 0 0 1 ]@ column vector

m43_to_m44 :: Num a => M43 a -> M44 a
m43_to_m44 :: forall a. Num a => M43 a -> M44 a
m43_to_m44
  (V4 (V3 a
a a
b a
c)
      (V3 a
d a
e a
f)
      (V3 a
g a
h a
i)
      (V3 a
j a
k a
l)) =
  forall a. a -> a -> a -> a -> V4 a
V4 (forall a. a -> a -> a -> a -> V4 a
V4 a
a a
b a
c a
0)
     (forall a. a -> a -> a -> a -> V4 a
V4 a
d a
e a
f a
0)
     (forall a. a -> a -> a -> a -> V4 a
V4 a
g a
h a
i a
0)
     (forall a. a -> a -> a -> a -> V4 a
V4 a
j a
k a
l a
1)

-- | Convert a 3x3 matrix to a 4x4 matrix extending it with 0's in the new row and column.

m33_to_m44 :: Num a => M33 a -> M44 a
m33_to_m44 :: forall a. Num a => M33 a -> M44 a
m33_to_m44 (V3 V3 a
r1 V3 a
r2 V3 a
r3) = forall a. a -> a -> a -> a -> V4 a
V4 (forall a. Num a => V3 a -> V4 a
vector V3 a
r1) (forall a. Num a => V3 a -> V4 a
vector V3 a
r2) (forall a. Num a => V3 a -> V4 a
vector V3 a
r3) (forall a. Num a => V3 a -> V4 a
point V3 a
0)

-- |The identity matrix for any dimension vector.

--

-- >>> identity :: M44 Int

-- V4 (V4 1 0 0 0) (V4 0 1 0 0) (V4 0 0 1 0) (V4 0 0 0 1)

-- >>> identity :: V3 (V3 Int)

-- V3 (V3 1 0 0) (V3 0 1 0) (V3 0 0 1)

identity :: (Num a, Traversable t, Applicative t) => t (t a)
identity :: forall a (t :: * -> *).
(Num a, Traversable t, Applicative t) =>
t (t a)
identity = forall (t :: * -> *) a. (Traversable t, Num a) => t a -> t (t a)
scaled (forall (f :: * -> *) a. Applicative f => a -> f a
pure a
1)

-- |Extract the translation vector (first three entries of the last

-- column) from a 3x4 or 4x4 matrix.

translation :: (Representable t, R3 t, R4 v) => Lens' (t (v a)) (V3 a)
translation :: forall (t :: * -> *) (v :: * -> *) a.
(Representable t, R3 t, R4 v) =>
Lens' (t (v a)) (V3 a)
translation = forall (f :: * -> *) a b s t.
Representable f =>
LensLike (Context a b) s t a b -> Lens (f s) (f t) (f a) (f b)
column forall (t :: * -> *) a. R4 t => Lens' (t a) a
_wforall b c a. (b -> c) -> (a -> b) -> a -> c
.forall (t :: * -> *) a. R3 t => Lens' (t a) (V3 a)
_xyz
{-
translation f rs = aux <$> f (view _w <$> view _xyz rs)
 where aux (V3 x y z) = (_x._w .~ x) . (_y._w .~ y) . (_z._w .~ z) $ rs

-- translation :: (R3 t, R4 v, Functor f, Functor t) => (V3 a -> f (V3 a)) -> t (v a) -> f (t a)
-- translation = (. fmap (^._w)) . _xyz where
--   x ^. l = getConst (l Const x)
-}

-- |Extract a 2x2 matrix from a matrix of higher dimensions by dropping excess

-- rows and columns.

_m22 :: (Representable t, R2 t, R2 v) => Lens' (t (v a)) (M22 a)
_m22 :: forall (t :: * -> *) (v :: * -> *) a.
(Representable t, R2 t, R2 v) =>
Lens' (t (v a)) (M22 a)
_m22 = forall (f :: * -> *) a b s t.
Representable f =>
LensLike (Context a b) s t a b -> Lens (f s) (f t) (f a) (f b)
column forall (t :: * -> *) a. R2 t => Lens' (t a) (V2 a)
_xyforall b c a. (b -> c) -> (a -> b) -> a -> c
.forall (t :: * -> *) a. R2 t => Lens' (t a) (V2 a)
_xy

-- |Extract a 2x3 matrix from a matrix of higher dimensions by dropping excess

-- rows and columns.

_m23 :: (Representable t, R2 t, R3 v) => Lens' (t (v a)) (M23 a)
_m23 :: forall (t :: * -> *) (v :: * -> *) a.
(Representable t, R2 t, R3 v) =>
Lens' (t (v a)) (M23 a)
_m23 = forall (f :: * -> *) a b s t.
Representable f =>
LensLike (Context a b) s t a b -> Lens (f s) (f t) (f a) (f b)
column forall (t :: * -> *) a. R3 t => Lens' (t a) (V3 a)
_xyzforall b c a. (b -> c) -> (a -> b) -> a -> c
.forall (t :: * -> *) a. R2 t => Lens' (t a) (V2 a)
_xy

-- |Extract a 2x4 matrix from a matrix of higher dimensions by dropping excess

-- rows and columns.

_m24 :: (Representable t, R2 t, R4 v) => Lens' (t (v a)) (M24 a)
_m24 :: forall (t :: * -> *) (v :: * -> *) a.
(Representable t, R2 t, R4 v) =>
Lens' (t (v a)) (M24 a)
_m24 = forall (f :: * -> *) a b s t.
Representable f =>
LensLike (Context a b) s t a b -> Lens (f s) (f t) (f a) (f b)
column forall (t :: * -> *) a. R4 t => Lens' (t a) (V4 a)
_xyzwforall b c a. (b -> c) -> (a -> b) -> a -> c
.forall (t :: * -> *) a. R2 t => Lens' (t a) (V2 a)
_xy

-- |Extract a 3x2 matrix from a matrix of higher dimensions by dropping excess

-- rows and columns.

_m32 :: (Representable t, R3 t, R2 v) => Lens' (t (v a)) (M32 a)
_m32 :: forall (t :: * -> *) (v :: * -> *) a.
(Representable t, R3 t, R2 v) =>
Lens' (t (v a)) (M32 a)
_m32 = forall (f :: * -> *) a b s t.
Representable f =>
LensLike (Context a b) s t a b -> Lens (f s) (f t) (f a) (f b)
column forall (t :: * -> *) a. R2 t => Lens' (t a) (V2 a)
_xyforall b c a. (b -> c) -> (a -> b) -> a -> c
.forall (t :: * -> *) a. R3 t => Lens' (t a) (V3 a)
_xyz

-- |Extract a 3x3 matrix from a matrix of higher dimensions by dropping excess

-- rows and columns.

_m33 :: (Representable t, R3 t, R3 v) => Lens' (t (v a)) (M33 a)
_m33 :: forall (t :: * -> *) (v :: * -> *) a.
(Representable t, R3 t, R3 v) =>
Lens' (t (v a)) (M33 a)
_m33 = forall (f :: * -> *) a b s t.
Representable f =>
LensLike (Context a b) s t a b -> Lens (f s) (f t) (f a) (f b)
column forall (t :: * -> *) a. R3 t => Lens' (t a) (V3 a)
_xyzforall b c a. (b -> c) -> (a -> b) -> a -> c
.forall (t :: * -> *) a. R3 t => Lens' (t a) (V3 a)
_xyz

-- |Extract a 3x4 matrix from a matrix of higher dimensions by dropping excess

-- rows and columns.

_m34 :: (Representable t, R3 t, R4 v) => Lens' (t (v a)) (M34 a)
_m34 :: forall (t :: * -> *) (v :: * -> *) a.
(Representable t, R3 t, R4 v) =>
Lens' (t (v a)) (M34 a)
_m34 = forall (f :: * -> *) a b s t.
Representable f =>
LensLike (Context a b) s t a b -> Lens (f s) (f t) (f a) (f b)
column forall (t :: * -> *) a. R4 t => Lens' (t a) (V4 a)
_xyzwforall b c a. (b -> c) -> (a -> b) -> a -> c
.forall (t :: * -> *) a. R3 t => Lens' (t a) (V3 a)
_xyz

-- |Extract a 4x2 matrix from a matrix of higher dimensions by dropping excess

-- rows and columns.

_m42 :: (Representable t, R4 t, R2 v) => Lens' (t (v a)) (M42 a)
_m42 :: forall (t :: * -> *) (v :: * -> *) a.
(Representable t, R4 t, R2 v) =>
Lens' (t (v a)) (M42 a)
_m42 = forall (f :: * -> *) a b s t.
Representable f =>
LensLike (Context a b) s t a b -> Lens (f s) (f t) (f a) (f b)
column forall (t :: * -> *) a. R2 t => Lens' (t a) (V2 a)
_xyforall b c a. (b -> c) -> (a -> b) -> a -> c
.forall (t :: * -> *) a. R4 t => Lens' (t a) (V4 a)
_xyzw

-- |Extract a 4x3 matrix from a matrix of higher dimensions by dropping excess

-- rows and columns.

_m43 :: (Representable t, R4 t, R3 v) => Lens' (t (v a)) (M43 a)
_m43 :: forall (t :: * -> *) (v :: * -> *) a.
(Representable t, R4 t, R3 v) =>
Lens' (t (v a)) (M43 a)
_m43 = forall (f :: * -> *) a b s t.
Representable f =>
LensLike (Context a b) s t a b -> Lens (f s) (f t) (f a) (f b)
column forall (t :: * -> *) a. R3 t => Lens' (t a) (V3 a)
_xyzforall b c a. (b -> c) -> (a -> b) -> a -> c
.forall (t :: * -> *) a. R4 t => Lens' (t a) (V4 a)
_xyzw

-- |Extract a 4x4 matrix from a matrix of higher dimensions by dropping excess

-- rows and columns.

_m44 :: (Representable t, R4 t, R4 v) => Lens' (t (v a)) (M44 a)
_m44 :: forall (t :: * -> *) (v :: * -> *) a.
(Representable t, R4 t, R4 v) =>
Lens' (t (v a)) (M44 a)
_m44 = forall (f :: * -> *) a b s t.
Representable f =>
LensLike (Context a b) s t a b -> Lens (f s) (f t) (f a) (f b)
column forall (t :: * -> *) a. R4 t => Lens' (t a) (V4 a)
_xyzwforall b c a. (b -> c) -> (a -> b) -> a -> c
.forall (t :: * -> *) a. R4 t => Lens' (t a) (V4 a)
_xyzw

-- |2x2 matrix determinant.

--

-- >>> det22 (V2 (V2 a b) (V2 c d))

-- a * d - b * c

det22 :: Num a => M22 a -> a
det22 :: forall a. Num a => M22 a -> a
det22 (V2 (V2 a
a a
b) (V2 a
c a
d)) = a
a forall a. Num a => a -> a -> a
* a
d forall a. Num a => a -> a -> a
- a
b forall a. Num a => a -> a -> a
* a
c
{-# INLINE det22 #-}

-- |3x3 matrix determinant.

--

-- >>> det33 (V3 (V3 a b c) (V3 d e f) (V3 g h i))

-- a * (e * i - f * h) - d * (b * i - c * h) + g * (b * f - c * e)

det33 :: Num a => M33 a -> a
det33 :: forall a. Num a => M33 a -> a
det33 (V3 (V3 a
a a
b a
c)
          (V3 a
d a
e a
f)
          (V3 a
g a
h a
i)) = a
a forall a. Num a => a -> a -> a
* (a
eforall a. Num a => a -> a -> a
*a
iforall a. Num a => a -> a -> a
-a
fforall a. Num a => a -> a -> a
*a
h) forall a. Num a => a -> a -> a
- a
d forall a. Num a => a -> a -> a
* (a
bforall a. Num a => a -> a -> a
*a
iforall a. Num a => a -> a -> a
-a
cforall a. Num a => a -> a -> a
*a
h) forall a. Num a => a -> a -> a
+ a
g forall a. Num a => a -> a -> a
* (a
bforall a. Num a => a -> a -> a
*a
fforall a. Num a => a -> a -> a
-a
cforall a. Num a => a -> a -> a
*a
e)
{-# INLINE det33 #-}

-- |4x4 matrix determinant.

det44 :: Num a => M44 a -> a
det44 :: forall a. Num a => M44 a -> a
det44 (V4 (V4 a
i00 a
i01 a
i02 a
i03)
          (V4 a
i10 a
i11 a
i12 a
i13)
          (V4 a
i20 a
i21 a
i22 a
i23)
          (V4 a
i30 a
i31 a
i32 a
i33)) =
  let
    s0 :: a
s0 = a
i00 forall a. Num a => a -> a -> a
* a
i11 forall a. Num a => a -> a -> a
- a
i10 forall a. Num a => a -> a -> a
* a
i01
    s1 :: a
s1 = a
i00 forall a. Num a => a -> a -> a
* a
i12 forall a. Num a => a -> a -> a
- a
i10 forall a. Num a => a -> a -> a
* a
i02
    s2 :: a
s2 = a
i00 forall a. Num a => a -> a -> a
* a
i13 forall a. Num a => a -> a -> a
- a
i10 forall a. Num a => a -> a -> a
* a
i03
    s3 :: a
s3 = a
i01 forall a. Num a => a -> a -> a
* a
i12 forall a. Num a => a -> a -> a
- a
i11 forall a. Num a => a -> a -> a
* a
i02
    s4 :: a
s4 = a
i01 forall a. Num a => a -> a -> a
* a
i13 forall a. Num a => a -> a -> a
- a
i11 forall a. Num a => a -> a -> a
* a
i03
    s5 :: a
s5 = a
i02 forall a. Num a => a -> a -> a
* a
i13 forall a. Num a => a -> a -> a
- a
i12 forall a. Num a => a -> a -> a
* a
i03

    c5 :: a
c5 = a
i22 forall a. Num a => a -> a -> a
* a
i33 forall a. Num a => a -> a -> a
- a
i32 forall a. Num a => a -> a -> a
* a
i23
    c4 :: a
c4 = a
i21 forall a. Num a => a -> a -> a
* a
i33 forall a. Num a => a -> a -> a
- a
i31 forall a. Num a => a -> a -> a
* a
i23
    c3 :: a
c3 = a
i21 forall a. Num a => a -> a -> a
* a
i32 forall a. Num a => a -> a -> a
- a
i31 forall a. Num a => a -> a -> a
* a
i22
    c2 :: a
c2 = a
i20 forall a. Num a => a -> a -> a
* a
i33 forall a. Num a => a -> a -> a
- a
i30 forall a. Num a => a -> a -> a
* a
i23
    c1 :: a
c1 = a
i20 forall a. Num a => a -> a -> a
* a
i32 forall a. Num a => a -> a -> a
- a
i30 forall a. Num a => a -> a -> a
* a
i22
    c0 :: a
c0 = a
i20 forall a. Num a => a -> a -> a
* a
i31 forall a. Num a => a -> a -> a
- a
i30 forall a. Num a => a -> a -> a
* a
i21
  in a
s0 forall a. Num a => a -> a -> a
* a
c5 forall a. Num a => a -> a -> a
- a
s1 forall a. Num a => a -> a -> a
* a
c4 forall a. Num a => a -> a -> a
+ a
s2 forall a. Num a => a -> a -> a
* a
c3 forall a. Num a => a -> a -> a
+ a
s3 forall a. Num a => a -> a -> a
* a
c2 forall a. Num a => a -> a -> a
- a
s4 forall a. Num a => a -> a -> a
* a
c1 forall a. Num a => a -> a -> a
+ a
s5 forall a. Num a => a -> a -> a
* a
c0
{-# INLINE det44 #-}

-- |2x2 matrix inverse.

--

-- >>> inv22 $ V2 (V2 1 2) (V2 3 4)

-- V2 (V2 (-2.0) 1.0) (V2 1.5 (-0.5))

inv22 :: Fractional a => M22 a -> M22 a
inv22 :: forall a. Fractional a => M22 a -> M22 a
inv22 m :: M22 a
m@(V2 (V2 a
a a
b) (V2 a
c a
d)) = (a
1 forall a. Fractional a => a -> a -> a
/ a
det) forall (m :: * -> *) (r :: * -> *) a.
(Functor m, Functor r, Num a) =>
a -> m (r a) -> m (r a)
*!! forall a. a -> a -> V2 a
V2 (forall a. a -> a -> V2 a
V2 a
d (-a
b)) (forall a. a -> a -> V2 a
V2 (-a
c) a
a)
  where det :: a
det = forall a. Num a => M22 a -> a
det22 M22 a
m
{-# INLINE inv22 #-}

-- |3x3 matrix inverse.

--

-- >>> inv33 $ V3 (V3 1 2 4) (V3 4 2 2) (V3 1 1 1)

-- V3 (V3 0.0 0.5 (-1.0)) (V3 (-0.5) (-0.75) 3.5) (V3 0.5 0.25 (-1.5))

inv33 :: Fractional a => M33 a -> M33 a
inv33 :: forall a. Fractional a => M33 a -> M33 a
inv33 m :: M33 a
m@(V3 (V3 a
a a
b a
c)
            (V3 a
d a
e a
f)
            (V3 a
g a
h a
i))
  = (a
1 forall a. Fractional a => a -> a -> a
/ a
det) forall (m :: * -> *) (r :: * -> *) a.
(Functor m, Functor r, Num a) =>
a -> m (r a) -> m (r a)
*!! forall a. a -> a -> a -> V3 a
V3 (forall a. a -> a -> a -> V3 a
V3 a
a' a
b' a
c')
                     (forall a. a -> a -> a -> V3 a
V3 a
d' a
e' a
f')
                     (forall a. a -> a -> a -> V3 a
V3 a
g' a
h' a
i')
  where a' :: a
a' = forall {a}. Num a => (a, a, a, a) -> a
cofactor (a
e,a
f,a
h,a
i)
        b' :: a
b' = forall {a}. Num a => (a, a, a, a) -> a
cofactor (a
c,a
b,a
i,a
h)
        c' :: a
c' = forall {a}. Num a => (a, a, a, a) -> a
cofactor (a
b,a
c,a
e,a
f)
        d' :: a
d' = forall {a}. Num a => (a, a, a, a) -> a
cofactor (a
f,a
d,a
i,a
g)
        e' :: a
e' = forall {a}. Num a => (a, a, a, a) -> a
cofactor (a
a,a
c,a
g,a
i)
        f' :: a
f' = forall {a}. Num a => (a, a, a, a) -> a
cofactor (a
c,a
a,a
f,a
d)
        g' :: a
g' = forall {a}. Num a => (a, a, a, a) -> a
cofactor (a
d,a
e,a
g,a
h)
        h' :: a
h' = forall {a}. Num a => (a, a, a, a) -> a
cofactor (a
b,a
a,a
h,a
g)
        i' :: a
i' = forall {a}. Num a => (a, a, a, a) -> a
cofactor (a
a,a
b,a
d,a
e)
        cofactor :: (a, a, a, a) -> a
cofactor (a
q,a
r,a
s,a
t) = forall a. Num a => M22 a -> a
det22 (forall a. a -> a -> V2 a
V2 (forall a. a -> a -> V2 a
V2 a
q a
r) (forall a. a -> a -> V2 a
V2 a
s a
t))
        det :: a
det = forall a. Num a => M33 a -> a
det33 M33 a
m
{-# INLINE inv33 #-}


-- | 'transpose' is just an alias for 'distribute'

--

-- > transpose (V3 (V2 1 2) (V2 3 4) (V2 5 6))

-- V2 (V3 1 3 5) (V3 2 4 6)

transpose :: (Distributive g, Functor f) => f (g a) -> g (f a)
transpose :: forall (g :: * -> *) (f :: * -> *) a.
(Distributive g, Functor f) =>
f (g a) -> g (f a)
transpose = forall (g :: * -> *) (f :: * -> *) a.
(Distributive g, Functor f) =>
f (g a) -> g (f a)
distribute
{-# INLINE transpose #-}

-- |4x4 matrix inverse.

inv44 :: Fractional a => M44 a -> M44 a
inv44 :: forall a. Fractional a => M44 a -> M44 a
inv44 (V4 (V4 a
i00 a
i01 a
i02 a
i03)
          (V4 a
i10 a
i11 a
i12 a
i13)
          (V4 a
i20 a
i21 a
i22 a
i23)
          (V4 a
i30 a
i31 a
i32 a
i33)) =
  let s0 :: a
s0 = a
i00 forall a. Num a => a -> a -> a
* a
i11 forall a. Num a => a -> a -> a
- a
i10 forall a. Num a => a -> a -> a
* a
i01
      s1 :: a
s1 = a
i00 forall a. Num a => a -> a -> a
* a
i12 forall a. Num a => a -> a -> a
- a
i10 forall a. Num a => a -> a -> a
* a
i02
      s2 :: a
s2 = a
i00 forall a. Num a => a -> a -> a
* a
i13 forall a. Num a => a -> a -> a
- a
i10 forall a. Num a => a -> a -> a
* a
i03
      s3 :: a
s3 = a
i01 forall a. Num a => a -> a -> a
* a
i12 forall a. Num a => a -> a -> a
- a
i11 forall a. Num a => a -> a -> a
* a
i02
      s4 :: a
s4 = a
i01 forall a. Num a => a -> a -> a
* a
i13 forall a. Num a => a -> a -> a
- a
i11 forall a. Num a => a -> a -> a
* a
i03
      s5 :: a
s5 = a
i02 forall a. Num a => a -> a -> a
* a
i13 forall a. Num a => a -> a -> a
- a
i12 forall a. Num a => a -> a -> a
* a
i03
      c5 :: a
c5 = a
i22 forall a. Num a => a -> a -> a
* a
i33 forall a. Num a => a -> a -> a
- a
i32 forall a. Num a => a -> a -> a
* a
i23
      c4 :: a
c4 = a
i21 forall a. Num a => a -> a -> a
* a
i33 forall a. Num a => a -> a -> a
- a
i31 forall a. Num a => a -> a -> a
* a
i23
      c3 :: a
c3 = a
i21 forall a. Num a => a -> a -> a
* a
i32 forall a. Num a => a -> a -> a
- a
i31 forall a. Num a => a -> a -> a
* a
i22
      c2 :: a
c2 = a
i20 forall a. Num a => a -> a -> a
* a
i33 forall a. Num a => a -> a -> a
- a
i30 forall a. Num a => a -> a -> a
* a
i23
      c1 :: a
c1 = a
i20 forall a. Num a => a -> a -> a
* a
i32 forall a. Num a => a -> a -> a
- a
i30 forall a. Num a => a -> a -> a
* a
i22
      c0 :: a
c0 = a
i20 forall a. Num a => a -> a -> a
* a
i31 forall a. Num a => a -> a -> a
- a
i30 forall a. Num a => a -> a -> a
* a
i21
      det :: a
det = a
s0 forall a. Num a => a -> a -> a
* a
c5 forall a. Num a => a -> a -> a
- a
s1 forall a. Num a => a -> a -> a
* a
c4 forall a. Num a => a -> a -> a
+ a
s2 forall a. Num a => a -> a -> a
* a
c3 forall a. Num a => a -> a -> a
+ a
s3 forall a. Num a => a -> a -> a
* a
c2 forall a. Num a => a -> a -> a
- a
s4 forall a. Num a => a -> a -> a
* a
c1 forall a. Num a => a -> a -> a
+ a
s5 forall a. Num a => a -> a -> a
* a
c0
      invDet :: a
invDet = forall a. Fractional a => a -> a
recip a
det
  in a
invDet forall (m :: * -> *) (r :: * -> *) a.
(Functor m, Functor r, Num a) =>
a -> m (r a) -> m (r a)
*!! forall a. a -> a -> a -> a -> V4 a
V4 (forall a. a -> a -> a -> a -> V4 a
V4 (a
i11 forall a. Num a => a -> a -> a
* a
c5 forall a. Num a => a -> a -> a
- a
i12 forall a. Num a => a -> a -> a
* a
c4 forall a. Num a => a -> a -> a
+ a
i13 forall a. Num a => a -> a -> a
* a
c3)
                       (-a
i01 forall a. Num a => a -> a -> a
* a
c5 forall a. Num a => a -> a -> a
+ a
i02 forall a. Num a => a -> a -> a
* a
c4 forall a. Num a => a -> a -> a
- a
i03 forall a. Num a => a -> a -> a
* a
c3)
                       (a
i31 forall a. Num a => a -> a -> a
* a
s5 forall a. Num a => a -> a -> a
- a
i32 forall a. Num a => a -> a -> a
* a
s4 forall a. Num a => a -> a -> a
+ a
i33 forall a. Num a => a -> a -> a
* a
s3)
                       (-a
i21 forall a. Num a => a -> a -> a
* a
s5 forall a. Num a => a -> a -> a
+ a
i22 forall a. Num a => a -> a -> a
* a
s4 forall a. Num a => a -> a -> a
- a
i23 forall a. Num a => a -> a -> a
* a
s3))
                   (forall a. a -> a -> a -> a -> V4 a
V4 (-a
i10 forall a. Num a => a -> a -> a
* a
c5 forall a. Num a => a -> a -> a
+ a
i12 forall a. Num a => a -> a -> a
* a
c2 forall a. Num a => a -> a -> a
- a
i13 forall a. Num a => a -> a -> a
* a
c1)
                       (a
i00 forall a. Num a => a -> a -> a
* a
c5 forall a. Num a => a -> a -> a
- a
i02 forall a. Num a => a -> a -> a
* a
c2 forall a. Num a => a -> a -> a
+ a
i03 forall a. Num a => a -> a -> a
* a
c1)
                       (-a
i30 forall a. Num a => a -> a -> a
* a
s5 forall a. Num a => a -> a -> a
+ a
i32 forall a. Num a => a -> a -> a
* a
s2 forall a. Num a => a -> a -> a
- a
i33 forall a. Num a => a -> a -> a
* a
s1)
                       (a
i20 forall a. Num a => a -> a -> a
* a
s5 forall a. Num a => a -> a -> a
- a
i22 forall a. Num a => a -> a -> a
* a
s2 forall a. Num a => a -> a -> a
+ a
i23 forall a. Num a => a -> a -> a
* a
s1))
                   (forall a. a -> a -> a -> a -> V4 a
V4 (a
i10 forall a. Num a => a -> a -> a
* a
c4 forall a. Num a => a -> a -> a
- a
i11 forall a. Num a => a -> a -> a
* a
c2 forall a. Num a => a -> a -> a
+ a
i13 forall a. Num a => a -> a -> a
* a
c0)
                       (-a
i00 forall a. Num a => a -> a -> a
* a
c4 forall a. Num a => a -> a -> a
+ a
i01 forall a. Num a => a -> a -> a
* a
c2 forall a. Num a => a -> a -> a
- a
i03 forall a. Num a => a -> a -> a
* a
c0)
                       (a
i30 forall a. Num a => a -> a -> a
* a
s4 forall a. Num a => a -> a -> a
- a
i31 forall a. Num a => a -> a -> a
* a
s2 forall a. Num a => a -> a -> a
+ a
i33 forall a. Num a => a -> a -> a
* a
s0)
                       (-a
i20 forall a. Num a => a -> a -> a
* a
s4 forall a. Num a => a -> a -> a
+ a
i21 forall a. Num a => a -> a -> a
* a
s2 forall a. Num a => a -> a -> a
- a
i23 forall a. Num a => a -> a -> a
* a
s0))
                   (forall a. a -> a -> a -> a -> V4 a
V4 (-a
i10 forall a. Num a => a -> a -> a
* a
c3 forall a. Num a => a -> a -> a
+ a
i11 forall a. Num a => a -> a -> a
* a
c1 forall a. Num a => a -> a -> a
- a
i12 forall a. Num a => a -> a -> a
* a
c0)
                       (a
i00 forall a. Num a => a -> a -> a
* a
c3 forall a. Num a => a -> a -> a
- a
i01 forall a. Num a => a -> a -> a
* a
c1 forall a. Num a => a -> a -> a
+ a
i02 forall a. Num a => a -> a -> a
* a
c0)
                       (-a
i30 forall a. Num a => a -> a -> a
* a
s3 forall a. Num a => a -> a -> a
+ a
i31 forall a. Num a => a -> a -> a
* a
s1 forall a. Num a => a -> a -> a
- a
i32 forall a. Num a => a -> a -> a
* a
s0)
                       (a
i20 forall a. Num a => a -> a -> a
* a
s3 forall a. Num a => a -> a -> a
- a
i21 forall a. Num a => a -> a -> a
* a
s1 forall a. Num a => a -> a -> a
+ a
i22 forall a. Num a => a -> a -> a
* a
s0))
{-# INLINE inv44 #-}

-- | Compute the (L, U) decomposition of a square matrix using Crout's

--   algorithm. The 'Index' of the vectors must be 'Integral'.

lu :: ( Num a
      , Fractional a
      , Foldable m
      , Traversable m
      , Applicative m
      , Additive m
      , Ixed (m a)
      , Ixed (m (m a))
      , i ~ Index (m a)
      , i ~ Index (m (m a))
      , Eq i
      , Integral i
      , a ~ IxValue (m a)
      , m a ~ IxValue (m (m a))
      , Num (m a)
      )
   => m (m a)
   -> (m (m a), m (m a))
lu :: forall a (m :: * -> *) i.
(Num a, Fractional a, Foldable m, Traversable m, Applicative m,
 Additive m, Ixed (m a), Ixed (m (m a)), i ~ Index (m a),
 i ~ Index (m (m a)), Eq i, Integral i, a ~ IxValue (m a),
 m a ~ IxValue (m (m a)), Num (m a)) =>
m (m a) -> (m (m a), m (m a))
lu m (m a)
a =
    let n :: i
n = forall a b. (Integral a, Num b) => a -> b
fromIntegral (forall (t :: * -> *) a. Foldable t => t a -> Int
length m (m a)
a)
        initU :: m (m a)
initU = forall a (t :: * -> *).
(Num a, Traversable t, Applicative t) =>
t (t a)
identity
        initL :: m (m a)
initL = forall (f :: * -> *) a. (Additive f, Num a) => f a
zero
        buildLVal :: i -> i -> m (m a) -> m (m a) -> m (m a)
buildLVal !i
i !i
j !m (m a)
l !m (m a)
u =
            let go :: i -> a -> a
go !i
k !a
s
                    | i
k forall a. Eq a => a -> a -> Bool
== i
j = a
s
                    | Bool
otherwise = i -> a -> a
go (i
kforall a. Num a => a -> a -> a
+i
1)
                                     ( a
s
                                      forall a. Num a => a -> a -> a
+ ( (m (m a)
l forall s a. HasCallStack => s -> Getting (Endo a) s a -> a
^?! forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
i forall s a. HasCallStack => s -> Getting (Endo a) s a -> a
^?! forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
k)
                                        forall a. Num a => a -> a -> a
* (m (m a)
u forall s a. HasCallStack => s -> Getting (Endo a) s a -> a
^?! forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
k forall s a. HasCallStack => s -> Getting (Endo a) s a -> a
^?! forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
j)
                                        )
                                      )
                s' :: a
s' = i -> a -> a
go i
0 a
0
            in m (m a)
l forall a b. a -> (a -> b) -> b
& (forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
i forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
j) forall s t a b. ASetter s t a b -> b -> s -> t
.~ ((m (m a)
a forall s a. HasCallStack => s -> Getting (Endo a) s a -> a
^?! forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
i forall s a. HasCallStack => s -> Getting (Endo a) s a -> a
^?! forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
j) forall a. Num a => a -> a -> a
- a
s')
        buildL :: i -> i -> m (m a) -> m (m a) -> m (m a)
buildL !i
i !i
j !m (m a)
l !m (m a)
u
            | i
i forall a. Eq a => a -> a -> Bool
== i
n = m (m a)
l
            | Bool
otherwise = i -> i -> m (m a) -> m (m a) -> m (m a)
buildL (i
iforall a. Num a => a -> a -> a
+i
1) i
j (i -> i -> m (m a) -> m (m a) -> m (m a)
buildLVal i
i i
j m (m a)
l m (m a)
u) m (m a)
u
        buildUVal :: i -> i -> m (m a) -> m (m a) -> m (m a)
buildUVal !i
i !i
j !m (m a)
l !m (m a)
u =
            let go :: i -> a -> a
go !i
k !a
s
                    | i
k forall a. Eq a => a -> a -> Bool
== i
j = a
s
                    | Bool
otherwise = i -> a -> a
go (i
kforall a. Num a => a -> a -> a
+i
1)
                                     ( a
s
                                     forall a. Num a => a -> a -> a
+ ( (m (m a)
l forall s a. HasCallStack => s -> Getting (Endo a) s a -> a
^?! forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
j forall s a. HasCallStack => s -> Getting (Endo a) s a -> a
^?! forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
k)
                                       forall a. Num a => a -> a -> a
* (m (m a)
u forall s a. HasCallStack => s -> Getting (Endo a) s a -> a
^?! forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
k forall s a. HasCallStack => s -> Getting (Endo a) s a -> a
^?! forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
i)
                                       )
                                     )
                s' :: a
s' = i -> a -> a
go i
0 a
0
            in m (m a)
u forall a b. a -> (a -> b) -> b
& (forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
j forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
i) forall s t a b. ASetter s t a b -> b -> s -> t
.~ ( ((m (m a)
a forall s a. HasCallStack => s -> Getting (Endo a) s a -> a
^?! forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
j forall s a. HasCallStack => s -> Getting (Endo a) s a -> a
^?! forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
i) forall a. Num a => a -> a -> a
- a
s')
                                    forall a. Fractional a => a -> a -> a
/ (m (m a)
l forall s a. HasCallStack => s -> Getting (Endo a) s a -> a
^?! forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
j forall s a. HasCallStack => s -> Getting (Endo a) s a -> a
^?! forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
j)
                                    )
        buildU :: i -> i -> m (m a) -> m (m a) -> m (m a)
buildU !i
i !i
j !m (m a)
l !m (m a)
u
            | i
i forall a. Eq a => a -> a -> Bool
== i
n = m (m a)
u
            | Bool
otherwise = i -> i -> m (m a) -> m (m a) -> m (m a)
buildU (i
iforall a. Num a => a -> a -> a
+i
1) i
j m (m a)
l (i -> i -> m (m a) -> m (m a) -> m (m a)
buildUVal i
i i
j m (m a)
l m (m a)
u)
        buildLU :: i -> m (m a) -> m (m a) -> (m (m a), m (m a))
buildLU !i
j !m (m a)
l !m (m a)
u
            | i
j forall a. Eq a => a -> a -> Bool
== i
n = (m (m a)
l, m (m a)
u)
            | Bool
otherwise =
                let l' :: m (m a)
l' = i -> i -> m (m a) -> m (m a) -> m (m a)
buildL i
j i
j m (m a)
l m (m a)
u
                    u' :: m (m a)
u' = i -> i -> m (m a) -> m (m a) -> m (m a)
buildU i
j i
j m (m a)
l' m (m a)
u
                in i -> m (m a) -> m (m a) -> (m (m a), m (m a))
buildLU (i
jforall a. Num a => a -> a -> a
+i
1) m (m a)
l' m (m a)
u'
    in i -> m (m a) -> m (m a) -> (m (m a), m (m a))
buildLU i
0 m (m a)
initL m (m a)
initU

-- | Compute the (L, U) decomposition of a square matrix using Crout's

--   algorithm, using the vector's 'Finite' instance to provide an index.

luFinite :: ( Num a
            , Fractional a
            , Functor m
            , Finite m
            , n ~ Size m
            , KnownNat n
            , Num (m a)
            )
         => m (m a)
         -> (m (m a), m (m a))
luFinite :: forall a (m :: * -> *) (n :: Nat).
(Num a, Fractional a, Functor m, Finite m, n ~ Size m, KnownNat n,
 Num (m a)) =>
m (m a) -> (m (m a), m (m a))
luFinite m (m a)
a =
    forall (p :: * -> * -> *) a b c d.
Bifunctor p =>
(a -> b) -> (c -> d) -> p a c -> p b d
bimap (forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap forall (v :: * -> *) a. Finite v => V (Size v) a -> v a
fromV forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall (v :: * -> *) a. Finite v => V (Size v) a -> v a
fromV)
          (forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap forall (v :: * -> *) a. Finite v => V (Size v) a -> v a
fromV forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall (v :: * -> *) a. Finite v => V (Size v) a -> v a
fromV)
          (forall a (m :: * -> *) i.
(Num a, Fractional a, Foldable m, Traversable m, Applicative m,
 Additive m, Ixed (m a), Ixed (m (m a)), i ~ Index (m a),
 i ~ Index (m (m a)), Eq i, Integral i, a ~ IxValue (m a),
 m a ~ IxValue (m (m a)), Num (m a)) =>
m (m a) -> (m (m a), m (m a))
lu (forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap forall (v :: * -> *) a. Finite v => v a -> V (Size v) a
toV (forall (v :: * -> *) a. Finite v => v a -> V (Size v) a
toV m (m a)
a)))

-- | Solve a linear system with a lower-triangular matrix of coefficients with

--   forwards substitution.

forwardSub :: ( Num a
              , Fractional a
              , Foldable m
              , Additive m
              , Ixed (m a)
              , Ixed (m (m a))
              , i ~ Index (m a)
              , i ~ Index (m (m a))
              , Eq i
              , Ord i
              , Integral i
              , a ~ IxValue (m a)
              , m a ~ IxValue (m (m a))
              )
           => m (m a)
           -> m a
           -> m a
forwardSub :: forall a (m :: * -> *) i.
(Num a, Fractional a, Foldable m, Additive m, Ixed (m a),
 Ixed (m (m a)), i ~ Index (m a), i ~ Index (m (m a)), Eq i, Ord i,
 Integral i, a ~ IxValue (m a), m a ~ IxValue (m (m a))) =>
m (m a) -> m a -> m a
forwardSub m (m a)
a m a
b =
    let n :: i
n = forall a b. (Integral a, Num b) => a -> b
fromIntegral (forall (t :: * -> *) a. Foldable t => t a -> Int
length m a
b)
        initX :: m a
initX = forall (f :: * -> *) a. (Additive f, Num a) => f a
zero
        coeff :: i -> i -> a -> m a -> a
coeff !i
i !i
j !a
s !m a
x
            | i
j forall a. Eq a => a -> a -> Bool
== i
i = a
s
            | Bool
otherwise = i -> i -> a -> m a -> a
coeff i
i (i
jforall a. Num a => a -> a -> a
+i
1) (a
s forall a. Num a => a -> a -> a
+ ((m (m a)
a forall s a. HasCallStack => s -> Getting (Endo a) s a -> a
^?! forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
i forall s a. HasCallStack => s -> Getting (Endo a) s a -> a
^?! forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
j) forall a. Num a => a -> a -> a
* (m a
x forall s a. HasCallStack => s -> Getting (Endo a) s a -> a
^?! forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
j))) m a
x
        go :: i -> m a -> m a
go !i
i !m a
x
            | i
i forall a. Eq a => a -> a -> Bool
== i
n = m a
x
            | Bool
otherwise = i -> m a -> m a
go (i
i forall a. Num a => a -> a -> a
+ i
1) (m a
x forall a b. a -> (a -> b) -> b
& forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
i forall s t a b. ASetter s t a b -> b -> s -> t
.~ ( ((m a
b forall s a. HasCallStack => s -> Getting (Endo a) s a -> a
^?! forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
i) forall a. Num a => a -> a -> a
- i -> i -> a -> m a -> a
coeff i
i i
0 a
0 m a
x)
                                                  forall a. Fractional a => a -> a -> a
/ (m (m a)
a forall s a. HasCallStack => s -> Getting (Endo a) s a -> a
^?! forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
i forall s a. HasCallStack => s -> Getting (Endo a) s a -> a
^?! forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
i)
                                                  ))
    in i -> m a -> m a
go i
0 m a
initX

-- | Solve a linear system with a lower-triangular matrix of coefficients with

--   forwards substitution, using the vector's 'Finite' instance to provide an

--   index.

forwardSubFinite :: ( Num a
                    , Fractional a
                    , Foldable m
                    , n ~ Size m
                    , KnownNat n
                    , Additive m
                    , Finite m
                    )
                 => m (m a)
                 -> m a
                 -> m a
forwardSubFinite :: forall a (m :: * -> *) (n :: Nat).
(Num a, Fractional a, Foldable m, n ~ Size m, KnownNat n,
 Additive m, Finite m) =>
m (m a) -> m a -> m a
forwardSubFinite m (m a)
a m a
b = forall (v :: * -> *) a. Finite v => V (Size v) a -> v a
fromV (forall a (m :: * -> *) i.
(Num a, Fractional a, Foldable m, Additive m, Ixed (m a),
 Ixed (m (m a)), i ~ Index (m a), i ~ Index (m (m a)), Eq i, Ord i,
 Integral i, a ~ IxValue (m a), m a ~ IxValue (m (m a))) =>
m (m a) -> m a -> m a
forwardSub (forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap forall (v :: * -> *) a. Finite v => v a -> V (Size v) a
toV (forall (v :: * -> *) a. Finite v => v a -> V (Size v) a
toV m (m a)
a)) (forall (v :: * -> *) a. Finite v => v a -> V (Size v) a
toV m a
b))

-- | Solve a linear system with an upper-triangular matrix of coefficients with

--   backwards substitution.

backwardSub :: ( Num a
               , Fractional a
               , Foldable m
               , Additive m
               , Ixed (m a)
               , Ixed (m (m a))
               , i ~ Index (m a)
               , i ~ Index (m (m a))
               , Eq i
               , Ord i
               , Integral i
               , a ~ IxValue (m a)
               , m a ~ IxValue (m (m a))
               )
            => m (m a)
            -> m a
            -> m a
backwardSub :: forall a (m :: * -> *) i.
(Num a, Fractional a, Foldable m, Additive m, Ixed (m a),
 Ixed (m (m a)), i ~ Index (m a), i ~ Index (m (m a)), Eq i, Ord i,
 Integral i, a ~ IxValue (m a), m a ~ IxValue (m (m a))) =>
m (m a) -> m a -> m a
backwardSub m (m a)
a m a
b =
    let n :: i
n = forall a b. (Integral a, Num b) => a -> b
fromIntegral (forall (t :: * -> *) a. Foldable t => t a -> Int
length m a
b)
        initX :: m a
initX = forall (f :: * -> *) a. (Additive f, Num a) => f a
zero
        coeff :: i -> i -> a -> m a -> a
coeff !i
i !i
j !a
s !m a
x
            | i
j forall a. Eq a => a -> a -> Bool
== i
n = a
s
            | Bool
otherwise = i -> i -> a -> m a -> a
coeff i
i
                                (i
jforall a. Num a => a -> a -> a
+i
1)
                                (a
s forall a. Num a => a -> a -> a
+ ((m (m a)
a forall s a. HasCallStack => s -> Getting (Endo a) s a -> a
^?! forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
i forall s a. HasCallStack => s -> Getting (Endo a) s a -> a
^?! forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
j) forall a. Num a => a -> a -> a
* (m a
x forall s a. HasCallStack => s -> Getting (Endo a) s a -> a
^?! forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
j)))
                                m a
x
        go :: i -> m a -> m a
go !i
i !m a
x
            | i
i forall a. Ord a => a -> a -> Bool
< i
0 = m a
x
            | Bool
otherwise = i -> m a -> m a
go (i
iforall a. Num a => a -> a -> a
-i
1)
                             (m a
x forall a b. a -> (a -> b) -> b
& forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
i forall s t a b. ASetter s t a b -> b -> s -> t
.~ ( ((m a
b forall s a. HasCallStack => s -> Getting (Endo a) s a -> a
^?! forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
i) forall a. Num a => a -> a -> a
- i -> i -> a -> m a -> a
coeff i
i (i
iforall a. Num a => a -> a -> a
+i
1) a
0 m a
x)
                                          forall a. Fractional a => a -> a -> a
/ (m (m a)
a forall s a. HasCallStack => s -> Getting (Endo a) s a -> a
^?! forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
i forall s a. HasCallStack => s -> Getting (Endo a) s a -> a
^?! forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
i)
                                          ))
    in i -> m a -> m a
go (i
nforall a. Num a => a -> a -> a
-i
1) m a
initX

-- | Solve a linear system with an upper-triangular matrix of coefficients with

--   backwards substitution, using the vector's 'Finite' instance to provide an

--   index.

backwardSubFinite :: ( Num a
                     , Fractional a
                     , Foldable m
                     , n ~ Size m
                     , KnownNat n
                     , Additive m
                     , Finite m
                     )
                  => m (m a)
                  -> m a
                  -> m a
backwardSubFinite :: forall a (m :: * -> *) (n :: Nat).
(Num a, Fractional a, Foldable m, n ~ Size m, KnownNat n,
 Additive m, Finite m) =>
m (m a) -> m a -> m a
backwardSubFinite m (m a)
a m a
b = forall (v :: * -> *) a. Finite v => V (Size v) a -> v a
fromV (forall a (m :: * -> *) i.
(Num a, Fractional a, Foldable m, Additive m, Ixed (m a),
 Ixed (m (m a)), i ~ Index (m a), i ~ Index (m (m a)), Eq i, Ord i,
 Integral i, a ~ IxValue (m a), m a ~ IxValue (m (m a))) =>
m (m a) -> m a -> m a
backwardSub (forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap forall (v :: * -> *) a. Finite v => v a -> V (Size v) a
toV (forall (v :: * -> *) a. Finite v => v a -> V (Size v) a
toV m (m a)
a)) (forall (v :: * -> *) a. Finite v => v a -> V (Size v) a
toV m a
b))

-- | Solve a linear system with LU decomposition.

luSolve :: ( Num a
           , Fractional a
           , Foldable m
           , Traversable m
           , Applicative m
           , Additive m
           , Ixed (m a)
           , Ixed (m (m a))
           , i ~ Index (m a)
           , i ~ Index (m (m a))
           , Eq i
           , Integral i
           , a ~ IxValue (m a)
           , m a ~ IxValue (m (m a))
           , Num (m a)
           )
        => m (m a)
        -> m a
        -> m a
luSolve :: forall a (m :: * -> *) i.
(Num a, Fractional a, Foldable m, Traversable m, Applicative m,
 Additive m, Ixed (m a), Ixed (m (m a)), i ~ Index (m a),
 i ~ Index (m (m a)), Eq i, Integral i, a ~ IxValue (m a),
 m a ~ IxValue (m (m a)), Num (m a)) =>
m (m a) -> m a -> m a
luSolve m (m a)
a m a
b =
    let (m (m a)
l, m (m a)
u) = forall a (m :: * -> *) i.
(Num a, Fractional a, Foldable m, Traversable m, Applicative m,
 Additive m, Ixed (m a), Ixed (m (m a)), i ~ Index (m a),
 i ~ Index (m (m a)), Eq i, Integral i, a ~ IxValue (m a),
 m a ~ IxValue (m (m a)), Num (m a)) =>
m (m a) -> (m (m a), m (m a))
lu m (m a)
a
    in forall a (m :: * -> *) i.
(Num a, Fractional a, Foldable m, Additive m, Ixed (m a),
 Ixed (m (m a)), i ~ Index (m a), i ~ Index (m (m a)), Eq i, Ord i,
 Integral i, a ~ IxValue (m a), m a ~ IxValue (m (m a))) =>
m (m a) -> m a -> m a
backwardSub m (m a)
u (forall a (m :: * -> *) i.
(Num a, Fractional a, Foldable m, Additive m, Ixed (m a),
 Ixed (m (m a)), i ~ Index (m a), i ~ Index (m (m a)), Eq i, Ord i,
 Integral i, a ~ IxValue (m a), m a ~ IxValue (m (m a))) =>
m (m a) -> m a -> m a
forwardSub m (m a)
l m a
b)

-- | Solve a linear system with LU decomposition, using the vector's 'Finite'

--   instance to provide an index.

luSolveFinite :: ( Num a
                 , Fractional a
                 , Functor m
                 , Finite m
                 , n ~ Size m
                 , KnownNat n
                 , Num (m a)
                 )
              => m (m a)
              -> m a
              -> m a
luSolveFinite :: forall a (m :: * -> *) (n :: Nat).
(Num a, Fractional a, Functor m, Finite m, n ~ Size m, KnownNat n,
 Num (m a)) =>
m (m a) -> m a -> m a
luSolveFinite m (m a)
a m a
b = forall (v :: * -> *) a. Finite v => V (Size v) a -> v a
fromV (forall a (m :: * -> *) i.
(Num a, Fractional a, Foldable m, Traversable m, Applicative m,
 Additive m, Ixed (m a), Ixed (m (m a)), i ~ Index (m a),
 i ~ Index (m (m a)), Eq i, Integral i, a ~ IxValue (m a),
 m a ~ IxValue (m (m a)), Num (m a)) =>
m (m a) -> m a -> m a
luSolve (forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap forall (v :: * -> *) a. Finite v => v a -> V (Size v) a
toV (forall (v :: * -> *) a. Finite v => v a -> V (Size v) a
toV m (m a)
a)) (forall (v :: * -> *) a. Finite v => v a -> V (Size v) a
toV m a
b))

-- | Invert a matrix with LU decomposition.

luInv :: ( Num a
         , Fractional a
         , Foldable m
         , Traversable m
         , Applicative m
         , Additive m
         , Distributive m
         , Ixed (m a)
         , Ixed (m (m a))
         , i ~ Index (m a)
         , i ~ Index (m (m a))
         , Eq i
         , Integral i
         , a ~ IxValue (m a)
         , m a ~ IxValue (m (m a))
         , Num (m a)
         )
      => m (m a)
      -> m (m a)
luInv :: forall a (m :: * -> *) i.
(Num a, Fractional a, Foldable m, Traversable m, Applicative m,
 Additive m, Distributive m, Ixed (m a), Ixed (m (m a)),
 i ~ Index (m a), i ~ Index (m (m a)), Eq i, Integral i,
 a ~ IxValue (m a), m a ~ IxValue (m (m a)), Num (m a)) =>
m (m a) -> m (m a)
luInv m (m a)
a =
    let n :: i
n = forall a b. (Integral a, Num b) => a -> b
fromIntegral (forall (t :: * -> *) a. Foldable t => t a -> Int
length m (m a)
a)
        initA' :: m (m a)
initA' = forall (f :: * -> *) a. (Additive f, Num a) => f a
zero
        (m (m a)
l, m (m a)
u) = forall a (m :: * -> *) i.
(Num a, Fractional a, Foldable m, Traversable m, Applicative m,
 Additive m, Ixed (m a), Ixed (m (m a)), i ~ Index (m a),
 i ~ Index (m (m a)), Eq i, Integral i, a ~ IxValue (m a),
 m a ~ IxValue (m (m a)), Num (m a)) =>
m (m a) -> (m (m a), m (m a))
lu m (m a)
a
        go :: i -> m (m a) -> m (m a)
go !i
i !m (m a)
a'
            | i
i forall a. Eq a => a -> a -> Bool
== i
n = m (m a)
a'
            | Bool
otherwise = let e :: m a
e   = forall (f :: * -> *) a. (Additive f, Num a) => f a
zero forall a b. a -> (a -> b) -> b
& forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
i forall s t a b. ASetter s t a b -> b -> s -> t
.~ a
1
                              a'r :: m a
a'r = forall a (m :: * -> *) i.
(Num a, Fractional a, Foldable m, Additive m, Ixed (m a),
 Ixed (m (m a)), i ~ Index (m a), i ~ Index (m (m a)), Eq i, Ord i,
 Integral i, a ~ IxValue (m a), m a ~ IxValue (m (m a))) =>
m (m a) -> m a -> m a
backwardSub m (m a)
u (forall a (m :: * -> *) i.
(Num a, Fractional a, Foldable m, Additive m, Ixed (m a),
 Ixed (m (m a)), i ~ Index (m a), i ~ Index (m (m a)), Eq i, Ord i,
 Integral i, a ~ IxValue (m a), m a ~ IxValue (m (m a))) =>
m (m a) -> m a -> m a
forwardSub m (m a)
l m a
e)
                          in i -> m (m a) -> m (m a)
go (i
iforall a. Num a => a -> a -> a
+i
1) (m (m a)
a' forall a b. a -> (a -> b) -> b
& forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
i forall s t a b. ASetter s t a b -> b -> s -> t
.~ m a
a'r)
    in forall (g :: * -> *) (f :: * -> *) a.
(Distributive g, Functor f) =>
f (g a) -> g (f a)
transpose (i -> m (m a) -> m (m a)
go i
0 m (m a)
initA')

-- | Invert a matrix with LU decomposition, using the vector's 'Finite' instance

--   to provide an index.

luInvFinite :: ( Num a
               , Fractional a
               , Functor m
               , Finite m
               , n ~ Size m
               , KnownNat n
               , Num (m a)
               )
            => m (m a)
            -> m (m a)
luInvFinite :: forall a (m :: * -> *) (n :: Nat).
(Num a, Fractional a, Functor m, Finite m, n ~ Size m, KnownNat n,
 Num (m a)) =>
m (m a) -> m (m a)
luInvFinite m (m a)
a = forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap forall (v :: * -> *) a. Finite v => V (Size v) a -> v a
fromV (forall (v :: * -> *) a. Finite v => V (Size v) a -> v a
fromV (forall a (m :: * -> *) i.
(Num a, Fractional a, Foldable m, Traversable m, Applicative m,
 Additive m, Distributive m, Ixed (m a), Ixed (m (m a)),
 i ~ Index (m a), i ~ Index (m (m a)), Eq i, Integral i,
 a ~ IxValue (m a), m a ~ IxValue (m (m a)), Num (m a)) =>
m (m a) -> m (m a)
luInv (forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap forall (v :: * -> *) a. Finite v => v a -> V (Size v) a
toV (forall (v :: * -> *) a. Finite v => v a -> V (Size v) a
toV m (m a)
a))))

-- | Compute the determinant of a matrix using LU decomposition.

luDet :: ( Num a
         , Fractional a
         , Foldable m
         , Traversable m
         , Applicative m
         , Additive m
         , Trace m
         , Ixed (m a)
         , Ixed (m (m a))
         , i ~ Index (m a)
         , i ~ Index (m (m a))
         , Eq i
         , Integral i
         , a ~ IxValue (m a)
         , m a ~ IxValue (m (m a))
         , Num (m a)
         )
      => m (m a)
      -> a
luDet :: forall a (m :: * -> *) i.
(Num a, Fractional a, Foldable m, Traversable m, Applicative m,
 Additive m, Trace m, Ixed (m a), Ixed (m (m a)), i ~ Index (m a),
 i ~ Index (m (m a)), Eq i, Integral i, a ~ IxValue (m a),
 m a ~ IxValue (m (m a)), Num (m a)) =>
m (m a) -> a
luDet m (m a)
a =
    let (m (m a)
l, m (m a)
u) = forall a (m :: * -> *) i.
(Num a, Fractional a, Foldable m, Traversable m, Applicative m,
 Additive m, Ixed (m a), Ixed (m (m a)), i ~ Index (m a),
 i ~ Index (m (m a)), Eq i, Integral i, a ~ IxValue (m a),
 m a ~ IxValue (m (m a)), Num (m a)) =>
m (m a) -> (m (m a), m (m a))
lu m (m a)
a
        p :: m a -> a
p      = forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
Foldable.foldl forall a. Num a => a -> a -> a
(*) a
1
    in m a -> a
p (forall (m :: * -> *) a. Trace m => m (m a) -> m a
diagonal m (m a)
l) forall a. Num a => a -> a -> a
* m a -> a
p (forall (m :: * -> *) a. Trace m => m (m a) -> m a
diagonal m (m a)
u)

-- | Compute the determinant of a matrix using LU decomposition, using the

--   vector's 'Finite' instance to provide an index.

luDetFinite :: ( Num a
               , Fractional a
               , Functor m
               , Finite m
               , n ~ Size m
               , KnownNat n
               , Num (m a)
               )
            => m (m a)
            -> a
luDetFinite :: forall a (m :: * -> *) (n :: Nat).
(Num a, Fractional a, Functor m, Finite m, n ~ Size m, KnownNat n,
 Num (m a)) =>
m (m a) -> a
luDetFinite = forall a (m :: * -> *) i.
(Num a, Fractional a, Foldable m, Traversable m, Applicative m,
 Additive m, Trace m, Ixed (m a), Ixed (m (m a)), i ~ Index (m a),
 i ~ Index (m (m a)), Eq i, Integral i, a ~ IxValue (m a),
 m a ~ IxValue (m (m a)), Num (m a)) =>
m (m a) -> a
luDet forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap forall (v :: * -> *) a. Finite v => v a -> V (Size v) a
toV forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall (v :: * -> *) a. Finite v => v a -> V (Size v) a
toV