{-# OPTIONS_GHC -fno-warn-name-shadowing #-}
{-# LANGUAGE BangPatterns, DataKinds, FlexibleContexts, FlexibleInstances  #-}
{-# LANGUAGE GADTs, GeneralizedNewtypeDeriving, KindSignatures, LambdaCase #-}
{-# LANGUAGE MultiParamTypeClasses, NamedFieldPuns, NoImplicitPrelude      #-}
{-# LANGUAGE NoMonomorphismRestriction, PolyKinds, RankNTypes              #-}
{-# LANGUAGE ScopedTypeVariables, StandaloneDeriving, TemplateHaskell      #-}
{-# LANGUAGE TupleSections, UndecidableInstances, ViewPatterns             #-}
{-# OPTIONS_GHC -funbox-strict-fields -Wno-type-defaults -Wno-redundant-constraints #-}
module Algebra.LinkedMatrix (Matrix, toLists, fromLists, fromList,
                             swapRows, identity,nonZeroRows,nonZeroCols,
                             swapCols, switchCols, switchRows, addRow,
                             addCol, ncols, nrows, getRow, getCol, triangulateModular,
                             scaleRow, combineRows, combineCols, transpose,
                             inBound, height, width, cmap, empty, rowVector,
                             colVector, rowCount, colCount, traverseRow,
                             traverseCol, Entry, idx, value, substMatrix,
                             catRow, catCol, (<||>), (<-->), toRows, toCols,
                             zeroMat, getDiag, trace, diagProd, diag,
                             scaleCol, clearRow, clearCol, index, (!),
                             nonZeroEntries, rankLM, splitIndependentDirs,
                             structuredGauss, multWithVector, solveWiedemann,
                             henselLift, solveHensel, structuredGauss', intDet) where
import Algebra.Algorithms.ChineseRemainder
import Algebra.Field.Prime
import Algebra.Instances                   ()
import Algebra.Prelude.Core                hiding (Vector, empty, fromList,
                                            generate, insert, transpose, (%),
                                            (<.>), (<>))

import           Control.Applicative          ((<|>))
import           Control.Arrow                ((&&&))
import           Control.Arrow                ((>>>))
import           Control.Lens                 hiding (index, (<.>))
import           Control.Monad                (replicateM)
import           Control.Monad.Loops          (iterateUntil)
import           Control.Monad.Random         hiding (fromList)
import           Control.Monad.ST.Strict      (ST, runST)
import           Control.Monad.State.Strict   (evalState, runState)
import           Control.Parallel.Strategies  (parMap)
import           Control.Parallel.Strategies  (rdeepseq)
import           Data.IntMap.Strict           (IntMap, alter, insert,
                                               mapMaybeWithKey, minViewWithKey)
import qualified Data.IntMap.Strict           as IM
import           Data.IntSet                  (IntSet)
import qualified Data.IntSet                  as IS
import           Data.List                    (minimumBy, sort)
import           Data.List                    (sortBy)
import           Data.Maybe                   (fromJust, fromMaybe, mapMaybe)
import           Data.Monoid                  (First (..))
import           Data.Numbers.Primes          (primes)
import           Data.Ord                     (comparing)
import           Data.Proxy                   (Proxy (..))
import           Data.Reflection              (Reifies (..), reify)
import           Data.Semigroup               hiding (First (..))
import           Data.Tuple                   (swap)
import           Data.Vector                  (Vector, create, generate, thaw,
                                               unsafeFreeze)
import qualified Data.Vector                  as V
import           Data.Vector.Mutable          (grow)
import qualified Data.Vector.Mutable          as MV
import           Numeric.Decidable.Zero       (isZero)
import           Numeric.Domain.GCD           (gcd, lcm)
import           Numeric.Field.Fraction       ((%))
import           Numeric.Semiring.ZeroProduct (ZeroProductSemiring)
import qualified Prelude                      as P

data Entry a = Entry { Entry a -> a
_value   :: !a
                     , Entry a -> (Int, Int)
_idx     :: !(Int, Int)
                     , Entry a -> Maybe Int
_rowNext :: !(Maybe Int)
                     , Entry a -> Maybe Int
_colNext :: !(Maybe Int)
                     } deriving (ReadPrec [Entry a]
ReadPrec (Entry a)
Int -> ReadS (Entry a)
ReadS [Entry a]
(Int -> ReadS (Entry a))
-> ReadS [Entry a]
-> ReadPrec (Entry a)
-> ReadPrec [Entry a]
-> Read (Entry a)
forall a. Read a => ReadPrec [Entry a]
forall a. Read a => ReadPrec (Entry a)
forall a. Read a => Int -> ReadS (Entry a)
forall a. Read a => ReadS [Entry a]
forall a.
(Int -> ReadS a)
-> ReadS [a] -> ReadPrec a -> ReadPrec [a] -> Read a
readListPrec :: ReadPrec [Entry a]
$creadListPrec :: forall a. Read a => ReadPrec [Entry a]
readPrec :: ReadPrec (Entry a)
$creadPrec :: forall a. Read a => ReadPrec (Entry a)
readList :: ReadS [Entry a]
$creadList :: forall a. Read a => ReadS [Entry a]
readsPrec :: Int -> ReadS (Entry a)
$creadsPrec :: forall a. Read a => Int -> ReadS (Entry a)
Read, Int -> Entry a -> ShowS
[Entry a] -> ShowS
Entry a -> String
(Int -> Entry a -> ShowS)
-> (Entry a -> String) -> ([Entry a] -> ShowS) -> Show (Entry a)
forall a. Show a => Int -> Entry a -> ShowS
forall a. Show a => [Entry a] -> ShowS
forall a. Show a => Entry a -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [Entry a] -> ShowS
$cshowList :: forall a. Show a => [Entry a] -> ShowS
show :: Entry a -> String
$cshow :: forall a. Show a => Entry a -> String
showsPrec :: Int -> Entry a -> ShowS
$cshowsPrec :: forall a. Show a => Int -> Entry a -> ShowS
Show, Entry a -> Entry a -> Bool
(Entry a -> Entry a -> Bool)
-> (Entry a -> Entry a -> Bool) -> Eq (Entry a)
forall a. Eq a => Entry a -> Entry a -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: Entry a -> Entry a -> Bool
$c/= :: forall a. Eq a => Entry a -> Entry a -> Bool
== :: Entry a -> Entry a -> Bool
$c== :: forall a. Eq a => Entry a -> Entry a -> Bool
Eq, Eq (Entry a)
Eq (Entry a)
-> (Entry a -> Entry a -> Ordering)
-> (Entry a -> Entry a -> Bool)
-> (Entry a -> Entry a -> Bool)
-> (Entry a -> Entry a -> Bool)
-> (Entry a -> Entry a -> Bool)
-> (Entry a -> Entry a -> Entry a)
-> (Entry a -> Entry a -> Entry a)
-> Ord (Entry a)
Entry a -> Entry a -> Bool
Entry a -> Entry a -> Ordering
Entry a -> Entry a -> Entry a
forall a.
Eq a
-> (a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
forall a. Ord a => Eq (Entry a)
forall a. Ord a => Entry a -> Entry a -> Bool
forall a. Ord a => Entry a -> Entry a -> Ordering
forall a. Ord a => Entry a -> Entry a -> Entry a
min :: Entry a -> Entry a -> Entry a
$cmin :: forall a. Ord a => Entry a -> Entry a -> Entry a
max :: Entry a -> Entry a -> Entry a
$cmax :: forall a. Ord a => Entry a -> Entry a -> Entry a
>= :: Entry a -> Entry a -> Bool
$c>= :: forall a. Ord a => Entry a -> Entry a -> Bool
> :: Entry a -> Entry a -> Bool
$c> :: forall a. Ord a => Entry a -> Entry a -> Bool
<= :: Entry a -> Entry a -> Bool
$c<= :: forall a. Ord a => Entry a -> Entry a -> Bool
< :: Entry a -> Entry a -> Bool
$c< :: forall a. Ord a => Entry a -> Entry a -> Bool
compare :: Entry a -> Entry a -> Ordering
$ccompare :: forall a. Ord a => Entry a -> Entry a -> Ordering
$cp1Ord :: forall a. Ord a => Eq (Entry a)
Ord)

makeLenses ''Entry

newEntry :: a -> Entry a
newEntry :: a -> Entry a
newEntry a
v = a -> (Int, Int) -> Maybe Int -> Maybe Int -> Entry a
forall a. a -> (Int, Int) -> Maybe Int -> Maybe Int -> Entry a
Entry a
v (-Int
1,-Int
1) Maybe Int
forall a. Maybe a
Nothing Maybe Int
forall a. Maybe a
Nothing

data Matrix a = Matrix { Matrix a -> Vector (Entry a)
_coefficients :: !(Vector (Entry a))
                       , Matrix a -> IntMap Int
_rowStart     :: !(IntMap Int)
                       , Matrix a -> IntMap Int
_colStart     :: !(IntMap Int)
                       , Matrix a -> Int
_height       :: !Int
                       , Matrix a -> Int
_width        :: !Int
                       } deriving (ReadPrec [Matrix a]
ReadPrec (Matrix a)
Int -> ReadS (Matrix a)
ReadS [Matrix a]
(Int -> ReadS (Matrix a))
-> ReadS [Matrix a]
-> ReadPrec (Matrix a)
-> ReadPrec [Matrix a]
-> Read (Matrix a)
forall a. Read a => ReadPrec [Matrix a]
forall a. Read a => ReadPrec (Matrix a)
forall a. Read a => Int -> ReadS (Matrix a)
forall a. Read a => ReadS [Matrix a]
forall a.
(Int -> ReadS a)
-> ReadS [a] -> ReadPrec a -> ReadPrec [a] -> Read a
readListPrec :: ReadPrec [Matrix a]
$creadListPrec :: forall a. Read a => ReadPrec [Matrix a]
readPrec :: ReadPrec (Matrix a)
$creadPrec :: forall a. Read a => ReadPrec (Matrix a)
readList :: ReadS [Matrix a]
$creadList :: forall a. Read a => ReadS [Matrix a]
readsPrec :: Int -> ReadS (Matrix a)
$creadsPrec :: forall a. Read a => Int -> ReadS (Matrix a)
Read, Int -> Matrix a -> ShowS
[Matrix a] -> ShowS
Matrix a -> String
(Int -> Matrix a -> ShowS)
-> (Matrix a -> String) -> ([Matrix a] -> ShowS) -> Show (Matrix a)
forall a. Show a => Int -> Matrix a -> ShowS
forall a. Show a => [Matrix a] -> ShowS
forall a. Show a => Matrix a -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [Matrix a] -> ShowS
$cshowList :: forall a. Show a => [Matrix a] -> ShowS
show :: Matrix a -> String
$cshow :: forall a. Show a => Matrix a -> String
showsPrec :: Int -> Matrix a -> ShowS
$cshowsPrec :: forall a. Show a => Int -> Matrix a -> ShowS
Show)

makeLenses ''Matrix

data BuildState = BuildState { BuildState -> IntMap Int
_colMap :: !(IntMap Int)
                             , BuildState -> IntMap Int
_rowMap :: !(IntMap Int)
                             , BuildState -> Int
_curIdx :: !Int
                             }
makeLenses ''BuildState

data GaussianState a = GaussianState { GaussianState a -> Matrix a
_input     :: !(Matrix a)
                                     , GaussianState a -> Matrix a
_output    :: !(Matrix a)
                                     , GaussianState a -> Int
_prevCol   :: !Int
                                     , GaussianState a -> IntSet
_heavyCols :: !IntSet
                                     , GaussianState a -> Int
_curRow    :: !Int
                                     , GaussianState a -> a
_detAcc    :: !a
                                     }
makeLenses ''GaussianState

instance Eq a => Eq (Matrix a) where
  Matrix a
n == :: Matrix a -> Matrix a -> Bool
== Matrix a
m =
    Matrix a
nMatrix a -> Getting Int (Matrix a) Int -> Int
forall s a. s -> Getting a s a -> a
^.Getting Int (Matrix a) Int
forall a. Lens' (Matrix a) Int
height Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Matrix a
mMatrix a -> Getting Int (Matrix a) Int -> Int
forall s a. s -> Getting a s a -> a
^.Getting Int (Matrix a) Int
forall a. Lens' (Matrix a) Int
height Bool -> Bool -> Bool
&& Matrix a
nMatrix a -> Getting Int (Matrix a) Int -> Int
forall s a. s -> Getting a s a -> a
^.Getting Int (Matrix a) Int
forall a. Lens' (Matrix a) Int
width Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Matrix a
mMatrix a -> Getting Int (Matrix a) Int -> Int
forall s a. s -> Getting a s a -> a
^.Getting Int (Matrix a) Int
forall a. Lens' (Matrix a) Int
width Bool -> Bool -> Bool
&& Matrix a -> [((Int, Int), a)]
forall b. Matrix b -> [((Int, Int), b)]
content Matrix a
n [((Int, Int), a)] -> [((Int, Int), a)] -> Bool
forall a. Eq a => a -> a -> Bool
== Matrix a -> [((Int, Int), a)]
forall b. Matrix b -> [((Int, Int), b)]
content Matrix a
m
    where
      content :: Matrix b -> [((Int, Int), b)]
content = (((Int, Int), b) -> ((Int, Int), b) -> Ordering)
-> [((Int, Int), b)] -> [((Int, Int), b)]
forall a. (a -> a -> Ordering) -> [a] -> [a]
sortBy ((((Int, Int), b) -> (Int, Int))
-> ((Int, Int), b) -> ((Int, Int), b) -> Ordering
forall a b. Ord a => (b -> a) -> b -> b -> Ordering
comparing ((Int, Int), b) -> (Int, Int)
forall a b. (a, b) -> a
fst) ([((Int, Int), b)] -> [((Int, Int), b)])
-> (Matrix b -> [((Int, Int), b)]) -> Matrix b -> [((Int, Int), b)]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector ((Int, Int), b) -> [((Int, Int), b)]
forall a. Vector a -> [a]
V.toList (Vector ((Int, Int), b) -> [((Int, Int), b)])
-> (Matrix b -> Vector ((Int, Int), b))
-> Matrix b
-> [((Int, Int), b)]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Matrix b -> Vector ((Int, Int), b)
forall a. Matrix a -> Vector ((Int, Int), a)
nonZeroEntries

data MaxEntry a b = MaxEntry { MaxEntry a b -> a
_weight :: !a
                             , MaxEntry a b -> b
entry   :: b
                             } deriving (ReadPrec [MaxEntry a b]
ReadPrec (MaxEntry a b)
Int -> ReadS (MaxEntry a b)
ReadS [MaxEntry a b]
(Int -> ReadS (MaxEntry a b))
-> ReadS [MaxEntry a b]
-> ReadPrec (MaxEntry a b)
-> ReadPrec [MaxEntry a b]
-> Read (MaxEntry a b)
forall a.
(Int -> ReadS a)
-> ReadS [a] -> ReadPrec a -> ReadPrec [a] -> Read a
forall a b. (Read a, Read b) => ReadPrec [MaxEntry a b]
forall a b. (Read a, Read b) => ReadPrec (MaxEntry a b)
forall a b. (Read a, Read b) => Int -> ReadS (MaxEntry a b)
forall a b. (Read a, Read b) => ReadS [MaxEntry a b]
readListPrec :: ReadPrec [MaxEntry a b]
$creadListPrec :: forall a b. (Read a, Read b) => ReadPrec [MaxEntry a b]
readPrec :: ReadPrec (MaxEntry a b)
$creadPrec :: forall a b. (Read a, Read b) => ReadPrec (MaxEntry a b)
readList :: ReadS [MaxEntry a b]
$creadList :: forall a b. (Read a, Read b) => ReadS [MaxEntry a b]
readsPrec :: Int -> ReadS (MaxEntry a b)
$creadsPrec :: forall a b. (Read a, Read b) => Int -> ReadS (MaxEntry a b)
Read, Int -> MaxEntry a b -> ShowS
[MaxEntry a b] -> ShowS
MaxEntry a b -> String
(Int -> MaxEntry a b -> ShowS)
-> (MaxEntry a b -> String)
-> ([MaxEntry a b] -> ShowS)
-> Show (MaxEntry a b)
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
forall a b. (Show a, Show b) => Int -> MaxEntry a b -> ShowS
forall a b. (Show a, Show b) => [MaxEntry a b] -> ShowS
forall a b. (Show a, Show b) => MaxEntry a b -> String
showList :: [MaxEntry a b] -> ShowS
$cshowList :: forall a b. (Show a, Show b) => [MaxEntry a b] -> ShowS
show :: MaxEntry a b -> String
$cshow :: forall a b. (Show a, Show b) => MaxEntry a b -> String
showsPrec :: Int -> MaxEntry a b -> ShowS
$cshowsPrec :: forall a b. (Show a, Show b) => Int -> MaxEntry a b -> ShowS
Show, MaxEntry a b -> MaxEntry a b -> Bool
(MaxEntry a b -> MaxEntry a b -> Bool)
-> (MaxEntry a b -> MaxEntry a b -> Bool) -> Eq (MaxEntry a b)
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
forall a b. (Eq a, Eq b) => MaxEntry a b -> MaxEntry a b -> Bool
/= :: MaxEntry a b -> MaxEntry a b -> Bool
$c/= :: forall a b. (Eq a, Eq b) => MaxEntry a b -> MaxEntry a b -> Bool
== :: MaxEntry a b -> MaxEntry a b -> Bool
$c== :: forall a b. (Eq a, Eq b) => MaxEntry a b -> MaxEntry a b -> Bool
Eq, Eq (MaxEntry a b)
Eq (MaxEntry a b)
-> (MaxEntry a b -> MaxEntry a b -> Ordering)
-> (MaxEntry a b -> MaxEntry a b -> Bool)
-> (MaxEntry a b -> MaxEntry a b -> Bool)
-> (MaxEntry a b -> MaxEntry a b -> Bool)
-> (MaxEntry a b -> MaxEntry a b -> Bool)
-> (MaxEntry a b -> MaxEntry a b -> MaxEntry a b)
-> (MaxEntry a b -> MaxEntry a b -> MaxEntry a b)
-> Ord (MaxEntry a b)
MaxEntry a b -> MaxEntry a b -> Bool
MaxEntry a b -> MaxEntry a b -> Ordering
MaxEntry a b -> MaxEntry a b -> MaxEntry a b
forall a.
Eq a
-> (a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
forall a b. (Ord a, Ord b) => Eq (MaxEntry a b)
forall a b. (Ord a, Ord b) => MaxEntry a b -> MaxEntry a b -> Bool
forall a b.
(Ord a, Ord b) =>
MaxEntry a b -> MaxEntry a b -> Ordering
forall a b.
(Ord a, Ord b) =>
MaxEntry a b -> MaxEntry a b -> MaxEntry a b
min :: MaxEntry a b -> MaxEntry a b -> MaxEntry a b
$cmin :: forall a b.
(Ord a, Ord b) =>
MaxEntry a b -> MaxEntry a b -> MaxEntry a b
max :: MaxEntry a b -> MaxEntry a b -> MaxEntry a b
$cmax :: forall a b.
(Ord a, Ord b) =>
MaxEntry a b -> MaxEntry a b -> MaxEntry a b
>= :: MaxEntry a b -> MaxEntry a b -> Bool
$c>= :: forall a b. (Ord a, Ord b) => MaxEntry a b -> MaxEntry a b -> Bool
> :: MaxEntry a b -> MaxEntry a b -> Bool
$c> :: forall a b. (Ord a, Ord b) => MaxEntry a b -> MaxEntry a b -> Bool
<= :: MaxEntry a b -> MaxEntry a b -> Bool
$c<= :: forall a b. (Ord a, Ord b) => MaxEntry a b -> MaxEntry a b -> Bool
< :: MaxEntry a b -> MaxEntry a b -> Bool
$c< :: forall a b. (Ord a, Ord b) => MaxEntry a b -> MaxEntry a b -> Bool
compare :: MaxEntry a b -> MaxEntry a b -> Ordering
$ccompare :: forall a b.
(Ord a, Ord b) =>
MaxEntry a b -> MaxEntry a b -> Ordering
$cp1Ord :: forall a b. (Ord a, Ord b) => Eq (MaxEntry a b)
Ord)

empty :: Matrix a
empty :: Matrix a
empty = Vector (Entry a)
-> IntMap Int -> IntMap Int -> Int -> Int -> Matrix a
forall a.
Vector (Entry a)
-> IntMap Int -> IntMap Int -> Int -> Int -> Matrix a
Matrix Vector (Entry a)
forall a. Vector a
V.empty IntMap Int
forall a. IntMap a
IM.empty IntMap Int
forall a. IntMap a
IM.empty Int
0 Int
0

fromLists :: DecidableZero a => [[a]] -> Matrix a
fromLists :: [[a]] -> Matrix a
fromLists [[a]]
xss =
  [((Int, Int), a)] -> Matrix a
forall a. DecidableZero a => [((Int, Int), a)] -> Matrix a
fromList ([[((Int, Int), a)]] -> [((Int, Int), a)]
forall w. Monoid w => [w] -> w
concat ([[((Int, Int), a)]] -> [((Int, Int), a)])
-> [[((Int, Int), a)]] -> [((Int, Int), a)]
forall a b. (a -> b) -> a -> b
$ (Int -> [a] -> [((Int, Int), a)])
-> [Int] -> [[a]] -> [[((Int, Int), a)]]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Int
i -> (Int -> a -> ((Int, Int), a)) -> [Int] -> [a] -> [((Int, Int), a)]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Int
j -> ((Int
i,Int
j),)) [Int
0..]) [Int
0..] [[a]]
xss)
  Matrix a -> (Matrix a -> Matrix a) -> Matrix a
forall a b. a -> (a -> b) -> b
& (Int -> Identity Int) -> Matrix a -> Identity (Matrix a)
forall a. Lens' (Matrix a) Int
width  ((Int -> Identity Int) -> Matrix a -> Identity (Matrix a))
-> Int -> Matrix a -> Matrix a
forall s t a b. ASetter s t a b -> b -> s -> t
.~ [Int] -> Int
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum (Int
0 Int -> [Int] -> [Int]
forall a. a -> [a] -> [a]
: ([a] -> Int) -> [[a]] -> [Int]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
map [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[a]]
xss)
  Matrix a -> (Matrix a -> Matrix a) -> Matrix a
forall a b. a -> (a -> b) -> b
& (Int -> Identity Int) -> Matrix a -> Identity (Matrix a)
forall a. Lens' (Matrix a) Int
height ((Int -> Identity Int) -> Matrix a -> Identity (Matrix a))
-> Int -> Matrix a -> Matrix a
forall s t a b. ASetter s t a b -> b -> s -> t
.~ [[a]] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[a]]
xss

fromCols :: DecidableZero a => [Vector a] -> Matrix a
fromCols :: [Vector a] -> Matrix a
fromCols [Vector a]
xss =
  let h :: Int
h = [Int] -> Int
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum ([Int] -> Int) -> [Int] -> Int
forall a b. (a -> b) -> a -> b
$ (Vector a -> Int) -> [Vector a] -> [Int]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
map Vector a -> Int
forall a. Vector a -> Int
V.length [Vector a]
xss
      w :: Int
w = [Vector a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Vector a]
xss
  in [((Int, Int), a)] -> Matrix a
forall a. DecidableZero a => [((Int, Int), a)] -> Matrix a
fromList ([[((Int, Int), a)]] -> [((Int, Int), a)]
forall w. Monoid w => [w] -> w
concat ([[((Int, Int), a)]] -> [((Int, Int), a)])
-> [[((Int, Int), a)]] -> [((Int, Int), a)]
forall a b. (a -> b) -> a -> b
$ (Int -> Vector a -> [((Int, Int), a)])
-> [Int] -> [Vector a] -> [[((Int, Int), a)]]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Int
i -> Vector ((Int, Int), a) -> [((Int, Int), a)]
forall a. Vector a -> [a]
V.toList (Vector ((Int, Int), a) -> [((Int, Int), a)])
-> (Vector a -> Vector ((Int, Int), a))
-> Vector a
-> [((Int, Int), a)]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Int -> a -> ((Int, Int), a)) -> Vector a -> Vector ((Int, Int), a)
forall a b. (Int -> a -> b) -> Vector a -> Vector b
V.imap (\Int
j -> ((Int
j,Int
i),))) [Int
0..] [Vector a]
xss)
     Matrix a -> (Matrix a -> Matrix a) -> Matrix a
forall a b. a -> (a -> b) -> b
& (Int -> Identity Int) -> Matrix a -> Identity (Matrix a)
forall a. Lens' (Matrix a) Int
width ((Int -> Identity Int) -> Matrix a -> Identity (Matrix a))
-> Int -> Matrix a -> Matrix a
forall s t a b. ASetter s t a b -> b -> s -> t
.~ Int
w
     Matrix a -> (Matrix a -> Matrix a) -> Matrix a
forall a b. a -> (a -> b) -> b
& (Int -> Identity Int) -> Matrix a -> Identity (Matrix a)
forall a. Lens' (Matrix a) Int
height ((Int -> Identity Int) -> Matrix a -> Identity (Matrix a))
-> Int -> Matrix a -> Matrix a
forall s t a b. ASetter s t a b -> b -> s -> t
.~ Int
h

fromList :: DecidableZero a => [((Int, Int), a)] -> Matrix a
fromList :: [((Int, Int), a)] -> Matrix a
fromList [((Int, Int), a)]
cs =
  let ([Entry a]
as, BuildState
bs) = State BuildState [Entry a] -> BuildState -> ([Entry a], BuildState)
forall s a. State s a -> s -> (a, s)
runState ((((Int, Int), a) -> StateT BuildState Identity (Entry a))
-> [((Int, Int), a)] -> State BuildState [Entry a]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
mapM ((Int, Int), a) -> StateT BuildState Identity (Entry a)
forall (m :: * -> *) a.
MonadState BuildState m =>
((Int, Int), a) -> m (Entry a)
initialize ([((Int, Int), a)] -> State BuildState [Entry a])
-> [((Int, Int), a)] -> State BuildState [Entry a]
forall a b. (a -> b) -> a -> b
$ (((Int, Int), a) -> Bool) -> [((Int, Int), a)] -> [((Int, Int), a)]
forall a. (a -> Bool) -> [a] -> [a]
filter (Getting Bool ((Int, Int), a) Bool -> ((Int, Int), a) -> Bool
forall s (m :: * -> *) a. MonadReader s m => Getting a s a -> m a
view (Getting Bool ((Int, Int), a) Bool -> ((Int, Int), a) -> Bool)
-> Getting Bool ((Int, Int), a) Bool -> ((Int, Int), a) -> Bool
forall a b. (a -> b) -> a -> b
$ (a -> Const Bool a)
-> ((Int, Int), a) -> Const Bool ((Int, Int), a)
forall s t a b. Field2 s t a b => Lens s t a b
_2 ((a -> Const Bool a)
 -> ((Int, Int), a) -> Const Bool ((Int, Int), a))
-> ((Bool -> Const Bool Bool) -> a -> Const Bool a)
-> Getting Bool ((Int, Int), a) Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (a -> Bool) -> (Bool -> Const Bool Bool) -> a -> Const Bool a
forall (p :: * -> * -> *) (f :: * -> *) s a.
(Profunctor p, Contravariant f) =>
(s -> a) -> Optic' p f s a
to (Bool -> Bool
not(Bool -> Bool) -> (a -> Bool) -> a -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
.a -> Bool
forall r. DecidableZero r => r -> Bool
isZero)) [((Int, Int), a)]
cs)
                 (IntMap Int -> IntMap Int -> Int -> BuildState
BuildState IntMap Int
forall a. IntMap a
IM.empty IntMap Int
forall a. IntMap a
IM.empty (-Int
1))
      vec :: Vector (Entry a)
vec = [Entry a] -> Vector (Entry a)
forall a. [a] -> Vector a
V.fromList [Entry a]
as
      h :: Int
h = [Int] -> Int
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum (Int
0Int -> [Int] -> [Int]
forall a. a -> [a] -> [a]
:(((Int, Int), a) -> Int) -> [((Int, Int), a)] -> [Int]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
map (Getting Int ((Int, Int), a) Int -> ((Int, Int), a) -> Int
forall s (m :: * -> *) a. MonadReader s m => Getting a s a -> m a
view (Getting Int ((Int, Int), a) Int -> ((Int, Int), a) -> Int)
-> Getting Int ((Int, Int), a) Int -> ((Int, Int), a) -> Int
forall a b. (a -> b) -> a -> b
$ ((Int, Int) -> Const Int (Int, Int))
-> ((Int, Int), a) -> Const Int ((Int, Int), a)
forall s t a b. Field1 s t a b => Lens s t a b
_1(((Int, Int) -> Const Int (Int, Int))
 -> ((Int, Int), a) -> Const Int ((Int, Int), a))
-> ((Int -> Const Int Int) -> (Int, Int) -> Const Int (Int, Int))
-> Getting Int ((Int, Int), a) Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
.(Int -> Const Int Int) -> (Int, Int) -> Const Int (Int, Int)
forall s t a b. Field1 s t a b => Lens s t a b
_1) [((Int, Int), a)]
cs) Int -> Int -> Int
forall r. Additive r => r -> r -> r
+ Int
1
      w :: Int
w = [Int] -> Int
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum (Int
0Int -> [Int] -> [Int]
forall a. a -> [a] -> [a]
:(((Int, Int), a) -> Int) -> [((Int, Int), a)] -> [Int]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
map (Getting Int ((Int, Int), a) Int -> ((Int, Int), a) -> Int
forall s (m :: * -> *) a. MonadReader s m => Getting a s a -> m a
view (Getting Int ((Int, Int), a) Int -> ((Int, Int), a) -> Int)
-> Getting Int ((Int, Int), a) Int -> ((Int, Int), a) -> Int
forall a b. (a -> b) -> a -> b
$ ((Int, Int) -> Const Int (Int, Int))
-> ((Int, Int), a) -> Const Int ((Int, Int), a)
forall s t a b. Field1 s t a b => Lens s t a b
_1(((Int, Int) -> Const Int (Int, Int))
 -> ((Int, Int), a) -> Const Int ((Int, Int), a))
-> ((Int -> Const Int Int) -> (Int, Int) -> Const Int (Int, Int))
-> Getting Int ((Int, Int), a) Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
.(Int -> Const Int Int) -> (Int, Int) -> Const Int (Int, Int)
forall s t a b. Field2 s t a b => Lens s t a b
_2) [((Int, Int), a)]
cs) Int -> Int -> Int
forall r. Additive r => r -> r -> r
+ Int
1
  in Vector (Entry a)
-> IntMap Int -> IntMap Int -> Int -> Int -> Matrix a
forall a.
Vector (Entry a)
-> IntMap Int -> IntMap Int -> Int -> Int -> Matrix a
Matrix Vector (Entry a)
vec (BuildState
bsBuildState
-> Getting (IntMap Int) BuildState (IntMap Int) -> IntMap Int
forall s a. s -> Getting a s a -> a
^.Getting (IntMap Int) BuildState (IntMap Int)
Lens' BuildState (IntMap Int)
rowMap) (BuildState
bsBuildState
-> Getting (IntMap Int) BuildState (IntMap Int) -> IntMap Int
forall s a. s -> Getting a s a -> a
^.Getting (IntMap Int) BuildState (IntMap Int)
Lens' BuildState (IntMap Int)
colMap) Int
h Int
w
  where
    initialize :: ((Int, Int), a) -> m (Entry a)
initialize ((Int
i, Int
j), a
c) =  do
        (Int -> Identity Int) -> BuildState -> Identity BuildState
Lens' BuildState Int
curIdx ((Int -> Identity Int) -> BuildState -> Identity BuildState)
-> Int -> m ()
forall s (m :: * -> *) a.
(MonadState s m, Num a) =>
ASetter' s a -> a -> m ()
+= Int
1
        Int
n <- Getting Int BuildState Int -> m Int
forall s (m :: * -> *) a. MonadState s m => Getting a s a -> m a
use Getting Int BuildState Int
Lens' BuildState Int
curIdx
        Maybe Int
nc <- Getting (Maybe Int) BuildState (Maybe Int) -> m (Maybe Int)
forall s (m :: * -> *) a. MonadState s m => Getting a s a -> m a
use (Getting (Maybe Int) BuildState (Maybe Int) -> m (Maybe Int))
-> Getting (Maybe Int) BuildState (Maybe Int) -> m (Maybe Int)
forall a b. (a -> b) -> a -> b
$ (IntMap Int -> Const (Maybe Int) (IntMap Int))
-> BuildState -> Const (Maybe Int) BuildState
Lens' BuildState (IntMap Int)
colMap((IntMap Int -> Const (Maybe Int) (IntMap Int))
 -> BuildState -> Const (Maybe Int) BuildState)
-> ((Maybe Int -> Const (Maybe Int) (Maybe Int))
    -> IntMap Int -> Const (Maybe Int) (IntMap Int))
-> Getting (Maybe Int) BuildState (Maybe Int)
forall b c a. (b -> c) -> (a -> b) -> a -> c
.Index (IntMap Int)
-> Lens' (IntMap Int) (Maybe (IxValue (IntMap Int)))
forall m. At m => Index m -> Lens' m (Maybe (IxValue m))
at Int
Index (IntMap Int)
j
        Maybe Int
nr <- Getting (Maybe Int) BuildState (Maybe Int) -> m (Maybe Int)
forall s (m :: * -> *) a. MonadState s m => Getting a s a -> m a
use (Getting (Maybe Int) BuildState (Maybe Int) -> m (Maybe Int))
-> Getting (Maybe Int) BuildState (Maybe Int) -> m (Maybe Int)
forall a b. (a -> b) -> a -> b
$ (IntMap Int -> Const (Maybe Int) (IntMap Int))
-> BuildState -> Const (Maybe Int) BuildState
Lens' BuildState (IntMap Int)
rowMap((IntMap Int -> Const (Maybe Int) (IntMap Int))
 -> BuildState -> Const (Maybe Int) BuildState)
-> ((Maybe Int -> Const (Maybe Int) (Maybe Int))
    -> IntMap Int -> Const (Maybe Int) (IntMap Int))
-> Getting (Maybe Int) BuildState (Maybe Int)
forall b c a. (b -> c) -> (a -> b) -> a -> c
.Index (IntMap Int)
-> Lens' (IntMap Int) (Maybe (IxValue (IntMap Int)))
forall m. At m => Index m -> Lens' m (Maybe (IxValue m))
at Int
Index (IntMap Int)
i
        (IntMap Int -> Identity (IntMap Int))
-> BuildState -> Identity BuildState
Lens' BuildState (IntMap Int)
colMap ((IntMap Int -> Identity (IntMap Int))
 -> BuildState -> Identity BuildState)
-> (IntMap Int -> IntMap Int) -> m ()
forall s (m :: * -> *) a b.
MonadState s m =>
ASetter s s a b -> (a -> b) -> m ()
%= Int -> Int -> IntMap Int -> IntMap Int
forall a. Int -> a -> IntMap a -> IntMap a
insert Int
j Int
n
        (IntMap Int -> Identity (IntMap Int))
-> BuildState -> Identity BuildState
Lens' BuildState (IntMap Int)
rowMap ((IntMap Int -> Identity (IntMap Int))
 -> BuildState -> Identity BuildState)
-> (IntMap Int -> IntMap Int) -> m ()
forall s (m :: * -> *) a b.
MonadState s m =>
ASetter s s a b -> (a -> b) -> m ()
%= Int -> Int -> IntMap Int -> IntMap Int
forall a. Int -> a -> IntMap a -> IntMap a
insert Int
i Int
n
        Entry a -> m (Entry a)
forall (m :: * -> *) a. Monad m => a -> m a
return (Entry a -> m (Entry a)) -> Entry a -> m (Entry a)
forall a b. (a -> b) -> a -> b
$ Entry :: forall a. a -> (Int, Int) -> Maybe Int -> Maybe Int -> Entry a
Entry { _value :: a
_value = a
c
                       , _idx :: (Int, Int)
_idx = (Int
i, Int
j)
                       , _rowNext :: Maybe Int
_rowNext = Maybe Int
nr
                       , _colNext :: Maybe Int
_colNext = Maybe Int
nc
                       }


getDiag :: Monoidal a => Matrix a -> Vector a
getDiag :: Matrix a -> Vector a
getDiag Matrix a
mat = Int -> (Int -> a) -> Vector a
forall a. Int -> (Int -> a) -> Vector a
V.generate (Int -> Int -> Int
forall a. Ord a => a -> a -> a
min (Matrix a
matMatrix a -> Getting Int (Matrix a) Int -> Int
forall s a. s -> Getting a s a -> a
^.Getting Int (Matrix a) Int
forall a. Lens' (Matrix a) Int
height) (Matrix a
matMatrix a -> Getting Int (Matrix a) Int -> Int
forall s a. s -> Getting a s a -> a
^.Getting Int (Matrix a) Int
forall a. Lens' (Matrix a) Int
width)) ((Int -> a) -> Vector a) -> (Int -> a) -> Vector a
forall a b. (a -> b) -> a -> b
$ \Int
i ->
  a -> Maybe a -> a
forall a. a -> Maybe a -> a
fromMaybe a
forall m. Monoidal m => m
zero (Maybe a -> a) -> Maybe a -> a
forall a b. (a -> b) -> a -> b
$ Maybe a
-> (Maybe a -> Int -> Entry a -> Maybe a)
-> Direction
-> Int
-> Matrix a
-> Maybe a
forall b a.
b
-> (b -> Int -> Entry a -> b) -> Direction -> Int -> Matrix a -> b
traverseDir Maybe a
forall a. Maybe a
Nothing (\Maybe a
a Int
_ Entry a
e -> Maybe a
a Maybe a -> Maybe a -> Maybe a
forall (f :: * -> *) a. Alternative f => f a -> f a -> f a
<|> if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Entry a
eEntry a -> Getting Int (Entry a) Int -> Int
forall s a. s -> Getting a s a -> a
^.Direction -> Lens' (Entry a) Int
forall a. Direction -> Lens' (Entry a) Int
nthL Direction
Row
                                                        then a -> Maybe a
forall a. a -> Maybe a
Just (Entry a
eEntry a -> Getting a (Entry a) a -> a
forall s a. s -> Getting a s a -> a
^.Getting a (Entry a) a
forall a a. Lens (Entry a) (Entry a) a a
value )
                                                        else Maybe a
forall a. Maybe a
Nothing) Direction
Row Int
i Matrix a
mat

diagProd :: (Unital c, Monoidal c) => Matrix c -> c
diagProd :: Matrix c -> c
diagProd = (c -> c -> c) -> c -> Vector c -> c
forall a b. (a -> b -> b) -> b -> Vector a -> b
V.foldr' c -> c -> c
forall r. Multiplicative r => r -> r -> r
(*) c
forall r. Unital r => r
one (Vector c -> c) -> (Matrix c -> Vector c) -> Matrix c -> c
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Matrix c -> Vector c
forall a. Monoidal a => Matrix a -> Vector a
getDiag

trace :: Monoidal c => Matrix c -> c
trace :: Matrix c -> c
trace = (c -> c -> c) -> c -> Vector c -> c
forall a b. (a -> b -> b) -> b -> Vector a -> b
V.foldr' c -> c -> c
forall r. Additive r => r -> r -> r
(+) c
forall m. Monoidal m => m
zero (Vector c -> c) -> (Matrix c -> Vector c) -> Matrix c -> c
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Matrix c -> Vector c
forall a. Monoidal a => Matrix a -> Vector a
getDiag

toLists :: forall a. Monoidal a => Matrix a -> [[a]]
toLists :: Matrix a -> [[a]]
toLists Matrix a
mat =
  let orig :: [[a]]
orig = Int -> [a] -> [[a]]
forall a. Int -> a -> [a]
replicate (Matrix a -> Int
forall a. Matrix a -> Int
_height Matrix a
mat) ([a] -> [[a]]) -> [a] -> [[a]]
forall a b. (a -> b) -> a -> b
$ Int -> a -> [a]
forall a. Int -> a -> [a]
replicate (Matrix a -> Int
forall a. Matrix a -> Int
_width Matrix a
mat) (a
forall m. Monoidal m => m
zero :: a)
  in [Entry (IxValue (IxValue [[a]]))] -> [[a]] -> [[a]]
forall t.
(Ixed t, Ixed (IxValue t), Index (IxValue t) ~ Int,
 Index t ~ Int) =>
[Entry (IxValue (IxValue t))] -> t -> t
go (Vector (Entry a) -> [Entry a]
forall a. Vector a -> [a]
V.toList (Vector (Entry a) -> [Entry a]) -> Vector (Entry a) -> [Entry a]
forall a b. (a -> b) -> a -> b
$ Matrix a -> Vector (Entry a)
forall a. Matrix a -> Vector (Entry a)
_coefficients Matrix a
mat) [[a]]
orig
  where
    go :: [Entry (IxValue (IxValue t))] -> t -> t
go [] t
m = t
m
    go (Entry{_value :: forall a. Entry a -> a
_value = IxValue (IxValue t)
v, _idx :: forall a. Entry a -> (Int, Int)
_idx = (Int
i,Int
j) }:[Entry (IxValue (IxValue t))]
es) t
m =
      [Entry (IxValue (IxValue t))] -> t -> t
go [Entry (IxValue (IxValue t))]
es (t
mt -> (t -> t) -> t
forall a b. a -> (a -> b) -> b
&Index t -> Traversal' t (IxValue t)
forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix Int
Index t
i((IxValue t -> Identity (IxValue t)) -> t -> Identity t)
-> ((IxValue (IxValue t) -> Identity (IxValue (IxValue t)))
    -> IxValue t -> Identity (IxValue t))
-> (IxValue (IxValue t) -> Identity (IxValue (IxValue t)))
-> t
-> Identity t
forall b c a. (b -> c) -> (a -> b) -> a -> c
.Index (IxValue t) -> Traversal' (IxValue t) (IxValue (IxValue t))
forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix Int
Index (IxValue t)
j ((IxValue (IxValue t) -> Identity (IxValue (IxValue t)))
 -> t -> Identity t)
-> IxValue (IxValue t) -> t -> t
forall s t a b. ASetter s t a b -> b -> s -> t
.~ IxValue (IxValue t)
v)

swapRows :: Int -> Int -> Matrix a -> Matrix a
swapRows :: Int -> Int -> Matrix a -> Matrix a
swapRows = Direction -> Int -> Int -> Matrix a -> Matrix a
forall a. Direction -> Int -> Int -> Matrix a -> Matrix a
swapper Direction
Row

swapCols :: Int -> Int -> Matrix a -> Matrix a
swapCols :: Int -> Int -> Matrix a -> Matrix a
swapCols = Direction -> Int -> Int -> Matrix a -> Matrix a
forall a. Direction -> Int -> Int -> Matrix a -> Matrix a
swapper Direction
Column

swapper :: Direction -> Int -> Int -> Matrix a -> Matrix a
swapper :: Direction -> Int -> Int -> Matrix a -> Matrix a
swapper Direction
dir Int
i Int
j Matrix a
mat =
  let ith :: Maybe Int
ith = Matrix a
matMatrix a -> Getting (Maybe Int) (Matrix a) (Maybe Int) -> Maybe Int
forall s a. s -> Getting a s a -> a
^.Direction -> Lens' (Matrix a) (IntMap Int)
forall a. Direction -> Lens' (Matrix a) (IntMap Int)
startL Direction
dir((IntMap Int -> Const (Maybe Int) (IntMap Int))
 -> Matrix a -> Const (Maybe Int) (Matrix a))
-> ((Maybe Int -> Const (Maybe Int) (Maybe Int))
    -> IntMap Int -> Const (Maybe Int) (IntMap Int))
-> Getting (Maybe Int) (Matrix a) (Maybe Int)
forall b c a. (b -> c) -> (a -> b) -> a -> c
.Index (IntMap Int)
-> Lens' (IntMap Int) (Maybe (IxValue (IntMap Int)))
forall m. At m => Index m -> Lens' m (Maybe (IxValue m))
at Int
Index (IntMap Int)
i
      jth :: Maybe Int
jth = Matrix a
matMatrix a -> Getting (Maybe Int) (Matrix a) (Maybe Int) -> Maybe Int
forall s a. s -> Getting a s a -> a
^.Direction -> Lens' (Matrix a) (IntMap Int)
forall a. Direction -> Lens' (Matrix a) (IntMap Int)
startL Direction
dir((IntMap Int -> Const (Maybe Int) (IntMap Int))
 -> Matrix a -> Const (Maybe Int) (Matrix a))
-> ((Maybe Int -> Const (Maybe Int) (Maybe Int))
    -> IntMap Int -> Const (Maybe Int) (IntMap Int))
-> Getting (Maybe Int) (Matrix a) (Maybe Int)
forall b c a. (b -> c) -> (a -> b) -> a -> c
.Index (IntMap Int)
-> Lens' (IntMap Int) (Maybe (IxValue (IntMap Int)))
forall m. At m => Index m -> Lens' m (Maybe (IxValue m))
at Int
Index (IntMap Int)
j
  in  Matrix a
mat Matrix a -> (Matrix a -> Matrix a) -> Matrix a
forall a b. a -> (a -> b) -> b
& Direction -> Lens' (Matrix a) (IntMap Int)
forall a. Direction -> Lens' (Matrix a) (IntMap Int)
startL Direction
dir   ((IntMap Int -> Identity (IntMap Int))
 -> Matrix a -> Identity (Matrix a))
-> (IntMap Int -> IntMap Int) -> Matrix a -> Matrix a
forall s t a b. ASetter s t a b -> (a -> b) -> s -> t
%~ (Maybe Int -> Maybe Int) -> Int -> IntMap Int -> IntMap Int
forall a. (Maybe a -> Maybe a) -> Int -> IntMap a -> IntMap a
alter (Maybe Int -> Maybe Int -> Maybe Int
forall a b. a -> b -> a
const Maybe Int
jth) Int
i (IntMap Int -> IntMap Int)
-> (IntMap Int -> IntMap Int) -> IntMap Int -> IntMap Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Maybe Int -> Maybe Int) -> Int -> IntMap Int -> IntMap Int
forall a. (Maybe a -> Maybe a) -> Int -> IntMap a -> IntMap a
alter (Maybe Int -> Maybe Int -> Maybe Int
forall a b. a -> b -> a
const Maybe Int
ith) Int
j
          Matrix a -> (Matrix a -> Matrix a) -> Matrix a
forall a b. a -> (a -> b) -> b
& (Vector (Entry a) -> Identity (Vector (Entry a)))
-> Matrix a -> Identity (Matrix a)
forall a a.
Lens (Matrix a) (Matrix a) (Vector (Entry a)) (Vector (Entry a))
coefficients ((Vector (Entry a) -> Identity (Vector (Entry a)))
 -> Matrix a -> Identity (Matrix a))
-> (Vector (Entry a) -> Vector (Entry a)) -> Matrix a -> Matrix a
forall s t a b. ASetter s t a b -> (a -> b) -> s -> t
%~ Maybe Int -> Vector (Entry a) -> Vector (Entry a)
go Maybe Int
ith (Vector (Entry a) -> Vector (Entry a))
-> (Vector (Entry a) -> Vector (Entry a))
-> Vector (Entry a)
-> Vector (Entry a)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Maybe Int -> Vector (Entry a) -> Vector (Entry a)
go Maybe Int
jth
  where
    go :: Maybe Int -> Vector (Entry a) -> Vector (Entry a)
go Maybe Int
Nothing Vector (Entry a)
v = Vector (Entry a)
v
    go (Just Int
k) Vector (Entry a)
vec =
      let !cur :: Entry a
cur = Vector (Entry a)
vec Vector (Entry a) -> Int -> Entry a
forall a. Vector a -> Int -> a
V.! Int
k
      in Maybe Int -> Vector (Entry a) -> Vector (Entry a)
go (Entry a
cur Entry a -> Getting (Maybe Int) (Entry a) (Maybe Int) -> Maybe Int
forall s a. s -> Getting a s a -> a
^. Direction -> Lens' (Entry a) (Maybe Int)
forall a. Direction -> Lens' (Entry a) (Maybe Int)
nextL Direction
dir) (Vector (Entry a)
vec Vector (Entry a)
-> (Vector (Entry a) -> Vector (Entry a)) -> Vector (Entry a)
forall a b. a -> (a -> b) -> b
& Index (Vector (Entry a))
-> Traversal' (Vector (Entry a)) (IxValue (Vector (Entry a)))
forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix Int
Index (Vector (Entry a))
k ((Entry a -> Identity (Entry a))
 -> Vector (Entry a) -> Identity (Vector (Entry a)))
-> ((Int -> Identity Int) -> Entry a -> Identity (Entry a))
-> (Int -> Identity Int)
-> Vector (Entry a)
-> Identity (Vector (Entry a))
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Direction -> Lens' (Entry a) Int
forall a. Direction -> Lens' (Entry a) Int
coordL Direction
dir ((Int -> Identity Int)
 -> Vector (Entry a) -> Identity (Vector (Entry a)))
-> (Int -> Int) -> Vector (Entry a) -> Vector (Entry a)
forall s t a b. ASetter s t a b -> (a -> b) -> s -> t
%~ Int -> Int
change)
    change :: Int -> Int
change Int
k | Int
k Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
i = Int
j
             | Int
k Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j = Int
i
             | Bool
otherwise = Int
k

scaleDir :: (DecidableZero a, Multiplicative a) => Direction -> a -> Int -> Matrix a -> Matrix a
scaleDir :: Direction -> a -> Int -> Matrix a -> Matrix a
scaleDir Direction
dir a
a Int
i Matrix a
mat
  | Bool
otherwise = (a -> a) -> Direction -> Int -> Matrix a -> Matrix a
forall a. (a -> a) -> Direction -> Int -> Matrix a -> Matrix a
mapDir (a -> a -> a
forall r. Multiplicative r => r -> r -> r
*a
a) Direction
dir Int
i Matrix a
mat
  | a -> Bool
forall r. DecidableZero r => r -> Bool
isZero a
a  = Direction -> Int -> Matrix a -> Matrix a
forall a. Direction -> Int -> Matrix a -> Matrix a
clearDir Direction
dir Int
i Matrix a
mat

clearAt :: Int -> Matrix a -> Matrix a
clearAt :: Int -> Matrix a -> Matrix a
clearAt Int
k Matrix a
mat = Matrix a
mat Matrix a -> (Matrix a -> Matrix a) -> Matrix a
forall a b. a -> (a -> b) -> b
& (Vector (Entry a) -> Identity (Vector (Entry a)))
-> Matrix a -> Identity (Matrix a)
forall a a.
Lens (Matrix a) (Matrix a) (Vector (Entry a)) (Vector (Entry a))
coefficients ((Vector (Entry a) -> Identity (Vector (Entry a)))
 -> Matrix a -> Identity (Matrix a))
-> (Vector (Entry a) -> Vector (Entry a)) -> Matrix a -> Matrix a
forall s t a b. ASetter s t a b -> (a -> b) -> s -> t
%~ Vector (Entry a) -> Vector (Entry a)
go
                    Matrix a -> (Matrix a -> Matrix a) -> Matrix a
forall a b. a -> (a -> b) -> b
& Direction -> Matrix a -> Matrix a
forwardStart Direction
Column
                    Matrix a -> (Matrix a -> Matrix a) -> Matrix a
forall a b. a -> (a -> b) -> b
& Direction -> Matrix a -> Matrix a
forwardStart Direction
Row
  where
    !old :: Entry a
old = ((Matrix a
mat Matrix a
-> Getting (Vector (Entry a)) (Matrix a) (Vector (Entry a))
-> Vector (Entry a)
forall s a. s -> Getting a s a -> a
^. Getting (Vector (Entry a)) (Matrix a) (Vector (Entry a))
forall a a.
Lens (Matrix a) (Matrix a) (Vector (Entry a)) (Vector (Entry a))
coefficients) Vector (Entry a) -> Int -> Entry a
forall a. Vector a -> Int -> a
V.! Int
k)
           Entry a -> (Entry a -> Entry a) -> Entry a
forall a b. a -> (a -> b) -> b
& (Maybe Int -> Identity (Maybe Int))
-> Entry a -> Identity (Entry a)
forall a. Lens' (Entry a) (Maybe Int)
colNext((Maybe Int -> Identity (Maybe Int))
 -> Entry a -> Identity (Entry a))
-> ((Int -> Identity Int) -> Maybe Int -> Identity (Maybe Int))
-> (Int -> Identity Int)
-> Entry a
-> Identity (Entry a)
forall b c a. (b -> c) -> (a -> b) -> a -> c
.(Int -> Identity Int) -> Maybe Int -> Identity (Maybe Int)
forall a b. Prism (Maybe a) (Maybe b) a b
_Just ((Int -> Identity Int) -> Entry a -> Identity (Entry a))
-> (Int -> Int) -> Entry a -> Entry a
forall s t a b. ASetter s t a b -> (a -> b) -> s -> t
%~ Int -> Int
shifter
           Entry a -> (Entry a -> Entry a) -> Entry a
forall a b. a -> (a -> b) -> b
& (Maybe Int -> Identity (Maybe Int))
-> Entry a -> Identity (Entry a)
forall a. Lens' (Entry a) (Maybe Int)
rowNext((Maybe Int -> Identity (Maybe Int))
 -> Entry a -> Identity (Entry a))
-> ((Int -> Identity Int) -> Maybe Int -> Identity (Maybe Int))
-> (Int -> Identity Int)
-> Entry a
-> Identity (Entry a)
forall b c a. (b -> c) -> (a -> b) -> a -> c
.(Int -> Identity Int) -> Maybe Int -> Identity (Maybe Int)
forall a b. Prism (Maybe a) (Maybe b) a b
_Just ((Int -> Identity Int) -> Entry a -> Identity (Entry a))
-> (Int -> Int) -> Entry a -> Entry a
forall s t a b. ASetter s t a b -> (a -> b) -> s -> t
%~ Int -> Int
shifter
    forwardStart :: Direction -> Matrix a -> Matrix a
forwardStart Direction
dir =
      let l :: Int
l = (Entry a
old Entry a -> Getting Int (Entry a) Int -> Int
forall s a. s -> Getting a s a -> a
^. Direction -> Lens' (Entry a) Int
forall a. Direction -> Lens' (Entry a) Int
coordL Direction
dir)
      in Direction -> Lens' (Matrix a) (IntMap Int)
forall a. Direction -> Lens' (Matrix a) (IntMap Int)
startL Direction
dir ((IntMap Int -> Identity (IntMap Int))
 -> Matrix a -> Identity (Matrix a))
-> (IntMap Int -> IntMap Int) -> Matrix a -> Matrix a
forall s t a b. ASetter s t a b -> (a -> b) -> s -> t
%~ (Int -> Int -> Maybe Int) -> IntMap Int -> IntMap Int
forall a b. (Int -> a -> Maybe b) -> IntMap a -> IntMap b
mapMaybeWithKey
                       (\Int
d Int
v -> if Int
d Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
l Bool -> Bool -> Bool
&& Int
v Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
k
                                then Entry a
old Entry a -> Getting (Maybe Int) (Entry a) (Maybe Int) -> Maybe Int
forall s a. s -> Getting a s a -> a
^. Direction -> Lens' (Entry a) (Maybe Int)
forall a. Direction -> Lens' (Entry a) (Maybe Int)
nextL Direction
dir
                                else Int -> Maybe Int
forall a. a -> Maybe a
Just (Int -> Maybe Int) -> Int -> Maybe Int
forall a b. (a -> b) -> a -> b
$ Int -> Int
shifter Int
v)
    shiftDir :: Direction -> Entry a -> Entry a
shiftDir Direction
sel = Direction -> Lens' (Entry a) (Maybe Int)
forall a. Direction -> Lens' (Entry a) (Maybe Int)
nextL Direction
sel ((Maybe Int -> Identity (Maybe Int))
 -> Entry a -> Identity (Entry a))
-> (Maybe Int -> Maybe Int) -> Entry a -> Entry a
forall s t a b. ASetter s t a b -> (a -> b) -> s -> t
%~ \case
      Maybe Int
Nothing -> Maybe Int
forall a. Maybe a
Nothing
      Just Int
l ->
        if Int
l Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
k
        then Entry a
old Entry a -> Getting (Maybe Int) (Entry a) (Maybe Int) -> Maybe Int
forall s a. s -> Getting a s a -> a
^. Direction -> Lens' (Entry a) (Maybe Int)
forall a. Direction -> Lens' (Entry a) (Maybe Int)
nextL Direction
sel
        else Int -> Maybe Int
forall a. a -> Maybe a
Just (Int -> Maybe Int) -> Int -> Maybe Int
forall a b. (a -> b) -> a -> b
$ Int -> Int
shifter Int
l
    shifter :: Int -> Int
shifter Int
n = if Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
k then Int
n else Int
n Int -> Int -> Int
forall r. Group r => r -> r -> r
- Int
1
    go :: Vector (Entry a) -> Vector (Entry a)
go Vector (Entry a)
vs = Int -> (Int -> Entry a) -> Vector (Entry a)
forall a. Int -> (Int -> a) -> Vector a
generate (Vector (Entry a) -> Int
forall a. Vector a -> Int
V.length Vector (Entry a)
vs Int -> Int -> Int
forall r. Group r => r -> r -> r
- Int
1) ((Int -> Entry a) -> Vector (Entry a))
-> (Int -> Entry a) -> Vector (Entry a)
forall a b. (a -> b) -> a -> b
$ \Int
n ->
      Vector (Entry a)
vs Vector (Entry a) -> Int -> Entry a
forall a. Vector a -> Int -> a
V.! (if Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
k then Int
n else Int
n Int -> Int -> Int
forall r. Additive r => r -> r -> r
+ Int
1) Entry a -> (Entry a -> Entry a) -> Entry a
forall a b. a -> (a -> b) -> b
& Direction -> Entry a -> Entry a
shiftDir Direction
Column Entry a -> (Entry a -> Entry a) -> Entry a
forall a b. a -> (a -> b) -> b
& Direction -> Entry a -> Entry a
shiftDir Direction
Row

clearDir :: Direction -> Int -> Matrix a -> Matrix a
clearDir :: Direction -> Int -> Matrix a -> Matrix a
clearDir Direction
dir Int
i Matrix a
mat = (Matrix a -> Int -> Matrix a) -> Matrix a -> [Int] -> Matrix a
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl ((Int -> Matrix a -> Matrix a) -> Matrix a -> Int -> Matrix a
forall a b c. (a -> b -> c) -> b -> a -> c
flip Int -> Matrix a -> Matrix a
forall a. Int -> Matrix a -> Matrix a
clearAt) Matrix a
mat ([Int] -> Matrix a) -> [Int] -> Matrix a
forall a b. (a -> b) -> a -> b
$ [Int] -> [Int]
forall a. Ord a => [a] -> [a]
sort ([Int] -> [Int]) -> [Int] -> [Int]
forall a b. (a -> b) -> a -> b
$ (Maybe (Int, Entry a) -> Maybe Int)
-> [Maybe (Int, Entry a)] -> [Int]
forall a b. (a -> Maybe b) -> [a] -> [b]
mapMaybe (((Int, Entry a) -> Int) -> Maybe (Int, Entry a) -> Maybe Int
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (Int, Entry a) -> Int
forall a b. (a, b) -> a
fst) ([Maybe (Int, Entry a)] -> [Int])
-> [Maybe (Int, Entry a)] -> [Int]
forall a b. (a -> b) -> a -> b
$ Vector (Maybe (Int, Entry a)) -> [Maybe (Int, Entry a)]
forall a. Vector a -> [a]
V.toList (Vector (Maybe (Int, Entry a)) -> [Maybe (Int, Entry a)])
-> Vector (Maybe (Int, Entry a)) -> [Maybe (Int, Entry a)]
forall a b. (a -> b) -> a -> b
$ Direction -> Int -> Matrix a -> Vector (Maybe (Int, Entry a))
forall a.
Direction -> Int -> Matrix a -> Vector (Maybe (Int, Entry a))
igetDir Direction
dir Int
i Matrix a
mat

clearRow :: Int -> Matrix a -> Matrix a
clearRow :: Int -> Matrix a -> Matrix a
clearRow = Direction -> Int -> Matrix a -> Matrix a
forall a. Direction -> Int -> Matrix a -> Matrix a
clearDir Direction
Row

clearCol :: Int -> Matrix a -> Matrix a
clearCol :: Int -> Matrix a -> Matrix a
clearCol = Direction -> Int -> Matrix a -> Matrix a
forall a. Direction -> Int -> Matrix a -> Matrix a
clearDir Direction
Column

scaleRow :: (DecidableZero a, Multiplicative a) => a -> Int -> Matrix a -> Matrix a
scaleRow :: a -> Int -> Matrix a -> Matrix a
scaleRow = Direction -> a -> Int -> Matrix a -> Matrix a
forall a.
(DecidableZero a, Multiplicative a) =>
Direction -> a -> Int -> Matrix a -> Matrix a
scaleDir Direction
Row

scaleCol :: (DecidableZero a, Multiplicative a) => a -> Int -> Matrix a -> Matrix a
scaleCol :: a -> Int -> Matrix a -> Matrix a
scaleCol = Direction -> a -> Int -> Matrix a -> Matrix a
forall a.
(DecidableZero a, Multiplicative a) =>
Direction -> a -> Int -> Matrix a -> Matrix a
scaleDir Direction
Column

mapDir :: (a -> a) -> Direction
       -> Int -> Matrix a -> Matrix a
mapDir :: (a -> a) -> Direction -> Int -> Matrix a -> Matrix a
mapDir a -> a
f Direction
dir Int
i Matrix a
mat = Matrix a
-> (Matrix a -> Int -> Entry a -> Matrix a)
-> Direction
-> Int
-> Matrix a
-> Matrix a
forall b a.
b
-> (b -> Int -> Entry a -> b) -> Direction -> Int -> Matrix a -> b
traverseDir Matrix a
mat Matrix a -> Int -> Entry a -> Matrix a
trv Direction
dir Int
i Matrix a
mat
  where
    trv :: Matrix a -> Int -> Entry a -> Matrix a
trv Matrix a
m Int
k Entry a
_ = Matrix a
m Matrix a -> (Matrix a -> Matrix a) -> Matrix a
forall a b. a -> (a -> b) -> b
& (Vector (Entry a) -> Identity (Vector (Entry a)))
-> Matrix a -> Identity (Matrix a)
forall a a.
Lens (Matrix a) (Matrix a) (Vector (Entry a)) (Vector (Entry a))
coefficients ((Vector (Entry a) -> Identity (Vector (Entry a)))
 -> Matrix a -> Identity (Matrix a))
-> ((a -> Identity a)
    -> Vector (Entry a) -> Identity (Vector (Entry a)))
-> (a -> Identity a)
-> Matrix a
-> Identity (Matrix a)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Index (Vector (Entry a))
-> Traversal' (Vector (Entry a)) (IxValue (Vector (Entry a)))
forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix Int
Index (Vector (Entry a))
k ((Entry a -> Identity (Entry a))
 -> Vector (Entry a) -> Identity (Vector (Entry a)))
-> ((a -> Identity a) -> Entry a -> Identity (Entry a))
-> (a -> Identity a)
-> Vector (Entry a)
-> Identity (Vector (Entry a))
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (a -> Identity a) -> Entry a -> Identity (Entry a)
forall a a. Lens (Entry a) (Entry a) a a
value ((a -> Identity a) -> Matrix a -> Identity (Matrix a))
-> (a -> a) -> Matrix a -> Matrix a
forall s t a b. ASetter s t a b -> (a -> b) -> s -> t
%~ a -> a
f

traverseDir :: b -> (b -> Int -> Entry a -> b)
               -> Direction
               -> Int -> Matrix a -> b
traverseDir :: b
-> (b -> Int -> Entry a -> b) -> Direction -> Int -> Matrix a -> b
traverseDir b
ini b -> Int -> Entry a -> b
f Direction
dir Int
i Matrix a
mat =
  Identity b -> b
forall a. Identity a -> a
runIdentity (Identity b -> b) -> Identity b -> b
forall a b. (a -> b) -> a -> b
$  b
-> (b -> Int -> Entry a -> Identity b)
-> Direction
-> Int
-> Matrix a
-> Identity b
forall (m :: * -> *) b a.
Monad m =>
b
-> (b -> Int -> Entry a -> m b)
-> Direction
-> Int
-> Matrix a
-> m b
traverseDirM b
ini (\b
b Int
j Entry a
e -> b -> Identity b
forall (m :: * -> *) a. Monad m => a -> m a
return (b -> Identity b) -> b -> Identity b
forall a b. (a -> b) -> a -> b
$ b -> Int -> Entry a -> b
f b
b Int
j Entry a
e) Direction
dir Int
i Matrix a
mat

traverseDirM :: Monad m => b -> (b -> Int -> Entry a -> m b)
                -> Direction
                -> Int -> Matrix a -> m b
traverseDirM :: b
-> (b -> Int -> Entry a -> m b)
-> Direction
-> Int
-> Matrix a
-> m b
traverseDirM b
ini b -> Int -> Entry a -> m b
f Direction
dir Int
i Matrix a
mat = Maybe Int -> b -> m b
go (Int -> IntMap Int -> Maybe Int
forall a. Int -> IntMap a -> Maybe a
IM.lookup Int
i (Matrix a
matMatrix a
-> Getting (IntMap Int) (Matrix a) (IntMap Int) -> IntMap Int
forall s a. s -> Getting a s a -> a
^.Direction -> Lens' (Matrix a) (IntMap Int)
forall a. Direction -> Lens' (Matrix a) (IntMap Int)
startL Direction
dir)) b
ini
  where
    vec :: Vector (Entry a)
vec = Matrix a
mat Matrix a
-> Getting (Vector (Entry a)) (Matrix a) (Vector (Entry a))
-> Vector (Entry a)
forall s a. s -> Getting a s a -> a
^. Getting (Vector (Entry a)) (Matrix a) (Vector (Entry a))
forall a a.
Lens (Matrix a) (Matrix a) (Vector (Entry a)) (Vector (Entry a))
coefficients
    go :: Maybe Int -> b -> m b
go Maybe Int
Nothing  b
b = b -> m b
forall (m :: * -> *) a. Monad m => a -> m a
return b
b
    go (Just Int
k) b
b = do
      let !cur :: Entry a
cur = Vector (Entry a)
vec Vector (Entry a) -> Int -> Entry a
forall a. Vector a -> Int -> a
V.! Int
k
      Maybe Int -> b -> m b
go (Entry a
cur Entry a -> Getting (Maybe Int) (Entry a) (Maybe Int) -> Maybe Int
forall s a. s -> Getting a s a -> a
^. Direction -> Lens' (Entry a) (Maybe Int)
forall a. Direction -> Lens' (Entry a) (Maybe Int)
nextL Direction
dir) (b -> m b) -> m b -> m b
forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
=<< b -> Int -> Entry a -> m b
f b
b Int
k Entry a
cur

getDir :: forall a. Monoidal a
       => Direction -> Int -> Matrix a -> Vector a
getDir :: Direction -> Int -> Matrix a -> Vector a
getDir Direction
dir Int
i Matrix a
mat =
  (forall s. ST s (MVector s a)) -> Vector a
forall a. (forall s. ST s (MVector s a)) -> Vector a
create ((forall s. ST s (MVector s a)) -> Vector a)
-> (forall s. ST s (MVector s a)) -> Vector a
forall a b. (a -> b) -> a -> b
$ do
    MVector s a
v <- Int -> a -> ST s (MVector (PrimState (ST s)) a)
forall (m :: * -> *) a.
PrimMonad m =>
Int -> a -> m (MVector (PrimState m) a)
MV.replicate (Matrix a
mat Matrix a -> Getting Int (Matrix a) Int -> Int
forall s a. s -> Getting a s a -> a
^. Direction -> Lens' (Matrix a) Int
forall a. Direction -> Lens' (Matrix a) Int
lenL Direction
dir) (a
forall m. Monoidal m => m
zero :: a)
    ()
-> (() -> Int -> Entry a -> ST s ())
-> Direction
-> Int
-> Matrix a
-> ST s ()
forall (m :: * -> *) b a.
Monad m =>
b
-> (b -> Int -> Entry a -> m b)
-> Direction
-> Int
-> Matrix a
-> m b
traverseDirM () (MVector s a -> () -> Int -> Entry a -> ST s ()
forall s. MVector s a -> () -> Int -> Entry a -> ST s ()
trav MVector s a
v) Direction
dir Int
i Matrix a
mat
    MVector s a -> ST s (MVector s a)
forall (m :: * -> *) a. Monad m => a -> m a
return MVector s a
v
  where
    trav :: forall s. MV.MVector s a -> () -> Int -> Entry a -> ST s ()
    trav :: MVector s a -> () -> Int -> Entry a -> ST s ()
trav MVector s a
v ()
_ Int
_ Entry a
ent = MVector (PrimState (ST s)) a -> Int -> a -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
MV.write MVector s a
MVector (PrimState (ST s)) a
v (Entry a
ent Entry a -> Getting Int (Entry a) Int -> Int
forall s a. s -> Getting a s a -> a
^. Direction -> Lens' (Entry a) Int
forall a. Direction -> Lens' (Entry a) Int
nthL Direction
dir) (Entry a
ent Entry a -> Getting a (Entry a) a -> a
forall s a. s -> Getting a s a -> a
^. Getting a (Entry a) a
forall a a. Lens (Entry a) (Entry a) a a
value)

igetDir :: forall a. Direction -> Int -> Matrix a -> Vector (Maybe (Int, Entry a))
igetDir :: Direction -> Int -> Matrix a -> Vector (Maybe (Int, Entry a))
igetDir Direction
dir Int
i Matrix a
mat =
  (forall s. ST s (MVector s (Maybe (Int, Entry a))))
-> Vector (Maybe (Int, Entry a))
forall a. (forall s. ST s (MVector s a)) -> Vector a
create ((forall s. ST s (MVector s (Maybe (Int, Entry a))))
 -> Vector (Maybe (Int, Entry a)))
-> (forall s. ST s (MVector s (Maybe (Int, Entry a))))
-> Vector (Maybe (Int, Entry a))
forall a b. (a -> b) -> a -> b
$ do
    MVector s (Maybe (Int, Entry a))
v <- Int
-> Maybe (Int, Entry a)
-> ST s (MVector (PrimState (ST s)) (Maybe (Int, Entry a)))
forall (m :: * -> *) a.
PrimMonad m =>
Int -> a -> m (MVector (PrimState m) a)
MV.replicate (Matrix a
mat Matrix a -> Getting Int (Matrix a) Int -> Int
forall s a. s -> Getting a s a -> a
^. Direction -> Lens' (Matrix a) Int
forall a. Direction -> Lens' (Matrix a) Int
lenL Direction
dir) Maybe (Int, Entry a)
forall a. Maybe a
Nothing
    ()
-> (() -> Int -> Entry a -> ST s ())
-> Direction
-> Int
-> Matrix a
-> ST s ()
forall (m :: * -> *) b a.
Monad m =>
b
-> (b -> Int -> Entry a -> m b)
-> Direction
-> Int
-> Matrix a
-> m b
traverseDirM () (MVector s (Maybe (Int, Entry a)) -> () -> Int -> Entry a -> ST s ()
forall s.
MVector s (Maybe (Int, Entry a)) -> () -> Int -> Entry a -> ST s ()
trav MVector s (Maybe (Int, Entry a))
v) Direction
dir Int
i Matrix a
mat
    MVector s (Maybe (Int, Entry a))
-> ST s (MVector s (Maybe (Int, Entry a)))
forall (m :: * -> *) a. Monad m => a -> m a
return MVector s (Maybe (Int, Entry a))
v
  where
    trav :: forall s. MV.MVector s (Maybe (Int, Entry a)) -> () -> Int -> Entry a -> ST s ()
    trav :: MVector s (Maybe (Int, Entry a)) -> () -> Int -> Entry a -> ST s ()
trav MVector s (Maybe (Int, Entry a))
v ()
_ Int
k Entry a
ent = MVector (PrimState (ST s)) (Maybe (Int, Entry a))
-> Int -> Maybe (Int, Entry a) -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
MV.write MVector s (Maybe (Int, Entry a))
MVector (PrimState (ST s)) (Maybe (Int, Entry a))
v (Entry a
ent Entry a -> Getting Int (Entry a) Int -> Int
forall s a. s -> Getting a s a -> a
^. Direction -> Lens' (Entry a) Int
forall a. Direction -> Lens' (Entry a) Int
nthL Direction
dir) ((Int, Entry a) -> Maybe (Int, Entry a)
forall a. a -> Maybe a
Just (Int
k, Entry a
ent))

getRow :: Monoidal a => Int -> Matrix a -> Vector a
getRow :: Int -> Matrix a -> Vector a
getRow = Direction -> Int -> Matrix a -> Vector a
forall a. Monoidal a => Direction -> Int -> Matrix a -> Vector a
getDir Direction
Row

getCol :: Monoidal a => Int -> Matrix a -> Vector a
getCol :: Int -> Matrix a -> Vector a
getCol = Direction -> Int -> Matrix a -> Vector a
forall a. Monoidal a => Direction -> Int -> Matrix a -> Vector a
getDir Direction
Column

data Direction = Row | Column deriving (ReadPrec [Direction]
ReadPrec Direction
Int -> ReadS Direction
ReadS [Direction]
(Int -> ReadS Direction)
-> ReadS [Direction]
-> ReadPrec Direction
-> ReadPrec [Direction]
-> Read Direction
forall a.
(Int -> ReadS a)
-> ReadS [a] -> ReadPrec a -> ReadPrec [a] -> Read a
readListPrec :: ReadPrec [Direction]
$creadListPrec :: ReadPrec [Direction]
readPrec :: ReadPrec Direction
$creadPrec :: ReadPrec Direction
readList :: ReadS [Direction]
$creadList :: ReadS [Direction]
readsPrec :: Int -> ReadS Direction
$creadsPrec :: Int -> ReadS Direction
Read, Int -> Direction -> ShowS
[Direction] -> ShowS
Direction -> String
(Int -> Direction -> ShowS)
-> (Direction -> String)
-> ([Direction] -> ShowS)
-> Show Direction
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [Direction] -> ShowS
$cshowList :: [Direction] -> ShowS
show :: Direction -> String
$cshow :: Direction -> String
showsPrec :: Int -> Direction -> ShowS
$cshowsPrec :: Int -> Direction -> ShowS
Show, Direction -> Direction -> Bool
(Direction -> Direction -> Bool)
-> (Direction -> Direction -> Bool) -> Eq Direction
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: Direction -> Direction -> Bool
$c/= :: Direction -> Direction -> Bool
== :: Direction -> Direction -> Bool
$c== :: Direction -> Direction -> Bool
Eq, Eq Direction
Eq Direction
-> (Direction -> Direction -> Ordering)
-> (Direction -> Direction -> Bool)
-> (Direction -> Direction -> Bool)
-> (Direction -> Direction -> Bool)
-> (Direction -> Direction -> Bool)
-> (Direction -> Direction -> Direction)
-> (Direction -> Direction -> Direction)
-> Ord Direction
Direction -> Direction -> Bool
Direction -> Direction -> Ordering
Direction -> Direction -> Direction
forall a.
Eq a
-> (a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
min :: Direction -> Direction -> Direction
$cmin :: Direction -> Direction -> Direction
max :: Direction -> Direction -> Direction
$cmax :: Direction -> Direction -> Direction
>= :: Direction -> Direction -> Bool
$c>= :: Direction -> Direction -> Bool
> :: Direction -> Direction -> Bool
$c> :: Direction -> Direction -> Bool
<= :: Direction -> Direction -> Bool
$c<= :: Direction -> Direction -> Bool
< :: Direction -> Direction -> Bool
$c< :: Direction -> Direction -> Bool
compare :: Direction -> Direction -> Ordering
$ccompare :: Direction -> Direction -> Ordering
$cp1Ord :: Eq Direction
Ord)

lenL, countL :: Direction -> Lens' (Matrix a) Int
lenL :: Direction -> Lens' (Matrix a) Int
lenL Direction
Row = (Int -> f Int) -> Matrix a -> f (Matrix a)
forall a. Lens' (Matrix a) Int
width
lenL Direction
Column = (Int -> f Int) -> Matrix a -> f (Matrix a)
forall a. Lens' (Matrix a) Int
height
countL :: Direction -> Lens' (Matrix a) Int
countL Direction
Row = (Int -> f Int) -> Matrix a -> f (Matrix a)
forall a. Lens' (Matrix a) Int
height
countL Direction
Column = (Int -> f Int) -> Matrix a -> f (Matrix a)
forall a. Lens' (Matrix a) Int
width

nthL, coordL :: Direction -> Lens' (Entry a) Int
coordL :: Direction -> Lens' (Entry a) Int
coordL Direction
Row = ((Int, Int) -> f (Int, Int)) -> Entry a -> f (Entry a)
forall a. Lens' (Entry a) (Int, Int)
idx (((Int, Int) -> f (Int, Int)) -> Entry a -> f (Entry a))
-> ((Int -> f Int) -> (Int, Int) -> f (Int, Int))
-> (Int -> f Int)
-> Entry a
-> f (Entry a)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Int -> f Int) -> (Int, Int) -> f (Int, Int)
forall s t a b. Field1 s t a b => Lens s t a b
_1
coordL Direction
Column = ((Int, Int) -> f (Int, Int)) -> Entry a -> f (Entry a)
forall a. Lens' (Entry a) (Int, Int)
idx (((Int, Int) -> f (Int, Int)) -> Entry a -> f (Entry a))
-> ((Int -> f Int) -> (Int, Int) -> f (Int, Int))
-> (Int -> f Int)
-> Entry a
-> f (Entry a)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Int -> f Int) -> (Int, Int) -> f (Int, Int)
forall s t a b. Field2 s t a b => Lens s t a b
_2

nthL :: Direction -> Lens' (Entry a) Int
nthL Direction
Row = ((Int, Int) -> f (Int, Int)) -> Entry a -> f (Entry a)
forall a. Lens' (Entry a) (Int, Int)
idx (((Int, Int) -> f (Int, Int)) -> Entry a -> f (Entry a))
-> ((Int -> f Int) -> (Int, Int) -> f (Int, Int))
-> (Int -> f Int)
-> Entry a
-> f (Entry a)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Int -> f Int) -> (Int, Int) -> f (Int, Int)
forall s t a b. Field2 s t a b => Lens s t a b
_2
nthL Direction
Column = ((Int, Int) -> f (Int, Int)) -> Entry a -> f (Entry a)
forall a. Lens' (Entry a) (Int, Int)
idx (((Int, Int) -> f (Int, Int)) -> Entry a -> f (Entry a))
-> ((Int -> f Int) -> (Int, Int) -> f (Int, Int))
-> (Int -> f Int)
-> Entry a
-> f (Entry a)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Int -> f Int) -> (Int, Int) -> f (Int, Int)
forall s t a b. Field1 s t a b => Lens s t a b
_1

startL :: Direction -> Lens' (Matrix a) (IntMap Int)
startL :: Direction -> Lens' (Matrix a) (IntMap Int)
startL Direction
Row = (IntMap Int -> f (IntMap Int)) -> Matrix a -> f (Matrix a)
forall a. Lens' (Matrix a) (IntMap Int)
rowStart
startL Direction
Column = (IntMap Int -> f (IntMap Int)) -> Matrix a -> f (Matrix a)
forall a. Lens' (Matrix a) (IntMap Int)
colStart

nextL :: Direction -> Lens' (Entry a) (Maybe Int)
nextL :: Direction -> Lens' (Entry a) (Maybe Int)
nextL Direction
Row    = (Maybe Int -> f (Maybe Int)) -> Entry a -> f (Entry a)
forall a. Lens' (Entry a) (Maybe Int)
rowNext
nextL Direction
Column = (Maybe Int -> f (Maybe Int)) -> Entry a -> f (Entry a)
forall a. Lens' (Entry a) (Maybe Int)
colNext

addDir :: forall a. (DecidableZero a)
       => Direction -> Vector a -> Int -> Matrix a -> Matrix a
addDir :: Direction -> Vector a -> Int -> Matrix a -> Matrix a
addDir Direction
dir Vector a
vec Int
i Matrix a
mat = (forall s. ST s (Matrix a)) -> Matrix a
forall a. (forall s. ST s a) -> a
runST ((forall s. ST s (Matrix a)) -> Matrix a)
-> (forall s. ST s (Matrix a)) -> Matrix a
forall a b. (a -> b) -> a -> b
$ do
    MVector s (Entry a)
mv <- Vector (Entry a) -> ST s (MVector (PrimState (ST s)) (Entry a))
forall (m :: * -> *) a.
PrimMonad m =>
Vector a -> m (MVector (PrimState m) a)
thaw (Vector (Entry a) -> ST s (MVector (PrimState (ST s)) (Entry a)))
-> Vector (Entry a) -> ST s (MVector (PrimState (ST s)) (Entry a))
forall a b. (a -> b) -> a -> b
$ Matrix a
mat Matrix a
-> Getting (Vector (Entry a)) (Matrix a) (Vector (Entry a))
-> Vector (Entry a)
forall s a. s -> Getting a s a -> a
^. Getting (Vector (Entry a)) (Matrix a) (Vector (Entry a))
forall a a.
Lens (Matrix a) (Matrix a) (Vector (Entry a)) (Vector (Entry a))
coefficients
    let n :: Int
n = MVector s (Entry a) -> Int
forall s a. MVector s a -> Int
MV.length MVector s (Entry a)
mv
        upd :: (IntMap a, [Int]) -> Int -> Entry a -> ST s (IntMap a, [Int])
upd (IntMap a
dic, [Int]
del) Int
k Entry a
e = do
          let v' :: a
v' = Entry a
e Entry a -> Getting a (Entry a) a -> a
forall s a. s -> Getting a s a -> a
^. Getting a (Entry a) a
forall a a. Lens (Entry a) (Entry a) a a
value a -> a -> a
forall r. Additive r => r -> r -> r
+ a -> Int -> IntMap a -> a
forall a. a -> Int -> IntMap a -> a
IM.findWithDefault a
forall m. Monoidal m => m
zero (Entry a
e Entry a -> Getting Int (Entry a) Int -> Int
forall s a. s -> Getting a s a -> a
^. Direction -> Lens' (Entry a) Int
forall a. Direction -> Lens' (Entry a) Int
nthL Direction
dir) IntMap a
mp
          [Int]
d' <- if a -> Bool
forall r. DecidableZero r => r -> Bool
isZero a
v'
                then [Int] -> ST s [Int]
forall (m :: * -> *) a. Monad m => a -> m a
return ([Int] -> ST s [Int]) -> [Int] -> ST s [Int]
forall a b. (a -> b) -> a -> b
$ Int
k Int -> [Int] -> [Int]
forall a. a -> [a] -> [a]
: [Int]
del
                else MVector (PrimState (ST s)) (Entry a) -> Int -> Entry a -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
MV.write MVector s (Entry a)
MVector (PrimState (ST s)) (Entry a)
mv Int
k (Entry a
e Entry a -> (Entry a -> Entry a) -> Entry a
forall a b. a -> (a -> b) -> b
& (a -> Identity a) -> Entry a -> Identity (Entry a)
forall a a. Lens (Entry a) (Entry a) a a
value ((a -> Identity a) -> Entry a -> Identity (Entry a))
-> a -> Entry a -> Entry a
forall s t a b. ASetter s t a b -> b -> s -> t
.~ a
v') ST s () -> ST s [Int] -> ST s [Int]
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> [Int] -> ST s [Int]
forall (m :: * -> *) a. Monad m => a -> m a
return [Int]
del
          (IntMap a, [Int]) -> ST s (IntMap a, [Int])
forall (m :: * -> *) a. Monad m => a -> m a
return (Int -> IntMap a -> IntMap a
forall a. Int -> IntMap a -> IntMap a
IM.delete (Entry a
e Entry a -> Getting Int (Entry a) Int -> Int
forall s a. s -> Getting a s a -> a
^. Direction -> Lens' (Entry a) Int
forall a. Direction -> Lens' (Entry a) Int
nthL Direction
dir) IntMap a
dic, [Int]
d')
    (IntMap a
rest, [Int]
dels) <- (IntMap a, [Int])
-> ((IntMap a, [Int]) -> Int -> Entry a -> ST s (IntMap a, [Int]))
-> Direction
-> Int
-> Matrix a
-> ST s (IntMap a, [Int])
forall (m :: * -> *) b a.
Monad m =>
b
-> (b -> Int -> Entry a -> m b)
-> Direction
-> Int
-> Matrix a
-> m b
traverseDirM (IntMap a
mp, []) (IntMap a, [Int]) -> Int -> Entry a -> ST s (IntMap a, [Int])
upd Direction
dir Int
i Matrix a
mat
    MVector s (Entry a)
mv' <- if IntMap a -> Bool
forall a. IntMap a -> Bool
IM.null IntMap a
rest
           then MVector s (Entry a) -> ST s (MVector s (Entry a))
forall (m :: * -> *) a. Monad m => a -> m a
return MVector s (Entry a)
mv
           else MVector (PrimState (ST s)) (Entry a)
-> Int -> ST s (MVector (PrimState (ST s)) (Entry a))
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> m (MVector (PrimState m) a)
grow MVector s (Entry a)
MVector (PrimState (ST s)) (Entry a)
mv (IntMap a -> Int
forall a. IntMap a -> Int
IM.size IntMap a
rest)
    let app :: Int
-> (Maybe Int, Int, IntMap Int)
-> a
-> ST s (Maybe Int, Int, IntMap Int)
app Int
j (Maybe Int
p, Int
k, IntMap Int
opo) a
v = do
          let preOpo :: Maybe Int
preOpo = Matrix a
mat Matrix a -> Getting (Maybe Int) (Matrix a) (Maybe Int) -> Maybe Int
forall s a. s -> Getting a s a -> a
^. Direction -> Lens' (Matrix a) (IntMap Int)
forall a. Direction -> Lens' (Matrix a) (IntMap Int)
startL (Direction -> Direction
perp Direction
dir) ((IntMap Int -> Const (Maybe Int) (IntMap Int))
 -> Matrix a -> Const (Maybe Int) (Matrix a))
-> ((Maybe Int -> Const (Maybe Int) (Maybe Int))
    -> IntMap Int -> Const (Maybe Int) (IntMap Int))
-> Getting (Maybe Int) (Matrix a) (Maybe Int)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Index (IntMap Int)
-> Lens' (IntMap Int) (Maybe (IxValue (IntMap Int)))
forall m. At m => Index m -> Lens' m (Maybe (IxValue m))
at Int
Index (IntMap Int)
j
          MVector (PrimState (ST s)) (Entry a) -> Int -> Entry a -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
MV.write MVector s (Entry a)
MVector (PrimState (ST s)) (Entry a)
mv' Int
k (Entry a -> ST s ()) -> Entry a -> ST s ()
forall a b. (a -> b) -> a -> b
$ a -> Entry a
forall a. a -> Entry a
newEntry a
v
                         Entry a -> (Entry a -> Entry a) -> Entry a
forall a b. a -> (a -> b) -> b
& Direction -> Lens' (Entry a) (Maybe Int)
forall a. Direction -> Lens' (Entry a) (Maybe Int)
nextL Direction
dir ((Maybe Int -> Identity (Maybe Int))
 -> Entry a -> Identity (Entry a))
-> Maybe Int -> Entry a -> Entry a
forall s t a b. ASetter s t a b -> b -> s -> t
.~ Maybe Int
p
                         Entry a -> (Entry a -> Entry a) -> Entry a
forall a b. a -> (a -> b) -> b
& Direction -> Lens' (Entry a) (Maybe Int)
forall a. Direction -> Lens' (Entry a) (Maybe Int)
nextL (Direction -> Direction
perp Direction
dir) ((Maybe Int -> Identity (Maybe Int))
 -> Entry a -> Identity (Entry a))
-> Maybe Int -> Entry a -> Entry a
forall s t a b. ASetter s t a b -> b -> s -> t
.~ Maybe Int
preOpo
                         Entry a -> (Entry a -> Entry a) -> Entry a
forall a b. a -> (a -> b) -> b
& Direction -> Lens' (Entry a) Int
forall a. Direction -> Lens' (Entry a) Int
coordL Direction
dir ((Int -> Identity Int) -> Entry a -> Identity (Entry a))
-> Int -> Entry a -> Entry a
forall s t a b. ASetter s t a b -> b -> s -> t
.~ Int
i
                         Entry a -> (Entry a -> Entry a) -> Entry a
forall a b. a -> (a -> b) -> b
& Direction -> Lens' (Entry a) Int
forall a. Direction -> Lens' (Entry a) Int
nthL Direction
dir ((Int -> Identity Int) -> Entry a -> Identity (Entry a))
-> Int -> Entry a -> Entry a
forall s t a b. ASetter s t a b -> b -> s -> t
.~ Int
j

          (Maybe Int, Int, IntMap Int) -> ST s (Maybe Int, Int, IntMap Int)
forall (m :: * -> *) a. Monad m => a -> m a
return (Int -> Maybe Int
forall a. a -> Maybe a
Just Int
k, Int
kInt -> Int -> Int
forall r. Additive r => r -> r -> r
+Int
1, (Maybe Int -> Maybe Int) -> Int -> IntMap Int -> IntMap Int
forall a. (Maybe a -> Maybe a) -> Int -> IntMap a -> IntMap a
alter (Maybe Int -> Maybe Int -> Maybe Int
forall a b. a -> b -> a
const (Maybe Int -> Maybe Int -> Maybe Int)
-> Maybe Int -> Maybe Int -> Maybe Int
forall a b. (a -> b) -> a -> b
$ Int -> Maybe Int
forall a. a -> Maybe a
Just Int
k) Int
j IntMap Int
opo)
    (Maybe Int
l, Int
_, IntMap Int
opoStart) <- (Int
 -> (Maybe Int, Int, IntMap Int)
 -> a
 -> ST s (Maybe Int, Int, IntMap Int))
-> (Maybe Int, Int, IntMap Int)
-> IntMap a
-> ST s (Maybe Int, Int, IntMap Int)
forall i (f :: * -> *) (m :: * -> *) b a.
(FoldableWithIndex i f, Monad m) =>
(i -> b -> a -> m b) -> b -> f a -> m b
ifoldlM Int
-> (Maybe Int, Int, IntMap Int)
-> a
-> ST s (Maybe Int, Int, IntMap Int)
app (Matrix a
mat Matrix a -> Getting (Maybe Int) (Matrix a) (Maybe Int) -> Maybe Int
forall s a. s -> Getting a s a -> a
^. Direction -> Lens' (Matrix a) (IntMap Int)
forall a. Direction -> Lens' (Matrix a) (IntMap Int)
startL Direction
dir ((IntMap Int -> Const (Maybe Int) (IntMap Int))
 -> Matrix a -> Const (Maybe Int) (Matrix a))
-> ((Maybe Int -> Const (Maybe Int) (Maybe Int))
    -> IntMap Int -> Const (Maybe Int) (IntMap Int))
-> Getting (Maybe Int) (Matrix a) (Maybe Int)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Index (IntMap Int)
-> Lens' (IntMap Int) (Maybe (IxValue (IntMap Int)))
forall m. At m => Index m -> Lens' m (Maybe (IxValue m))
at Int
Index (IntMap Int)
i, Int
n, Matrix a
mat Matrix a
-> Getting (IntMap Int) (Matrix a) (IntMap Int) -> IntMap Int
forall s a. s -> Getting a s a -> a
^. Direction -> Lens' (Matrix a) (IntMap Int)
forall a. Direction -> Lens' (Matrix a) (IntMap Int)
startL (Direction -> Direction
perp Direction
dir)) IntMap a
rest
    Vector (Entry a)
v' <- MVector (PrimState (ST s)) (Entry a) -> ST s (Vector (Entry a))
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> m (Vector a)
unsafeFreeze MVector s (Entry a)
MVector (PrimState (ST s)) (Entry a)
mv'
    let mat' :: Matrix a
mat' = Matrix a
mat Matrix a -> (Matrix a -> Matrix a) -> Matrix a
forall a b. a -> (a -> b) -> b
& (Vector (Entry a) -> Identity (Vector (Entry a)))
-> Matrix a -> Identity (Matrix a)
forall a a.
Lens (Matrix a) (Matrix a) (Vector (Entry a)) (Vector (Entry a))
coefficients ((Vector (Entry a) -> Identity (Vector (Entry a)))
 -> Matrix a -> Identity (Matrix a))
-> Vector (Entry a) -> Matrix a -> Matrix a
forall s t a b. ASetter s t a b -> b -> s -> t
.~ Vector (Entry a)
v'
                   Matrix a -> (Matrix a -> Matrix a) -> Matrix a
forall a b. a -> (a -> b) -> b
& Direction -> Lens' (Matrix a) (IntMap Int)
forall a. Direction -> Lens' (Matrix a) (IntMap Int)
startL Direction
dir ((IntMap Int -> Identity (IntMap Int))
 -> Matrix a -> Identity (Matrix a))
-> (IntMap Int -> IntMap Int) -> Matrix a -> Matrix a
forall s t a b. ASetter s t a b -> (a -> b) -> s -> t
%~ (Maybe Int -> Maybe Int) -> Int -> IntMap Int -> IntMap Int
forall a. (Maybe a -> Maybe a) -> Int -> IntMap a -> IntMap a
alter (Maybe Int -> Maybe Int -> Maybe Int
forall a b. a -> b -> a
const Maybe Int
l) Int
i
                   Matrix a -> (Matrix a -> Matrix a) -> Matrix a
forall a b. a -> (a -> b) -> b
& Direction -> Lens' (Matrix a) (IntMap Int)
forall a. Direction -> Lens' (Matrix a) (IntMap Int)
startL (Direction -> Direction
perp Direction
dir) ((IntMap Int -> Identity (IntMap Int))
 -> Matrix a -> Identity (Matrix a))
-> IntMap Int -> Matrix a -> Matrix a
forall s t a b. ASetter s t a b -> b -> s -> t
.~ IntMap Int
opoStart
    Matrix a -> ST s (Matrix a)
forall (m :: * -> *) a. Monad m => a -> m a
return (Matrix a -> ST s (Matrix a)) -> Matrix a -> ST s (Matrix a)
forall a b. (a -> b) -> a -> b
$ (Int -> Matrix a -> Matrix a) -> Matrix a -> [Int] -> Matrix a
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr Int -> Matrix a -> Matrix a
forall a. Int -> Matrix a -> Matrix a
clearAt Matrix a
mat' [Int]
dels
  where
    mp :: IntMap a
    mp :: IntMap a
mp = (Int -> a -> IntMap a -> IntMap a)
-> IntMap a -> Vector a -> IntMap a
forall a b. (Int -> a -> b -> b) -> b -> Vector a -> b
V.ifoldr (\Int
k a
v IntMap a
d -> if a -> Bool
forall r. DecidableZero r => r -> Bool
isZero a
v then IntMap a
d else Int -> a -> IntMap a -> IntMap a
forall a. Int -> a -> IntMap a -> IntMap a
IM.insert Int
k a
v IntMap a
d) IntMap a
forall a. IntMap a
IM.empty Vector a
vec

perp :: Direction -> Direction
perp :: Direction -> Direction
perp Direction
Row = Direction
Column
perp Direction
Column = Direction
Row

addRow :: (DecidableZero a) => Vector a -> Int -> Matrix a -> Matrix a
addRow :: Vector a -> Int -> Matrix a -> Matrix a
addRow = Direction -> Vector a -> Int -> Matrix a -> Matrix a
forall a.
DecidableZero a =>
Direction -> Vector a -> Int -> Matrix a -> Matrix a
addDir Direction
Row

addCol :: (DecidableZero a) => Vector a -> Int -> Matrix a -> Matrix a
addCol :: Vector a -> Int -> Matrix a -> Matrix a
addCol = Direction -> Vector a -> Int -> Matrix a -> Matrix a
forall a.
DecidableZero a =>
Direction -> Vector a -> Int -> Matrix a -> Matrix a
addDir Direction
Column

inBound :: (Int, Int) -> Matrix a -> Bool
inBound :: (Int, Int) -> Matrix a -> Bool
inBound (Int
i, Int
j) Matrix a
mat = Int
0 Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
i Bool -> Bool -> Bool
&& Int
i Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Matrix a
mat Matrix a -> Getting Int (Matrix a) Int -> Int
forall s a. s -> Getting a s a -> a
^. Getting Int (Matrix a) Int
forall a. Lens' (Matrix a) Int
height Bool -> Bool -> Bool
&& Int
0 Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
j Bool -> Bool -> Bool
&& Int
j Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Matrix a
mat Matrix a -> Getting Int (Matrix a) Int -> Int
forall s a. s -> Getting a s a -> a
^. Getting Int (Matrix a) Int
forall a. Lens' (Matrix a) Int
width

index :: Monoidal a => IM.Key -> Int -> Matrix a -> Maybe a
index :: Int -> Int -> Matrix a -> Maybe a
index Int
i Int
j Matrix a
mat
  | Bool -> Bool
not (Bool -> Bool) -> Bool -> Bool
forall a b. (a -> b) -> a -> b
$ (Int, Int) -> Matrix a -> Bool
forall a. (Int, Int) -> Matrix a -> Bool
inBound (Int
i, Int
j) Matrix a
mat = Maybe a
forall a. Maybe a
Nothing
  | Bool
otherwise = a -> Maybe a
forall a. a -> Maybe a
Just (a -> Maybe a) -> a -> Maybe a
forall a b. (a -> b) -> a -> b
$ Maybe Int -> a
go (Int -> IntMap Int -> Maybe Int
forall a. Int -> IntMap a -> Maybe a
IM.lookup Int
i (IntMap Int -> Maybe Int) -> IntMap Int -> Maybe Int
forall a b. (a -> b) -> a -> b
$ Matrix a
mat Matrix a
-> Getting (IntMap Int) (Matrix a) (IntMap Int) -> IntMap Int
forall s a. s -> Getting a s a -> a
^. Getting (IntMap Int) (Matrix a) (IntMap Int)
forall a. Lens' (Matrix a) (IntMap Int)
rowStart)
  where
    go :: Maybe Int -> a
go Maybe Int
Nothing  = a
forall m. Monoidal m => m
zero
    go (Just Int
k) =
      let e :: Entry a
e = (Matrix a
mat Matrix a
-> Getting (Vector (Entry a)) (Matrix a) (Vector (Entry a))
-> Vector (Entry a)
forall s a. s -> Getting a s a -> a
^. Getting (Vector (Entry a)) (Matrix a) (Vector (Entry a))
forall a a.
Lens (Matrix a) (Matrix a) (Vector (Entry a)) (Vector (Entry a))
coefficients) Vector (Entry a) -> Int -> Entry a
forall a. Vector a -> Int -> a
V.! Int
k
      in if Entry a
eEntry a -> Getting Int (Entry a) Int -> Int
forall s a. s -> Getting a s a -> a
^.((Int, Int) -> Const Int (Int, Int))
-> Entry a -> Const Int (Entry a)
forall a. Lens' (Entry a) (Int, Int)
idx(((Int, Int) -> Const Int (Int, Int))
 -> Entry a -> Const Int (Entry a))
-> ((Int -> Const Int Int) -> (Int, Int) -> Const Int (Int, Int))
-> Getting Int (Entry a) Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
.(Int -> Const Int Int) -> (Int, Int) -> Const Int (Int, Int)
forall s t a b. Field2 s t a b => Lens s t a b
_2 Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j
         then Entry a
e Entry a -> Getting a (Entry a) a -> a
forall s a. s -> Getting a s a -> a
^. Getting a (Entry a) a
forall a a. Lens (Entry a) (Entry a) a a
value
         else Maybe Int -> a
go (Entry a
eEntry a -> Getting (Maybe Int) (Entry a) (Maybe Int) -> Maybe Int
forall s a. s -> Getting a s a -> a
^.Getting (Maybe Int) (Entry a) (Maybe Int)
forall a. Lens' (Entry a) (Maybe Int)
rowNext)

(!) :: Monoidal a => Matrix a -> (Int, Int) -> a
(!) Matrix a
a (Int
i, Int
j) = Maybe a -> a
forall a. HasCallStack => Maybe a -> a
fromJust (Maybe a -> a) -> Maybe a -> a
forall a b. (a -> b) -> a -> b
$ Int -> Int -> Matrix a -> Maybe a
forall a. Monoidal a => Int -> Int -> Matrix a -> Maybe a
index Int
i Int
j Matrix a
a

combineDir :: (DecidableZero a, Multiplicative a) => Direction -> a -> Int -> Int -> Matrix a -> Matrix a
combineDir :: Direction -> a -> Int -> Int -> Matrix a -> Matrix a
combineDir Direction
dir a
alpha Int
i Int
j Matrix a
mat = Direction -> Vector a -> Int -> Matrix a -> Matrix a
forall a.
DecidableZero a =>
Direction -> Vector a -> Int -> Matrix a -> Matrix a
addDir Direction
dir ((a -> a) -> Vector a -> Vector a
forall a b. (a -> b) -> Vector a -> Vector b
V.map (a
alpha a -> a -> a
forall r. Multiplicative r => r -> r -> r
*) (Vector a -> Vector a) -> Vector a -> Vector a
forall a b. (a -> b) -> a -> b
$ Direction -> Int -> Matrix a -> Vector a
forall a. Monoidal a => Direction -> Int -> Matrix a -> Vector a
getDir Direction
dir Int
i Matrix a
mat) Int
j Matrix a
mat

combineRows :: (DecidableZero a, Multiplicative a) => a -> Int -> Int -> Matrix a -> Matrix a
combineRows :: a -> Int -> Int -> Matrix a -> Matrix a
combineRows = Direction -> a -> Int -> Int -> Matrix a -> Matrix a
forall a.
(DecidableZero a, Multiplicative a) =>
Direction -> a -> Int -> Int -> Matrix a -> Matrix a
combineDir Direction
Row

combineCols :: (DecidableZero a, Multiplicative a) => a -> Int -> Int -> Matrix a -> Matrix a
combineCols :: a -> Int -> Int -> Matrix a -> Matrix a
combineCols = Direction -> a -> Int -> Int -> Matrix a -> Matrix a
forall a.
(DecidableZero a, Multiplicative a) =>
Direction -> a -> Int -> Int -> Matrix a -> Matrix a
combineDir Direction
Column

nrows, ncols :: Matrix a -> Int
ncols :: Matrix a -> Int
ncols = Getting Int (Matrix a) Int -> Matrix a -> Int
forall s (m :: * -> *) a. MonadReader s m => Getting a s a -> m a
view Getting Int (Matrix a) Int
forall a. Lens' (Matrix a) Int
width
nrows :: Matrix a -> Int
nrows = Getting Int (Matrix a) Int -> Matrix a -> Int
forall s (m :: * -> *) a. MonadReader s m => Getting a s a -> m a
view Getting Int (Matrix a) Int
forall a. Lens' (Matrix a) Int
height

identity :: Unital a => Int -> Matrix a
identity :: Int -> Matrix a
identity Int
n =
  let idMap :: IntMap Int
idMap = [(Int, Int)] -> IntMap Int
forall a. [(Int, a)] -> IntMap a
IM.fromList [(Int
i,Int
i) | Int
i <- [Int
0..Int
nInt -> Int -> Int
forall r. Group r => r -> r -> r
-Int
1]]
  in Vector (Entry a)
-> IntMap Int -> IntMap Int -> Int -> Int -> Matrix a
forall a.
Vector (Entry a)
-> IntMap Int -> IntMap Int -> Int -> Int -> Matrix a
Matrix ([Entry a] -> Vector (Entry a)
forall a. [a] -> Vector a
V.fromList [a -> Entry a
forall a. a -> Entry a
newEntry a
forall r. Unital r => r
one Entry a -> (Entry a -> Entry a) -> Entry a
forall a b. a -> (a -> b) -> b
& ((Int, Int) -> Identity (Int, Int))
-> Entry a -> Identity (Entry a)
forall a. Lens' (Entry a) (Int, Int)
idx (((Int, Int) -> Identity (Int, Int))
 -> Entry a -> Identity (Entry a))
-> (Int, Int) -> Entry a -> Entry a
forall s t a b. ASetter s t a b -> b -> s -> t
.~ (Int
i,Int
i) | Int
i <- [Int
0..Int
nInt -> Int -> Int
forall r. Group r => r -> r -> r
-Int
1]])
            IntMap Int
idMap IntMap Int
idMap Int
n Int
n

diag :: DecidableZero a => Vector a -> Matrix a
diag :: Vector a -> Matrix a
diag Vector a
v =
  let n :: Int
n = Vector a -> Int
forall a. Vector a -> Int
V.length Vector a
v
      idMap :: IntMap Int
idMap = [(Int, Int)] -> IntMap Int
forall a. [(Int, a)] -> IntMap a
IM.fromList [(Int
i,Int
i) | Int
i <- [Int
0..Int
nInt -> Int -> Int
forall r. Group r => r -> r -> r
-Int
1]]
  in Matrix a -> Matrix a
forall a. DecidableZero a => Matrix a -> Matrix a
clearZero (Matrix a -> Matrix a) -> Matrix a -> Matrix a
forall a b. (a -> b) -> a -> b
$ Vector (Entry a)
-> IntMap Int -> IntMap Int -> Int -> Int -> Matrix a
forall a.
Vector (Entry a)
-> IntMap Int -> IntMap Int -> Int -> Int -> Matrix a
Matrix ((Int -> a -> Entry a) -> Vector a -> Vector (Entry a)
forall a b. (Int -> a -> b) -> Vector a -> Vector b
V.imap (\Int
i a
a -> a -> Entry a
forall a. a -> Entry a
newEntry a
a Entry a -> (Entry a -> Entry a) -> Entry a
forall a b. a -> (a -> b) -> b
& ((Int, Int) -> Identity (Int, Int))
-> Entry a -> Identity (Entry a)
forall a. Lens' (Entry a) (Int, Int)
idx (((Int, Int) -> Identity (Int, Int))
 -> Entry a -> Identity (Entry a))
-> (Int, Int) -> Entry a -> Entry a
forall s t a b. ASetter s t a b -> b -> s -> t
.~ (Int
i,Int
i)) Vector a
v)
                 IntMap Int
idMap IntMap Int
idMap Int
n Int
n

catDir :: DecidableZero b => Direction -> Matrix b -> Vector b -> Matrix b
catDir :: Direction -> Matrix b -> Vector b -> Matrix b
catDir Direction
dir Matrix b
mat Vector b
vec = (forall s. ST s (Matrix b)) -> Matrix b
forall a. (forall s. ST s a) -> a
runST ((forall s. ST s (Matrix b)) -> Matrix b)
-> (forall s. ST s (Matrix b)) -> Matrix b
forall a b. (a -> b) -> a -> b
$ do
  let seed :: Vector (Int, b)
seed = ((Int, b) -> Bool) -> Vector (Int, b) -> Vector (Int, b)
forall a. (a -> Bool) -> Vector a -> Vector a
V.filter (Bool -> Bool
not (Bool -> Bool) -> ((Int, b) -> Bool) -> (Int, b) -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. b -> Bool
forall r. DecidableZero r => r -> Bool
isZero (b -> Bool) -> ((Int, b) -> b) -> (Int, b) -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Int, b) -> b
forall a b. (a, b) -> b
snd) (Vector (Int, b) -> Vector (Int, b))
-> Vector (Int, b) -> Vector (Int, b)
forall a b. (a -> b) -> a -> b
$ Int -> Vector (Int, b) -> Vector (Int, b)
forall a. Int -> Vector a -> Vector a
V.take (Matrix b
mat Matrix b -> Getting Int (Matrix b) Int -> Int
forall s a. s -> Getting a s a -> a
^. Direction -> Lens' (Matrix b) Int
forall a. Direction -> Lens' (Matrix a) Int
lenL Direction
dir) (Vector (Int, b) -> Vector (Int, b))
-> Vector (Int, b) -> Vector (Int, b)
forall a b. (a -> b) -> a -> b
$ Vector b -> Vector (Int, b)
forall a. Vector a -> Vector (Int, a)
V.indexed Vector b
vec
      n :: Int
n    = Vector (Entry b) -> Int
forall a. Vector a -> Int
V.length (Vector (Entry b) -> Int) -> Vector (Entry b) -> Int
forall a b. (a -> b) -> a -> b
$ Matrix b
mat Matrix b
-> Getting (Vector (Entry b)) (Matrix b) (Vector (Entry b))
-> Vector (Entry b)
forall s a. s -> Getting a s a -> a
^. Getting (Vector (Entry b)) (Matrix b) (Vector (Entry b))
forall a a.
Lens (Matrix a) (Matrix a) (Vector (Entry a)) (Vector (Entry a))
coefficients
      curD :: Int
curD = Matrix b
mat Matrix b -> Getting Int (Matrix b) Int -> Int
forall s a. s -> Getting a s a -> a
^. Direction -> Lens' (Matrix b) Int
forall a. Direction -> Lens' (Matrix a) Int
countL Direction
dir
      getNextIdx :: Int -> Maybe Int
getNextIdx Int
l | Int
l Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 = Maybe Int
forall a. Maybe a
Nothing
                   | Bool
otherwise = Int -> Maybe Int
forall a. a -> Maybe a
Just (Int
nInt -> Int -> Int
forall r. Additive r => r -> r -> r
+Int
lInt -> Int -> Int
forall r. Group r => r -> r -> r
-Int
1)
  MVector s (Entry b)
mv <- (MVector s (Entry b) -> Int -> ST s (MVector s (Entry b)))
-> Int -> MVector s (Entry b) -> ST s (MVector s (Entry b))
forall a b c. (a -> b -> c) -> b -> a -> c
flip MVector s (Entry b) -> Int -> ST s (MVector s (Entry b))
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> m (MVector (PrimState m) a)
grow (Vector (Int, b) -> Int
forall a. Vector a -> Int
V.length Vector (Int, b)
seed) (MVector s (Entry b) -> ST s (MVector s (Entry b)))
-> ST s (MVector s (Entry b)) -> ST s (MVector s (Entry b))
forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
=<< Vector (Entry b) -> ST s (MVector (PrimState (ST s)) (Entry b))
forall (m :: * -> *) a.
PrimMonad m =>
Vector a -> m (MVector (PrimState m) a)
thaw (Matrix b
matMatrix b
-> Getting (Vector (Entry b)) (Matrix b) (Vector (Entry b))
-> Vector (Entry b)
forall s a. s -> Getting a s a -> a
^.Getting (Vector (Entry b)) (Matrix b) (Vector (Entry b))
forall a a.
Lens (Matrix a) (Matrix a) (Vector (Entry a)) (Vector (Entry a))
coefficients)
  let upd :: (Int, b) -> (Int, IntMap Int) -> ST s (Int, IntMap Int)
upd (Int
k, b
v) (Int
l, IntMap Int
opdic) = do
        MVector (PrimState (ST s)) (Entry b) -> Int -> Entry b -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
MV.write MVector s (Entry b)
MVector (PrimState (ST s)) (Entry b)
mv (Int
nInt -> Int -> Int
forall r. Additive r => r -> r -> r
+Int
l) (Entry b -> ST s ()) -> Entry b -> ST s ()
forall a b. (a -> b) -> a -> b
$ b -> Entry b
forall a. a -> Entry a
newEntry b
v
                          Entry b -> (Entry b -> Entry b) -> Entry b
forall a b. a -> (a -> b) -> b
& Direction -> Lens' (Entry b) Int
forall a. Direction -> Lens' (Entry a) Int
nthL Direction
dir ((Int -> Identity Int) -> Entry b -> Identity (Entry b))
-> Int -> Entry b -> Entry b
forall s t a b. ASetter s t a b -> b -> s -> t
.~ Int
k
                          Entry b -> (Entry b -> Entry b) -> Entry b
forall a b. a -> (a -> b) -> b
& Direction -> Lens' (Entry b) Int
forall a. Direction -> Lens' (Entry a) Int
coordL Direction
dir ((Int -> Identity Int) -> Entry b -> Identity (Entry b))
-> Int -> Entry b -> Entry b
forall s t a b. ASetter s t a b -> b -> s -> t
.~ Int
curD
                          Entry b -> (Entry b -> Entry b) -> Entry b
forall a b. a -> (a -> b) -> b
& Direction -> Lens' (Entry b) (Maybe Int)
forall a. Direction -> Lens' (Entry a) (Maybe Int)
nextL Direction
dir ((Maybe Int -> Identity (Maybe Int))
 -> Entry b -> Identity (Entry b))
-> Maybe Int -> Entry b -> Entry b
forall s t a b. ASetter s t a b -> b -> s -> t
.~ Int -> Maybe Int
getNextIdx Int
l
                          Entry b -> (Entry b -> Entry b) -> Entry b
forall a b. a -> (a -> b) -> b
& Direction -> Lens' (Entry b) (Maybe Int)
forall a. Direction -> Lens' (Entry a) (Maybe Int)
nextL (Direction -> Direction
perp Direction
dir) ((Maybe Int -> Identity (Maybe Int))
 -> Entry b -> Identity (Entry b))
-> Maybe Int -> Entry b -> Entry b
forall s t a b. ASetter s t a b -> b -> s -> t
.~ Int -> IntMap Int -> Maybe Int
forall a. Int -> IntMap a -> Maybe a
IM.lookup Int
k IntMap Int
opdic
        (Int, IntMap Int) -> ST s (Int, IntMap Int)
forall (m :: * -> *) a. Monad m => a -> m a
return (Int
lInt -> Int -> Int
forall r. Additive r => r -> r -> r
+Int
1, (Maybe Int -> Maybe Int) -> Int -> IntMap Int -> IntMap Int
forall a. (Maybe a -> Maybe a) -> Int -> IntMap a -> IntMap a
alter (Maybe Int -> Maybe Int -> Maybe Int
forall a b. a -> b -> a
const (Maybe Int -> Maybe Int -> Maybe Int)
-> Maybe Int -> Maybe Int -> Maybe Int
forall a b. (a -> b) -> a -> b
$ Int -> Maybe Int
forall a. a -> Maybe a
Just (Int -> Maybe Int) -> Int -> Maybe Int
forall a b. (a -> b) -> a -> b
$ Int
nInt -> Int -> Int
forall r. Additive r => r -> r -> r
+Int
l) Int
k IntMap Int
opdic)
  (Int
l, IntMap Int
op') <- Getting
  (Endo ((Int, IntMap Int) -> ST s (Int, IntMap Int)))
  (Vector (Int, b))
  (Int, b)
-> ((Int, IntMap Int) -> (Int, b) -> ST s (Int, IntMap Int))
-> (Int, IntMap Int)
-> Vector (Int, b)
-> ST s (Int, IntMap Int)
forall (m :: * -> *) r s a.
Monad m =>
Getting (Endo (r -> m r)) s a -> (r -> a -> m r) -> r -> s -> m r
foldlMOf Getting
  (Endo ((Int, IntMap Int) -> ST s (Int, IntMap Int)))
  (Vector (Int, b))
  (Int, b)
forall (f :: * -> *) a. Foldable f => IndexedFold Int (f a) a
folded (((Int, b) -> (Int, IntMap Int) -> ST s (Int, IntMap Int))
-> (Int, IntMap Int) -> (Int, b) -> ST s (Int, IntMap Int)
forall a b c. (a -> b -> c) -> b -> a -> c
flip (Int, b) -> (Int, IntMap Int) -> ST s (Int, IntMap Int)
upd) (Int
0, Matrix b
mat Matrix b
-> Getting (IntMap Int) (Matrix b) (IntMap Int) -> IntMap Int
forall s a. s -> Getting a s a -> a
^. Direction -> Lens' (Matrix b) (IntMap Int)
forall a. Direction -> Lens' (Matrix a) (IntMap Int)
startL (Direction -> Direction
perp Direction
dir)) Vector (Int, b)
seed
  Vector (Entry b)
v <- MVector (PrimState (ST s)) (Entry b) -> ST s (Vector (Entry b))
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> m (Vector a)
unsafeFreeze MVector s (Entry b)
MVector (PrimState (ST s)) (Entry b)
mv
  Matrix b -> ST s (Matrix b)
forall (m :: * -> *) a. Monad m => a -> m a
return (Matrix b -> ST s (Matrix b)) -> Matrix b -> ST s (Matrix b)
forall a b. (a -> b) -> a -> b
$ Matrix b
mat Matrix b -> (Matrix b -> Matrix b) -> Matrix b
forall a b. a -> (a -> b) -> b
& Direction -> Lens' (Matrix b) Int
forall a. Direction -> Lens' (Matrix a) Int
countL Direction
dir ((Int -> Identity Int) -> Matrix b -> Identity (Matrix b))
-> Int -> Matrix b -> Matrix b
forall a s t. Num a => ASetter s t a a -> a -> s -> t
+~ Int
1
               Matrix b -> (Matrix b -> Matrix b) -> Matrix b
forall a b. a -> (a -> b) -> b
& Direction -> Lens' (Matrix b) (IntMap Int)
forall a. Direction -> Lens' (Matrix a) (IntMap Int)
startL Direction
dir ((IntMap Int -> Identity (IntMap Int))
 -> Matrix b -> Identity (Matrix b))
-> (IntMap Int -> IntMap Int) -> Matrix b -> Matrix b
forall s t a b. ASetter s t a b -> (a -> b) -> s -> t
%~ (Maybe Int -> Maybe Int) -> Int -> IntMap Int -> IntMap Int
forall a. (Maybe a -> Maybe a) -> Int -> IntMap a -> IntMap a
alter (Maybe Int -> Maybe Int -> Maybe Int
forall a b. a -> b -> a
const (Maybe Int -> Maybe Int -> Maybe Int)
-> Maybe Int -> Maybe Int -> Maybe Int
forall a b. (a -> b) -> a -> b
$ Int -> Maybe Int
getNextIdx Int
l) Int
curD
               Matrix b -> (Matrix b -> Matrix b) -> Matrix b
forall a b. a -> (a -> b) -> b
& Direction -> Lens' (Matrix b) (IntMap Int)
forall a. Direction -> Lens' (Matrix a) (IntMap Int)
startL (Direction -> Direction
perp Direction
dir) ((IntMap Int -> Identity (IntMap Int))
 -> Matrix b -> Identity (Matrix b))
-> IntMap Int -> Matrix b -> Matrix b
forall s t a b. ASetter s t a b -> b -> s -> t
.~ IntMap Int
op'
               Matrix b -> (Matrix b -> Matrix b) -> Matrix b
forall a b. a -> (a -> b) -> b
& (Vector (Entry b) -> Identity (Vector (Entry b)))
-> Matrix b -> Identity (Matrix b)
forall a a.
Lens (Matrix a) (Matrix a) (Vector (Entry a)) (Vector (Entry a))
coefficients ((Vector (Entry b) -> Identity (Vector (Entry b)))
 -> Matrix b -> Identity (Matrix b))
-> Vector (Entry b) -> Matrix b -> Matrix b
forall s t a b. ASetter s t a b -> b -> s -> t
.~ Vector (Entry b)
v

dirVector :: DecidableZero a => Direction -> Vector a -> Matrix a
dirVector :: Direction -> Vector a -> Matrix a
dirVector Direction
Row = Vector a -> Matrix a
forall a. DecidableZero a => Vector a -> Matrix a
rowVector
dirVector Direction
Column = Vector a -> Matrix a
forall a. DecidableZero a => Vector a -> Matrix a
colVector

rowVector :: DecidableZero a => Vector a -> Matrix a
rowVector :: Vector a -> Matrix a
rowVector = [[a]] -> Matrix a
forall a. DecidableZero a => [[a]] -> Matrix a
fromLists([[a]] -> Matrix a) -> (Vector a -> [[a]]) -> Vector a -> Matrix a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ([a] -> [[a]] -> [[a]]
forall a. a -> [a] -> [a]
:[]) ([a] -> [[a]]) -> (Vector a -> [a]) -> Vector a -> [[a]]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector a -> [a]
forall a. Vector a -> [a]
V.toList

colVector :: DecidableZero a => Vector a -> Matrix a
colVector :: Vector a -> Matrix a
colVector = [[a]] -> Matrix a
forall a. DecidableZero a => [[a]] -> Matrix a
fromLists ([[a]] -> Matrix a) -> (Vector a -> [[a]]) -> Vector a -> Matrix a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (a -> [a]) -> [a] -> [[a]]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
map (a -> [a] -> [a]
forall a. a -> [a] -> [a]
:[]) ([a] -> [[a]]) -> (Vector a -> [a]) -> Vector a -> [[a]]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector a -> [a]
forall a. Vector a -> [a]
V.toList

toDirs :: Monoidal a => Direction -> Matrix a -> [Vector a]
toDirs :: Direction -> Matrix a -> [Vector a]
toDirs Direction
dir Matrix a
mat = [ Direction -> Int -> Matrix a -> Vector a
forall a. Monoidal a => Direction -> Int -> Matrix a -> Vector a
getDir Direction
dir Int
i Matrix a
mat | Int
i <- [Int
0..Matrix a
matMatrix a -> Getting Int (Matrix a) Int -> Int
forall s a. s -> Getting a s a -> a
^.Direction -> Lens' (Matrix a) Int
forall a. Direction -> Lens' (Matrix a) Int
countL Direction
dirInt -> Int -> Int
forall r. Group r => r -> r -> r
-Int
1]]

toRows :: Monoidal a => Matrix a -> [Vector a]
toRows :: Matrix a -> [Vector a]
toRows = Direction -> Matrix a -> [Vector a]
forall a. Monoidal a => Direction -> Matrix a -> [Vector a]
toDirs Direction
Row

toCols :: Monoidal a => Matrix a -> [Vector a]
toCols :: Matrix a -> [Vector a]
toCols = Direction -> Matrix a -> [Vector a]
forall a. Monoidal a => Direction -> Matrix a -> [Vector a]
toDirs Direction
Column

appendDir :: DecidableZero b => Direction -> Matrix b -> Matrix b -> Matrix b
appendDir :: Direction -> Matrix b -> Matrix b -> Matrix b
appendDir Direction
dir Matrix b
m = (Matrix b -> Vector b -> Matrix b)
-> Matrix b -> [Vector b] -> Matrix b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl (Direction -> Matrix b -> Vector b -> Matrix b
forall b.
DecidableZero b =>
Direction -> Matrix b -> Vector b -> Matrix b
catDir Direction
dir) Matrix b
m ([Vector b] -> Matrix b)
-> (Matrix b -> [Vector b]) -> Matrix b -> Matrix b
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Direction -> Matrix b -> [Vector b]
forall a. Monoidal a => Direction -> Matrix a -> [Vector a]
toDirs Direction
dir

(<-->) :: DecidableZero b => Matrix b -> Matrix b -> Matrix b
<--> :: Matrix b -> Matrix b -> Matrix b
(<-->) = Direction -> Matrix b -> Matrix b -> Matrix b
forall b.
DecidableZero b =>
Direction -> Matrix b -> Matrix b -> Matrix b
appendDir Direction
Row

(<||>) :: DecidableZero b => Matrix b -> Matrix b -> Matrix b
<||> :: Matrix b -> Matrix b -> Matrix b
(<||>) = Direction -> Matrix b -> Matrix b -> Matrix b
forall b.
DecidableZero b =>
Direction -> Matrix b -> Matrix b -> Matrix b
appendDir Direction
Column

catRow :: DecidableZero b => Matrix b -> Vector b -> Matrix b
catRow :: Matrix b -> Vector b -> Matrix b
catRow = Direction -> Matrix b -> Vector b -> Matrix b
forall b.
DecidableZero b =>
Direction -> Matrix b -> Vector b -> Matrix b
catDir Direction
Row

catCol :: DecidableZero b => Matrix b -> Vector b -> Matrix b
catCol :: Matrix b -> Vector b -> Matrix b
catCol = Direction -> Matrix b -> Vector b -> Matrix b
forall b.
DecidableZero b =>
Direction -> Matrix b -> Vector b -> Matrix b
catDir Direction
Column

switchRows :: Int -> Int -> Matrix a -> Matrix a
switchRows :: Int -> Int -> Matrix a -> Matrix a
switchRows = Int -> Int -> Matrix a -> Matrix a
forall a. Int -> Int -> Matrix a -> Matrix a
swapRows

switchCols :: Int -> Int -> Matrix a -> Matrix a
switchCols :: Int -> Int -> Matrix a -> Matrix a
switchCols = Int -> Int -> Matrix a -> Matrix a
forall a. Int -> Int -> Matrix a -> Matrix a
swapCols

cmap :: DecidableZero a => (a1 -> a) -> Matrix a1 -> Matrix a
cmap :: (a1 -> a) -> Matrix a1 -> Matrix a
cmap a1 -> a
f = Matrix a -> Matrix a
forall a. DecidableZero a => Matrix a -> Matrix a
clearZero (Matrix a -> Matrix a)
-> (Matrix a1 -> Matrix a) -> Matrix a1 -> Matrix a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((Vector (Entry a1) -> Identity (Vector (Entry a)))
-> Matrix a1 -> Identity (Matrix a)
forall a a.
Lens (Matrix a) (Matrix a) (Vector (Entry a)) (Vector (Entry a))
coefficients ((Vector (Entry a1) -> Identity (Vector (Entry a)))
 -> Matrix a1 -> Identity (Matrix a))
-> ((a1 -> Identity a)
    -> Vector (Entry a1) -> Identity (Vector (Entry a)))
-> (a1 -> Identity a)
-> Matrix a1
-> Identity (Matrix a)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Entry a1 -> Identity (Entry a))
-> Vector (Entry a1) -> Identity (Vector (Entry a))
forall (f :: * -> *) a b. Functor f => Setter (f a) (f b) a b
mapped ((Entry a1 -> Identity (Entry a))
 -> Vector (Entry a1) -> Identity (Vector (Entry a)))
-> ((a1 -> Identity a) -> Entry a1 -> Identity (Entry a))
-> (a1 -> Identity a)
-> Vector (Entry a1)
-> Identity (Vector (Entry a))
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (a1 -> Identity a) -> Entry a1 -> Identity (Entry a)
forall a a. Lens (Entry a) (Entry a) a a
value ((a1 -> Identity a) -> Matrix a1 -> Identity (Matrix a))
-> (a1 -> a) -> Matrix a1 -> Matrix a
forall s t a b. ASetter s t a b -> (a -> b) -> s -> t
%~ a1 -> a
f)

clearZero :: DecidableZero a => Matrix a -> Matrix a
clearZero :: Matrix a -> Matrix a
clearZero Matrix a
mat = (Int -> Entry a -> Matrix a -> Matrix a)
-> Matrix a -> Vector (Entry a) -> Matrix a
forall a b. (Int -> a -> b -> b) -> b -> Vector a -> b
V.ifoldr (\Int
i Entry a
v Matrix a
m -> if a -> Bool
forall r. DecidableZero r => r -> Bool
isZero (Entry a
vEntry a -> Getting a (Entry a) a -> a
forall s a. s -> Getting a s a -> a
^.Getting a (Entry a) a
forall a a. Lens (Entry a) (Entry a) a a
value) then Int -> Matrix a -> Matrix a
forall a. Int -> Matrix a -> Matrix a
clearAt Int
i Matrix a
m else Matrix a
m)
                Matrix a
mat (Matrix a
mat Matrix a
-> Getting (Vector (Entry a)) (Matrix a) (Vector (Entry a))
-> Vector (Entry a)
forall s a. s -> Getting a s a -> a
^. Getting (Vector (Entry a)) (Matrix a) (Vector (Entry a))
forall a a.
Lens (Matrix a) (Matrix a) (Vector (Entry a)) (Vector (Entry a))
coefficients)

transpose :: Matrix a -> Matrix a
transpose :: Matrix a -> Matrix a
transpose Matrix a
mat = Matrix a
mat Matrix a -> (Matrix a -> Matrix a) -> Matrix a
forall a b. a -> (a -> b) -> b
& (IntMap Int -> Identity (IntMap Int))
-> Matrix a -> Identity (Matrix a)
forall a. Lens' (Matrix a) (IntMap Int)
rowStart ((IntMap Int -> Identity (IntMap Int))
 -> Matrix a -> Identity (Matrix a))
-> IntMap Int -> Matrix a -> Matrix a
forall s t a b. ASetter s t a b -> b -> s -> t
.~ Matrix a
matMatrix a
-> Getting (IntMap Int) (Matrix a) (IntMap Int) -> IntMap Int
forall s a. s -> Getting a s a -> a
^.Getting (IntMap Int) (Matrix a) (IntMap Int)
forall a. Lens' (Matrix a) (IntMap Int)
colStart
                    Matrix a -> (Matrix a -> Matrix a) -> Matrix a
forall a b. a -> (a -> b) -> b
& (IntMap Int -> Identity (IntMap Int))
-> Matrix a -> Identity (Matrix a)
forall a. Lens' (Matrix a) (IntMap Int)
colStart ((IntMap Int -> Identity (IntMap Int))
 -> Matrix a -> Identity (Matrix a))
-> IntMap Int -> Matrix a -> Matrix a
forall s t a b. ASetter s t a b -> b -> s -> t
.~ Matrix a
matMatrix a
-> Getting (IntMap Int) (Matrix a) (IntMap Int) -> IntMap Int
forall s a. s -> Getting a s a -> a
^.Getting (IntMap Int) (Matrix a) (IntMap Int)
forall a. Lens' (Matrix a) (IntMap Int)
rowStart
                    Matrix a -> (Matrix a -> Matrix a) -> Matrix a
forall a b. a -> (a -> b) -> b
& (Int -> Identity Int) -> Matrix a -> Identity (Matrix a)
forall a. Lens' (Matrix a) Int
height   ((Int -> Identity Int) -> Matrix a -> Identity (Matrix a))
-> Int -> Matrix a -> Matrix a
forall s t a b. ASetter s t a b -> b -> s -> t
.~ Matrix a
matMatrix a -> Getting Int (Matrix a) Int -> Int
forall s a. s -> Getting a s a -> a
^.Getting Int (Matrix a) Int
forall a. Lens' (Matrix a) Int
width
                    Matrix a -> (Matrix a -> Matrix a) -> Matrix a
forall a b. a -> (a -> b) -> b
& (Int -> Identity Int) -> Matrix a -> Identity (Matrix a)
forall a. Lens' (Matrix a) Int
width    ((Int -> Identity Int) -> Matrix a -> Identity (Matrix a))
-> Int -> Matrix a -> Matrix a
forall s t a b. ASetter s t a b -> b -> s -> t
.~ Matrix a
matMatrix a -> Getting Int (Matrix a) Int -> Int
forall s a. s -> Getting a s a -> a
^.Getting Int (Matrix a) Int
forall a. Lens' (Matrix a) Int
height
                    Matrix a -> (Matrix a -> Matrix a) -> Matrix a
forall a b. a -> (a -> b) -> b
& (Vector (Entry a) -> Identity (Vector (Entry a)))
-> Matrix a -> Identity (Matrix a)
forall a a.
Lens (Matrix a) (Matrix a) (Vector (Entry a)) (Vector (Entry a))
coefficients ((Vector (Entry a) -> Identity (Vector (Entry a)))
 -> Matrix a -> Identity (Matrix a))
-> ((Entry a -> Identity (Entry a))
    -> Vector (Entry a) -> Identity (Vector (Entry a)))
-> (Entry a -> Identity (Entry a))
-> Matrix a
-> Identity (Matrix a)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Entry a -> Identity (Entry a))
-> Vector (Entry a) -> Identity (Vector (Entry a))
forall s t a b. Each s t a b => Traversal s t a b
each ((Entry a -> Identity (Entry a))
 -> Matrix a -> Identity (Matrix a))
-> (Entry a -> Entry a) -> Matrix a -> Matrix a
forall s t a b. ASetter s t a b -> (a -> b) -> s -> t
%~ Entry a -> Entry a
forall a. Entry a -> Entry a
swapEntry
  where
    swapEntry :: Entry a -> Entry a
swapEntry Entry a
ent = Entry a
ent Entry a -> (Entry a -> Entry a) -> Entry a
forall a b. a -> (a -> b) -> b
& ((Int, Int) -> Identity (Int, Int))
-> Entry a -> Identity (Entry a)
forall a. Lens' (Entry a) (Int, Int)
idx     (((Int, Int) -> Identity (Int, Int))
 -> Entry a -> Identity (Entry a))
-> ((Int, Int) -> (Int, Int)) -> Entry a -> Entry a
forall s t a b. ASetter s t a b -> (a -> b) -> s -> t
%~ (Int, Int) -> (Int, Int)
forall a b. (a, b) -> (b, a)
swap
                        Entry a -> (Entry a -> Entry a) -> Entry a
forall a b. a -> (a -> b) -> b
& (Maybe Int -> Identity (Maybe Int))
-> Entry a -> Identity (Entry a)
forall a. Lens' (Entry a) (Maybe Int)
rowNext ((Maybe Int -> Identity (Maybe Int))
 -> Entry a -> Identity (Entry a))
-> Maybe Int -> Entry a -> Entry a
forall s t a b. ASetter s t a b -> b -> s -> t
.~ Entry a
ent Entry a -> Getting (Maybe Int) (Entry a) (Maybe Int) -> Maybe Int
forall s a. s -> Getting a s a -> a
^. Getting (Maybe Int) (Entry a) (Maybe Int)
forall a. Lens' (Entry a) (Maybe Int)
colNext
                        Entry a -> (Entry a -> Entry a) -> Entry a
forall a b. a -> (a -> b) -> b
& (Maybe Int -> Identity (Maybe Int))
-> Entry a -> Identity (Entry a)
forall a. Lens' (Entry a) (Maybe Int)
colNext ((Maybe Int -> Identity (Maybe Int))
 -> Entry a -> Identity (Entry a))
-> Maybe Int -> Entry a -> Entry a
forall s t a b. ASetter s t a b -> b -> s -> t
.~ Entry a
ent Entry a -> Getting (Maybe Int) (Entry a) (Maybe Int) -> Maybe Int
forall s a. s -> Getting a s a -> a
^. Getting (Maybe Int) (Entry a) (Maybe Int)
forall a. Lens' (Entry a) (Maybe Int)
rowNext

zeroMat :: Int -> Int -> Matrix a
zeroMat :: Int -> Int -> Matrix a
zeroMat = Vector (Entry a)
-> IntMap Int -> IntMap Int -> Int -> Int -> Matrix a
forall a.
Vector (Entry a)
-> IntMap Int -> IntMap Int -> Int -> Int -> Matrix a
Matrix Vector (Entry a)
forall a. Vector a
V.empty IntMap Int
forall a. IntMap a
IM.empty IntMap Int
forall a. IntMap a
IM.empty

dirCount :: Direction -> Int -> Matrix a -> Int
dirCount :: Direction -> Int -> Matrix a -> Int
dirCount = Int
-> (Int -> Int -> Entry a -> Int)
-> Direction
-> Int
-> Matrix a
-> Int
forall b a.
b
-> (b -> Int -> Entry a -> b) -> Direction -> Int -> Matrix a -> b
traverseDir Int
0 (\Int
a Int
_ Entry a
_ -> Int -> Int
forall a. Enum a => a -> a
succ Int
a)

rowCount :: Int -> Matrix a -> Int
rowCount :: Int -> Matrix a -> Int
rowCount = Direction -> Int -> Matrix a -> Int
forall a. Direction -> Int -> Matrix a -> Int
dirCount Direction
Row

colCount :: Int -> Matrix a -> Int
colCount :: Int -> Matrix a -> Int
colCount = Direction -> Int -> Matrix a -> Int
forall a. Direction -> Int -> Matrix a -> Int
dirCount Direction
Column

instance (Ord a, Semigroup b) => Semigroup (MaxEntry a b) where
  MaxEntry a
a b
as <> :: MaxEntry a b -> MaxEntry a b -> MaxEntry a b
<> MaxEntry a
b b
bs =
    case a -> a -> Ordering
forall a. Ord a => a -> a -> Ordering
compare a
a a
b of
      Ordering
EQ -> a -> b -> MaxEntry a b
forall a b. a -> b -> MaxEntry a b
MaxEntry a
a (b
as b -> b -> b
forall a. Semigroup a => a -> a -> a
<> b
bs)
      Ordering
LT -> a -> b -> MaxEntry a b
forall a b. a -> b -> MaxEntry a b
MaxEntry a
b b
bs
      Ordering
GT -> a -> b -> MaxEntry a b
forall a b. a -> b -> MaxEntry a b
MaxEntry a
a b
as

instance (Ord a, Bounded a, Monoid b) => Monoid (MaxEntry a b) where
  mappend :: MaxEntry a b -> MaxEntry a b -> MaxEntry a b
mappend (MaxEntry a
a b
as) (MaxEntry a
b b
bs) =
    case a -> a -> Ordering
forall a. Ord a => a -> a -> Ordering
compare a
a a
b of
      Ordering
EQ -> a -> b -> MaxEntry a b
forall a b. a -> b -> MaxEntry a b
MaxEntry a
a (b
as b -> b -> b
forall a. Monoid a => a -> a -> a
`mappend` b
bs)
      Ordering
LT -> a -> b -> MaxEntry a b
forall a b. a -> b -> MaxEntry a b
MaxEntry a
b b
bs
      Ordering
GT -> a -> b -> MaxEntry a b
forall a b. a -> b -> MaxEntry a b
MaxEntry a
a b
as
  mempty :: MaxEntry a b
mempty = a -> b -> MaxEntry a b
forall a b. a -> b -> MaxEntry a b
MaxEntry a
forall a. Bounded a => a
minBound b
forall a. Monoid a => a
mempty


newGaussianState :: Unital a => Matrix a -> GaussianState a
newGaussianState :: Matrix a -> GaussianState a
newGaussianState Matrix a
inp =
  Matrix a
-> Matrix a -> Int -> IntSet -> Int -> a -> GaussianState a
forall a.
Matrix a
-> Matrix a -> Int -> IntSet -> Int -> a -> GaussianState a
GaussianState Matrix a
inp (Int -> Matrix a
forall a. Unital a => Int -> Matrix a
identity (Int -> Matrix a) -> Int -> Matrix a
forall a b. (a -> b) -> a -> b
$ Matrix a
inp Matrix a -> Getting Int (Matrix a) Int -> Int
forall s a. s -> Getting a s a -> a
^. Getting Int (Matrix a) Int
forall a. Lens' (Matrix a) Int
height) (-Int
1) (IntSet -> Matrix a -> IntSet
forall a. IntSet -> Matrix a -> IntSet
getHeaviest IntSet
IS.empty Matrix a
inp) Int
0 a
forall r. Unital r => r
one

getHeaviest :: IntSet -> Matrix a -> IntSet
getHeaviest :: IntSet -> Matrix a -> IntSet
getHeaviest IntSet
old Matrix a
inp =
  if IntSet -> Int
IS.size IntSet
old Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Matrix a
inpMatrix a -> Getting Int (Matrix a) Int -> Int
forall s a. s -> Getting a s a -> a
^.Getting Int (Matrix a) Int
forall a. Lens' (Matrix a) Int
widthInt -> Int -> Int
forall r. Multiplicative r => r -> r -> r
*Int
5Int -> Int -> Int
forall a. Integral a => a -> a -> a
`P.div`Int
100
  then IntSet
old
  else  let news :: IntSet
news = MaxEntry Int IntSet -> IntSet
forall a b. MaxEntry a b -> b
entry (MaxEntry Int IntSet -> IntSet) -> MaxEntry Int IntSet -> IntSet
forall a b. (a -> b) -> a -> b
$ [MaxEntry Int IntSet] -> MaxEntry Int IntSet
forall w. Monoid w => [w] -> w
mconcat ([MaxEntry Int IntSet] -> MaxEntry Int IntSet)
-> [MaxEntry Int IntSet] -> MaxEntry Int IntSet
forall a b. (a -> b) -> a -> b
$ (Int -> MaxEntry Int IntSet) -> [Int] -> [MaxEntry Int IntSet]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
map (\Int
k -> Int -> IntSet -> MaxEntry Int IntSet
forall a b. a -> b -> MaxEntry a b
MaxEntry (Int -> Matrix a -> Int
forall a. Int -> Matrix a -> Int
colCount Int
k Matrix a
inp) (Int -> IntSet
IS.singleton Int
k)) ([Int] -> [MaxEntry Int IntSet]) -> [Int] -> [MaxEntry Int IntSet]
forall a b. (a -> b) -> a -> b
$
                   IntSet -> [Int]
IS.toList (IntSet -> [Int]) -> IntSet -> [Int]
forall a b. (a -> b) -> a -> b
$ IntMap Int -> IntSet
forall a. IntMap a -> IntSet
IM.keysSet (Matrix a
inpMatrix a
-> Getting (IntMap Int) (Matrix a) (IntMap Int) -> IntMap Int
forall s a. s -> Getting a s a -> a
^.Getting (IntMap Int) (Matrix a) (IntMap Int)
forall a. Lens' (Matrix a) (IntMap Int)
colStart) IntSet -> IntSet -> IntSet
IS.\\ IntSet
old
        in IntSet
news IntSet -> IntSet -> IntSet
`IS.union` IntSet
old

traverseRow :: b -> (b -> Int -> Entry a -> b) -> Int -> Matrix a -> b
traverseRow :: b -> (b -> Int -> Entry a -> b) -> Int -> Matrix a -> b
traverseRow b
a b -> Int -> Entry a -> b
f = b
-> (b -> Int -> Entry a -> b) -> Direction -> Int -> Matrix a -> b
forall b a.
b
-> (b -> Int -> Entry a -> b) -> Direction -> Int -> Matrix a -> b
traverseDir b
a b -> Int -> Entry a -> b
f Direction
Row

traverseCol :: b -> (b -> Int -> Entry a -> b) -> Int -> Matrix a -> b
traverseCol :: b -> (b -> Int -> Entry a -> b) -> Int -> Matrix a -> b
traverseCol b
a b -> Int -> Entry a -> b
f = b
-> (b -> Int -> Entry a -> b) -> Direction -> Int -> Matrix a -> b
forall b a.
b
-> (b -> Int -> Entry a -> b) -> Direction -> Int -> Matrix a -> b
traverseDir b
a b -> Int -> Entry a -> b
f Direction
Column

structuredGauss :: (DecidableZero a, Division a, Group a)
                => Matrix a -> (Matrix a, Matrix a)
structuredGauss :: Matrix a -> (Matrix a, Matrix a)
structuredGauss = Matrix a -> (Matrix a, Matrix a, a)
forall a.
(DecidableZero a, Division a, Group a) =>
Matrix a -> (Matrix a, Matrix a, a)
structuredGauss' (Matrix a -> (Matrix a, Matrix a, a))
-> ((Matrix a, Matrix a, a) -> (Matrix a, Matrix a))
-> Matrix a
-> (Matrix a, Matrix a)
forall k (cat :: k -> k -> *) (a :: k) (b :: k) (c :: k).
Category cat =>
cat a b -> cat b c -> cat a c
>>> Getting (Matrix a) (Matrix a, Matrix a, a) (Matrix a)
-> (Matrix a, Matrix a, a) -> Matrix a
forall s (m :: * -> *) a. MonadReader s m => Getting a s a -> m a
view Getting (Matrix a) (Matrix a, Matrix a, a) (Matrix a)
forall s t a b. Field1 s t a b => Lens s t a b
_1 ((Matrix a, Matrix a, a) -> Matrix a)
-> ((Matrix a, Matrix a, a) -> Matrix a)
-> (Matrix a, Matrix a, a)
-> (Matrix a, Matrix a)
forall (a :: * -> * -> *) b c c'.
Arrow a =>
a b c -> a b c' -> a b (c, c')
&&& Getting (Matrix a) (Matrix a, Matrix a, a) (Matrix a)
-> (Matrix a, Matrix a, a) -> Matrix a
forall s (m :: * -> *) a. MonadReader s m => Getting a s a -> m a
view Getting (Matrix a) (Matrix a, Matrix a, a) (Matrix a)
forall s t a b. Field2 s t a b => Lens s t a b
_2

structuredGauss' :: (DecidableZero a, Division a, Group a)
                 => Matrix a -> (Matrix a, Matrix a, a)
structuredGauss' :: Matrix a -> (Matrix a, Matrix a, a)
structuredGauss' = State (GaussianState a) (Matrix a, Matrix a, a)
-> GaussianState a -> (Matrix a, Matrix a, a)
forall s a. State s a -> s -> a
evalState State (GaussianState a) (Matrix a, Matrix a, a)
forall a (m :: * -> *).
(MonadState (GaussianState a) m, Group a, Division a,
 DecidableZero a) =>
m (Matrix a, Matrix a, a)
go (GaussianState a -> (Matrix a, Matrix a, a))
-> (Matrix a -> GaussianState a)
-> Matrix a
-> (Matrix a, Matrix a, a)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Matrix a -> GaussianState a
forall a. Unital a => Matrix a -> GaussianState a
newGaussianState
  where
    countLight :: IntSet -> Int -> Matrix a -> Int
countLight IntSet
heavys = Int -> (Int -> Int -> Entry a -> Int) -> Int -> Matrix a -> Int
forall b a. b -> (b -> Int -> Entry a -> b) -> Int -> Matrix a -> b
traverseRow (Int
0 :: Int)
                           (\(!Int
c) Int
_ Entry a
ent -> if (Entry a
entEntry a -> Getting Int (Entry a) Int -> Int
forall s a. s -> Getting a s a -> a
^.Direction -> Lens' (Entry a) Int
forall a. Direction -> Lens' (Entry a) Int
coordL Direction
Column) Int -> IntSet -> Bool
`IS.member` IntSet
heavys
                                        then Int
c
                                        else Int
cInt -> Int -> Int
forall r. Additive r => r -> r -> r
+Int
1)
    go :: m (Matrix a, Matrix a, a)
go = do
      Matrix a
old <- Getting (Matrix a) (GaussianState a) (Matrix a) -> m (Matrix a)
forall s (m :: * -> *) a. MonadState s m => Getting a s a -> m a
use Getting (Matrix a) (GaussianState a) (Matrix a)
forall a. Lens' (GaussianState a) (Matrix a)
input
      Int
destRow <- Getting Int (GaussianState a) Int -> m Int
forall s (m :: * -> *) a. MonadState s m => Getting a s a -> m a
use Getting Int (GaussianState a) Int
forall a. Lens' (GaussianState a) Int
curRow
      Int
pcol <- Getting Int (GaussianState a) Int -> m Int
forall s (m :: * -> *) a. MonadState s m => Getting a s a -> m a
use Getting Int (GaussianState a) Int
forall a. Lens' (GaussianState a) Int
prevCol
      (IntMap Int
_, IntMap Int
rest) <- LensLike'
  (Const (IntMap Int, IntMap Int)) (GaussianState a) (IntMap Int)
-> (IntMap Int -> (IntMap Int, IntMap Int))
-> m (IntMap Int, IntMap Int)
forall s (m :: * -> *) r a.
MonadState s m =>
LensLike' (Const r) s a -> (a -> r) -> m r
uses ((Matrix a -> Const (IntMap Int, IntMap Int) (Matrix a))
-> GaussianState a
-> Const (IntMap Int, IntMap Int) (GaussianState a)
forall a. Lens' (GaussianState a) (Matrix a)
input((Matrix a -> Const (IntMap Int, IntMap Int) (Matrix a))
 -> GaussianState a
 -> Const (IntMap Int, IntMap Int) (GaussianState a))
-> ((IntMap Int -> Const (IntMap Int, IntMap Int) (IntMap Int))
    -> Matrix a -> Const (IntMap Int, IntMap Int) (Matrix a))
-> LensLike'
     (Const (IntMap Int, IntMap Int)) (GaussianState a) (IntMap Int)
forall b c a. (b -> c) -> (a -> b) -> a -> c
.(IntMap Int -> Const (IntMap Int, IntMap Int) (IntMap Int))
-> Matrix a -> Const (IntMap Int, IntMap Int) (Matrix a)
forall a. Lens' (Matrix a) (IntMap Int)
colStart) (Int -> IntMap Int -> (IntMap Int, IntMap Int)
forall a. Int -> IntMap a -> (IntMap a, IntMap a)
IM.split Int
pcol)
      case IntMap Int -> Maybe ((Int, Int), IntMap Int)
forall a. IntMap a -> Maybe ((Int, a), IntMap a)
minViewWithKey IntMap Int
rest of
        Maybe ((Int, Int), IntMap Int)
_ | Int
destRow Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Matrix a
old Matrix a -> Getting Int (Matrix a) Int -> Int
forall s a. s -> Getting a s a -> a
^. Getting Int (Matrix a) Int
forall a. Lens' (Matrix a) Int
height -> (,,) (Matrix a -> Matrix a -> a -> (Matrix a, Matrix a, a))
-> m (Matrix a) -> m (Matrix a -> a -> (Matrix a, Matrix a, a))
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Getting (Matrix a) (GaussianState a) (Matrix a) -> m (Matrix a)
forall s (m :: * -> *) a. MonadState s m => Getting a s a -> m a
use Getting (Matrix a) (GaussianState a) (Matrix a)
forall a. Lens' (GaussianState a) (Matrix a)
input m (Matrix a -> a -> (Matrix a, Matrix a, a))
-> m (Matrix a) -> m (a -> (Matrix a, Matrix a, a))
forall (f :: * -> *) a b. Applicative f => f (a -> b) -> f a -> f b
<*> Getting (Matrix a) (GaussianState a) (Matrix a) -> m (Matrix a)
forall s (m :: * -> *) a. MonadState s m => Getting a s a -> m a
use Getting (Matrix a) (GaussianState a) (Matrix a)
forall a. Lens' (GaussianState a) (Matrix a)
output m (a -> (Matrix a, Matrix a, a))
-> m a -> m (Matrix a, Matrix a, a)
forall (f :: * -> *) a b. Applicative f => f (a -> b) -> f a -> f b
<*> Getting a (GaussianState a) a -> m a
forall s (m :: * -> *) a. MonadState s m => Getting a s a -> m a
use Getting a (GaussianState a) a
forall a. Lens' (GaussianState a) a
detAcc
        Maybe ((Int, Int), IntMap Int)
Nothing -> (,,) (Matrix a -> Matrix a -> a -> (Matrix a, Matrix a, a))
-> m (Matrix a) -> m (Matrix a -> a -> (Matrix a, Matrix a, a))
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Getting (Matrix a) (GaussianState a) (Matrix a) -> m (Matrix a)
forall s (m :: * -> *) a. MonadState s m => Getting a s a -> m a
use Getting (Matrix a) (GaussianState a) (Matrix a)
forall a. Lens' (GaussianState a) (Matrix a)
input m (Matrix a -> a -> (Matrix a, Matrix a, a))
-> m (Matrix a) -> m (a -> (Matrix a, Matrix a, a))
forall (f :: * -> *) a b. Applicative f => f (a -> b) -> f a -> f b
<*> Getting (Matrix a) (GaussianState a) (Matrix a) -> m (Matrix a)
forall s (m :: * -> *) a. MonadState s m => Getting a s a -> m a
use Getting (Matrix a) (GaussianState a) (Matrix a)
forall a. Lens' (GaussianState a) (Matrix a)
output m (a -> (Matrix a, Matrix a, a))
-> m a -> m (Matrix a, Matrix a, a)
forall (f :: * -> *) a b. Applicative f => f (a -> b) -> f a -> f b
<*> Getting a (GaussianState a) a -> m a
forall s (m :: * -> *) a. MonadState s m => Getting a s a -> m a
use Getting a (GaussianState a) a
forall a. Lens' (GaussianState a) a
detAcc
        Just ((Int
pivCol, Int
_), IntMap Int
_) -> do
          IntSet
heavys <- Getting IntSet (GaussianState a) IntSet -> m IntSet
forall s (m :: * -> *) a. MonadState s m => Getting a s a -> m a
use Getting IntSet (GaussianState a) IntSet
forall a. Lens' (GaussianState a) IntSet
heavyCols
          (Int -> Identity Int)
-> GaussianState a -> Identity (GaussianState a)
forall a. Lens' (GaussianState a) Int
prevCol ((Int -> Identity Int)
 -> GaussianState a -> Identity (GaussianState a))
-> Int -> m ()
forall s (m :: * -> *) a b.
MonadState s m =>
ASetter s s a b -> b -> m ()
.= Int
pivCol
          let trav :: Maybe (Entry a, Int) -> Int -> Entry a -> m (Maybe (Entry a, Int))
trav Maybe (Entry a, Int)
b Int
_ Entry a
ent = do
                if (Entry a
ent Entry a -> Getting Int (Entry a) Int -> Int
forall s a. s -> Getting a s a -> a
^. Direction -> Lens' (Entry a) Int
forall a. Direction -> Lens' (Entry a) Int
coordL Direction
Row) Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
destRow
                  then do
                  Maybe (Entry a, Int) -> m (Maybe (Entry a, Int))
forall (m :: * -> *) a. Monad m => a -> m a
return Maybe (Entry a, Int)
b
                  else do
                  let lc :: Int
lc = IntSet -> Int -> Matrix a -> Int
forall a. IntSet -> Int -> Matrix a -> Int
countLight IntSet
heavys (Entry a
ent Entry a -> Getting Int (Entry a) Int -> Int
forall s a. s -> Getting a s a -> a
^. Direction -> Lens' (Entry a) Int
forall a. Direction -> Lens' (Entry a) Int
coordL Direction
Row) Matrix a
old
                  Maybe (Entry a, Int) -> m (Maybe (Entry a, Int))
forall (m :: * -> *) a. Monad m => a -> m a
return (Maybe (Entry a, Int) -> m (Maybe (Entry a, Int)))
-> Maybe (Entry a, Int) -> m (Maybe (Entry a, Int))
forall a b. (a -> b) -> a -> b
$ case Maybe (Entry a, Int)
b of
                    Maybe (Entry a, Int)
Nothing -> (Entry a, Int) -> Maybe (Entry a, Int)
forall a. a -> Maybe a
Just (Entry a
ent, Int
lc)
                    Just (Entry a
p, Int
l0)
                      | Int
l0 Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
lc  -> (Entry a, Int) -> Maybe (Entry a, Int)
forall a. a -> Maybe a
Just (Entry a
p, Int
l0)
                      | Bool
otherwise -> (Entry a, Int) -> Maybe (Entry a, Int)
forall a. a -> Maybe a
Just (Entry a
ent, Int
lc)
          Maybe (Entry a, Int)
mans <- Maybe (Entry a, Int)
-> (Maybe (Entry a, Int)
    -> Int -> Entry a -> m (Maybe (Entry a, Int)))
-> Direction
-> Int
-> Matrix a
-> m (Maybe (Entry a, Int))
forall (m :: * -> *) b a.
Monad m =>
b
-> (b -> Int -> Entry a -> m b)
-> Direction
-> Int
-> Matrix a
-> m b
traverseDirM Maybe (Entry a, Int)
forall a. Maybe a
Nothing Maybe (Entry a, Int) -> Int -> Entry a -> m (Maybe (Entry a, Int))
trav Direction
Column Int
pivCol Matrix a
old
          case Maybe (Entry a, Int)
mans of
            Maybe (Entry a, Int)
Nothing -> m (Matrix a, Matrix a, a)
nextElim
            Just (Entry a
pivot, Int
_) -> do
              let pivRow :: Int
pivRow = Entry a
pivot Entry a -> Getting Int (Entry a) Int -> Int
forall s a. s -> Getting a s a -> a
^. Direction -> Lens' (Entry a) Int
forall a. Direction -> Lens' (Entry a) Int
coordL Direction
Row
                  pivCoe :: a
pivCoe = Entry a
pivot Entry a -> Getting a (Entry a) a -> a
forall s a. s -> Getting a s a -> a
^. Getting a (Entry a) a
forall a a. Lens (Entry a) (Entry a) a a
value
                  sgn :: a
sgn    = if Int
pivRow Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
destRow then a
forall r. Unital r => r
one else a -> a
forall r. Group r => r -> r
negate a
forall r. Unital r => r
one
              Matrix a
p0 <- Getting (Matrix a) (GaussianState a) (Matrix a) -> m (Matrix a)
forall s (m :: * -> *) a. MonadState s m => Getting a s a -> m a
use Getting (Matrix a) (GaussianState a) (Matrix a)
forall a. Lens' (GaussianState a) (Matrix a)
output
              let elim :: (Matrix a, Matrix a) -> Int -> Entry a -> m (Matrix a, Matrix a)
elim (Matrix a
m, Matrix a
p) Int
_ Entry a
ent = do
                    if Entry a
entEntry a -> Getting Int (Entry a) Int -> Int
forall s a. s -> Getting a s a -> a
^.Direction -> Lens' (Entry a) Int
forall a. Direction -> Lens' (Entry a) Int
coordL Direction
Row Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
pivRow
                      then do
                        let coe :: a
coe = a -> a
forall r. Group r => r -> r
negate (Entry a
ent Entry a -> Getting a (Entry a) a -> a
forall s a. s -> Getting a s a -> a
^. Getting a (Entry a) a
forall a a. Lens (Entry a) (Entry a) a a
value) a -> a -> a
forall r. Division r => r -> r -> r
/ a
pivCoe
                        (Matrix a, Matrix a) -> m (Matrix a, Matrix a)
forall (m :: * -> *) a. Monad m => a -> m a
return ((Matrix a, Matrix a) -> m (Matrix a, Matrix a))
-> (Matrix a, Matrix a) -> m (Matrix a, Matrix a)
forall a b. (a -> b) -> a -> b
$ (Matrix a
m, Matrix a
p) (Matrix a, Matrix a)
-> ((Matrix a, Matrix a) -> (Matrix a, Matrix a))
-> (Matrix a, Matrix a)
forall a b. a -> (a -> b) -> b
& (Matrix a -> Identity (Matrix a))
-> (Matrix a, Matrix a) -> Identity (Matrix a, Matrix a)
forall (r :: * -> * -> *) a b.
Bitraversable r =>
Traversal (r a a) (r b b) a b
both ((Matrix a -> Identity (Matrix a))
 -> (Matrix a, Matrix a) -> Identity (Matrix a, Matrix a))
-> (Matrix a -> Matrix a)
-> (Matrix a, Matrix a)
-> (Matrix a, Matrix a)
forall s t a b. ASetter s t a b -> (a -> b) -> s -> t
%~ a -> Int -> Int -> Matrix a -> Matrix a
forall a.
(DecidableZero a, Multiplicative a) =>
a -> Int -> Int -> Matrix a -> Matrix a
combineRows a
coe Int
pivRow (Entry a
ent Entry a -> Getting Int (Entry a) Int -> Int
forall s a. s -> Getting a s a -> a
^. Direction -> Lens' (Entry a) Int
forall a. Direction -> Lens' (Entry a) Int
coordL Direction
Row)
                      else do
                        (Matrix a, Matrix a) -> m (Matrix a, Matrix a)
forall (m :: * -> *) a. Monad m => a -> m a
return (Matrix a
m, Matrix a
p)
              (Matrix a
input', Matrix a
output') <- (Matrix a, Matrix a)
-> ((Matrix a, Matrix a)
    -> Int -> Entry a -> m (Matrix a, Matrix a))
-> Direction
-> Int
-> Matrix a
-> m (Matrix a, Matrix a)
forall (m :: * -> *) b a.
Monad m =>
b
-> (b -> Int -> Entry a -> m b)
-> Direction
-> Int
-> Matrix a
-> m b
traverseDirM (Matrix a
old, Matrix a
p0) (Matrix a, Matrix a) -> Int -> Entry a -> m (Matrix a, Matrix a)
elim Direction
Column Int
pivCol Matrix a
old
                    m (Matrix a, Matrix a)
-> ((Matrix a, Matrix a) -> (Matrix a, Matrix a))
-> m (Matrix a, Matrix a)
forall (f :: * -> *) a b. Functor f => f a -> (a -> b) -> f b
<&> (Matrix a -> Identity (Matrix a))
-> (Matrix a, Matrix a) -> Identity (Matrix a, Matrix a)
forall (r :: * -> * -> *) a b.
Bitraversable r =>
Traversal (r a a) (r b b) a b
both ((Matrix a -> Identity (Matrix a))
 -> (Matrix a, Matrix a) -> Identity (Matrix a, Matrix a))
-> (Matrix a -> Matrix a)
-> (Matrix a, Matrix a)
-> (Matrix a, Matrix a)
forall s t a b. ASetter s t a b -> (a -> b) -> s -> t
%~ a -> Int -> Matrix a -> Matrix a
forall a.
(DecidableZero a, Multiplicative a) =>
a -> Int -> Matrix a -> Matrix a
scaleRow (a -> a
forall r. Division r => r -> r
recip a
pivCoe) Int
destRow (Matrix a -> Matrix a)
-> (Matrix a -> Matrix a) -> Matrix a -> Matrix a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> Int -> Matrix a -> Matrix a
forall a. Int -> Int -> Matrix a -> Matrix a
switchRows Int
destRow Int
pivRow
              (Matrix a -> Identity (Matrix a))
-> GaussianState a -> Identity (GaussianState a)
forall a. Lens' (GaussianState a) (Matrix a)
input ((Matrix a -> Identity (Matrix a))
 -> GaussianState a -> Identity (GaussianState a))
-> Matrix a -> m ()
forall s (m :: * -> *) a b.
MonadState s m =>
ASetter s s a b -> b -> m ()
.= Matrix a
input'
              (Matrix a -> Identity (Matrix a))
-> GaussianState a -> Identity (GaussianState a)
forall a. Lens' (GaussianState a) (Matrix a)
output ((Matrix a -> Identity (Matrix a))
 -> GaussianState a -> Identity (GaussianState a))
-> Matrix a -> m ()
forall s (m :: * -> *) a b.
MonadState s m =>
ASetter s s a b -> b -> m ()
.= Matrix a
output'
              (Int -> Identity Int)
-> GaussianState a -> Identity (GaussianState a)
forall a. Lens' (GaussianState a) Int
curRow ((Int -> Identity Int)
 -> GaussianState a -> Identity (GaussianState a))
-> Int -> m ()
forall s (m :: * -> *) a.
(MonadState s m, Num a) =>
ASetter' s a -> a -> m ()
+= Int
1
              (a -> Identity a) -> GaussianState a -> Identity (GaussianState a)
forall a. Lens' (GaussianState a) a
detAcc ((a -> Identity a)
 -> GaussianState a -> Identity (GaussianState a))
-> (a -> a) -> m ()
forall s (m :: * -> *) a b.
MonadState s m =>
ASetter s s a b -> (a -> b) -> m ()
%= (a -> a -> a
forall r. Multiplicative r => r -> r -> r
*(a
pivCoea -> a -> a
forall r. Multiplicative r => r -> r -> r
*a
sgn))
              m (Matrix a, Matrix a, a)
nextElim
    nextElim :: m (Matrix a, Matrix a, a)
nextElim = do
      IntSet
oldHeavys <- Getting IntSet (GaussianState a) IntSet -> m IntSet
forall s (m :: * -> *) a. MonadState s m => Getting a s a -> m a
use Getting IntSet (GaussianState a) IntSet
forall a. Lens' (GaussianState a) IntSet
heavyCols
      IntSet
newHeavyCols <- LensLike' (Const IntSet) (GaussianState a) (Matrix a)
-> (Matrix a -> IntSet) -> m IntSet
forall s (m :: * -> *) r a.
MonadState s m =>
LensLike' (Const r) s a -> (a -> r) -> m r
uses LensLike' (Const IntSet) (GaussianState a) (Matrix a)
forall a. Lens' (GaussianState a) (Matrix a)
input (IntSet -> Matrix a -> IntSet
forall a. IntSet -> Matrix a -> IntSet
getHeaviest IntSet
oldHeavys)
      (IntSet -> Identity IntSet)
-> GaussianState a -> Identity (GaussianState a)
forall a. Lens' (GaussianState a) IntSet
heavyCols ((IntSet -> Identity IntSet)
 -> GaussianState a -> Identity (GaussianState a))
-> (IntSet -> IntSet) -> m ()
forall s (m :: * -> *) a b.
MonadState s m =>
ASetter s s a b -> (a -> b) -> m ()
%= IntSet -> IntSet -> IntSet
IS.union IntSet
newHeavyCols
      m (Matrix a, Matrix a, a)
go

nonZeroEntries :: Matrix a -> Vector ((Int, Int), a)
nonZeroEntries :: Matrix a -> Vector ((Int, Int), a)
nonZeroEntries Matrix a
mat = (Entry a -> ((Int, Int), a))
-> Vector (Entry a) -> Vector ((Int, Int), a)
forall a b. (a -> b) -> Vector a -> Vector b
V.map (Getting (Int, Int) (Entry a) (Int, Int) -> Entry a -> (Int, Int)
forall s (m :: * -> *) a. MonadReader s m => Getting a s a -> m a
view Getting (Int, Int) (Entry a) (Int, Int)
forall a. Lens' (Entry a) (Int, Int)
idx (Entry a -> (Int, Int))
-> (Entry a -> a) -> Entry a -> ((Int, Int), a)
forall (a :: * -> * -> *) b c c'.
Arrow a =>
a b c -> a b c' -> a b (c, c')
&&& Getting a (Entry a) a -> Entry a -> a
forall s (m :: * -> *) a. MonadReader s m => Getting a s a -> m a
view Getting a (Entry a) a
forall a a. Lens (Entry a) (Entry a) a a
value) (Vector (Entry a) -> Vector ((Int, Int), a))
-> Vector (Entry a) -> Vector ((Int, Int), a)
forall a b. (a -> b) -> a -> b
$ Matrix a
mat Matrix a
-> Getting (Vector (Entry a)) (Matrix a) (Vector (Entry a))
-> Vector (Entry a)
forall s a. s -> Getting a s a -> a
^. Getting (Vector (Entry a)) (Matrix a) (Vector (Entry a))
forall a a.
Lens (Matrix a) (Matrix a) (Vector (Entry a)) (Vector (Entry a))
coefficients

multWithVector :: (Multiplicative a, Monoidal a)
               => Matrix a -> Vector a -> Vector a
multWithVector :: Matrix a -> Vector a -> Vector a
multWithVector Matrix a
mat Vector a
v =
  Int -> (Int -> a) -> Vector a
forall a. Int -> (Int -> a) -> Vector a
V.generate (Matrix a
matMatrix a -> Getting Int (Matrix a) Int -> Int
forall s a. s -> Getting a s a -> a
^.Getting Int (Matrix a) Int
forall a. Lens' (Matrix a) Int
height) ((Int -> a) -> Vector a) -> (Int -> a) -> Vector a
forall a b. (a -> b) -> a -> b
$ \Int
i ->
  a -> (a -> Int -> Entry a -> a) -> Int -> Matrix a -> a
forall b a. b -> (b -> Int -> Entry a -> b) -> Int -> Matrix a -> b
traverseRow a
forall m. Monoidal m => m
zero (\a
acc Int
_ Entry a
ent -> a
acc a -> a -> a
forall r. Additive r => r -> r -> r
+ (Entry a
entEntry a -> Getting a (Entry a) a -> a
forall s a. s -> Getting a s a -> a
^.Getting a (Entry a) a
forall a a. Lens (Entry a) (Entry a) a a
value)a -> a -> a
forall r. Multiplicative r => r -> r -> r
*(Vector a
v Vector a -> Int -> a
forall a. Vector a -> Int -> a
V.! (Entry a
entEntry a -> Getting Int (Entry a) Int -> Int
forall s a. s -> Getting a s a -> a
^.Direction -> Lens' (Entry a) Int
forall a. Direction -> Lens' (Entry a) Int
nthL Direction
Row))) Int
i Matrix a
mat

nonZeroDirs :: Direction -> Matrix r -> [Int]
nonZeroDirs :: Direction -> Matrix r -> [Int]
nonZeroDirs Direction
dir = Getting [Int] (Matrix r) [Int] -> Matrix r -> [Int]
forall s (m :: * -> *) a. MonadReader s m => Getting a s a -> m a
view (Getting [Int] (Matrix r) [Int] -> Matrix r -> [Int])
-> Getting [Int] (Matrix r) [Int] -> Matrix r -> [Int]
forall a b. (a -> b) -> a -> b
$ Direction -> Lens' (Matrix r) (IntMap Int)
forall a. Direction -> Lens' (Matrix a) (IntMap Int)
startL Direction
dir ((IntMap Int -> Const [Int] (IntMap Int))
 -> Matrix r -> Const [Int] (Matrix r))
-> (([Int] -> Const [Int] [Int])
    -> IntMap Int -> Const [Int] (IntMap Int))
-> Getting [Int] (Matrix r) [Int]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (IntMap Int -> [Int])
-> ([Int] -> Const [Int] [Int])
-> IntMap Int
-> Const [Int] (IntMap Int)
forall (p :: * -> * -> *) (f :: * -> *) s a.
(Profunctor p, Contravariant f) =>
(s -> a) -> Optic' p f s a
to IntMap Int -> [Int]
forall a. IntMap a -> [Int]
IM.keys

nonZeroRows :: Matrix r -> [Int]
nonZeroRows :: Matrix r -> [Int]
nonZeroRows = Direction -> Matrix r -> [Int]
forall r. Direction -> Matrix r -> [Int]
nonZeroDirs Direction
Row

nonZeroCols :: Matrix r -> [Int]
nonZeroCols :: Matrix r -> [Int]
nonZeroCols = Direction -> Matrix r -> [Int]
forall r. Direction -> Matrix r -> [Int]
nonZeroDirs Direction
Column

{-
testCase :: Matrix (Fraction Integer)
testCase = fromLists [[0,0,0,0,0,0,2,-3,-1,0]
                      ,[0,0,0,2,-3,-1,0,0,0,0]
                      ,[0,2,-3,0,-1,0,0,0,0,0]
                      ,[1,0,1,0,0,1,-2,0,0,0]
                      ,[2,-3,0,-1,0,0,0,0,0,0]
                      ,[1,0,1,0,0,1,0,0,0,-1]
                      ,[1,0,1,0,0,1,-2,0,0,0]]
-}

newtype Square n r = Square { Square n r -> Matrix r
runSquare :: Matrix r
                            } deriving (Int -> Square n r -> ShowS
[Square n r] -> ShowS
Square n r -> String
(Int -> Square n r -> ShowS)
-> (Square n r -> String)
-> ([Square n r] -> ShowS)
-> Show (Square n r)
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
forall k (n :: k) r. Show r => Int -> Square n r -> ShowS
forall k (n :: k) r. Show r => [Square n r] -> ShowS
forall k (n :: k) r. Show r => Square n r -> String
showList :: [Square n r] -> ShowS
$cshowList :: forall k (n :: k) r. Show r => [Square n r] -> ShowS
show :: Square n r -> String
$cshow :: forall k (n :: k) r. Show r => Square n r -> String
showsPrec :: Int -> Square n r -> ShowS
$cshowsPrec :: forall k (n :: k) r. Show r => Int -> Square n r -> ShowS
Show, Square n r -> Square n r -> Bool
(Square n r -> Square n r -> Bool)
-> (Square n r -> Square n r -> Bool) -> Eq (Square n r)
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
forall k (n :: k) r. Eq r => Square n r -> Square n r -> Bool
/= :: Square n r -> Square n r -> Bool
$c/= :: forall k (n :: k) r. Eq r => Square n r -> Square n r -> Bool
== :: Square n r -> Square n r -> Bool
$c== :: forall k (n :: k) r. Eq r => Square n r -> Square n r -> Bool
Eq, Natural -> Square n r -> Square n r
Square n r -> Square n r -> Square n r
(a -> Square n r) -> f a -> Square n r
(Square n r -> Square n r -> Square n r)
-> (Natural -> Square n r -> Square n r)
-> (forall (f :: * -> *) a.
    Foldable1 f =>
    (a -> Square n r) -> f a -> Square n r)
-> Additive (Square n r)
forall r.
(r -> r -> r)
-> (Natural -> r -> r)
-> (forall (f :: * -> *) a. Foldable1 f => (a -> r) -> f a -> r)
-> Additive r
forall k (n :: k) r.
DecidableZero r =>
Natural -> Square n r -> Square n r
forall k (n :: k) r.
DecidableZero r =>
Square n r -> Square n r -> Square n r
forall k (n :: k) r (f :: * -> *) a.
(DecidableZero r, Foldable1 f) =>
(a -> Square n r) -> f a -> Square n r
forall (f :: * -> *) a.
Foldable1 f =>
(a -> Square n r) -> f a -> Square n r
sumWith1 :: (a -> Square n r) -> f a -> Square n r
$csumWith1 :: forall k (n :: k) r (f :: * -> *) a.
(DecidableZero r, Foldable1 f) =>
(a -> Square n r) -> f a -> Square n r
sinnum1p :: Natural -> Square n r -> Square n r
$csinnum1p :: forall k (n :: k) r.
DecidableZero r =>
Natural -> Square n r -> Square n r
+ :: Square n r -> Square n r -> Square n r
$c+ :: forall k (n :: k) r.
DecidableZero r =>
Square n r -> Square n r -> Square n r
Additive, Square n r -> Natural -> Square n r
Square n r -> Square n r -> Square n r
(a -> Square n r) -> f a -> Square n r
(Square n r -> Square n r -> Square n r)
-> (Square n r -> Natural -> Square n r)
-> (forall (f :: * -> *) a.
    Foldable1 f =>
    (a -> Square n r) -> f a -> Square n r)
-> Multiplicative (Square n r)
forall r.
(r -> r -> r)
-> (r -> Natural -> r)
-> (forall (f :: * -> *) a. Foldable1 f => (a -> r) -> f a -> r)
-> Multiplicative r
forall k (n :: k) r.
(DecidableZero r, Multiplicative r) =>
Square n r -> Natural -> Square n r
forall k (n :: k) r.
(DecidableZero r, Multiplicative r) =>
Square n r -> Square n r -> Square n r
forall k (n :: k) r (f :: * -> *) a.
(DecidableZero r, Multiplicative r, Foldable1 f) =>
(a -> Square n r) -> f a -> Square n r
forall (f :: * -> *) a.
Foldable1 f =>
(a -> Square n r) -> f a -> Square n r
productWith1 :: (a -> Square n r) -> f a -> Square n r
$cproductWith1 :: forall k (n :: k) r (f :: * -> *) a.
(DecidableZero r, Multiplicative r, Foldable1 f) =>
(a -> Square n r) -> f a -> Square n r
pow1p :: Square n r -> Natural -> Square n r
$cpow1p :: forall k (n :: k) r.
(DecidableZero r, Multiplicative r) =>
Square n r -> Natural -> Square n r
* :: Square n r -> Square n r -> Square n r
$c* :: forall k (n :: k) r.
(DecidableZero r, Multiplicative r) =>
Square n r -> Square n r -> Square n r
Multiplicative)

deriving instance (DecidableZero r, Semiring r, Multiplicative r)
               => LeftModule (Scalar r) (Square n r)
deriving instance (DecidableZero r, Semiring r, Multiplicative r)
               => RightModule (Scalar r) (Square n r)

instance (Unital r, Multiplicative r, Reifies n Integer, DecidableZero r) => Unital (Square n r) where
  one :: Square n r
one = Matrix r -> Square n r
forall k (n :: k) r. Matrix r -> Square n r
Square (Matrix r -> Square n r) -> Matrix r -> Square n r
forall a b. (a -> b) -> a -> b
$ Int -> Matrix r
forall a. Unital a => Int -> Matrix a
identity (Int -> Matrix r) -> Int -> Matrix r
forall a b. (a -> b) -> a -> b
$ Integer -> Int
forall r. Num r => Integer -> r
fromInteger (Integer -> Int) -> Integer -> Int
forall a b. (a -> b) -> a -> b
$ Proxy n -> Integer
forall k (s :: k) a (proxy :: k -> *). Reifies s a => proxy s -> a
reflect (Proxy n
forall k (t :: k). Proxy t
Proxy :: Proxy n)

instance (DecidableZero r, Multiplicative r) => Multiplicative (Matrix r) where
  Matrix r
m * :: Matrix r -> Matrix r -> Matrix r
* Matrix r
n = [((Int, Int), r)] -> Matrix r
forall a. DecidableZero a => [((Int, Int), a)] -> Matrix a
fromList [ ((Int
i,Int
j),Vector r -> r
forall (f :: * -> *) m. (Foldable f, Monoidal m) => f m -> m
sum (Vector r -> r) -> Vector r -> r
forall a b. (a -> b) -> a -> b
$ (r -> r -> r) -> Vector r -> Vector r -> Vector r
forall a b c. (a -> b -> c) -> Vector a -> Vector b -> Vector c
V.zipWith r -> r -> r
forall r. Multiplicative r => r -> r -> r
(*) (Int -> Matrix r -> Vector r
forall a. Monoidal a => Int -> Matrix a -> Vector a
getRow Int
i Matrix r
m) (Int -> Matrix r -> Vector r
forall a. Monoidal a => Int -> Matrix a -> Vector a
getCol Int
j Matrix r
n))
                   | Int
i <- Matrix r -> [Int]
forall r. Matrix r -> [Int]
nonZeroRows Matrix r
m
                   , Int
j <- Matrix r -> [Int]
forall r. Matrix r -> [Int]
nonZeroCols Matrix r
n
                   ] Matrix r -> (Matrix r -> Matrix r) -> Matrix r
forall a b. a -> (a -> b) -> b
& (Int -> Identity Int) -> Matrix r -> Identity (Matrix r)
forall a. Lens' (Matrix a) Int
width ((Int -> Identity Int) -> Matrix r -> Identity (Matrix r))
-> Int -> Matrix r -> Matrix r
forall s t a b. ASetter s t a b -> b -> s -> t
.~ Matrix r
nMatrix r -> Getting Int (Matrix r) Int -> Int
forall s a. s -> Getting a s a -> a
^.Getting Int (Matrix r) Int
forall a. Lens' (Matrix a) Int
width
                     Matrix r -> (Matrix r -> Matrix r) -> Matrix r
forall a b. a -> (a -> b) -> b
& (Int -> Identity Int) -> Matrix r -> Identity (Matrix r)
forall a. Lens' (Matrix a) Int
height ((Int -> Identity Int) -> Matrix r -> Identity (Matrix r))
-> Int -> Matrix r -> Matrix r
forall s t a b. ASetter s t a b -> b -> s -> t
.~ Matrix r
mMatrix r -> Getting Int (Matrix r) Int -> Int
forall s a. s -> Getting a s a -> a
^.Getting Int (Matrix r) Int
forall a. Lens' (Matrix a) Int
height


instance (DecidableZero r, RightModule Natural r) => RightModule Natural (Matrix r) where
  Matrix r
m *. :: Matrix r -> Natural -> Matrix r
*. Natural
n = (r -> r) -> Matrix r -> Matrix r
forall a a1. DecidableZero a => (a1 -> a) -> Matrix a1 -> Matrix a
cmap (r -> Natural -> r
forall r m. RightModule r m => m -> r -> m
*. Natural
n) Matrix r
m

instance (DecidableZero r, LeftModule Natural r) => LeftModule Natural (Matrix r) where
  Natural
n .* :: Natural -> Matrix r -> Matrix r
.* Matrix r
m = (r -> r) -> Matrix r -> Matrix r
forall a a1. DecidableZero a => (a1 -> a) -> Matrix a1 -> Matrix a
cmap (Natural
n Natural -> r -> r
forall r m. LeftModule r m => r -> m -> m
.*) Matrix r
m

instance (DecidableZero r, RightModule Integer r) => RightModule Integer (Matrix r) where
  Matrix r
m *. :: Matrix r -> Integer -> Matrix r
*. Integer
n = (r -> r) -> Matrix r -> Matrix r
forall a a1. DecidableZero a => (a1 -> a) -> Matrix a1 -> Matrix a
cmap (r -> Integer -> r
forall r m. RightModule r m => m -> r -> m
*. Integer
n) Matrix r
m

instance (DecidableZero r, LeftModule Integer r) => LeftModule Integer (Matrix r) where
  Integer
n .* :: Integer -> Matrix r -> Matrix r
.* Matrix r
m = (r -> r) -> Matrix r -> Matrix r
forall a a1. DecidableZero a => (a1 -> a) -> Matrix a1 -> Matrix a
cmap (Integer
n Integer -> r -> r
forall r m. LeftModule r m => r -> m -> m
.*) Matrix r
m

instance (DecidableZero r)
      => Monoidal (Matrix r) where
  zero :: Matrix r
zero = Int -> Int -> Matrix r
forall a. Int -> Int -> Matrix a
zeroMat Int
0 Int
0

instance (DecidableZero r) => Additive (Matrix r) where
  Matrix r
m + :: Matrix r -> Matrix r -> Matrix r
+ Matrix r
n =
    let dir :: Direction
dir = (Direction -> Direction -> Ordering) -> [Direction] -> Direction
forall (t :: * -> *) a.
Foldable t =>
(a -> a -> Ordering) -> t a -> a
minimumBy ((Direction -> Int) -> Direction -> Direction -> Ordering
forall a b. Ord a => (b -> a) -> b -> b -> Ordering
comparing ((Direction -> Int) -> Direction -> Direction -> Ordering)
-> (Direction -> Int) -> Direction -> Direction -> Ordering
forall a b. (a -> b) -> a -> b
$ [Int] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ([Int] -> Int) -> (Direction -> [Int]) -> Direction -> Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Direction -> Matrix r -> [Int]) -> Matrix r -> Direction -> [Int]
forall a b c. (a -> b -> c) -> b -> a -> c
flip Direction -> Matrix r -> [Int]
forall r. Direction -> Matrix r -> [Int]
nonZeroDirs Matrix r
n)
              [Direction
Row, Direction
Column]
    in (Int -> Matrix r -> Matrix r) -> Matrix r -> [Int] -> Matrix r
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr (\Int
i Matrix r
l -> Direction -> Vector r -> Int -> Matrix r -> Matrix r
forall a.
DecidableZero a =>
Direction -> Vector a -> Int -> Matrix a -> Matrix a
addDir Direction
dir (Direction -> Int -> Matrix r -> Vector r
forall a. Monoidal a => Direction -> Int -> Matrix a -> Vector a
getDir Direction
dir Int
i Matrix r
n) Int
i Matrix r
l) Matrix r
m (Direction -> Matrix r -> [Int]
forall r. Direction -> Matrix r -> [Int]
nonZeroDirs Direction
dir Matrix r
n)

