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

---------------------------------------------------------------------------
-- |
-- 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 :: 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 (f b -> f t) -> f (f b) -> f (f t)
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 ((b -> b) -> a -> Context a b b
forall a b t. (b -> t) -> a -> Context a b t
Context b -> b
forall a. a -> a
id)
   i :: f a
i = (Rep f -> a) -> f a
forall (f :: * -> *) a. Representable f => (Rep f -> a) -> f a
tabulate ((Rep f -> a) -> f a) -> (Rep f -> a) -> f a
forall a b. (a -> b) -> a -> b
$ \ Rep f
e -> Context a b t -> a
forall (w :: * -> * -> * -> *) a c t.
IndexedComonadStore w =>
w a c t -> a
ipos (Context a b t -> a) -> Context a b t -> a
forall a b. (a -> b) -> a -> b
$ s -> Context a b t
go (f s -> Rep f -> s
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 = (Rep f -> t) -> f t
forall (f :: * -> *) a. Representable f => (Rep f -> a) -> f a
tabulate ((Rep f -> t) -> f t) -> (Rep f -> t) -> f t
forall a b. (a -> b) -> a -> b
$ \ Rep f
e -> b -> Context a b t -> t
forall (w :: * -> * -> * -> *) c a t.
IndexedComonadStore w =>
c -> w a c t -> t
ipeek (f b -> Rep f -> b
forall (f :: * -> *) a. Representable f => f a -> Rep f -> a
index f b
eb Rep f
e) (s -> Context a b t
go (f s -> Rep f -> s
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 !*! :: m (t a) -> t (n a) -> m (n a)
!*! t (n a)
g = (t a -> n a) -> m (t a) -> m (n a)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (\ t a
f' -> (n a -> n a -> n a) -> n a -> t (n a) -> n a
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
Foldable.foldl' n a -> n a -> n a
forall (f :: * -> *) a. (Additive f, Num a) => f a -> f a -> f a
(^+^) n a
forall (f :: * -> *) a. (Additive f, Num a) => f a
zero (t (n a) -> n a) -> t (n a) -> n a
forall a b. (a -> b) -> a -> b
$ (a -> n a -> n a) -> t a -> t (n a) -> t (n a)
forall (f :: * -> *) a b c.
Additive f =>
(a -> b -> c) -> f a -> f b -> f c
liftI2 a -> n a -> n a
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 !+! :: m (n a) -> m (n a) -> m (n a)
!+! m (n a)
bs = (n a -> n a -> n a) -> m (n a) -> m (n a) -> m (n a)
forall (f :: * -> *) a.
Additive f =>
(a -> a -> a) -> f a -> f a -> f a
liftU2 n a -> n a -> n a
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 !-! :: m (n a) -> m (n a) -> m (n a)
!-! m (n a)
bs = (n a -> n a -> n a) -> m (n a) -> m (n a) -> m (n a)
forall (f :: * -> *) a.
Additive f =>
(a -> a -> a) -> f a -> f a -> f a
liftU2 n a -> n a -> n a
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 !* :: m (r a) -> r a -> m a
!* r a
v = (r a -> a) -> m (r a) -> m a
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (\r a
r -> r a -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
Foldable.sum (r a -> a) -> r a -> a
forall a b. (a -> b) -> a -> b
$ (a -> a -> a) -> r a -> r a -> r a
forall (f :: * -> *) a b c.
Additive f =>
(a -> b -> c) -> f a -> f b -> f c
liftI2 a -> a -> a
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 *! :: t a -> t (f a) -> f a
*! t (f a)
g = t (f a) -> f a
forall (f :: * -> *) (v :: * -> *) a.
(Foldable f, Additive v, Num a) =>
f (v a) -> v a
sumV (t (f a) -> f a) -> t (f a) -> f a
forall a b. (a -> b) -> a -> b
$ (a -> f a -> f a) -> t a -> t (f a) -> t (f a)
forall (f :: * -> *) a b c.
Additive f =>
(a -> b -> c) -> f a -> f b -> f c
liftI2 a -> f a -> f a
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 *!! :: a -> m (r a) -> m (r a)
*!! m (r a)
m = (r a -> r a) -> m (r a) -> m (r a)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (a
s a -> r a -> r a
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)
!!* :: m (r a) -> a -> m (r a)
(!!*) = (a -> m (r a) -> m (r a)) -> m (r a) -> a -> m (r a)
forall a b c. (a -> b -> c) -> b -> a -> c
flip a -> m (r a) -> m (r a)
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 !!/ :: m (r a) -> a -> m (r a)
!!/ a
s = (r a -> r a) -> m (r a) -> m (r a)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (r a -> a -> r a
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 :: m (n a) -> n (m a)
adjoint = (n a -> n a) -> m (n a) -> n (m a)
forall (g :: * -> *) (f :: * -> *) a b.
(Distributive g, Functor f) =>
(a -> g b) -> f a -> g (f b)
collect ((a -> a) -> n a -> n a
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap a -> a
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 :: Quaternion a -> M33 a
fromQuaternion (Quaternion a
w (V3 a
x a
y a
z)) =
  V3 a -> V3 a -> V3 a -> M33 a
forall a. a -> a -> a -> V3 a
V3 (a -> a -> a -> V3 a
forall a. a -> a -> a -> V3 a
V3 (a
1a -> a -> a
forall a. Num a => a -> a -> a
-a
2a -> a -> a
forall a. Num a => a -> a -> a
*(a
y2a -> a -> a
forall a. Num a => a -> a -> a
+a
z2)) (a
2a -> a -> a
forall a. Num a => a -> a -> a
*(a
xya -> a -> a
forall a. Num a => a -> a -> a
-a
zw)) (a
2a -> a -> a
forall a. Num a => a -> a -> a
*(a
xza -> a -> a
forall a. Num a => a -> a -> a
+a
yw)))
     (a -> a -> a -> V3 a
forall a. a -> a -> a -> V3 a
V3 (a
2a -> a -> a
forall a. Num a => a -> a -> a
*(a
xya -> a -> a
forall a. Num a => a -> a -> a
+a
zw)) (a
1a -> a -> a
forall a. Num a => a -> a -> a
-a
2a -> a -> a
forall a. Num a => a -> a -> a
*(a
x2a -> a -> a
forall a. Num a => a -> a -> a
+a
z2)) (a
2a -> a -> a
forall a. Num a => a -> a -> a
*(a
yza -> a -> a
forall a. Num a => a -> a -> a
-a
xw)))
     (a -> a -> a -> V3 a
forall a. a -> a -> a -> V3 a
V3 (a
2a -> a -> a
forall a. Num a => a -> a -> a
*(a
xza -> a -> a
forall a. Num a => a -> a -> a
-a
yw)) (a
2a -> a -> a
forall a. Num a => a -> a -> a
*(a
yza -> a -> a
forall a. Num a => a -> a -> a
+a
xw)) (a
1a -> a -> a
forall a. Num a => a -> a -> a
-a
2a -> a -> a
forall a. Num a => a -> a -> a
*(a
x2a -> a -> a
forall a. Num a => a -> a -> a
+a
y2)))
  where x2 :: a
x2 = a
xa -> a -> a
forall a. Num a => a -> a -> a
*a
x
        y2 :: a
y2 = a
ya -> a -> a
forall a. Num a => a -> a -> a
*a
y
        z2 :: a
z2 = a
za -> a -> a
forall a. Num a => a -> a -> a
*a
z
        xy :: a
xy = a
xa -> a -> a
forall a. Num a => a -> a -> a
*a
y
        xz :: a
xz = a
xa -> a -> a
forall a. Num a => a -> a -> a
*a
z
        xw :: a
xw = a
xa -> a -> a
forall a. Num a => a -> a -> a
*a
w
        yz :: a
yz = a
ya -> a -> a
forall a. Num a => a -> a -> a
*a
z
        yw :: a
yw = a
ya -> a -> a
forall a. Num a => a -> a -> a
*a
w
        zw :: a
zw = a
za -> a -> a
forall 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 :: M33 a -> V3 a -> M44 a
mkTransformationMat (V3 V3 a
r1 V3 a
r2 V3 a
r3) (V3 a
tx a
ty a
tz) =
  V4 a -> V4 a -> V4 a -> V4 a -> M44 a
forall a. a -> a -> a -> a -> V4 a
V4 (V3 a -> a -> V4 a
forall a. V3 a -> a -> V4 a
snoc3 V3 a
r1 a
tx) (V3 a -> a -> V4 a
forall a. V3 a -> a -> V4 a
snoc3 V3 a
r2 a
ty) (V3 a -> a -> V4 a
forall a. V3 a -> a -> V4 a
snoc3 V3 a
r3 a
tz) (a -> a -> a -> a -> V4 a
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) = a -> a -> a -> a -> V4 a
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 :: Quaternion a -> V3 a -> M44 a
mkTransformation = M33 a -> V3 a -> M44 a
forall a. Num a => M33 a -> V3 a -> M44 a
mkTransformationMat (M33 a -> V3 a -> M44 a)
-> (Quaternion a -> M33 a) -> Quaternion a -> V3 a -> M44 a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Quaternion a -> M33 a
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 :: 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)) =
  V4 a -> V4 a -> V4 a -> V4 a -> M44 a
forall a. a -> a -> a -> a -> V4 a
V4 (a -> a -> a -> a -> V4 a
forall a. a -> a -> a -> a -> V4 a
V4 a
a a
b a
c a
0)
     (a -> a -> a -> a -> V4 a
forall a. a -> a -> a -> a -> V4 a
V4 a
d a
e a
f a
0)
     (a -> a -> a -> a -> V4 a
forall a. a -> a -> a -> a -> V4 a
V4 a
g a
h a
i a
0)
     (a -> a -> a -> a -> V4 a
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 :: M33 a -> M44 a
m33_to_m44 (V3 V3 a
r1 V3 a
r2 V3 a
r3) = V4 a -> V4 a -> V4 a -> V4 a -> M44 a
forall a. a -> a -> a -> a -> V4 a
V4 (V3 a -> V4 a
forall a. Num a => V3 a -> V4 a
vector V3 a
r1) (V3 a -> V4 a
forall a. Num a => V3 a -> V4 a
vector V3 a
r2) (V3 a -> V4 a
forall a. Num a => V3 a -> V4 a
vector V3 a
r3) (V3 a -> V4 a
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 :: t (t a)
identity = t a -> t (t a)
forall (t :: * -> *) a. (Traversable t, Num a) => t a -> t (t a)
scaled (a -> t a
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 :: Lens' (t (v a)) (V3 a)
translation = LensLike (Context a a) (v a) (v a) a a
-> Lens (t (v a)) (t (v a)) (t a) (t a)
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 a) (v a) (v a) a a
forall (t :: * -> *) a. R4 t => Lens' (t a) a
_w((t a -> f (t a)) -> t (v a) -> f (t (v a)))
-> ((V3 a -> f (V3 a)) -> t a -> f (t a))
-> (V3 a -> f (V3 a))
-> t (v a)
-> f (t (v a))
forall b c a. (b -> c) -> (a -> b) -> a -> c
.(V3 a -> f (V3 a)) -> t a -> f (t a)
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 :: Lens' (t (v a)) (M22 a)
_m22 = LensLike (Context (V2 a) (V2 a)) (v a) (v a) (V2 a) (V2 a)
-> Lens (t (v a)) (t (v a)) (t (V2 a)) (t (V2 a))
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 (V2 a) (V2 a)) (v a) (v a) (V2 a) (V2 a)
forall (t :: * -> *) a. R2 t => Lens' (t a) (V2 a)
_xy((t (V2 a) -> f (t (V2 a))) -> t (v a) -> f (t (v a)))
-> ((M22 a -> f (M22 a)) -> t (V2 a) -> f (t (V2 a)))
-> (M22 a -> f (M22 a))
-> t (v a)
-> f (t (v a))
forall b c a. (b -> c) -> (a -> b) -> a -> c
.(M22 a -> f (M22 a)) -> t (V2 a) -> f (t (V2 a))
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 :: Lens' (t (v a)) (M23 a)
_m23 = LensLike (Context (V3 a) (V3 a)) (v a) (v a) (V3 a) (V3 a)
-> Lens (t (v a)) (t (v a)) (t (V3 a)) (t (V3 a))
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 (V3 a) (V3 a)) (v a) (v a) (V3 a) (V3 a)
forall (t :: * -> *) a. R3 t => Lens' (t a) (V3 a)
_xyz((t (V3 a) -> f (t (V3 a))) -> t (v a) -> f (t (v a)))
-> ((M23 a -> f (M23 a)) -> t (V3 a) -> f (t (V3 a)))
-> (M23 a -> f (M23 a))
-> t (v a)
-> f (t (v a))
forall b c a. (b -> c) -> (a -> b) -> a -> c
.(M23 a -> f (M23 a)) -> t (V3 a) -> f (t (V3 a))
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 :: Lens' (t (v a)) (M24 a)
_m24 = LensLike (Context (V4 a) (V4 a)) (v a) (v a) (V4 a) (V4 a)
-> Lens (t (v a)) (t (v a)) (t (V4 a)) (t (V4 a))
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 (V4 a) (V4 a)) (v a) (v a) (V4 a) (V4 a)
forall (t :: * -> *) a. R4 t => Lens' (t a) (V4 a)
_xyzw((t (V4 a) -> f (t (V4 a))) -> t (v a) -> f (t (v a)))
-> ((M24 a -> f (M24 a)) -> t (V4 a) -> f (t (V4 a)))
-> (M24 a -> f (M24 a))
-> t (v a)
-> f (t (v a))
forall b c a. (b -> c) -> (a -> b) -> a -> c
.(M24 a -> f (M24 a)) -> t (V4 a) -> f (t (V4 a))
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 :: Lens' (t (v a)) (M32 a)
_m32 = LensLike (Context (V2 a) (V2 a)) (v a) (v a) (V2 a) (V2 a)
-> Lens (t (v a)) (t (v a)) (t (V2 a)) (t (V2 a))
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 (V2 a) (V2 a)) (v a) (v a) (V2 a) (V2 a)
forall (t :: * -> *) a. R2 t => Lens' (t a) (V2 a)
_xy((t (V2 a) -> f (t (V2 a))) -> t (v a) -> f (t (v a)))
-> ((M32 a -> f (M32 a)) -> t (V2 a) -> f (t (V2 a)))
-> (M32 a -> f (M32 a))
-> t (v a)
-> f (t (v a))
forall b c a. (b -> c) -> (a -> b) -> a -> c
.(M32 a -> f (M32 a)) -> t (V2 a) -> f (t (V2 a))
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 :: Lens' (t (v a)) (M33 a)
_m33 = LensLike (Context (V3 a) (V3 a)) (v a) (v a) (V3 a) (V3 a)
-> Lens (t (v a)) (t (v a)) (t (V3 a)) (t (V3 a))
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 (V3 a) (V3 a)) (v a) (v a) (V3 a) (V3 a)
forall (t :: * -> *) a. R3 t => Lens' (t a) (V3 a)
_xyz((t (V3 a) -> f (t (V3 a))) -> t (v a) -> f (t (v a)))
-> ((M33 a -> f (M33 a)) -> t (V3 a) -> f (t (V3 a)))
-> (M33 a -> f (M33 a))
-> t (v a)
-> f (t (v a))
forall b c a. (b -> c) -> (a -> b) -> a -> c
.(M33 a -> f (M33 a)) -> t (V3 a) -> f (t (V3 a))
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 :: Lens' (t (v a)) (M34 a)
_m34 = LensLike (Context (V4 a) (V4 a)) (v a) (v a) (V4 a) (V4 a)
-> Lens (t (v a)) (t (v a)) (t (V4 a)) (t (V4 a))
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 (V4 a) (V4 a)) (v a) (v a) (V4 a) (V4 a)
forall (t :: * -> *) a. R4 t => Lens' (t a) (V4 a)
_xyzw((t (V4 a) -> f (t (V4 a))) -> t (v a) -> f (t (v a)))
-> ((M34 a -> f (M34 a)) -> t (V4 a) -> f (t (V4 a)))
-> (M34 a -> f (M34 a))
-> t (v a)
-> f (t (v a))
forall b c a. (b -> c) -> (a -> b) -> a -> c
.(M34 a -> f (M34 a)) -> t (V4 a) -> f (t (V4 a))
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 :: Lens' (t (v a)) (M42 a)
_m42 = LensLike (Context (V2 a) (V2 a)) (v a) (v a) (V2 a) (V2 a)
-> Lens (t (v a)) (t (v a)) (t (V2 a)) (t (V2 a))
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 (V2 a) (V2 a)) (v a) (v a) (V2 a) (V2 a)
forall (t :: * -> *) a. R2 t => Lens' (t a) (V2 a)
_xy((t (V2 a) -> f (t (V2 a))) -> t (v a) -> f (t (v a)))
-> ((M42 a -> f (M42 a)) -> t (V2 a) -> f (t (V2 a)))
-> (M42 a -> f (M42 a))
-> t (v a)
-> f (t (v a))
forall b c a. (b -> c) -> (a -> b) -> a -> c
.(M42 a -> f (M42 a)) -> t (V2 a) -> f (t (V2 a))
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 :: Lens' (t (v a)) (M43 a)
_m43 = LensLike (Context (V3 a) (V3 a)) (v a) (v a) (V3 a) (V3 a)
-> Lens (t (v a)) (t (v a)) (t (V3 a)) (t (V3 a))
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 (V3 a) (V3 a)) (v a) (v a) (V3 a) (V3 a)
forall (t :: * -> *) a. R3 t => Lens' (t a) (V3 a)
_xyz((t (V3 a) -> f (t (V3 a))) -> t (v a) -> f (t (v a)))
-> ((M43 a -> f (M43 a)) -> t (V3 a) -> f (t (V3 a)))
-> (M43 a -> f (M43 a))
-> t (v a)
-> f (t (v a))
forall b c a. (b -> c) -> (a -> b) -> a -> c
.(M43 a -> f (M43 a)) -> t (V3 a) -> f (t (V3 a))
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 :: Lens' (t (v a)) (M44 a)
_m44 = LensLike (Context (V4 a) (V4 a)) (v a) (v a) (V4 a) (V4 a)
-> Lens (t (v a)) (t (v a)) (t (V4 a)) (t (V4 a))
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 (V4 a) (V4 a)) (v a) (v a) (V4 a) (V4 a)
forall (t :: * -> *) a. R4 t => Lens' (t a) (V4 a)
_xyzw((t (V4 a) -> f (t (V4 a))) -> t (v a) -> f (t (v a)))
-> ((M44 a -> f (M44 a)) -> t (V4 a) -> f (t (V4 a)))
-> (M44 a -> f (M44 a))
-> t (v a)
-> f (t (v a))
forall b c a. (b -> c) -> (a -> b) -> a -> c
.(M44 a -> f (M44 a)) -> t (V4 a) -> f (t (V4 a))
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 :: M22 a -> a
det22 (V2 (V2 a
a a
b) (V2 a
c a
d)) = a
a a -> a -> a
forall a. Num a => a -> a -> a
* a
d a -> a -> a
forall a. Num a => a -> a -> a
- a
b a -> a -> a
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 :: 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 a -> a -> a
forall a. Num a => a -> a -> a
* (a
ea -> a -> a
forall a. Num a => a -> a -> a
*a
ia -> a -> a
forall a. Num a => a -> a -> a
-a
fa -> a -> a
forall a. Num a => a -> a -> a
*a
h) a -> a -> a
forall a. Num a => a -> a -> a
- a
d a -> a -> a
forall a. Num a => a -> a -> a
* (a
ba -> a -> a
forall a. Num a => a -> a -> a
*a
ia -> a -> a
forall a. Num a => a -> a -> a
-a
ca -> a -> a
forall a. Num a => a -> a -> a
*a
h) a -> a -> a
forall a. Num a => a -> a -> a
+ a
g a -> a -> a
forall a. Num a => a -> a -> a
* (a
ba -> a -> a
forall a. Num a => a -> a -> a
*a
fa -> a -> a
forall a. Num a => a -> a -> a
-a
ca -> a -> a
forall a. Num a => a -> a -> a
*a
e)
{-# INLINE det33 #-}

-- |4x4 matrix determinant.
det44 :: Num a => M44 a -> a
det44 :: 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 a -> a -> a
forall a. Num a => a -> a -> a
* a
i11 a -> a -> a
forall a. Num a => a -> a -> a
- a
i10 a -> a -> a
forall a. Num a => a -> a -> a
* a
i01
    s1 :: a
s1 = a
i00 a -> a -> a
forall a. Num a => a -> a -> a
* a
i12 a -> a -> a
forall a. Num a => a -> a -> a
- a
i10 a -> a -> a
forall a. Num a => a -> a -> a
* a
i02
    s2 :: a
s2 = a
i00 a -> a -> a
forall a. Num a => a -> a -> a
* a
i13 a -> a -> a
forall a. Num a => a -> a -> a
- a
i10 a -> a -> a
forall a. Num a => a -> a -> a
* a
i03
    s3 :: a
s3 = a
i01 a -> a -> a
forall a. Num a => a -> a -> a
* a
i12 a -> a -> a
forall a. Num a => a -> a -> a
- a
i11 a -> a -> a
forall a. Num a => a -> a -> a
* a
i02
    s4 :: a
s4 = a
i01 a -> a -> a
forall a. Num a => a -> a -> a
* a
i13 a -> a -> a
forall a. Num a => a -> a -> a
- a
i11 a -> a -> a
forall a. Num a => a -> a -> a
* a
i03
    s5 :: a
s5 = a
i02 a -> a -> a
forall a. Num a => a -> a -> a
* a
i13 a -> a -> a
forall a. Num a => a -> a -> a
- a
i12 a -> a -> a
forall a. Num a => a -> a -> a
* a
i03

    c5 :: a
c5 = a
i22 a -> a -> a
forall a. Num a => a -> a -> a
* a
i33 a -> a -> a
forall a. Num a => a -> a -> a
- a
i32 a -> a -> a
forall a. Num a => a -> a -> a
* a
i23
    c4 :: a
c4 = a
i21 a -> a -> a
forall a. Num a => a -> a -> a
* a
i33 a -> a -> a
forall a. Num a => a -> a -> a
- a
i31 a -> a -> a
forall a. Num a => a -> a -> a
* a
i23
    c3 :: a
c3 = a
i21 a -> a -> a
forall a. Num a => a -> a -> a
* a
i32 a -> a -> a
forall a. Num a => a -> a -> a
- a
i31 a -> a -> a
forall a. Num a => a -> a -> a
* a
i22
    c2 :: a
c2 = a
i20 a -> a -> a
forall a. Num a => a -> a -> a
* a
i33 a -> a -> a
forall a. Num a => a -> a -> a
- a
i30 a -> a -> a
forall a. Num a => a -> a -> a
* a
i23
    c1 :: a
c1 = a
i20 a -> a -> a
forall a. Num a => a -> a -> a
* a
i32 a -> a -> a
forall a. Num a => a -> a -> a
- a
i30 a -> a -> a
forall a. Num a => a -> a -> a
* a
i22
    c0 :: a
c0 = a
i20 a -> a -> a
forall a. Num a => a -> a -> a
* a
i31 a -> a -> a
forall a. Num a => a -> a -> a
- a
i30 a -> a -> a
forall a. Num a => a -> a -> a
* a
i21
  in a
s0 a -> a -> a
forall a. Num a => a -> a -> a
* a
c5 a -> a -> a
forall a. Num a => a -> a -> a
- a
s1 a -> a -> a
forall a. Num a => a -> a -> a
* a
c4 a -> a -> a
forall a. Num a => a -> a -> a
+ a
s2 a -> a -> a
forall a. Num a => a -> a -> a
* a
c3 a -> a -> a
forall a. Num a => a -> a -> a
+ a
s3 a -> a -> a
forall a. Num a => a -> a -> a
* a
c2 a -> a -> a
forall a. Num a => a -> a -> a
- a
s4 a -> a -> a
forall a. Num a => a -> a -> a
* a
c1 a -> a -> a
forall a. Num a => a -> a -> a
+ a
s5 a -> a -> a
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 :: M22 a -> M22 a
inv22 m :: M22 a
m@(V2 (V2 a
a a
b) (V2 a
c a
d)) = (a
1 a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
det) a -> M22 a -> M22 a
forall (m :: * -> *) (r :: * -> *) a.
(Functor m, Functor r, Num a) =>
a -> m (r a) -> m (r a)
*!! V2 a -> V2 a -> M22 a
forall a. a -> a -> V2 a
V2 (a -> a -> V2 a
forall a. a -> a -> V2 a
V2 a
d (-a
b)) (a -> a -> V2 a
forall a. a -> a -> V2 a
V2 (-a
c) a
a)
  where det :: a
det = M22 a -> a
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 :: 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 a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
det) a -> M33 a -> M33 a
forall (m :: * -> *) (r :: * -> *) a.
(Functor m, Functor r, Num a) =>
a -> m (r a) -> m (r a)
*!! V3 a -> V3 a -> V3 a -> M33 a
forall a. a -> a -> a -> V3 a
V3 (a -> a -> a -> V3 a
forall a. a -> a -> a -> V3 a
V3 a
a' a
b' a
c')
                     (a -> a -> a -> V3 a
forall a. a -> a -> a -> V3 a
V3 a
d' a
e' a
f')
                     (a -> a -> a -> V3 a
forall a. a -> a -> a -> V3 a
V3 a
g' a
h' a
i')
  where a' :: a
a' = (a, a, a, a) -> a
forall a. Num a => (a, a, a, a) -> a
cofactor (a
e,a
f,a
h,a
i)
        b' :: a
b' = (a, a, a, a) -> a
forall a. Num a => (a, a, a, a) -> a
cofactor (a
c,a
b,a
i,a
h)
        c' :: a
c' = (a, a, a, a) -> a
forall a. Num a => (a, a, a, a) -> a
cofactor (a
b,a
c,a
e,a
f)
        d' :: a
d' = (a, a, a, a) -> a
forall a. Num a => (a, a, a, a) -> a
cofactor (a
f,a
d,a
i,a
g)
        e' :: a
e' = (a, a, a, a) -> a
forall a. Num a => (a, a, a, a) -> a
cofactor (a
a,a
c,a
g,a
i)
        f' :: a
f' = (a, a, a, a) -> a
forall a. Num a => (a, a, a, a) -> a
cofactor (a
c,a
a,a
f,a
d)
        g' :: a
g' = (a, a, a, a) -> a
forall a. Num a => (a, a, a, a) -> a
cofactor (a
d,a
e,a
g,a
h)
        h' :: a
h' = (a, a, a, a) -> a
forall a. Num a => (a, a, a, a) -> a
cofactor (a
b,a
a,a
h,a
g)
        i' :: a
i' = (a, a, a, a) -> a
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) = M22 a -> a
forall a. Num a => M22 a -> a
det22 (V2 a -> V2 a -> M22 a
forall a. a -> a -> V2 a
V2 (a -> a -> V2 a
forall a. a -> a -> V2 a
V2 a
q a
r) (a -> a -> V2 a
forall a. a -> a -> V2 a
V2 a
s a
t))
        det :: a
det = M33 a -> a
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 :: f (g a) -> g (f a)
transpose = f (g a) -> g (f a)
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 :: 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 a -> a -> a
forall a. Num a => a -> a -> a
* a
i11 a -> a -> a
forall a. Num a => a -> a -> a
- a
i10 a -> a -> a
forall a. Num a => a -> a -> a
* a
i01
      s1 :: a
s1 = a
i00 a -> a -> a
forall a. Num a => a -> a -> a
* a
i12 a -> a -> a
forall a. Num a => a -> a -> a
- a
i10 a -> a -> a
forall a. Num a => a -> a -> a
* a
i02
      s2 :: a
s2 = a
i00 a -> a -> a
forall a. Num a => a -> a -> a
* a
i13 a -> a -> a
forall a. Num a => a -> a -> a
- a
i10 a -> a -> a
forall a. Num a => a -> a -> a
* a
i03
      s3 :: a
s3 = a
i01 a -> a -> a
forall a. Num a => a -> a -> a
* a
i12 a -> a -> a
forall a. Num a => a -> a -> a
- a
i11 a -> a -> a
forall a. Num a => a -> a -> a
* a
i02
      s4 :: a
s4 = a
i01 a -> a -> a
forall a. Num a => a -> a -> a
* a
i13 a -> a -> a
forall a. Num a => a -> a -> a
- a
i11 a -> a -> a
forall a. Num a => a -> a -> a
* a
i03
      s5 :: a
s5 = a
i02 a -> a -> a
forall a. Num a => a -> a -> a
* a
i13 a -> a -> a
forall a. Num a => a -> a -> a
- a
i12 a -> a -> a
forall a. Num a => a -> a -> a
* a
i03
      c5 :: a
c5 = a
i22 a -> a -> a
forall a. Num a => a -> a -> a
* a
i33 a -> a -> a
forall a. Num a => a -> a -> a
- a
i32 a -> a -> a
forall a. Num a => a -> a -> a
* a
i23
      c4 :: a
c4 = a
i21 a -> a -> a
forall a. Num a => a -> a -> a
* a
i33 a -> a -> a
forall a. Num a => a -> a -> a
- a
i31 a -> a -> a
forall a. Num a => a -> a -> a
* a
i23
      c3 :: a
c3 = a
i21 a -> a -> a
forall a. Num a => a -> a -> a
* a
i32 a -> a -> a
forall a. Num a => a -> a -> a
- a
i31 a -> a -> a
forall a. Num a => a -> a -> a
* a
i22
      c2 :: a
c2 = a
i20 a -> a -> a
forall a. Num a => a -> a -> a
* a
i33 a -> a -> a
forall a. Num a => a -> a -> a
- a
i30 a -> a -> a
forall a. Num a => a -> a -> a
* a
i23
      c1 :: a
c1 = a
i20 a -> a -> a
forall a. Num a => a -> a -> a
* a
i32 a -> a -> a
forall a. Num a => a -> a -> a
- a
i30 a -> a -> a
forall a. Num a => a -> a -> a
* a
i22
      c0 :: a
c0 = a
i20 a -> a -> a
forall a. Num a => a -> a -> a
* a
i31 a -> a -> a
forall a. Num a => a -> a -> a
- a
i30 a -> a -> a
forall a. Num a => a -> a -> a
* a
i21
      det :: a
det = a
s0 a -> a -> a
forall a. Num a => a -> a -> a
* a
c5 a -> a -> a
forall a. Num a => a -> a -> a
- a
s1 a -> a -> a
forall a. Num a => a -> a -> a
* a
c4 a -> a -> a
forall a. Num a => a -> a -> a
+ a
s2 a -> a -> a
forall a. Num a => a -> a -> a
* a
c3 a -> a -> a
forall a. Num a => a -> a -> a
+ a
s3 a -> a -> a
forall a. Num a => a -> a -> a
* a
c2 a -> a -> a
forall a. Num a => a -> a -> a
- a
s4 a -> a -> a
forall a. Num a => a -> a -> a
* a
c1 a -> a -> a
forall a. Num a => a -> a -> a
+ a
s5 a -> a -> a
forall a. Num a => a -> a -> a
* a
c0
      invDet :: a
invDet = a -> a
forall a. Fractional a => a -> a
recip a
det
  in a
invDet a -> M44 a -> M44 a
forall (m :: * -> *) (r :: * -> *) a.
(Functor m, Functor r, Num a) =>
a -> m (r a) -> m (r a)
*!! V4 a -> V4 a -> V4 a -> V4 a -> M44 a
forall a. a -> a -> a -> a -> V4 a
V4 (a -> a -> a -> a -> V4 a
forall a. a -> a -> a -> a -> V4 a
V4 (a
i11 a -> a -> a
forall a. Num a => a -> a -> a
* a
c5 a -> a -> a
forall a. Num a => a -> a -> a
- a
i12 a -> a -> a
forall a. Num a => a -> a -> a
* a
c4 a -> a -> a
forall a. Num a => a -> a -> a
+ a
i13 a -> a -> a
forall a. Num a => a -> a -> a
* a
c3)
                       (-a
i01 a -> a -> a
forall a. Num a => a -> a -> a
* a
c5 a -> a -> a
forall a. Num a => a -> a -> a
+ a
i02 a -> a -> a
forall a. Num a => a -> a -> a
* a
c4 a -> a -> a
forall a. Num a => a -> a -> a
- a
i03 a -> a -> a
forall a. Num a => a -> a -> a
* a
c3)
                       (a
i31 a -> a -> a
forall a. Num a => a -> a -> a
* a
s5 a -> a -> a
forall a. Num a => a -> a -> a
- a
i32 a -> a -> a
forall a. Num a => a -> a -> a
* a
s4 a -> a -> a
forall a. Num a => a -> a -> a
+ a
i33 a -> a -> a
forall a. Num a => a -> a -> a
* a
s3)
                       (-a
i21 a -> a -> a
forall a. Num a => a -> a -> a
* a
s5 a -> a -> a
forall a. Num a => a -> a -> a
+ a
i22 a -> a -> a
forall a. Num a => a -> a -> a
* a
s4 a -> a -> a
forall a. Num a => a -> a -> a
- a
i23 a -> a -> a
forall a. Num a => a -> a -> a
* a
s3))
                   (a -> a -> a -> a -> V4 a
forall a. a -> a -> a -> a -> V4 a
V4 (-a
i10 a -> a -> a
forall a. Num a => a -> a -> a
* a
c5 a -> a -> a
forall a. Num a => a -> a -> a
+ a
i12 a -> a -> a
forall a. Num a => a -> a -> a
* a
c2 a -> a -> a
forall a. Num a => a -> a -> a
- a
i13 a -> a -> a
forall a. Num a => a -> a -> a
* a
c1)
                       (a
i00 a -> a -> a
forall a. Num a => a -> a -> a
* a
c5 a -> a -> a
forall a. Num a => a -> a -> a
- a
i02 a -> a -> a
forall a. Num a => a -> a -> a
* a
c2 a -> a -> a
forall a. Num a => a -> a -> a
+ a
i03 a -> a -> a
forall a. Num a => a -> a -> a
* a
c1)
                       (-a
i30 a -> a -> a
forall a. Num a => a -> a -> a
* a
s5 a -> a -> a
forall a. Num a => a -> a -> a
+ a
i32 a -> a -> a
forall a. Num a => a -> a -> a
* a
s2 a -> a -> a
forall a. Num a => a -> a -> a
- a
i33 a -> a -> a
forall a. Num a => a -> a -> a
* a
s1)
                       (a
i20 a -> a -> a
forall a. Num a => a -> a -> a
* a
s5 a -> a -> a
forall a. Num a => a -> a -> a
- a
i22 a -> a -> a
forall a. Num a => a -> a -> a
* a
s2 a -> a -> a
forall a. Num a => a -> a -> a
+ a
i23 a -> a -> a
forall a. Num a => a -> a -> a
* a
s1))
                   (a -> a -> a -> a -> V4 a
forall a. a -> a -> a -> a -> V4 a
V4 (a
i10 a -> a -> a
forall a. Num a => a -> a -> a
* a
c4 a -> a -> a
forall a. Num a => a -> a -> a
- a
i11 a -> a -> a
forall a. Num a => a -> a -> a
* a
c2 a -> a -> a
forall a. Num a => a -> a -> a
+ a
i13 a -> a -> a
forall a. Num a => a -> a -> a
* a
c0)
                       (-a
i00 a -> a -> a
forall a. Num a => a -> a -> a
* a
c4 a -> a -> a
forall a. Num a => a -> a -> a
+ a
i01 a -> a -> a
forall a. Num a => a -> a -> a
* a
c2 a -> a -> a
forall a. Num a => a -> a -> a
- a
i03 a -> a -> a
forall a. Num a => a -> a -> a
* a
c0)
                       (a
i30 a -> a -> a
forall a. Num a => a -> a -> a
* a
s4 a -> a -> a
forall a. Num a => a -> a -> a
- a
i31 a -> a -> a
forall a. Num a => a -> a -> a
* a
s2 a -> a -> a
forall a. Num a => a -> a -> a
+ a
i33 a -> a -> a
forall a. Num a => a -> a -> a
* a
s0)
                       (-a
i20 a -> a -> a
forall a. Num a => a -> a -> a
* a
s4 a -> a -> a
forall a. Num a => a -> a -> a
+ a
i21 a -> a -> a
forall a. Num a => a -> a -> a
* a
s2 a -> a -> a
forall a. Num a => a -> a -> a
- a
i23 a -> a -> a
forall a. Num a => a -> a -> a
* a
s0))
                   (a -> a -> a -> a -> V4 a
forall a. a -> a -> a -> a -> V4 a
V4 (-a
i10 a -> a -> a
forall a. Num a => a -> a -> a
* a
c3 a -> a -> a
forall a. Num a => a -> a -> a
+ a
i11 a -> a -> a
forall a. Num a => a -> a -> a
* a
c1 a -> a -> a
forall a. Num a => a -> a -> a
- a
i12 a -> a -> a
forall a. Num a => a -> a -> a
* a
c0)
                       (a
i00 a -> a -> a
forall a. Num a => a -> a -> a
* a
c3 a -> a -> a
forall a. Num a => a -> a -> a
- a
i01 a -> a -> a
forall a. Num a => a -> a -> a
* a
c1 a -> a -> a
forall a. Num a => a -> a -> a
+ a
i02 a -> a -> a
forall a. Num a => a -> a -> a
* a
c0)
                       (-a
i30 a -> a -> a
forall a. Num a => a -> a -> a
* a
s3 a -> a -> a
forall a. Num a => a -> a -> a
+ a
i31 a -> a -> a
forall a. Num a => a -> a -> a
* a
s1 a -> a -> a
forall a. Num a => a -> a -> a
- a
i32 a -> a -> a
forall a. Num a => a -> a -> a
* a
s0)
                       (a
i20 a -> a -> a
forall a. Num a => a -> a -> a
* a
s3 a -> a -> a
forall a. Num a => a -> a -> a
- a
i21 a -> a -> a
forall a. Num a => a -> a -> a
* a
s1 a -> a -> a
forall a. Num a => a -> a -> a
+ a
i22 a -> a -> a
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 :: m (m a) -> (m (m a), m (m a))
lu m (m a)
a =
    let n :: i
n = Int -> i
forall a b. (Integral a, Num b) => a -> b
fromIntegral (m (m a) -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length m (m a)
a)
        initU :: m (m a)
initU = m (m a)
forall a (t :: * -> *).
(Num a, Traversable t, Applicative t) =>
t (t a)
identity
        initL :: m (m a)
initL = m (m a)
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 i -> i -> Bool
forall a. Eq a => a -> a -> Bool
== i
j = a
s
                    | Bool
otherwise = i -> a -> a
go (i
ki -> i -> i
forall a. Num a => a -> a -> a
+i
1)
                                     ( a
s
                                      a -> a -> a
forall a. Num a => a -> a -> a
+ ( (m (m a)
l m (m a) -> Getting (Endo (m a)) (m (m a)) (m a) -> m a
forall s a. HasCallStack => s -> Getting (Endo a) s a -> a
^?! Index (m (m a)) -> Traversal' (m (m a)) (IxValue (m (m a)))
forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
Index (m (m a))
i m a -> Getting (Endo a) (m a) a -> a
forall s a. HasCallStack => s -> Getting (Endo a) s a -> a
^?! Index (m a) -> Traversal' (m a) (IxValue (m a))
forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
Index (m a)
k)
                                        a -> a -> a
forall a. Num a => a -> a -> a
* (m (m a)
u m (m a) -> Getting (Endo (m a)) (m (m a)) (m a) -> m a
forall s a. HasCallStack => s -> Getting (Endo a) s a -> a
^?! Index (m (m a)) -> Traversal' (m (m a)) (IxValue (m (m a)))
forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
Index (m (m a))
k m a -> Getting (Endo a) (m a) a -> a
forall s a. HasCallStack => s -> Getting (Endo a) s a -> a
^?! Index (m a) -> Traversal' (m a) (IxValue (m a))
forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
Index (m a)
j)
                                        )
                                      )
                s' :: a
s' = i -> a -> a
go i
0 a
0
            in m (m a)
l m (m a) -> (m (m a) -> m (m a)) -> m (m a)
forall a b. a -> (a -> b) -> b
& (Index (m (m a)) -> Traversal' (m (m a)) (IxValue (m (m a)))
forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
Index (m (m a))
i ((m a -> Identity (m a)) -> m (m a) -> Identity (m (m a)))
-> ((a -> Identity a) -> m a -> Identity (m a))
-> (a -> Identity a)
-> m (m a)
-> Identity (m (m a))
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Index (m a) -> Traversal' (m a) (IxValue (m a))
forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
Index (m a)
j) ((a -> Identity a) -> m (m a) -> Identity (m (m a)))
-> a -> m (m a) -> m (m a)
forall s t a b. ASetter s t a b -> b -> s -> t
.~ ((m (m a)
a m (m a) -> Getting (Endo (m a)) (m (m a)) (m a) -> m a
forall s a. HasCallStack => s -> Getting (Endo a) s a -> a
^?! Index (m (m a)) -> Traversal' (m (m a)) (IxValue (m (m a)))
forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
Index (m (m a))
i m a -> Getting (Endo a) (m a) a -> a
forall s a. HasCallStack => s -> Getting (Endo a) s a -> a
^?! Index (m a) -> Traversal' (m a) (IxValue (m a))
forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
Index (m a)
j) a -> a -> a
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 i -> i -> Bool
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
ii -> i -> i
forall 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 i -> i -> Bool
forall a. Eq a => a -> a -> Bool
== i
j = a
s
                    | Bool
otherwise = i -> a -> a
go (i
ki -> i -> i
forall a. Num a => a -> a -> a
+i
1)
                                     ( a
s
                                     a -> a -> a
forall a. Num a => a -> a -> a
+ ( (m (m a)
l m (m a) -> Getting (Endo (m a)) (m (m a)) (m a) -> m a
forall s a. HasCallStack => s -> Getting (Endo a) s a -> a
^?! Index (m (m a)) -> Traversal' (m (m a)) (IxValue (m (m a)))
forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
Index (m (m a))
j m a -> Getting (Endo a) (m a) a -> a
forall s a. HasCallStack => s -> Getting (Endo a) s a -> a
^?! Index (m a) -> Traversal' (m a) (IxValue (m a))
forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
Index (m a)
k)
                                       a -> a -> a
forall a. Num a => a -> a -> a
* (m (m a)
u m (m a) -> Getting (Endo (m a)) (m (m a)) (m a) -> m a
forall s a. HasCallStack => s -> Getting (Endo a) s a -> a
^?! Index (m (m a)) -> Traversal' (m (m a)) (IxValue (m (m a)))
forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
Index (m (m a))
k m a -> Getting (Endo a) (m a) a -> a
forall s a. HasCallStack => s -> Getting (Endo a) s a -> a
^?! Index (m a) -> Traversal' (m a) (IxValue (m a))
forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
Index (m a)
i)
                                       )
                                     )
                s' :: a
s' = i -> a -> a
go i
0 a
0
            in m (m a)
u m (m a) -> (m (m a) -> m (m a)) -> m (m a)
forall a b. a -> (a -> b) -> b
& (Index (m (m a)) -> Traversal' (m (m a)) (IxValue (m (m a)))
forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
Index (m (m a))
j ((m a -> Identity (m a)) -> m (m a) -> Identity (m (m a)))
-> ((a -> Identity a) -> m a -> Identity (m a))
-> (a -> Identity a)
-> m (m a)
-> Identity (m (m a))
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Index (m a) -> Traversal' (m a) (IxValue (m a))
forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
Index (m a)
i) ((a -> Identity a) -> m (m a) -> Identity (m (m a)))
-> a -> m (m a) -> m (m a)
forall s t a b. ASetter s t a b -> b -> s -> t
.~ ( ((m (m a)
a m (m a) -> Getting (Endo (m a)) (m (m a)) (m a) -> m a
forall s a. HasCallStack => s -> Getting (Endo a) s a -> a
^?! Index (m (m a)) -> Traversal' (m (m a)) (IxValue (m (m a)))
forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
Index (m (m a))
j m a -> Getting (Endo a) (m a) a -> a
forall s a. HasCallStack => s -> Getting (Endo a) s a -> a
^?! Index (m a) -> Traversal' (m a) (IxValue (m a))
forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
Index (m a)
i) a -> a -> a
forall a. Num a => a -> a -> a
- a
s')
                                    a -> a -> a
forall a. Fractional a => a -> a -> a
/ (m (m a)
l m (m a) -> Getting (Endo (m a)) (m (m a)) (m a) -> m a
forall s a. HasCallStack => s -> Getting (Endo a) s a -> a
^?! Index (m (m a)) -> Traversal' (m (m a)) (IxValue (m (m a)))
forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
Index (m (m a))
j m a -> Getting (Endo a) (m a) a -> a
forall s a. HasCallStack => s -> Getting (Endo a) s a -> a
^?! Index (m a) -> Traversal' (m a) (IxValue (m a))
forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
Index (m a)
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 i -> i -> Bool
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
ii -> i -> i
forall 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 i -> i -> Bool
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
ji -> i -> i
forall 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 :: m (m a) -> (m (m a), m (m a))
luFinite m (m a)
a =
    (V n (V n a) -> m (m a))
-> (V n (V n a) -> m (m a))
-> (V n (V n a), V n (V n a))
-> (m (m a), m (m a))
forall (p :: * -> * -> *) a b c d.
Bifunctor p =>
(a -> b) -> (c -> d) -> p a c -> p b d
bimap ((V n a -> m a) -> m (V n a) -> m (m a)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap V n a -> m a
forall (v :: * -> *) a. Finite v => V (Size v) a -> v a
fromV (m (V n a) -> m (m a))
-> (V n (V n a) -> m (V n a)) -> V n (V n a) -> m (m a)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. V n (V n a) -> m (V n a)
forall (v :: * -> *) a. Finite v => V (Size v) a -> v a
fromV)
          ((V n a -> m a) -> m (V n a) -> m (m a)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap V n a -> m a
forall (v :: * -> *) a. Finite v => V (Size v) a -> v a
fromV (m (V n a) -> m (m a))
-> (V n (V n a) -> m (V n a)) -> V n (V n a) -> m (m a)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. V n (V n a) -> m (V n a)
forall (v :: * -> *) a. Finite v => V (Size v) a -> v a
fromV)
          (V n (V n a) -> (V n (V n a), V n (V n a))
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 a -> V n a) -> V n (m a) -> V n (V n a)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap m a -> V n a
forall (v :: * -> *) a. Finite v => v a -> V (Size v) a
toV (m (m a) -> V (Size m) (m a)
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 :: m (m a) -> m a -> m a
forwardSub m (m a)
a m a
b =
    let n :: i
n = Int -> i
forall a b. (Integral a, Num b) => a -> b
fromIntegral (m a -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length m a
b)
        initX :: m a
initX = m a
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 i -> i -> Bool
forall a. Eq a => a -> a -> Bool
== i
i = a
s
            | Bool
otherwise = i -> i -> a -> m a -> a
coeff i
i (i
ji -> i -> i
forall a. Num a => a -> a -> a
+i
1) (a
s a -> a -> a
forall a. Num a => a -> a -> a
+ ((m (m a)
a m (m a) -> Getting (Endo (m a)) (m (m a)) (m a) -> m a
forall s a. HasCallStack => s -> Getting (Endo a) s a -> a
^?! Index (m (m a)) -> Traversal' (m (m a)) (IxValue (m (m a)))
forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
Index (m (m a))
i m a -> Getting (Endo a) (m a) a -> a
forall s a. HasCallStack => s -> Getting (Endo a) s a -> a
^?! Index (m a) -> Traversal' (m a) (IxValue (m a))
forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
Index (m a)
j) a -> a -> a
forall a. Num a => a -> a -> a
* (m a
x m a -> Getting (Endo a) (m a) a -> a
forall s a. HasCallStack => s -> Getting (Endo a) s a -> a
^?! Index (m a) -> Traversal' (m a) (IxValue (m a))
forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
Index (m a)
j))) m a
x
        go :: i -> m a -> m a
go !i
i !m a
x
            | i
i i -> i -> Bool
forall a. Eq a => a -> a -> Bool
== i
n = m a
x
            | Bool
otherwise = i -> m a -> m a
go (i
i i -> i -> i
forall a. Num a => a -> a -> a
+ i
1) (m a
x m a -> (m a -> m a) -> m a
forall a b. a -> (a -> b) -> b
& Index (m a) -> Traversal' (m a) (IxValue (m a))
forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
Index (m a)
i ((a -> Identity a) -> m a -> Identity (m a)) -> a -> m a -> m a
forall s t a b. ASetter s t a b -> b -> s -> t
.~ ( ((m a
b m a -> Getting (Endo a) (m a) a -> a
forall s a. HasCallStack => s -> Getting (Endo a) s a -> a
^?! Index (m a) -> Traversal' (m a) (IxValue (m a))
forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
Index (m a)
i) a -> a -> a
forall a. Num a => a -> a -> a
- i -> i -> a -> m a -> a
coeff i
i i
0 a
0 m a
x)
                                                  a -> a -> a
forall a. Fractional a => a -> a -> a
/ (m (m a)
a m (m a) -> Getting (Endo (m a)) (m (m a)) (m a) -> m a
forall s a. HasCallStack => s -> Getting (Endo a) s a -> a
^?! Index (m (m a)) -> Traversal' (m (m a)) (IxValue (m (m a)))
forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
Index (m (m a))
i m a -> Getting (Endo a) (m a) a -> a
forall s a. HasCallStack => s -> Getting (Endo a) s a -> a
^?! Index (m a) -> Traversal' (m a) (IxValue (m a))
forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
Index (m a)
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 :: m (m a) -> m a -> m a
forwardSubFinite m (m a)
a m a
b = V (Size m) a -> m a
forall (v :: * -> *) a. Finite v => V (Size v) a -> v a
fromV (V n (V n a) -> V n a -> V n a
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 a -> V n a) -> V n (m a) -> V n (V n a)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap m a -> V n a
forall (v :: * -> *) a. Finite v => v a -> V (Size v) a
toV (m (m a) -> V (Size m) (m a)
forall (v :: * -> *) a. Finite v => v a -> V (Size v) a
toV m (m a)
a)) (m a -> V (Size m) 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 :: m (m a) -> m a -> m a
backwardSub m (m a)
a m a
b =
    let n :: i
n = Int -> i
forall a b. (Integral a, Num b) => a -> b
fromIntegral (m a -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length m a
b)
        initX :: m a
initX = m a
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 i -> i -> Bool
forall a. Eq a => a -> a -> Bool
== i
n = a
s
            | Bool
otherwise = i -> i -> a -> m a -> a
coeff i
i
                                (i
ji -> i -> i
forall a. Num a => a -> a -> a
+i
1)
                                (a
s a -> a -> a
forall a. Num a => a -> a -> a
+ ((m (m a)
a m (m a) -> Getting (Endo (m a)) (m (m a)) (m a) -> m a
forall s a. HasCallStack => s -> Getting (Endo a) s a -> a
^?! Index (m (m a)) -> Traversal' (m (m a)) (IxValue (m (m a)))
forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
Index (m (m a))
i m a -> Getting (Endo a) (m a) a -> a
forall s a. HasCallStack => s -> Getting (Endo a) s a -> a
^?! Index (m a) -> Traversal' (m a) (IxValue (m a))
forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
Index (m a)
j) a -> a -> a
forall a. Num a => a -> a -> a
* (m a
x m a -> Getting (Endo a) (m a) a -> a
forall s a. HasCallStack => s -> Getting (Endo a) s a -> a
^?! Index (m a) -> Traversal' (m a) (IxValue (m a))
forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
Index (m a)
j)))
                                m a
x
        go :: i -> m a -> m a
go !i
i !m a
x
            | i
i i -> i -> Bool
forall a. Ord a => a -> a -> Bool
< i
0 = m a
x
            | Bool
otherwise = i -> m a -> m a
go (i
ii -> i -> i
forall a. Num a => a -> a -> a
-i
1)
                             (m a
x m a -> (m a -> m a) -> m a
forall a b. a -> (a -> b) -> b
& Index (m a) -> Traversal' (m a) (IxValue (m a))
forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
Index (m a)
i ((a -> Identity a) -> m a -> Identity (m a)) -> a -> m a -> m a
forall s t a b. ASetter s t a b -> b -> s -> t
.~ ( ((m a
b m a -> Getting (Endo a) (m a) a -> a
forall s a. HasCallStack => s -> Getting (Endo a) s a -> a
^?! Index (m a) -> Traversal' (m a) (IxValue (m a))
forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
Index (m a)
i) a -> a -> a
forall a. Num a => a -> a -> a
- i -> i -> a -> m a -> a
coeff i
i (i
ii -> i -> i
forall a. Num a => a -> a -> a
+i
1) a
0 m a
x)
                                          a -> a -> a
forall a. Fractional a => a -> a -> a
/ (m (m a)
a m (m a) -> Getting (Endo (m a)) (m (m a)) (m a) -> m a
forall s a. HasCallStack => s -> Getting (Endo a) s a -> a
^?! Index (m (m a)) -> Traversal' (m (m a)) (IxValue (m (m a)))
forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
Index (m (m a))
i m a -> Getting (Endo a) (m a) a -> a
forall s a. HasCallStack => s -> Getting (Endo a) s a -> a
^?! Index (m a) -> Traversal' (m a) (IxValue (m a))
forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
Index (m a)
i)
                                          ))
    in i -> m a -> m a
go (i
ni -> i -> i
forall 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 :: m (m a) -> m a -> m a
backwardSubFinite m (m a)
a m a
b = V (Size m) a -> m a
forall (v :: * -> *) a. Finite v => V (Size v) a -> v a
fromV (V n (V n a) -> V n a -> V n a
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 a -> V n a) -> V n (m a) -> V n (V n a)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap m a -> V n a
forall (v :: * -> *) a. Finite v => v a -> V (Size v) a
toV (m (m a) -> V (Size m) (m a)
forall (v :: * -> *) a. Finite v => v a -> V (Size v) a
toV m (m a)
a)) (m a -> V (Size m) 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 :: m (m a) -> m a -> m a
luSolve m (m a)
a m a
b =
    let (m (m a)
l, m (m a)
u) = m (m a) -> (m (m a), m (m a))
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 m (m a) -> m a -> m a
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 (m (m a) -> m a -> m a
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 :: m (m a) -> m a -> m a
luSolveFinite m (m a)
a m a
b = V (Size m) a -> m a
forall (v :: * -> *) a. Finite v => V (Size v) a -> v a
fromV (V n (V n a) -> V n a -> V n a
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 a -> V n a) -> V n (m a) -> V n (V n a)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap m a -> V n a
forall (v :: * -> *) a. Finite v => v a -> V (Size v) a
toV (m (m a) -> V (Size m) (m a)
forall (v :: * -> *) a. Finite v => v a -> V (Size v) a
toV m (m a)
a)) (m a -> V (Size m) 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 :: m (m a) -> m (m a)
luInv m (m a)
a =
    let n :: i
n = Int -> i
forall a b. (Integral a, Num b) => a -> b
fromIntegral (m (m a) -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length m (m a)
a)
        initA' :: m (m a)
initA' = m (m a)
forall (f :: * -> *) a. (Additive f, Num a) => f a
zero
        (m (m a)
l, m (m a)
u) = m (m a) -> (m (m a), m (m a))
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 i -> i -> Bool
forall a. Eq a => a -> a -> Bool
== i
n = m (m a)
a'
            | Bool
otherwise = let e :: m a
e   = m a
forall (f :: * -> *) a. (Additive f, Num a) => f a
zero m a -> (m a -> m a) -> m a
forall a b. a -> (a -> b) -> b
& Index (m a) -> Traversal' (m a) (IxValue (m a))
forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
Index (m a)
i ((a -> Identity a) -> m a -> Identity (m a)) -> a -> m a -> m a
forall s t a b. ASetter s t a b -> b -> s -> t
.~ a
1
                              a'r :: m a
a'r = m (m a) -> m a -> m a
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 (m (m a) -> m a -> m a
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
ii -> i -> i
forall a. Num a => a -> a -> a
+i
1) (m (m a)
a' m (m a) -> (m (m a) -> m (m a)) -> m (m a)
forall a b. a -> (a -> b) -> b
& Index (m (m a)) -> Traversal' (m (m a)) (IxValue (m (m a)))
forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix i
Index (m (m a))
i ((m a -> Identity (m a)) -> m (m a) -> Identity (m (m a)))
-> m a -> m (m a) -> m (m a)
forall s t a b. ASetter s t a b -> b -> s -> t
.~ m a
a'r)
    in m (m a) -> m (m a)
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 :: m (m a) -> m (m a)
luInvFinite m (m a)
a = (V n a -> m a) -> m (V n a) -> m (m a)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap V n a -> m a
forall (v :: * -> *) a. Finite v => V (Size v) a -> v a
fromV (V (Size m) (V n a) -> m (V n a)
forall (v :: * -> *) a. Finite v => V (Size v) a -> v a
fromV (V n (V n a) -> V n (V n a)
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 a -> V n a) -> V n (m a) -> V n (V n a)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap m a -> V n a
forall (v :: * -> *) a. Finite v => v a -> V (Size v) a
toV (m (m a) -> V (Size m) (m a)
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 :: m (m a) -> a
luDet m (m a)
a =
    let (m (m a)
l, m (m a)
u) = m (m a) -> (m (m a), m (m a))
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      = (a -> a -> a) -> a -> m a -> a
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
Foldable.foldl a -> a -> a
forall a. Num a => a -> a -> a
(*) a
1
    in m a -> a
p (m (m a) -> m a
forall (m :: * -> *) a. Trace m => m (m a) -> m a
diagonal m (m a)
l) a -> a -> a
forall a. Num a => a -> a -> a
* m a -> a
p (m (m a) -> m a
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 :: m (m a) -> a
luDetFinite = V n (V n a) -> a
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 (V n (V n a) -> a) -> (m (m a) -> V n (V n a)) -> m (m a) -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (m a -> V n a) -> V n (m a) -> V n (V n a)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap m a -> V n a
forall (v :: * -> *) a. Finite v => v a -> V (Size v) a
toV (V n (m a) -> V n (V n a))
-> (m (m a) -> V n (m a)) -> m (m a) -> V n (V n a)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. m (m a) -> V n (m a)
forall (v :: * -> *) a. Finite v => v a -> V (Size v) a
toV