instance (DecidableZero r, Semiring r, Multiplicative r)
      => LeftModule (Scalar r) (Matrix r) where
  Scalar r
r .* :: Scalar r -> Matrix r -> Matrix r
.* Matrix r
mat = (r -> r) -> Matrix r -> Matrix r
forall a a1. DecidableZero a => (a1 -> a) -> Matrix a1 -> Matrix a
cmap (r
rr -> r -> r
forall r. Multiplicative r => r -> r -> r
*) Matrix r
mat

instance (DecidableZero r, Semiring r, Multiplicative r)
      => RightModule (Scalar r) (Matrix r) where
  Matrix r
mat *. :: Matrix r -> Scalar r -> Matrix r
*. Scalar r
r = (r -> r) -> Matrix r -> Matrix r
forall a a1. DecidableZero a => (a1 -> a) -> Matrix a1 -> Matrix a
cmap (r -> r -> r
forall r. Multiplicative r => r -> r -> r
*r
r) Matrix r
mat

instance (DecidableZero r, Group r) => Group (Matrix r) where
  negate :: Matrix r -> Matrix r
negate = (r -> r) -> Matrix r -> Matrix r
forall a a1. DecidableZero a => (a1 -> a) -> Matrix a1 -> Matrix a
cmap r -> r
forall r. Group r => r -> r
negate

instance (DecidableZero r, Abelian r) => Abelian (Matrix r)

instance (DecidableZero r, Semiring r) => Semiring (Matrix r)

substMatrix :: (CoeffRing r)
            => Matrix r -> Polynomial r 1 -> Matrix r
substMatrix :: Matrix r -> Polynomial r 1 -> Matrix r
substMatrix Matrix r
m Polynomial r 1
f =
  let n :: Int
n = Matrix r -> Int
forall a. Matrix a -> Int
ncols Matrix r
m
  in if Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Matrix r -> Int
forall a. Matrix a -> Int
nrows Matrix r
m
     then Integer
-> (forall s. Reifies s Integer => Proxy s -> Matrix r) -> Matrix r
forall a r. a -> (forall s. Reifies s a => Proxy s -> r) -> r
reify (Int -> Integer
forall a. Integral a => a -> Integer
P.toInteger Int
n) ((forall s. Reifies s Integer => Proxy s -> Matrix r) -> Matrix r)
-> (forall s. Reifies s Integer => Proxy s -> Matrix r) -> Matrix r
forall a b. (a -> b) -> a -> b
$ \Proxy s
pxy -> Square s r -> Matrix r
forall k (n :: k) r. Square n r -> Matrix r
runSquare (Square s r -> Matrix r) -> Square s r -> Matrix r
forall a b. (a -> b) -> a -> b
$ Square s r -> Polynomial r 1 -> Square s r
forall r b order.
(Module (Scalar r) b, Unital b, CoeffRing r,
 IsMonomialOrder 1 order) =>
b -> OrderedPolynomial r order 1 -> b
substUnivariate (Proxy s -> Matrix r -> Square s r
forall k (proxy :: k -> *) (n :: k) r.
proxy n -> Matrix r -> Square n r
toSquare Proxy s
pxy Matrix r
m) Polynomial r 1
f
     else String -> Matrix r
forall a. HasCallStack => String -> a
error String
"Matrix must be square"

toSquare :: proxy n -> Matrix r -> Square n r
toSquare :: proxy n -> Matrix r -> Square n r
toSquare proxy n
_ = Matrix r -> Square n r
forall k (n :: k) r. Matrix r -> Square n r
Square

(<.>) :: (Multiplicative m, Monoidal m) => Vector m -> Vector m -> m
Vector m
v <.> :: Vector m -> Vector m -> m
<.> Vector m
u = Vector m -> m
forall (f :: * -> *) m. (Foldable f, Monoidal m) => f m -> m
sum (Vector m -> m) -> Vector m -> m
forall a b. (a -> b) -> a -> b
$ (m -> m -> m) -> Vector m -> Vector m -> Vector m
forall a b c. (a -> b -> c) -> Vector a -> Vector b -> Vector c
V.zipWith m -> m -> m
forall r. Multiplicative r => r -> r -> r
(*) Vector m
v Vector m
u

krylovMinpol :: (Eq a, Ring a, DecidableZero a, DecidableUnits a,
                 Field a, ZeroProductSemiring a,
                 Random a, MonadRandom m)
             => Matrix a -> Vector a -> m (Polynomial a 1)
krylovMinpol :: Matrix a -> Vector a -> m (Polynomial a 1)
krylovMinpol Matrix a
m Vector a
b
  | (a -> Bool) -> Vector a -> Bool
forall a. (a -> Bool) -> Vector a -> Bool
V.all a -> Bool
forall r. DecidableZero r => r -> Bool
isZero Vector a
b = Polynomial a 1 -> m (Polynomial a 1)
forall (m :: * -> *) a. Monad m => a -> m a
return Polynomial a 1
forall r. Unital r => r
one
  | Bool
otherwise = Integer
-> (forall s. Reifies s Integer => Proxy s -> m (Polynomial a 1))
-> m (Polynomial a 1)
forall a r. a -> (forall s. Reifies s a => Proxy s -> r) -> r
reify (Int -> Integer
forall a. Integral a => a -> Integer
P.toInteger Int
n) ((forall s. Reifies s Integer => Proxy s -> m (Polynomial a 1))
 -> m (Polynomial a 1))
-> (forall s. Reifies s Integer => Proxy s -> m (Polynomial a 1))
-> m (Polynomial a 1)
forall a b. (a -> b) -> a -> b
$ \Proxy s
pxy -> do
    (Polynomial a 1 -> Bool)
-> m (Polynomial a 1) -> m (Polynomial a 1)
forall (m :: * -> *) a. Monad m => (a -> Bool) -> m a -> m a
iterateUntil (\Polynomial a 1
h -> (a -> Bool) -> Vector a -> Bool
forall a. (a -> Bool) -> Vector a -> Bool
V.all a -> Bool
forall r. DecidableZero r => r -> Bool
isZero (Vector a -> Bool) -> Vector a -> Bool
forall a b. (a -> b) -> a -> b
$ Matrix a -> Vector a -> Vector a
forall a.
(Multiplicative a, Monoidal a) =>
Matrix a -> Vector a -> Vector a
multWithVector (Matrix a -> Polynomial a 1 -> Matrix a
forall r. CoeffRing r => Matrix r -> Polynomial r 1 -> Matrix r
substMatrix Matrix a
m Polynomial a 1
h) Vector a
b) (m (Polynomial a 1) -> m (Polynomial a 1))
-> m (Polynomial a 1) -> m (Polynomial a 1)
forall a b. (a -> b) -> a -> b
$ do
      [a]
u <- Int -> m a -> m [a]
forall (m :: * -> *) a. Applicative m => Int -> m a -> m [a]
replicateM Int
n m a
forall (m :: * -> *) a. (MonadRandom m, Random a) => m a
getRandom
      Polynomial a 1 -> m (Polynomial a 1)
forall (m :: * -> *) a. Monad m => a -> m a
return (Polynomial a 1 -> m (Polynomial a 1))
-> Polynomial a 1 -> m (Polynomial a 1)
forall a b. (a -> b) -> a -> b
$ Natural -> [a] -> Polynomial a 1
forall k.
(Eq k, ZeroProductSemiring k, DecidableUnits k, DecidableZero k,
 Field k) =>
Natural -> [k] -> Polynomial k 1
minpolRecurrent (Int -> Natural
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n)
        [ [a] -> Vector a
forall a. [a] -> Vector a
V.fromList [a]
u Vector a -> Vector a -> a
forall m.
(Multiplicative m, Monoidal m) =>
Vector m -> Vector m -> m
<.> Matrix a -> Vector a -> Vector a
forall a.
(Multiplicative a, Monoidal a) =>
Matrix a -> Vector a -> Vector a
multWithVector (Square s a -> Matrix a
forall k (n :: k) r. Square n r -> Matrix r
runSquare (Square s a -> Matrix a) -> Square s a -> Matrix a
forall a b. (a -> b) -> a -> b
$ Proxy s -> Matrix a -> Square s a
forall k (proxy :: k -> *) (n :: k) r.
proxy n -> Matrix r -> Square n r
toSquare Proxy s
pxy Matrix a
m Square s a -> Natural -> Square s a
forall r. Unital r => r -> Natural -> r
^ Int -> Natural
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
i) Vector a
b
        | Int
i <- [Int
0..Int
2Int -> Int -> Int
forall r. Multiplicative r => r -> r -> r
*Int
nInt -> Int -> Int
forall r. Group r => r -> r -> r
-Int
1]]
    where
      n :: Int
n = Matrix a -> Int
forall a. Matrix a -> Int
ncols Matrix a
m

-- | Solving linear equation using linearly recurrent sequence (Wiedemann algorithm).
solveWiedemann :: (Eq a, Field a, DecidableZero a, DecidableUnits a,
                ZeroProductSemiring a, Random a, MonadRandom m)
            => Matrix a -> Vector a -> m (Either (Vector a) (Vector a))
solveWiedemann :: Matrix a -> Vector a -> m (Either (Vector a) (Vector a))
solveWiedemann Matrix a
a Vector a
b = do
  Polynomial a 1
m <- Matrix a -> Vector a -> m (Polynomial a 1)
forall a (m :: * -> *).
(Eq a, Ring a, DecidableZero a, DecidableUnits a, Field a,
 ZeroProductSemiring a, Random a, MonadRandom m) =>
Matrix a -> Vector a -> m (Polynomial a 1)
krylovMinpol Matrix a
a Vector a
b
  Either (Vector a) (Vector a) -> m (Either (Vector a) (Vector a))
forall (m :: * -> *) a. Monad m => a -> m a
return (Either (Vector a) (Vector a) -> m (Either (Vector a) (Vector a)))
-> Either (Vector a) (Vector a) -> m (Either (Vector a) (Vector a))
forall a b. (a -> b) -> a -> b
$
    let m0 :: Polynomial a 1
m0 = Coefficient (Polynomial a 1) -> Polynomial a 1
forall poly. IsPolynomial poly => Coefficient poly -> poly
injectCoeff (OrderedMonomial (MOrder (Polynomial a 1)) (Arity (Polynomial a 1))
-> Polynomial a 1 -> Coefficient (Polynomial a 1)
forall poly.
IsOrderedPolynomial poly =>
OrderedMonomial (MOrder poly) (Arity poly)
-> poly -> Coefficient poly
coeff OrderedMonomial (MOrder (Polynomial a 1)) (Arity (Polynomial a 1))
forall r. Unital r => r
one Polynomial a 1
m)
        g :: Polynomial a 1
g = (Polynomial a 1
m Polynomial a 1 -> Polynomial a 1 -> Polynomial a 1
forall r. Group r => r -> r -> r
- Polynomial a 1
m0) Polynomial a 1 -> Polynomial a 1 -> Polynomial a 1
forall d. Euclidean d => d -> d -> d
`quot` Polynomial a 1
forall r (n :: Nat) order.
(CoeffRing r, KnownNat n, IsMonomialOrder n order, 0 < n) =>
OrderedPolynomial r order n
varX
    in if a -> Bool
forall r. DecidableZero r => r -> Bool
isZero (OrderedMonomial (MOrder (Polynomial a 1)) (Arity (Polynomial a 1))
-> Polynomial a 1 -> Coefficient (Polynomial a 1)
forall poly.
IsOrderedPolynomial poly =>
OrderedMonomial (MOrder poly) (Arity poly)
-> poly -> Coefficient poly
coeff OrderedMonomial (MOrder (Polynomial a 1)) (Arity (Polynomial a 1))
forall r. Unital r => r
one Polynomial a 1
m)
       then Vector a -> Either (Vector a) (Vector a)
forall a b. a -> Either a b
Left (Vector a -> Either (Vector a) (Vector a))
-> Vector a -> Either (Vector a) (Vector a)
forall a b. (a -> b) -> a -> b
$ Matrix a -> Polynomial a 1 -> Matrix a
forall r. CoeffRing r => Matrix r -> Polynomial r 1 -> Matrix r
substMatrix Matrix a
a Polynomial a 1
g Matrix a -> Vector a -> Vector a
forall a.
(Multiplicative a, Monoidal a) =>
Matrix a -> Vector a -> Vector a
`multWithVector` Vector a
b
       else let h :: Polynomial a 1
h = Polynomial a 1 -> Polynomial a 1
forall r. Group r => r -> r
negate Polynomial a 1
g Polynomial a 1 -> Polynomial a 1 -> Polynomial a 1
forall d. Euclidean d => d -> d -> d
`quot` Polynomial a 1
m0
            in Vector a -> Either (Vector a) (Vector a)
forall a b. b -> Either a b
Right (Vector a -> Either (Vector a) (Vector a))
-> Vector a -> Either (Vector a) (Vector a)
forall a b. (a -> b) -> a -> b
$ Matrix a -> Polynomial a 1 -> Matrix a
forall r. CoeffRing r => Matrix r -> Polynomial r 1 -> Matrix r
substMatrix Matrix a
a Polynomial a 1
h Matrix a -> Vector a -> Vector a
forall a.
(Multiplicative a, Monoidal a) =>
Matrix a -> Vector a -> Vector a
`multWithVector` Vector a
b

rankLM :: (DecidableZero r, Division r, Group r) => Matrix r -> Int
rankLM :: Matrix r -> Int
rankLM Matrix r
mat =
  let m' :: Matrix r
m' = (Matrix r, Matrix r) -> Matrix r
forall a b. (a, b) -> a
fst ((Matrix r, Matrix r) -> Matrix r)
-> (Matrix r, Matrix r) -> Matrix r
forall a b. (a -> b) -> a -> b
$ Matrix r -> (Matrix r, Matrix r)
forall a.
(DecidableZero a, Division a, Group a) =>
Matrix a -> (Matrix a, Matrix a)
structuredGauss Matrix r
mat
  in Int -> Int -> Int
forall a. Ord a => a -> a -> a
min ([Int] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ([Int] -> Int) -> [Int] -> Int
forall a b. (a -> b) -> a -> b
$ Matrix r -> [Int]
forall r. Matrix r -> [Int]
nonZeroRows Matrix r
m') ([Int] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ([Int] -> Int) -> [Int] -> Int
forall a b. (a -> b) -> a -> b
$ Matrix r -> [Int]
forall r. Matrix r -> [Int]
nonZeroCols Matrix r
m')

splitIndependentDirs :: (DecidableZero a, Field a)
                     => Direction -> Matrix a
                     -> (Matrix a, [Int], [Int])
                     -- ^ @(m', bs, as)@ with @m@ is full-rank submatrix,
                     --   @bs@ are independent and @as@ are dependent.
splitIndependentDirs :: Direction -> Matrix a -> (Matrix a, [Int], [Int])
splitIndependentDirs Direction
dir Matrix a
mat =
  case Direction -> Matrix a -> [Int]
forall r. Direction -> Matrix r -> [Int]
nonZeroDirs Direction
dir Matrix a
mat of
    []  -> (Matrix a
forall m. Monoidal m => m
zero, [], [])
    [Int
a] -> (Direction -> Vector a -> Matrix a
forall a. DecidableZero a => Direction -> Vector a -> Matrix a
dirVector Direction
dir (Vector a -> Matrix a) -> Vector a -> Matrix a
forall a b. (a -> b) -> a -> b
$ Direction -> Int -> Matrix a -> Vector a
forall a. Monoidal a => Direction -> Int -> Matrix a -> Vector a
getDir Direction
dir Int
a Matrix a
mat, [Int
a], [])
    (Int
x:[Int]
xs)  -> Int
-> [Int] -> Matrix a -> [Int] -> [Int] -> (Matrix a, [Int], [Int])
go Int
1 [Int]
xs (Direction -> Vector a -> Matrix a
forall a. DecidableZero a => Direction -> Vector a -> Matrix a
dirVector Direction
dir (Vector a -> Matrix a) -> Vector a -> Matrix a
forall a b. (a -> b) -> a -> b
$ Direction -> Int -> Matrix a -> Vector a
forall a. Monoidal a => Direction -> Int -> Matrix a -> Vector a
getDir Direction
dir Int
x Matrix a
mat) [Int
x] []
  where
    n :: Int
n = Int -> Int -> Int
forall a. Ord a => a -> a -> a
min (Matrix a -> Int
forall a. Matrix a -> Int
nrows Matrix a
mat) (Matrix a -> Int
forall a. Matrix a -> Int
ncols Matrix a
mat)
    go :: Int
-> [Int] -> Matrix a -> [Int] -> [Int] -> (Matrix a, [Int], [Int])
go Int
_ []     Matrix a
nat [Int]
ok [Int]
bad = (Matrix a
nat, [Int]
ok, [Int]
bad)
    go Int
i (Int
k:[Int]
ks) Matrix a
nat [Int]
ok [Int]
bad
      | Int
i Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
n = (Matrix a
nat, [Int]
ok, [Int]
bad)
      | Bool
otherwise =
        let nat' :: Matrix a
nat' = Direction -> Matrix a -> Vector a -> Matrix a
forall b.
DecidableZero b =>
Direction -> Matrix b -> Vector b -> Matrix b
catDir Direction
dir Matrix a
nat (Vector a -> Matrix a) -> Vector a -> Matrix a
forall a b. (a -> b) -> a -> b
$ Direction -> Int -> Matrix a -> Vector a
forall a. Monoidal a => Direction -> Int -> Matrix a -> Vector a
getDir Direction
dir Int
k Matrix a
mat
        in if ({-# SCC "rankLM" #-} Matrix a -> Int
forall r. (DecidableZero r, Division r, Group r) => Matrix r -> Int
rankLM Matrix a
nat') Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
i
           then Int
-> [Int] -> Matrix a -> [Int] -> [Int] -> (Matrix a, [Int], [Int])
go Int
i     [Int]
ks Matrix a
nat  [Int]
ok     (Int
kInt -> [Int] -> [Int]
forall a. a -> [a] -> [a]
:[Int]
bad)
           else Int
-> [Int] -> Matrix a -> [Int] -> [Int] -> (Matrix a, [Int], [Int])
go (Int
iInt -> Int -> Int
forall r. Additive r => r -> r -> r
+Int
1) [Int]
ks Matrix a
nat' (Int
kInt -> [Int] -> [Int]
forall a. a -> [a] -> [a]
:[Int]
ok) [Int]
bad

intDet :: Matrix Integer -> Integer
intDet :: Matrix Integer -> Integer
intDet Matrix Integer
mat =
  let b :: Natural
b = Vector Natural -> Natural
forall a. Ord a => Vector a -> a
V.maximum (Vector Natural -> Natural) -> Vector Natural -> Natural
forall a b. (a -> b) -> a -> b
$ (((Int, Int), Integer) -> Natural)
-> Vector ((Int, Int), Integer) -> Vector Natural
forall a b. (a -> b) -> Vector a -> Vector b
V.map (Integer -> Natural
forall r. Num r => Integer -> r
P.fromInteger (Integer -> Natural)
-> (((Int, Int), Integer) -> Integer)
-> ((Int, Int), Integer)
-> Natural
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Integer -> Integer
forall a. Num a => a -> a
abs (Integer -> Integer)
-> (((Int, Int), Integer) -> Integer)
-> ((Int, Int), Integer)
-> Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((Int, Int), Integer) -> Integer
forall a b. (a, b) -> b
snd) (Vector ((Int, Int), Integer) -> Vector Natural)
-> Vector ((Int, Int), Integer) -> Vector Natural
forall a b. (a -> b) -> a -> b
$ Matrix Integer -> Vector ((Int, Int), Integer)
forall a. Matrix a -> Vector ((Int, Int), a)
nonZeroEntries Matrix Integer
mat
      n :: Natural
n = Int -> Natural
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> Natural) -> Int -> Natural
forall a b. (a -> b) -> a -> b
$ Matrix Integer -> Int
forall a. Matrix a -> Int
ncols Matrix Integer
mat
      c :: Natural
c = Natural
nNatural -> Natural -> Natural
forall r. Unital r => r -> Natural -> r
^(Natural
n Natural -> Natural -> Natural
forall a. Integral a => a -> a -> a
`P.div` Natural
2) Natural -> Natural -> Natural
forall r. Multiplicative r => r -> r -> r
* Natural
bNatural -> Natural -> Natural
forall r. Unital r => r -> Natural -> r
^Natural
n
      r :: Int
r = Int -> Int
ceilingLogBase2 (Int
2Int -> Int -> Int
forall r. Multiplicative r => r -> r -> r
*Natural -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral Natural
c Int -> Int -> Int
forall r. Additive r => r -> r -> r
+ Int
1)
      ps :: [Integer]
ps = Int -> [Integer] -> [Integer]
forall a. Int -> [a] -> [a]
take Int
r [Integer]
forall int. Integral int => [int]
primes
      m :: Integer
m  = [Integer] -> Integer
forall (f :: * -> *) r. (Foldable f, Unital r) => f r -> r
product [Integer]
ps
      d :: Integer
d  = [(Integer, Integer)] -> Integer
forall r. Euclidean r => [(r, r)] -> r
chineseRemainder [ (Integer
p,
                               Integer
-> (forall (p :: Nat). KnownNat p => Proxy (F p) -> Integer)
-> Integer
forall a.
Integer -> (forall (p :: Nat). KnownNat p => Proxy (F p) -> a) -> a
reifyPrimeField Integer
p ((forall (p :: Nat). KnownNat p => Proxy (F p) -> Integer)
 -> Integer)
-> (forall (p :: Nat). KnownNat p => Proxy (F p) -> Integer)
-> Integer
forall a b. (a -> b) -> a -> b
$ \Proxy (F p)
pxy ->
                               Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
shiftHalf Integer
p (Integer -> Integer) -> Integer -> Integer
forall a b. (a -> b) -> a -> b
$ F p -> Integer
forall k (p :: k). F p -> Integer
naturalRepr (F p -> Integer) -> F p -> Integer
forall a b. (a -> b) -> a -> b
$ Getting (F p) (Matrix (F p), Matrix (F p), F p) (F p)
-> (Matrix (F p), Matrix (F p), F p) -> F p
forall s (m :: * -> *) a. MonadReader s m => Getting a s a -> m a
view Getting (F p) (Matrix (F p), Matrix (F p), F p) (F p)
forall s t a b. Field3 s t a b => Lens s t a b
_3 ((Matrix (F p), Matrix (F p), F p) -> F p)
-> (Matrix (F p), Matrix (F p), F p) -> F p
forall a b. (a -> b) -> a -> b
$
                               Matrix (F p) -> (Matrix (F p), Matrix (F p), F p)
forall a.
(DecidableZero a, Division a, Group a) =>
Matrix a -> (Matrix a, Matrix a, a)
structuredGauss' ((Integer -> F p) -> Matrix Integer -> Matrix (F p)
forall a a1. DecidableZero a => (a1 -> a) -> Matrix a1 -> Matrix a
cmap (Proxy (F p) -> Integer -> F p
forall k (proxy :: * -> *) (p :: k).
Reifies p Integer =>
proxy (F p) -> Integer -> F p
modNat' Proxy (F p)
pxy) Matrix Integer
mat))
                            | Integer
p <- [Integer]
ps]
      off :: Integer
off = Integer
d Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`div` Integer
m
  in if Integer
d Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
0
     then Integer
0
     else (Integer -> Integer -> Ordering) -> [Integer] -> Integer
forall (t :: * -> *) a.
Foldable t =>
(a -> a -> Ordering) -> t a -> a
minimumBy ((Integer -> Integer) -> Integer -> Integer -> Ordering
forall a b. Ord a => (b -> a) -> b -> b -> Ordering
comparing Integer -> Integer
forall a. Num a => a -> a
abs) [Integer
d Integer -> Integer -> Integer
forall r. Group r => r -> r -> r
- Integer
m Integer -> Integer -> Integer
forall r. Multiplicative r => r -> r -> r
* Integer
off, Integer
d Integer -> Integer -> Integer
forall r. Group r => r -> r -> r
- Integer
m Integer -> Integer -> Integer
forall r. Multiplicative r => r -> r -> r
* (Integer
off Integer -> Integer -> Integer
forall r. Additive r => r -> r -> r
+ Integer
1)]

shiftHalf :: P.Integral a => a -> a -> a
shiftHalf :: a -> a -> a
shiftHalf a
p a
n =
  let s :: a
s = a
p a -> a -> a
forall a. Integral a => a -> a -> a
`P.div` a
2
  in (a
n a -> a -> a
forall a. Num a => a -> a -> a
P.+ a
s) a -> a -> a
forall a. Integral a => a -> a -> a
`P.mod` a
p a -> a -> a
forall a. Num a => a -> a -> a
P.- a
s

triangulateModular :: Matrix (Fraction Integer)
                   -> (Matrix (Fraction Integer),
                       Matrix (Fraction Integer))
triangulateModular :: Matrix (Fraction Integer)
-> (Matrix (Fraction Integer), Matrix (Fraction Integer))
triangulateModular Matrix (Fraction Integer)
mat0 =
  let ps :: [Integer]
ps = (Integer -> Bool) -> [Integer] -> [Integer]
forall a. (a -> Bool) -> [a] -> [a]
filter ((Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
/= Integer
0) (Integer -> Bool) -> (Integer -> Integer) -> Integer -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
P.mod Integer
l)) [Integer]
forall int. Integral int => [int]
primes
  in [Integer] -> (Matrix (Fraction Integer), Matrix (Fraction Integer))
go [Integer]
ps
  where
    ds :: Integer
ds = (((Int, Int), Fraction Integer) -> Integer -> Integer)
-> Integer -> Vector ((Int, Int), Fraction Integer) -> Integer
forall a b. (a -> b -> b) -> b -> Vector a -> b
V.foldr (Integer -> Integer -> Integer
forall d. Euclidean d => d -> d -> d
lcm' (Integer -> Integer -> Integer)
-> (((Int, Int), Fraction Integer) -> Integer)
-> ((Int, Int), Fraction Integer)
-> Integer
-> Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Integer -> Integer
forall a. Num a => a -> a
abs (Integer -> Integer)
-> (((Int, Int), Fraction Integer) -> Integer)
-> ((Int, Int), Fraction Integer)
-> Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Fraction Integer -> Integer
forall t. Fraction t -> t
denominator(Fraction Integer -> Integer)
-> (((Int, Int), Fraction Integer) -> Fraction Integer)
-> ((Int, Int), Fraction Integer)
-> Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
.((Int, Int), Fraction Integer) -> Fraction Integer
forall a b. (a, b) -> b
snd) Integer
1 (Vector ((Int, Int), Fraction Integer) -> Integer)
-> Vector ((Int, Int), Fraction Integer) -> Integer
forall a b. (a -> b) -> a -> b
$ Matrix (Fraction Integer) -> Vector ((Int, Int), Fraction Integer)
forall a. Matrix a -> Vector ((Int, Int), a)
nonZeroEntries Matrix (Fraction Integer)
mat0
    mN :: Integer
mN = (((Int, Int), Fraction Integer) -> Integer -> Integer)
-> Integer -> Vector ((Int, Int), Fraction Integer) -> Integer
forall a b. (a -> b -> b) -> b -> Vector a -> b
V.foldr (Integer -> Integer -> Integer
forall d. Euclidean d => d -> d -> d
lcm' (Integer -> Integer -> Integer)
-> (((Int, Int), Fraction Integer) -> Integer)
-> ((Int, Int), Fraction Integer)
-> Integer
-> Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Integer -> Integer
forall a. Num a => a -> a
abs (Integer -> Integer)
-> (((Int, Int), Fraction Integer) -> Integer)
-> ((Int, Int), Fraction Integer)
-> Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Fraction Integer -> Integer
forall t. Fraction t -> t
numerator (Fraction Integer -> Integer)
-> (((Int, Int), Fraction Integer) -> Fraction Integer)
-> ((Int, Int), Fraction Integer)
-> Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((Int, Int), Fraction Integer) -> Fraction Integer
forall a b. (a, b) -> b
snd) Integer
1 (Vector ((Int, Int), Fraction Integer) -> Integer)
-> Vector ((Int, Int), Fraction Integer) -> Integer
forall a b. (a -> b) -> a -> b
$ Matrix (Fraction Integer) -> Vector ((Int, Int), Fraction Integer)
forall a. Matrix a -> Vector ((Int, Int), a)
nonZeroEntries Matrix (Fraction Integer)
mat0
    l :: Integer
l  = Integer -> Integer -> Integer
forall d. Euclidean d => d -> d -> d
lcm' Integer
ds Integer
mN
    go :: [Integer] -> (Matrix (Fraction Integer), Matrix (Fraction Integer))
go (Integer
p:[Integer]
ps) =
      let (IntSet
indepRows, IntSet
_, IntSet
indepCols, IntSet
depCols) = Integer
-> (forall (p :: Nat).
    KnownNat p =>
    Proxy (F p) -> (IntSet, IntSet, IntSet, IntSet))
-> (IntSet, IntSet, IntSet, IntSet)
forall a.
Integer -> (forall (p :: Nat). KnownNat p => Proxy (F p) -> a) -> a
reifyPrimeField Integer
p ((forall (p :: Nat).
  KnownNat p =>
  Proxy (F p) -> (IntSet, IntSet, IntSet, IntSet))
 -> (IntSet, IntSet, IntSet, IntSet))
-> (forall (p :: Nat).
    KnownNat p =>
    Proxy (F p) -> (IntSet, IntSet, IntSet, IntSet))
-> (IntSet, IntSet, IntSet, IntSet)
forall a b. (a -> b) -> a -> b
$ \Proxy (F p)
pxy ->
            let mat :: Matrix (F p)
mat = (Fraction Integer -> F p)
-> Matrix (Fraction Integer) -> Matrix (F p)
forall a a1. DecidableZero a => (a1 -> a) -> Matrix a1 -> Matrix a
cmap (Proxy (F p) -> Fraction Integer -> F p
forall k. FiniteField k => Proxy k -> Fraction Integer -> k
modRat Proxy (F p)
pxy) Matrix (Fraction Integer)
mat0
                (Matrix (F p)
koho, [Int] -> IntSet
IS.fromList -> IntSet
irs, [Int] -> IntSet
IS.fromList -> IntSet
drs) =
                  {-# SCC "splitRow" #-} Direction -> Matrix (F p) -> (Matrix (F p), [Int], [Int])
forall a.
(DecidableZero a, Field a) =>
Direction -> Matrix a -> (Matrix a, [Int], [Int])
splitIndependentDirs Direction
Row Matrix (F p)
mat
                (Matrix (F p)
_, [Int] -> IntSet
IS.fromList -> IntSet
ics, [Int] -> IntSet
IS.fromList -> IntSet
dcs) =
                  {-# SCC "splitCol" #-} Direction -> Matrix (F p) -> (Matrix (F p), [Int], [Int])
forall a.
(DecidableZero a, Field a) =>
Direction -> Matrix a -> (Matrix a, [Int], [Int])
splitIndependentDirs Direction
Column Matrix (F p)
koho
            in (IntSet
irs,
                IntSet
drs IntSet -> IntSet -> IntSet
`IS.union` ([Int] -> IntSet
IS.fromList (Matrix (Fraction Integer) -> [Int]
forall r. Matrix r -> [Int]
nonZeroRows Matrix (Fraction Integer)
mat0) IntSet -> IntSet -> IntSet
IS.\\ IntSet
irs),
                IntSet
ics,
                IntSet
dcs IntSet -> IntSet -> IntSet
`IS.union` ([Int] -> IntSet
IS.fromList (Matrix (Fraction Integer) -> [Int]
forall r. Matrix r -> [Int]
nonZeroCols (Matrix (Fraction Integer) -> [Int])
-> Matrix (Fraction Integer) -> [Int]
forall a b. (a -> b) -> a -> b
$ IntMap Int -> IntMap Int -> Matrix (Fraction Integer)
extract IntMap Int
irdic IntMap Int
colIdentDic) IntSet -> IntSet -> IntSet
IS.\\ IntSet
ics))
          colIdentDic :: IntMap Int
colIdentDic = [(Int, Int)] -> IntMap Int
forall a. [(Int, a)] -> IntMap a
IM.fromList ([(Int, Int)] -> IntMap Int) -> [(Int, Int)] -> IntMap Int
forall a b. (a -> b) -> a -> b
$ [Int] -> [Int] -> [(Int, Int)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0..Matrix (Fraction Integer) -> Int
forall a. Matrix a -> Int
ncols Matrix (Fraction Integer)
mat0 Int -> Int -> Int
forall r. Group r => r -> r -> r
- Int
1] [Int
0..]
          irdic :: IntMap Int
irdic = [(Int, Int)] -> IntMap Int
forall a. [(Int, a)] -> IntMap a
IM.fromList ([(Int, Int)] -> IntMap Int) -> [(Int, Int)] -> IntMap Int
forall a b. (a -> b) -> a -> b
$ [Int] -> [Int] -> [(Int, Int)]
forall a b. [a] -> [b] -> [(a, b)]
zip (IntSet -> [Int]
IS.toAscList IntSet
indepRows) [Int
0..]
          icdic :: IntMap Int
icdic = [(Int, Int)] -> IntMap Int
forall a. [(Int, a)] -> IntMap a
IM.fromList ([(Int, Int)] -> IntMap Int) -> [(Int, Int)] -> IntMap Int
forall a b. (a -> b) -> a -> b
$ [Int] -> [Int] -> [(Int, Int)]
forall a b. [a] -> [b] -> [(a, b)]
zip (IntSet -> [Int]
IS.toAscList IntSet
indepCols) [Int
0..]
          dcdic :: IntMap Int
dcdic = [(Int, Int)] -> IntMap Int
forall a. [(Int, a)] -> IntMap a
IM.fromList ([(Int, Int)] -> IntMap Int) -> [(Int, Int)] -> IntMap Int
forall a b. (a -> b) -> a -> b
$ [Int] -> [Int] -> [(Int, Int)]
forall a b. [a] -> [b] -> [(a, b)]
zip (IntSet -> [Int]
IS.toDescList IntSet
depCols) [Int
0..]
          newIdx :: IntMap a -> IntMap a -> (Int, Int) -> Maybe (a, a)
newIdx IntMap a
rd IntMap a
cd (Int
i, Int
j) = (,) (a -> a -> (a, a)) -> Maybe a -> Maybe (a -> (a, a))
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Int -> IntMap a -> Maybe a
forall a. Int -> IntMap a -> Maybe a
IM.lookup Int
i IntMap a
rd Maybe (a -> (a, a)) -> Maybe a -> Maybe (a, a)
forall (f :: * -> *) a b. Applicative f => f (a -> b) -> f a -> f b
<*> Int -> IntMap a -> Maybe a
forall a. Int -> IntMap a -> Maybe a
IM.lookup Int
j IntMap a
cd
          extract :: IntMap Int -> IntMap Int -> Matrix (Fraction Integer)
extract IntMap Int
rd IntMap Int
cd = [((Int, Int), Fraction Integer)] -> Matrix (Fraction Integer)
forall a. DecidableZero a => [((Int, Int), a)] -> Matrix a
fromList ((((Int, Int), Fraction Integer)
 -> Maybe ((Int, Int), Fraction Integer))
-> [((Int, Int), Fraction Integer)]
-> [((Int, Int), Fraction Integer)]
forall a b. (a -> Maybe b) -> [a] -> [b]
mapMaybe (\((Int, Int)
ind, Fraction Integer
c) -> (,Fraction Integer
c) ((Int, Int) -> ((Int, Int), Fraction Integer))
-> Maybe (Int, Int) -> Maybe ((Int, Int), Fraction Integer)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> IntMap Int -> IntMap Int -> (Int, Int) -> Maybe (Int, Int)
forall a a. IntMap a -> IntMap a -> (Int, Int) -> Maybe (a, a)
newIdx IntMap Int
rd IntMap Int
cd (Int, Int)
ind) ([((Int, Int), Fraction Integer)]
 -> [((Int, Int), Fraction Integer)])
-> [((Int, Int), Fraction Integer)]
-> [((Int, Int), Fraction Integer)]
forall a b. (a -> b) -> a -> b
$
                          Vector ((Int, Int), Fraction Integer)
-> [((Int, Int), Fraction Integer)]
forall a. Vector a -> [a]
V.toList (Vector ((Int, Int), Fraction Integer)
 -> [((Int, Int), Fraction Integer)])
-> Vector ((Int, Int), Fraction Integer)
-> [((Int, Int), Fraction Integer)]
forall a b. (a -> b) -> a -> b
$ Matrix (Fraction Integer) -> Vector ((Int, Int), Fraction Integer)
forall a. Matrix a -> Vector ((Int, Int), a)
nonZeroEntries Matrix (Fraction Integer)
mat0)
                          Matrix (Fraction Integer)
-> (Matrix (Fraction Integer) -> Matrix (Fraction Integer))
-> Matrix (Fraction Integer)
forall a b. a -> (a -> b) -> b
& (Int -> Identity Int)
-> Matrix (Fraction Integer)
-> Identity (Matrix (Fraction Integer))
forall a. Lens' (Matrix a) Int
height ((Int -> Identity Int)
 -> Matrix (Fraction Integer)
 -> Identity (Matrix (Fraction Integer)))
-> Int -> Matrix (Fraction Integer) -> Matrix (Fraction Integer)
forall s t a b. ASetter s t a b -> b -> s -> t
.~ IntMap Int -> Int
forall a. IntMap a -> Int
IM.size IntMap Int
rd
                          Matrix (Fraction Integer)
-> (Matrix (Fraction Integer) -> Matrix (Fraction Integer))
-> Matrix (Fraction Integer)
forall a b. a -> (a -> b) -> b
& (Int -> Identity Int)
-> Matrix (Fraction Integer)
-> Identity (Matrix (Fraction Integer))
forall a. Lens' (Matrix a) Int
width  ((Int -> Identity Int)
 -> Matrix (Fraction Integer)
 -> Identity (Matrix (Fraction Integer)))
-> Int -> Matrix (Fraction Integer) -> Matrix (Fraction Integer)
forall s t a b. ASetter s t a b -> b -> s -> t
.~ IntMap Int -> Int
forall a. IntMap a -> Int
IM.size IntMap Int
cd
          spec :: Matrix (Fraction Integer)
spec = IntMap Int -> IntMap Int -> Matrix (Fraction Integer)
extract IntMap Int
irdic IntMap Int
icdic
          qs :: [Integer]
qs = (Integer -> Bool) -> [Integer] -> [Integer]
forall a. (a -> Bool) -> [a] -> [a]
filter (\Integer
r -> Integer
ds Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`mod` Integer
r Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
/= Integer
0 Bool -> Bool -> Bool
&& ({-# SCC "checkDet" #-} Integer
-> (forall (p :: Nat). KnownNat p => Proxy (F p) -> Bool) -> Bool
forall a.
Integer -> (forall (p :: Nat). KnownNat p => Proxy (F p) -> a) -> a
reifyPrimeField Integer
r ((forall (p :: Nat). KnownNat p => Proxy (F p) -> Bool) -> Bool)
-> (forall (p :: Nat). KnownNat p => Proxy (F p) -> Bool) -> Bool
forall a b. (a -> b) -> a -> b
$ \Proxy (F p)
pxy ->
                          Bool -> Bool
not (Bool -> Bool) -> Bool -> Bool
forall a b. (a -> b) -> a -> b
$ F p -> Bool
forall r. DecidableZero r => r -> Bool
isZero (F p -> Bool) -> F p -> Bool
forall a b. (a -> b) -> a -> b
$ Getting (F p) (Matrix (F p), Matrix (F p), F p) (F p)
-> (Matrix (F p), Matrix (F p), F p) -> F p
forall s (m :: * -> *) a. MonadReader s m => Getting a s a -> m a
view Getting (F p) (Matrix (F p), Matrix (F p), F p) (F p)
forall s t a b. Field3 s t a b => Lens s t a b
_3 ((Matrix (F p), Matrix (F p), F p) -> F p)
-> (Matrix (F p), Matrix (F p), F p) -> F p
forall a b. (a -> b) -> a -> b
$
                          Matrix (F p) -> (Matrix (F p), Matrix (F p), F p)
forall a.
(DecidableZero a, Division a, Group a) =>
Matrix a -> (Matrix a, Matrix a, a)
structuredGauss' (Matrix (F p) -> (Matrix (F p), Matrix (F p), F p))
-> Matrix (F p) -> (Matrix (F p), Matrix (F p), F p)
forall a b. (a -> b) -> a -> b
$ (Fraction Integer -> F p)
-> Matrix (Fraction Integer) -> Matrix (F p)
forall a a1. DecidableZero a => (a1 -> a) -> Matrix a1 -> Matrix a
cmap (Proxy (F p) -> Fraction Integer -> F p
forall k. FiniteField k => Proxy k -> Fraction Integer -> k
modRat Proxy (F p)
pxy) (Matrix (Fraction Integer) -> Matrix (F p))
-> Matrix (Fraction Integer) -> Matrix (F p)
forall a b. (a -> b) -> a -> b
$ Matrix (Fraction Integer)
spec)) [Integer]
forall int. Integral int => [int]
primes
          anss :: [Vector (Fraction Integer)]
anss = Strategy (Vector (Fraction Integer))
-> (Vector (Fraction Integer) -> Vector (Fraction Integer))
-> [Vector (Fraction Integer)]
-> [Vector (Fraction Integer)]
forall b a. Strategy b -> (a -> b) -> [a] -> [b]
parMap Strategy (Vector (Fraction Integer))
forall a. NFData a => Strategy a
rdeepseq (\Vector (Fraction Integer)
xs -> Maybe (Vector (Fraction Integer)) -> Vector (Fraction Integer)
forall a. HasCallStack => Maybe a -> a
fromJust (Maybe (Vector (Fraction Integer)) -> Vector (Fraction Integer))
-> Maybe (Vector (Fraction Integer)) -> Vector (Fraction Integer)
forall a b. (a -> b) -> a -> b
$ (Unwrapped (First (Vector (Fraction Integer)))
 -> First (Vector (Fraction Integer)))
-> ((Unwrapped (First (Vector (Fraction Integer)))
     -> First (Vector (Fraction Integer)))
    -> [Maybe (Vector (Fraction Integer))]
    -> First (Vector (Fraction Integer)))
-> [Maybe (Vector (Fraction Integer))]
-> Unwrapped (First (Vector (Fraction Integer)))
forall (f :: * -> *) s t.
(Functor f, Rewrapping s t) =>
(Unwrapped s -> s)
-> ((Unwrapped t -> t) -> f s) -> f (Unwrapped s)
ala Unwrapped (First (Vector (Fraction Integer)))
-> First (Vector (Fraction Integer))
forall a. Maybe a -> First a
First (Unwrapped (First (Vector (Fraction Integer)))
 -> First (Vector (Fraction Integer)))
-> [Maybe (Vector (Fraction Integer))]
-> First (Vector (Fraction Integer))
forall (t :: * -> *) m a.
(Foldable t, Monoid m) =>
(a -> m) -> t a -> m
foldMap ([Maybe (Vector (Fraction Integer))]
 -> Maybe (Vector (Fraction Integer)))
-> [Maybe (Vector (Fraction Integer))]
-> Maybe (Vector (Fraction Integer))
forall a b. (a -> b) -> a -> b
$
                                         (Integer -> Maybe (Vector (Fraction Integer)))
-> [Integer] -> [Maybe (Vector (Fraction Integer))]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
map (\Integer
q -> Int
-> Integer
-> Matrix (Fraction Integer)
-> Vector (Fraction Integer)
-> Maybe (Vector (Fraction Integer))
solveHensel Int
10 Integer
q Matrix (Fraction Integer)
spec Vector (Fraction Integer)
xs) [Integer]
qs) ([Vector (Fraction Integer)] -> [Vector (Fraction Integer)])
-> [Vector (Fraction Integer)] -> [Vector (Fraction Integer)]
forall a b. (a -> b) -> a -> b
$
                 Matrix (Fraction Integer) -> [Vector (Fraction Integer)]
forall a. Monoidal a => Matrix a -> [Vector a]
toCols (Matrix (Fraction Integer) -> [Vector (Fraction Integer)])
-> Matrix (Fraction Integer) -> [Vector (Fraction Integer)]
forall a b. (a -> b) -> a -> b
$ IntMap Int -> IntMap Int -> Matrix (Fraction Integer)
extract IntMap Int
irdic IntMap Int
dcdic
          permMat :: Matrix (Fraction Integer)
permMat = [Vector (Fraction Integer)]
-> Int
-> [(Int, Vector (Fraction Integer))]
-> [(Int, Vector (Fraction Integer))]
-> Matrix (Fraction Integer)
forall r a.
(Ord r, Num r, DecidableZero a, Group r) =>
[Vector a] -> r -> [(r, Vector a)] -> [(r, Vector a)] -> Matrix a
build [] (Matrix (Fraction Integer) -> Int
forall a. Matrix a -> Int
ncols Matrix (Fraction Integer)
mat0 Int -> Int -> Int
forall r. Group r => r -> r -> r
- Int
1)
                    ([Int]
-> [Vector (Fraction Integer)]
-> [(Int, Vector (Fraction Integer))]
forall a b. [a] -> [b] -> [(a, b)]
zip (IntSet -> [Int]
IS.toDescList IntSet
indepCols) ([Vector (Fraction Integer)] -> [(Int, Vector (Fraction Integer))])
-> [Vector (Fraction Integer)]
-> [(Int, Vector (Fraction Integer))]
forall a b. (a -> b) -> a -> b
$ [Vector (Fraction Integer)] -> [Vector (Fraction Integer)]
forall a. [a] -> [a]
reverse ([Vector (Fraction Integer)] -> [Vector (Fraction Integer)])
-> [Vector (Fraction Integer)] -> [Vector (Fraction Integer)]
forall a b. (a -> b) -> a -> b
$ Matrix (Fraction Integer) -> [Vector (Fraction Integer)]
forall a. Monoidal a => Matrix a -> [Vector a]
toCols (Matrix (Fraction Integer) -> [Vector (Fraction Integer)])
-> Matrix (Fraction Integer) -> [Vector (Fraction Integer)]
forall a b. (a -> b) -> a -> b
$ Int -> Matrix (Fraction Integer)
forall a. Unital a => Int -> Matrix a
identity (Int -> Matrix (Fraction Integer))
-> Int -> Matrix (Fraction Integer)
forall a b. (a -> b) -> a -> b
$ Matrix (Fraction Integer) -> Int
forall a. Matrix a -> Int
nrows Matrix (Fraction Integer)
spec) ([(Int, Vector (Fraction Integer))] -> Matrix (Fraction Integer))
-> [(Int, Vector (Fraction Integer))] -> Matrix (Fraction Integer)
forall a b. (a -> b) -> a -> b
$
                    [Int]
-> [Vector (Fraction Integer)]
-> [(Int, Vector (Fraction Integer))]
forall a b. [a] -> [b] -> [(a, b)]
zip (IntSet -> [Int]
IS.toDescList IntSet
depCols) [Vector (Fraction Integer)]
anss
          origDeled :: Matrix (Fraction Integer)
origDeled = IntMap Int -> IntMap Int -> Matrix (Fraction Integer)
extract IntMap Int
irdic IntMap Int
colIdentDic Matrix (Fraction Integer)
-> (Matrix (Fraction Integer) -> Matrix (Fraction Integer))
-> Matrix (Fraction Integer)
forall a b. a -> (a -> b) -> b
& (Int -> Identity Int)
-> Matrix (Fraction Integer)
-> Identity (Matrix (Fraction Integer))
forall a. Lens' (Matrix a) Int
width ((Int -> Identity Int)
 -> Matrix (Fraction Integer)
 -> Identity (Matrix (Fraction Integer)))
-> Int -> Matrix (Fraction Integer) -> Matrix (Fraction Integer)
forall s t a b. ASetter s t a b -> b -> s -> t
.~ Matrix (Fraction Integer)
mat0Matrix (Fraction Integer)
-> Getting Int (Matrix (Fraction Integer)) Int -> Int
forall s a. s -> Getting a s a -> a
^.Getting Int (Matrix (Fraction Integer)) Int
forall a. Lens' (Matrix a) Int
width
     in if (Matrix (Fraction Integer)
spec Matrix (Fraction Integer)
-> Matrix (Fraction Integer) -> Matrix (Fraction Integer)
forall r. Multiplicative r => r -> r -> r
* Matrix (Fraction Integer)
permMat) Matrix (Fraction Integer) -> Matrix (Fraction Integer) -> Bool
forall a. Eq a => a -> a -> Bool
== Matrix (Fraction Integer)
origDeled
        then (Matrix (Fraction Integer)
permMat, Matrix (Fraction Integer)
spec)
        else [Integer] -> (Matrix (Fraction Integer), Matrix (Fraction Integer))
go [Integer]
ps
    go [Integer]
_ = String -> (Matrix (Fraction Integer), Matrix (Fraction Integer))
forall a. HasCallStack => String -> a
error String
"Cannot happen!"
    build :: [Vector a] -> r -> [(r, Vector a)] -> [(r, Vector a)] -> Matrix a
build [Vector a]
ans r
i [(r, Vector a)]
mns [(r, Vector a)]
vecs
      | r
i r -> r -> Bool
forall a. Ord a => a -> a -> Bool
< r
0 = [Vector a] -> Matrix a
forall a. DecidableZero a => [Vector a] -> Matrix a
fromCols [Vector a]
ans
      | Bool
otherwise = {-# SCC "building" #-}
        case [(r, Vector a)]
vecs of
          ((r
k, Vector a
v) : [(r, Vector a)]
vs) | r
i r -> r -> Bool
forall a. Eq a => a -> a -> Bool
== r
k -> [Vector a] -> r -> [(r, Vector a)] -> [(r, Vector a)] -> Matrix a
build (Vector a
v Vector a -> [Vector a] -> [Vector a]
forall a. a -> [a] -> [a]
: [Vector a]
ans) (r
ir -> r -> r
forall r. Group r => r -> r -> r
-r
1) [(r, Vector a)]
mns [(r, Vector a)]
vs
          [(r, Vector a)]
_ ->
            case [(r, Vector a)]
mns of
              ((r
l,Vector a
m):[(r, Vector a)]
mn) | r
i r -> r -> Bool
forall a. Eq a => a -> a -> Bool
== r
l -> [Vector a] -> r -> [(r, Vector a)] -> [(r, Vector a)] -> Matrix a
build (Vector a
m Vector a -> [Vector a] -> [Vector a]
forall a. a -> [a] -> [a]
: [Vector a]
ans) (r
ir -> r -> r
forall r. Group r => r -> r -> r
-r
1) [(r, Vector a)]
mn [(r, Vector a)]
vecs
              [(r, Vector a)]
_ -> [Vector a] -> r -> [(r, Vector a)] -> [(r, Vector a)] -> Matrix a
build (Vector a
forall a. Vector a
V.empty Vector a -> [Vector a] -> [Vector a]
forall a. a -> [a] -> [a]
: [Vector a]
ans) (r
ir -> r -> r
forall r. Group r => r -> r -> r
-r
1) [(r, Vector a)]
mns [(r, Vector a)]
vecs

lcm' :: Euclidean r => r -> r -> r
lcm' :: r -> r -> r
lcm' r
n r
m = r
n r -> r -> r
forall r. Multiplicative r => r -> r -> r
* r
m r -> r -> r
forall d. Euclidean d => d -> d -> d
`quot` r -> r -> r
forall d. GCDDomain d => d -> d -> d
gcd r
n r
m

henselLift :: Integer           -- ^ prime number @p@
           -> Matrix Integer -- ^ original matrix @M@
           -> Matrix Integer -- ^ inverse matrix of @M@ mod @p@
           -> V.Vector Integer  -- ^ coefficient vector @v@
           -> [V.Vector Integer]  -- ^ vector @x@ with @Mx = b mod p@
henselLift :: Integer
-> Matrix Integer
-> Matrix Integer
-> Vector Integer
-> [Vector Integer]
henselLift Integer
p Matrix Integer
m Matrix Integer
q Vector Integer
b =
  ((Integer, Vector Integer, Vector Integer) -> Vector Integer)
-> [(Integer, Vector Integer, Vector Integer)] -> [Vector Integer]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
map (Getting
  (Vector Integer)
  (Integer, Vector Integer, Vector Integer)
  (Vector Integer)
-> (Integer, Vector Integer, Vector Integer) -> Vector Integer
forall s (m :: * -> *) a. MonadReader s m => Getting a s a -> m a
view Getting
  (Vector Integer)
  (Integer, Vector Integer, Vector Integer)
  (Vector Integer)
forall s t a b. Field2 s t a b => Lens s t a b
_2) ([(Integer, Vector Integer, Vector Integer)] -> [Vector Integer])
-> [(Integer, Vector Integer, Vector Integer)] -> [Vector Integer]
forall a b. (a -> b) -> a -> b
$ ((Integer, Vector Integer, Vector Integer)
 -> (Integer, Vector Integer, Vector Integer))
-> (Integer, Vector Integer, Vector Integer)
-> [(Integer, Vector Integer, Vector Integer)]
forall a. (a -> a) -> a -> [a]
iterate (Integer, Vector Integer, Vector Integer)
-> (Integer, Vector Integer, Vector Integer)
step (Integer
1, Int -> Integer -> Vector Integer
forall a. Int -> a -> Vector a
V.replicate (Vector Integer -> Int
forall a. Vector a -> Int
V.length Vector Integer
b) Integer
0, Vector Integer
b)
  where
    step :: (Integer, Vector Integer, Vector Integer)
-> (Integer, Vector Integer, Vector Integer)
step (Integer
s, Vector Integer
acc, Vector Integer
r)
      | Bool
otherwise =
        let u :: Vector Integer
u = Integer
-> (forall (p :: Nat). KnownNat p => Proxy (F p) -> Vector Integer)
-> Vector Integer
forall a.
Integer -> (forall (p :: Nat). KnownNat p => Proxy (F p) -> a) -> a
reifyPrimeField Integer
p ((forall (p :: Nat). KnownNat p => Proxy (F p) -> Vector Integer)
 -> Vector Integer)
-> (forall (p :: Nat). KnownNat p => Proxy (F p) -> Vector Integer)
-> Vector Integer
forall a b. (a -> b) -> a -> b
$ \Proxy (F p)
pxy ->
              (Integer -> Integer) -> Vector Integer -> Vector Integer
forall a b. (a -> b) -> Vector a -> Vector b
V.map (F p -> Integer
forall k (p :: k). F p -> Integer
naturalRepr (F p -> Integer) -> (Integer -> F p) -> Integer -> Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Proxy (F p) -> Integer -> F p
forall k (proxy :: * -> *) (p :: k).
Reifies p Integer =>
proxy (F p) -> Integer -> F p
modNat' Proxy (F p)
pxy) (Vector Integer -> Vector Integer)
-> Vector Integer -> Vector Integer
forall a b. (a -> b) -> a -> b
$ Matrix Integer
q Matrix Integer -> Vector Integer -> Vector Integer
forall a.
(Multiplicative a, Monoidal a) =>
Matrix a -> Vector a -> Vector a
`multWithVector` Vector Integer
r
            r' :: Vector Integer
r' = (Integer -> Integer) -> Vector Integer -> Vector Integer
forall a b. (a -> b) -> Vector a -> Vector b
V.map (Integer -> Integer -> Integer
forall d. Euclidean d => d -> d -> d
`quot` Integer
p) (Vector Integer -> Vector Integer)
-> Vector Integer -> Vector Integer
forall a b. (a -> b) -> a -> b
$ (Integer -> Integer -> Integer)
-> Vector Integer -> Vector Integer -> Vector Integer
forall a b c. (a -> b -> c) -> Vector a -> Vector b -> Vector c
V.zipWith (-) Vector Integer
r (Matrix Integer
m Matrix Integer -> Vector Integer -> Vector Integer
forall a.
(Multiplicative a, Monoidal a) =>
Matrix a -> Vector a -> Vector a
`multWithVector` Vector Integer
u)
        in (Integer
sInteger -> Integer -> Integer
forall r. Multiplicative r => r -> r -> r
*Integer
p, Vector Integer
acc Vector Integer -> Vector Integer -> Vector Integer
forall r. Additive r => r -> r -> r
+ (Integer -> Integer) -> Vector Integer -> Vector Integer
forall a b. (a -> b) -> Vector a -> Vector b
V.map (Integer
sInteger -> Integer -> Integer
forall r. Multiplicative r => r -> r -> r
*) Vector Integer
u, Vector Integer
r')


solveHensel :: Int -> Integer
            -> Matrix (Fraction Integer)
            -> Vector (Fraction Integer)
            -> Maybe (Vector (Fraction Integer))
solveHensel :: Int
-> Integer
-> Matrix (Fraction Integer)
-> Vector (Fraction Integer)
-> Maybe (Vector (Fraction Integer))
solveHensel Int
cyc Integer
p Matrix (Fraction Integer)
mat Vector (Fraction Integer)
b = {-# SCC "solveHensel" #-}
  let g0 :: Integer
g0 = (((Int, Int), Fraction Integer) -> Integer -> Integer)
-> Integer -> Vector ((Int, Int), Fraction Integer) -> Integer
forall a b. (a -> b -> b) -> b -> Vector a -> b
V.foldr (Integer -> Integer -> Integer
forall d. GCDDomain d => d -> d -> d
lcm (Integer -> Integer -> Integer)
-> (((Int, Int), Fraction Integer) -> Integer)
-> ((Int, Int), Fraction Integer)
-> Integer
-> Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Fraction Integer -> Integer
forall t. Fraction t -> t
denominator (Fraction Integer -> Integer)
-> (((Int, Int), Fraction Integer) -> Fraction Integer)
-> ((Int, Int), Fraction Integer)
-> Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Getting
  (Fraction Integer)
  ((Int, Int), Fraction Integer)
  (Fraction Integer)
-> ((Int, Int), Fraction Integer) -> Fraction Integer
forall s (m :: * -> *) a. MonadReader s m => Getting a s a -> m a
view Getting
  (Fraction Integer)
  ((Int, Int), Fraction Integer)
  (Fraction Integer)
forall s t a b. Field2 s t a b => Lens s t a b
_2) Integer
forall r. Unital r => r
one (Vector ((Int, Int), Fraction Integer) -> Integer)
-> Vector ((Int, Int), Fraction Integer) -> Integer
forall a b. (a -> b) -> a -> b
$ Matrix (Fraction Integer) -> Vector ((Int, Int), Fraction Integer)
forall a. Matrix a -> Vector ((Int, Int), a)
nonZeroEntries Matrix (Fraction Integer)
mat
      g1 :: Integer
g1 = (Fraction Integer -> Integer -> Integer)
-> Integer -> Vector (Fraction Integer) -> Integer
forall a b. (a -> b -> b) -> b -> Vector a -> b
V.foldr (Integer -> Integer -> Integer
forall d. GCDDomain d => d -> d -> d
lcm (Integer -> Integer -> Integer)
-> (Fraction Integer -> Integer)
-> Fraction Integer
-> Integer
-> Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Fraction Integer -> Integer
forall t. Fraction t -> t
denominator) Integer
forall r. Unital r => r
one Vector (Fraction Integer)
b
      g :: Fraction Integer
g  = Integer -> Integer -> Integer
forall d. GCDDomain d => d -> d -> d
lcm Integer
g0 Integer
g1 Integer -> Integer -> Fraction Integer
forall d. GCDDomain d => d -> d -> Fraction d
% Integer
1
      mat' :: Matrix Integer
mat' = (Fraction Integer -> Integer)
-> Matrix (Fraction Integer) -> Matrix Integer
forall a a1. DecidableZero a => (a1 -> a) -> Matrix a1 -> Matrix a
cmap (Fraction Integer -> Integer
forall t. Fraction t -> t
numerator (Fraction Integer -> Integer)
-> (Fraction Integer -> Fraction Integer)
-> Fraction Integer
-> Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Fraction Integer -> Fraction Integer -> Fraction Integer
forall r. Multiplicative r => r -> r -> r
*Fraction Integer
g)) Matrix (Fraction Integer)
mat
      b' :: Vector Integer
b'   = (Fraction Integer -> Integer)
-> Vector (Fraction Integer) -> Vector Integer
forall a b. (a -> b) -> Vector a -> Vector b
V.map (Fraction Integer -> Integer
forall t. Fraction t -> t
numerator (Fraction Integer -> Integer)
-> (Fraction Integer -> Fraction Integer)
-> Fraction Integer
-> Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Fraction Integer -> Fraction Integer -> Fraction Integer
forall r. Multiplicative r => r -> r -> r
*Fraction Integer
g)) Vector (Fraction Integer)
b
      q :: Matrix Integer
q = Integer
-> (forall (p :: Nat). KnownNat p => Proxy (F p) -> Matrix Integer)
-> Matrix Integer
forall a.
Integer -> (forall (p :: Nat). KnownNat p => Proxy (F p) -> a) -> a
reifyPrimeField Integer
p ((forall (p :: Nat). KnownNat p => Proxy (F p) -> Matrix Integer)
 -> Matrix Integer)
-> (forall (p :: Nat). KnownNat p => Proxy (F p) -> Matrix Integer)
-> Matrix Integer
forall a b. (a -> b) -> a -> b
$ \Proxy (F p)
pxy ->
        (F p -> Integer) -> Matrix (F p) -> Matrix Integer
forall a a1. DecidableZero a => (a1 -> a) -> Matrix a1 -> Matrix a
cmap F p -> Integer
forall k (p :: k). F p -> Integer
naturalRepr (Matrix (F p) -> Matrix Integer) -> Matrix (F p) -> Matrix Integer
forall a b. (a -> b) -> a -> b
$ (Matrix (F p), Matrix (F p)) -> Matrix (F p)
forall a b. (a, b) -> b
snd ((Matrix (F p), Matrix (F p)) -> Matrix (F p))
-> (Matrix (F p), Matrix (F p)) -> Matrix (F p)
forall a b. (a -> b) -> a -> b
$ Matrix (F p) -> (Matrix (F p), Matrix (F p))
forall a.
(DecidableZero a, Division a, Group a) =>
Matrix a -> (Matrix a, Matrix a)
structuredGauss (Matrix (F p) -> (Matrix (F p), Matrix (F p)))
-> Matrix (F p) -> (Matrix (F p), Matrix (F p))
forall a b. (a -> b) -> a -> b
$ (Integer -> F p) -> Matrix Integer -> Matrix (F p)
forall a a1. DecidableZero a => (a1 -> a) -> Matrix a1 -> Matrix a
cmap (Proxy (F p) -> Integer -> F p
forall k (proxy :: * -> *) (p :: k).
Reifies p Integer =>
proxy (F p) -> Integer -> F p
modNat' Proxy (F p)
pxy) Matrix Integer
mat'
      hls :: [Vector Integer]
hls = Integer
-> Matrix Integer
-> Matrix Integer
-> Vector Integer
-> [Vector Integer]
henselLift Integer
p Matrix Integer
mat' Matrix Integer
q Vector Integer
b'
  in Maybe (Vector (Fraction Integer))
-> [(Integer, Vector Integer)] -> Maybe (Vector (Fraction Integer))
go Maybe (Vector (Fraction Integer))
forall a. Maybe a
Nothing ([(Integer, Vector Integer)] -> Maybe (Vector (Fraction Integer)))
-> [(Integer, Vector Integer)] -> Maybe (Vector (Fraction Integer))
forall a b. (a -> b) -> a -> b
$ Int -> [(Integer, Vector Integer)] -> [(Integer, Vector Integer)]
forall a. Int -> [a] -> [a]
drop Int
cyc ([(Integer, Vector Integer)] -> [(Integer, Vector Integer)])
-> [(Integer, Vector Integer)] -> [(Integer, Vector Integer)]
forall a b. (a -> b) -> a -> b
$ [Integer] -> [Vector Integer] -> [(Integer, Vector Integer)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Integer
pInteger -> Natural -> Integer
forall r. Unital r => r -> Natural -> r
^Natural
i | Natural
i <- [Natural
0..]] [Vector Integer]
hls
  where
    go :: Maybe (Vector (Fraction Integer))
-> [(Integer, Vector Integer)] -> Maybe (Vector (Fraction Integer))
go Maybe (Vector (Fraction Integer))
_ [] = Maybe (Vector (Fraction Integer))
forall a. Maybe a
Nothing
    go Maybe (Vector (Fraction Integer))
prev ((Integer
q,Vector Integer
x):[(Integer, Vector Integer)]
xs) =
      let mans :: Maybe (Vector (Fraction Integer))
mans = (Integer -> Maybe (Fraction Integer))
-> Vector Integer -> Maybe (Vector (Fraction Integer))
forall (m :: * -> *) a b.
Monad m =>
(a -> m b) -> Vector a -> m (Vector b)
V.mapM (Integer -> Integer -> Integer -> Maybe (Fraction Integer)
recoverRat (Double -> Integer
forall a b. (RealFrac a, Integral b) => a -> b
P.floor (Double -> Integer) -> Double -> Integer
forall a b. (a -> b) -> a -> b
$ Double -> Double
forall a. Floating a => a -> a
P.sqrt (Integer -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Integer
q Double -> Double -> Double
forall a. Fractional a => a -> a -> a
P./ Double
2 :: Double)) Integer
q) Vector Integer
x
      in case Maybe (Vector (Fraction Integer))
mans of
        Just Vector (Fraction Integer)
x' | Matrix (Fraction Integer)
mat Matrix (Fraction Integer)
-> Vector (Fraction Integer) -> Vector (Fraction Integer)
forall a.
(Multiplicative a, Monoidal a) =>
Matrix a -> Vector a -> Vector a
`multWithVector` Vector (Fraction Integer)
x' Vector (Fraction Integer) -> Vector (Fraction Integer) -> Bool
forall a. Eq a => a -> a -> Bool
== Vector (Fraction Integer)
b -> Vector (Fraction Integer) -> Maybe (Vector (Fraction Integer))
forall a. a -> Maybe a
Just Vector (Fraction Integer)
x'
                | Maybe (Vector (Fraction Integer))
mans Maybe (Vector (Fraction Integer))
-> Maybe (Vector (Fraction Integer)) -> Bool
forall a. Eq a => a -> a -> Bool
== Maybe (Vector (Fraction Integer))
prev -> Maybe (Vector (Fraction Integer))
forall a. Maybe a
Nothing
        Maybe (Vector (Fraction Integer))
_ -> Maybe (Vector (Fraction Integer))
-> [(Integer, Vector Integer)] -> Maybe (Vector (Fraction Integer))
go Maybe (Vector (Fraction Integer))
mans (Int -> [(Integer, Vector Integer)] -> [(Integer, Vector Integer)]
forall a. Int -> [a] -> [a]
drop Int
cyc [(Integer, Vector Integer)]
xs)