{-# LANGUAGE ConstraintKinds #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE ExplicitNamespaces #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE ParallelListComp #-}
{-# LANGUAGE PatternSynonyms #-}
{-# LANGUAGE PolyKinds #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeApplications #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE NoMonomorphismRestriction #-}
{-# OPTIONS_GHC -fno-warn-name-shadowing #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.Presburger #-}
module Algebra.Algorithms.ZeroDim
(
solveM,
solve',
solveViaCompanion,
solveLinear,
radical,
isRadical,
fglm,
fglmMap,
solveWith,
univPoly,
reduction,
matrixRep,
subspMatrix,
vectorRep,
)
where
import Algebra.Algorithms.FGLM
import Algebra.Algorithms.Groebner
import Algebra.Instances ()
import Algebra.Internal hiding (OLt)
import qualified Algebra.Matrix as AM
import Algebra.Prelude.Core
import Algebra.Ring.Polynomial.Quotient
import Control.Lens
import Control.Monad.Loops (whileM_)
import Control.Monad.Random hiding (next)
import Control.Monad.Reader (runReaderT)
import Control.Monad.ST (runST)
import Data.Complex (Complex (..), magnitude)
import Data.Convertible (Convertible, convert)
import qualified Data.Matrix as M
import Data.Maybe (fromJust)
import Data.Reflection (Reifies)
import Data.STRef.Strict (newSTRef)
import qualified Data.Sized as SV
import qualified Data.Vector as V
import qualified Data.Vector.Mutable as MV
import qualified Numeric.Algebra as NA
import qualified Numeric.LinearAlgebra as LA
import qualified Prelude as P
solveM ::
forall m r ord n.
( Normed r
, Ord r
, MonadRandom m
, Field r
, CoeffRing r
, KnownNat n
, IsMonomialOrder n ord
, Convertible r Double
, 0 < n
) =>
Ideal (OrderedPolynomial r ord n) ->
m [Sized n (Complex Double)]
solveM :: Ideal (OrderedPolynomial r ord n) -> m [Sized n (Complex Double)]
solveM Ideal (OrderedPolynomial r ord n)
ideal =
{-# SCC "solveM" #-}
SNat (n + 1)
-> (KnownNat (n + 1) => m [Sized n (Complex Double)])
-> m [Sized n (Complex Double)]
forall (n :: Nat) r. SNat n -> (KnownNat n => r) -> r
withKnownNat (SNat n -> SNat (n + 1)
forall (n :: Nat). SNat n -> SNat (Succ n)
sSucc (SNat n
forall (n :: Nat). KnownNat n => SNat n
sNat :: SNat n)) ((KnownNat (n + 1) => m [Sized n (Complex Double)])
-> m [Sized n (Complex Double)])
-> (KnownNat (n + 1) => m [Sized n (Complex Double)])
-> m [Sized n (Complex Double)]
forall a b. (a -> b) -> a -> b
$
Ideal (OrderedPolynomial r ord n)
-> (forall ideal.
Reifies ideal (QIdeal (OrderedPolynomial r ord n)) =>
Proxy ideal -> m [Sized n (Complex Double)])
-> m [Sized n (Complex Double)]
forall poly a.
(IsOrderedPolynomial poly, Field (Coefficient poly)) =>
Ideal poly
-> (forall ideal. Reifies ideal (QIdeal poly) => Proxy ideal -> a)
-> a
reifyQuotient (Ideal (OrderedPolynomial r ord n)
-> Ideal (OrderedPolynomial r ord n)
forall r ord (n :: Nat).
(Ord r, CoeffRing r, KnownNat n, Field r, IsMonomialOrder n ord) =>
Ideal (OrderedPolynomial r ord n)
-> Ideal (OrderedPolynomial r ord n)
radical Ideal (OrderedPolynomial r ord n)
ideal) ((forall ideal.
Reifies ideal (QIdeal (OrderedPolynomial r ord n)) =>
Proxy ideal -> m [Sized n (Complex Double)])
-> m [Sized n (Complex Double)])
-> (forall ideal.
Reifies ideal (QIdeal (OrderedPolynomial r ord n)) =>
Proxy ideal -> m [Sized n (Complex Double)])
-> m [Sized n (Complex Double)]
forall a b. (a -> b) -> a -> b
$ \Proxy ideal
pxy ->
case Proxy ideal -> Maybe [Quotient (OrderedPolynomial r ord n) ideal]
forall k poly (ideal :: k).
(IsOrderedPolynomial poly, Field (Coefficient poly),
Reifies ideal (QIdeal poly)) =>
Proxy ideal -> Maybe [Quotient poly ideal]
standardMonomials' Proxy ideal
pxy of
Just [Quotient (OrderedPolynomial r ord n) ideal]
bs -> Integer -> Int -> m [Sized n (Complex Double)]
step Integer
10 ([Quotient (OrderedPolynomial r ord n) ideal] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Quotient (OrderedPolynomial r ord n) ideal]
bs)
Maybe [Quotient (OrderedPolynomial r ord n) ideal]
Nothing -> [Char] -> m [Sized n (Complex Double)]
forall a. HasCallStack => [Char] -> a
error [Char]
"Given ideal is not zero-dimensional"
where
step :: Integer -> Int -> m [Sized n (Complex Double)]
step Integer
bd Int
len =
{-# SCC "solveM/step" #-}
do
[Integer]
coeffs <-
{-# SCC "solveM/coeff-gen" #-}
Int -> m Integer -> m [Integer]
forall (m :: * -> *) a. Applicative m => Int -> m a -> m [a]
replicateM (SNat (n + 1) -> Int
forall (n :: Nat). SNat n -> Int
sNatToInt (SNat n -> SNat (n + 1)
forall (n :: Nat). SNat n -> SNat (Succ n)
sSucc (SNat n
forall (n :: Nat). KnownNat n => SNat n
sNat :: SNat n))) (m Integer -> m [Integer]) -> m Integer -> m [Integer]
forall a b. (a -> b) -> a -> b
$ (Integer, Integer) -> m Integer
forall (m :: * -> *) a. (MonadRandom m, Random a) => (a, a) -> m a
getRandomR (- Integer
bd, Integer
bd)
let vars :: [OrderedPolynomial k ord n]
vars = OrderedPolynomial k ord n
forall r. Unital r => r
one OrderedPolynomial k ord n
-> [OrderedPolynomial k ord n] -> [OrderedPolynomial k ord n]
forall a. a -> [a] -> [a]
: Sized Vector n (OrderedPolynomial k ord n)
-> [OrderedPolynomial k ord n]
forall (f :: * -> *) (n :: Nat) a.
(CFoldable f, Dom f a) =>
Sized f n a -> [a]
SV.toList Sized Vector n (OrderedPolynomial k ord n)
forall k ord (n :: Nat).
(IsMonomialOrder n ord, CoeffRing k, KnownNat n) =>
Sized n (OrderedPolynomial k ord n)
allVars
f :: OrderedPolynomial r ord n
f = [OrderedPolynomial r ord n] -> OrderedPolynomial r ord n
forall (f :: * -> *) m. (Foldable f, Monoidal m) => f m -> m
sum ([OrderedPolynomial r ord n] -> OrderedPolynomial r ord n)
-> [OrderedPolynomial r ord n] -> OrderedPolynomial r ord n
forall a b. (a -> b) -> a -> b
$ (r -> OrderedPolynomial r ord n -> OrderedPolynomial r ord n)
-> [r]
-> [OrderedPolynomial r ord n]
-> [OrderedPolynomial r ord n]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith r -> OrderedPolynomial r ord n -> OrderedPolynomial r ord n
forall r m. Module (Scalar r) m => r -> m -> m
(.*.) ((Integer -> r) -> [Integer] -> [r]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
map (Integer -> r
forall r. Ring r => Integer -> r
NA.fromInteger :: Integer -> r) [Integer]
coeffs) [OrderedPolynomial r ord n]
forall (n :: Nat) ord k.
(IsMonomialOrder n ord, DecidableZero k, Ring k, Commutative k,
Eq k, KnownNat n) =>
[OrderedPolynomial k ord n]
vars
case OrderedPolynomial r ord n
-> Ideal (OrderedPolynomial r ord n)
-> Maybe [Sized n (Complex Double)]
forall r (n :: Nat) ord.
(DecidableZero r, Normed r, Ord r, Field r, CoeffRing r, 0 < n,
IsMonomialOrder n ord, KnownNat n, Convertible r Double) =>
OrderedPolynomial r ord n
-> Ideal (OrderedPolynomial r ord n)
-> Maybe [Sized n (Complex Double)]
solveWith OrderedPolynomial r ord n
f Ideal (OrderedPolynomial r ord n)
ideal of
Maybe [Sized n (Complex Double)]
Nothing -> Integer -> Int -> m [Sized n (Complex Double)]
step (Integer
bd Integer -> Integer -> Integer
forall r. Multiplicative r => r -> r -> r
* Integer
2) Int
len
Just [Sized n (Complex Double)]
sols -> [Sized n (Complex Double)] -> m [Sized n (Complex Double)]
forall (m :: * -> *) a. Monad m => a -> m a
return [Sized n (Complex Double)]
sols
solveWith ::
forall r n ord.
( DecidableZero r
, Normed r
, Ord r
, Field r
, CoeffRing r
, 0 < n
, IsMonomialOrder n ord
, KnownNat n
, Convertible r Double
) =>
OrderedPolynomial r ord n ->
Ideal (OrderedPolynomial r ord n) ->
Maybe [Sized n (Complex Double)]
solveWith :: OrderedPolynomial r ord n
-> Ideal (OrderedPolynomial r ord n)
-> Maybe [Sized n (Complex Double)]
solveWith OrderedPolynomial r ord n
f0 Ideal (OrderedPolynomial r ord n)
i0 =
{-# SCC "solveWith" #-}
SNat (n + 1)
-> (KnownNat (n + 1) => Maybe [Sized n (Complex Double)])
-> Maybe [Sized n (Complex Double)]
forall (n :: Nat) r. SNat n -> (KnownNat n => r) -> r
withKnownNat (SNat n -> SNat (n + 1)
forall (n :: Nat). SNat n -> SNat (Succ n)
sSucc (SNat n
forall (n :: Nat). KnownNat n => SNat n
sNat :: SNat n)) ((KnownNat (n + 1) => Maybe [Sized n (Complex Double)])
-> Maybe [Sized n (Complex Double)])
-> (KnownNat (n + 1) => Maybe [Sized n (Complex Double)])
-> Maybe [Sized n (Complex Double)]
forall a b. (a -> b) -> a -> b
$
Ideal (OrderedPolynomial r ord n)
-> (forall ideal.
Reifies ideal (QIdeal (OrderedPolynomial r ord n)) =>
Proxy ideal -> Maybe [Sized n (Complex Double)])
-> Maybe [Sized n (Complex Double)]
forall poly a.
(IsOrderedPolynomial poly, Field (Coefficient poly)) =>
Ideal poly
-> (forall ideal. Reifies ideal (QIdeal poly) => Proxy ideal -> a)
-> a
reifyQuotient (Ideal (OrderedPolynomial r ord n)
-> Ideal (OrderedPolynomial r ord n)
forall r ord (n :: Nat).
(Ord r, CoeffRing r, KnownNat n, Field r, IsMonomialOrder n ord) =>
Ideal (OrderedPolynomial r ord n)
-> Ideal (OrderedPolynomial r ord n)
radical Ideal (OrderedPolynomial r ord n)
i0) ((forall ideal.
Reifies ideal (QIdeal (OrderedPolynomial r ord n)) =>
Proxy ideal -> Maybe [Sized n (Complex Double)])
-> Maybe [Sized n (Complex Double)])
-> (forall ideal.
Reifies ideal (QIdeal (OrderedPolynomial r ord n)) =>
Proxy ideal -> Maybe [Sized n (Complex Double)])
-> Maybe [Sized n (Complex Double)]
forall a b. (a -> b) -> a -> b
$ \Proxy ideal
pxy ->
let ideal :: [OrderedPolynomial r ord n]
ideal = Proxy ideal -> [OrderedPolynomial r ord n]
forall k (ideal :: k) poly.
Reifies ideal (QIdeal poly) =>
Proxy ideal -> [poly]
gBasis' Proxy ideal
pxy
Just [OrderedMonomial ord n]
base = (Quotient (OrderedPolynomial r ord n) ideal
-> OrderedMonomial ord n)
-> [Quotient (OrderedPolynomial r ord n) ideal]
-> [OrderedMonomial ord n]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
map (OrderedPolynomial r ord n -> OrderedMonomial ord n
forall poly.
IsOrderedPolynomial poly =>
poly -> OrderedMonomial (MOrder poly) (Arity poly)
leadingMonomial (OrderedPolynomial r ord n -> OrderedMonomial ord n)
-> (Quotient (OrderedPolynomial r ord n) ideal
-> OrderedPolynomial r ord n)
-> Quotient (OrderedPolynomial r ord n) ideal
-> OrderedMonomial ord n
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Quotient (OrderedPolynomial r ord n) ideal
-> OrderedPolynomial r ord n
forall k poly (ideal :: k). Quotient poly ideal -> poly
quotRepr) ([Quotient (OrderedPolynomial r ord n) ideal]
-> [OrderedMonomial ord n])
-> Maybe [Quotient (OrderedPolynomial r ord n) ideal]
-> Maybe [OrderedMonomial ord n]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Proxy ideal -> Maybe [Quotient (OrderedPolynomial r ord n) ideal]
forall k poly (ideal :: k).
(IsOrderedPolynomial poly, Field (Coefficient poly),
Reifies ideal (QIdeal poly)) =>
Proxy ideal -> Maybe [Quotient poly ideal]
standardMonomials' Proxy ideal
pxy
in case {-# SCC "findOne" #-} OrderedMonomial ord n -> [OrderedMonomial ord n] -> Maybe Int
forall a. Eq a => a -> [a] -> Maybe Int
elemIndex OrderedMonomial ord n
forall r. Unital r => r
one [OrderedMonomial ord n]
base of
Maybe Int
Nothing -> [Sized n (Complex Double)] -> Maybe [Sized n (Complex Double)]
forall a. a -> Maybe a
Just []
Just Int
cind ->
let f :: Quotient (OrderedPolynomial r ord n) ideal
f = Proxy ideal
-> OrderedPolynomial r ord n
-> Quotient (OrderedPolynomial r ord n) ideal
forall k poly (ideal :: k).
(IsOrderedPolynomial poly, Field (Coefficient poly),
Reifies ideal (QIdeal poly)) =>
Proxy ideal -> poly -> Quotient poly ideal
modIdeal' Proxy ideal
pxy OrderedPolynomial r ord n
f0
vars :: [(Ordinal n, OrderedMonomial ord n)]
vars =
((Ordinal n, OrderedMonomial ord n)
-> (Ordinal n, OrderedMonomial ord n) -> Ordering)
-> [(Ordinal n, OrderedMonomial ord n)]
-> [(Ordinal n, OrderedMonomial ord n)]
forall a. (a -> a -> Ordering) -> [a] -> [a]
sortBy (((Ordinal n, OrderedMonomial ord n)
-> (Ordinal n, OrderedMonomial ord n) -> Ordering)
-> (Ordinal n, OrderedMonomial ord n)
-> (Ordinal n, OrderedMonomial ord n)
-> Ordering
forall a b c. (a -> b -> c) -> b -> a -> c
flip (((Ordinal n, OrderedMonomial ord n)
-> (Ordinal n, OrderedMonomial ord n) -> Ordering)
-> (Ordinal n, OrderedMonomial ord n)
-> (Ordinal n, OrderedMonomial ord n)
-> Ordering)
-> ((Ordinal n, OrderedMonomial ord n)
-> (Ordinal n, OrderedMonomial ord n) -> Ordering)
-> (Ordinal n, OrderedMonomial ord n)
-> (Ordinal n, OrderedMonomial ord n)
-> Ordering
forall a b. (a -> b) -> a -> b
$ ((Ordinal n, OrderedMonomial ord n) -> OrderedMonomial ord n)
-> (Ordinal n, OrderedMonomial ord n)
-> (Ordinal n, OrderedMonomial ord n)
-> Ordering
forall a b. Ord a => (b -> a) -> b -> b -> Ordering
comparing (Ordinal n, OrderedMonomial ord n) -> OrderedMonomial ord n
forall a b. (a, b) -> b
snd) ([(Ordinal n, OrderedMonomial ord n)]
-> [(Ordinal n, OrderedMonomial ord n)])
-> [(Ordinal n, OrderedMonomial ord n)]
-> [(Ordinal n, OrderedMonomial ord n)]
forall a b. (a -> b) -> a -> b
$
(Ordinal n -> (Ordinal n, OrderedMonomial ord n))
-> [Ordinal n] -> [(Ordinal n, OrderedMonomial ord n)]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
map (\Ordinal n
on -> (Ordinal n
on, OrderedPolynomial r ord n
-> OrderedMonomial
(MOrder (OrderedPolynomial r ord n))
(Arity (OrderedPolynomial r ord n))
forall poly.
IsOrderedPolynomial poly =>
poly -> OrderedMonomial (MOrder poly) (Arity poly)
leadingMonomial (OrderedPolynomial r ord n
-> OrderedMonomial
(MOrder (OrderedPolynomial r ord n))
(Arity (OrderedPolynomial r ord n)))
-> OrderedPolynomial r ord n
-> OrderedMonomial
(MOrder (OrderedPolynomial r ord n))
(Arity (OrderedPolynomial r ord n))
forall a b. (a -> b) -> a -> b
$ Ordinal (Arity (OrderedPolynomial r ord n))
-> OrderedPolynomial r ord n
forall poly. IsPolynomial poly => Ordinal (Arity poly) -> poly
var Ordinal n
Ordinal (Arity (OrderedPolynomial r ord n))
on OrderedPolynomial r ord n
-> OrderedPolynomial r ord n -> OrderedPolynomial r ord n
forall a. a -> a -> a
`asTypeOf` OrderedPolynomial r ord n
f0)) ([Ordinal n] -> [(Ordinal n, OrderedMonomial ord n)])
-> [Ordinal n] -> [(Ordinal n, OrderedMonomial ord n)]
forall a b. (a -> b) -> a -> b
$
SNat n -> [Ordinal n]
forall (n :: Nat). SNat n -> [Ordinal n]
enumOrdinal (SNat n -> [Ordinal n]) -> SNat n -> [Ordinal n]
forall a b. (a -> b) -> a -> b
$ OrderedPolynomial r ord n
-> SNat (Arity (OrderedPolynomial r ord n))
forall poly. IsPolynomial poly => poly -> SNat (Arity poly)
sArity' OrderedPolynomial r ord n
f0
inds :: [(Ordinal n,
Either (OrderedPolynomial (Complex Double) ord n) Int)]
inds = (((Ordinal n, OrderedMonomial ord n)
-> (Ordinal n,
Either (OrderedPolynomial (Complex Double) ord n) Int))
-> [(Ordinal n, OrderedMonomial ord n)]
-> [(Ordinal n,
Either (OrderedPolynomial (Complex Double) ord n) Int)])
-> [(Ordinal n, OrderedMonomial ord n)]
-> ((Ordinal n, OrderedMonomial ord n)
-> (Ordinal n,
Either (OrderedPolynomial (Complex Double) ord n) Int))
-> [(Ordinal n,
Either (OrderedPolynomial (Complex Double) ord n) Int)]
forall a b c. (a -> b -> c) -> b -> a -> c
flip ((Ordinal n, OrderedMonomial ord n)
-> (Ordinal n,
Either (OrderedPolynomial (Complex Double) ord n) Int))
-> [(Ordinal n, OrderedMonomial ord n)]
-> [(Ordinal n,
Either (OrderedPolynomial (Complex Double) ord n) Int)]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
map [(Ordinal n, OrderedMonomial ord n)]
vars (((Ordinal n, OrderedMonomial ord n)
-> (Ordinal n,
Either (OrderedPolynomial (Complex Double) ord n) Int))
-> [(Ordinal n,
Either (OrderedPolynomial (Complex Double) ord n) Int)])
-> ((Ordinal n, OrderedMonomial ord n)
-> (Ordinal n,
Either (OrderedPolynomial (Complex Double) ord n) Int))
-> [(Ordinal n,
Either (OrderedPolynomial (Complex Double) ord n) Int)]
forall a b. (a -> b) -> a -> b
$
(OrderedMonomial ord n
-> Either (OrderedPolynomial (Complex Double) ord n) Int)
-> (Ordinal n, OrderedMonomial ord n)
-> (Ordinal n,
Either (OrderedPolynomial (Complex Double) ord n) Int)
forall (a :: * -> * -> *) b c d.
Arrow a =>
a b c -> a (d, b) (d, c)
second ((OrderedMonomial ord n
-> Either (OrderedPolynomial (Complex Double) ord n) Int)
-> (Ordinal n, OrderedMonomial ord n)
-> (Ordinal n,
Either (OrderedPolynomial (Complex Double) ord n) Int))
-> (OrderedMonomial ord n
-> Either (OrderedPolynomial (Complex Double) ord n) Int)
-> (Ordinal n, OrderedMonomial ord n)
-> (Ordinal n,
Either (OrderedPolynomial (Complex Double) ord n) Int)
forall a b. (a -> b) -> a -> b
$ \OrderedMonomial ord n
b ->
case (OrderedMonomial ord n -> Bool)
-> [OrderedMonomial ord n] -> Maybe Int
forall a. (a -> Bool) -> [a] -> Maybe Int
findIndex (OrderedMonomial ord n -> OrderedMonomial ord n -> Bool
forall a. Eq a => a -> a -> Bool
== OrderedMonomial ord n
b) [OrderedMonomial ord n]
base of
Just Int
ind -> Int -> Either (OrderedPolynomial (Complex Double) ord n) Int
forall a b. b -> Either a b
Right Int
ind
Maybe Int
Nothing ->
let Just OrderedPolynomial r ord n
g = (OrderedPolynomial r ord n -> Bool)
-> [OrderedPolynomial r ord n] -> Maybe (OrderedPolynomial r ord n)
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Maybe a
find ((OrderedMonomial ord n -> OrderedMonomial ord n -> Bool
forall a. Eq a => a -> a -> Bool
== OrderedMonomial ord n
b) (OrderedMonomial ord n -> Bool)
-> (OrderedPolynomial r ord n -> OrderedMonomial ord n)
-> OrderedPolynomial r ord n
-> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. OrderedPolynomial r ord n -> OrderedMonomial ord n
forall poly.
IsOrderedPolynomial poly =>
poly -> OrderedMonomial (MOrder poly) (Arity poly)
leadingMonomial) [OrderedPolynomial r ord n]
ideal
r :: Coefficient (OrderedPolynomial r ord n)
r = OrderedPolynomial r ord n
-> Coefficient (OrderedPolynomial r ord n)
forall poly. IsOrderedPolynomial poly => poly -> Coefficient poly
leadingCoeff OrderedPolynomial r ord n
g
answer :: OrderedPolynomial (Complex Double) ord n
answer = (r -> Complex Double)
-> OrderedPolynomial r ord n
-> OrderedPolynomial (Complex Double) ord n
forall (n :: Nat) b ord a.
(KnownNat n, CoeffRing b, IsMonomialOrder n ord) =>
(a -> b) -> OrderedPolynomial a ord n -> OrderedPolynomial b ord n
mapCoeff r -> Complex Double
forall a. Convertible a Double => a -> Complex Double
toComplex (OrderedPolynomial r ord n
-> OrderedPolynomial (Complex Double) ord n)
-> OrderedPolynomial r ord n
-> OrderedPolynomial (Complex Double) ord n
forall a b. (a -> b) -> a -> b
$ Coefficient (OrderedPolynomial r ord n)
-> OrderedPolynomial r ord n
forall poly. IsPolynomial poly => Coefficient poly -> poly
injectCoeff (r -> r
forall r. Division r => r -> r
recip r
Coefficient (OrderedPolynomial r ord n)
r) OrderedPolynomial r ord n
-> OrderedPolynomial r ord n -> OrderedPolynomial r ord n
forall r. Multiplicative r => r -> r -> r
* ((Coefficient (OrderedPolynomial r ord n),
OrderedMonomial
(MOrder (OrderedPolynomial r ord n))
(Arity (OrderedPolynomial r ord n)))
-> OrderedPolynomial r ord n
forall poly.
IsOrderedPolynomial poly =>
(Coefficient poly, OrderedMonomial (MOrder poly) (Arity poly))
-> poly
toPolynomial (OrderedPolynomial r ord n
-> (Coefficient (OrderedPolynomial r ord n),
OrderedMonomial
(MOrder (OrderedPolynomial r ord n))
(Arity (OrderedPolynomial r ord n)))
forall poly.
IsOrderedPolynomial poly =>
poly
-> (Coefficient poly, OrderedMonomial (MOrder poly) (Arity poly))
leadingTerm OrderedPolynomial r ord n
g) OrderedPolynomial r ord n
-> OrderedPolynomial r ord n -> OrderedPolynomial r ord n
forall r. Group r => r -> r -> r
- OrderedPolynomial r ord n
g)
in OrderedPolynomial (Complex Double) ord n
-> Either (OrderedPolynomial (Complex Double) ord n) Int
forall a b. a -> Either a b
Left OrderedPolynomial (Complex Double) ord n
answer
mf :: Matrix (Complex Double)
mf = [[Complex Double]] -> Matrix (Complex Double)
forall (mat :: * -> *) a.
(Matrix mat, Elem mat a) =>
[[a]] -> mat a
AM.fromLists ([[Complex Double]] -> Matrix (Complex Double))
-> [[Complex Double]] -> Matrix (Complex Double)
forall a b. (a -> b) -> a -> b
$ ([r] -> [Complex Double]) -> [[r]] -> [[Complex Double]]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
map ((r -> Complex Double) -> [r] -> [Complex Double]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
map r -> Complex Double
forall a. Convertible a Double => a -> Complex Double
toComplex) ([[r]] -> [[Complex Double]]) -> [[r]] -> [[Complex Double]]
forall a b. (a -> b) -> a -> b
$ Quotient (OrderedPolynomial r ord n) ideal -> [[r]]
forall k t (n :: Nat) order (ideal :: k).
(DecidableZero t, Eq t, Field t, KnownNat n,
IsMonomialOrder n order,
Reifies ideal (QIdeal (OrderedPolynomial t order n))) =>
Quotient (OrderedPolynomial t order n) ideal -> [[t]]
matrixRep Quotient (OrderedPolynomial r ord n) ideal
f
(Vector (Complex Double)
_, Matrix (Complex Double)
evecs) = Matrix (Complex Double)
-> (Vector (Complex Double), Matrix (Complex Double))
forall t.
Field t =>
Matrix t -> (Vector (Complex Double), Matrix (Complex Double))
LA.eig (Matrix (Complex Double)
-> (Vector (Complex Double), Matrix (Complex Double)))
-> Matrix (Complex Double)
-> (Vector (Complex Double), Matrix (Complex Double))
forall a b. (a -> b) -> a -> b
$ Matrix (Complex Double) -> Matrix (Complex Double)
forall m mt. Transposable m mt => m -> mt
LA.tr Matrix (Complex Double)
mf
calc :: Vector (Complex Double) -> Maybe (Sized n (Complex Double))
calc Vector (Complex Double)
vec =
{-# SCC "calc" #-}
let c :: Complex Double
c = Vector (Complex Double)
vec Vector (Complex Double) -> Int -> Complex Double
forall c t. Indexable c t => c -> Int -> t
LA.! Int
cind
phi :: (Ordinal n, Either (OrderedPolynomial (Complex Double) ord n) Int)
-> Sized n (Complex Double) -> Sized n (Complex Double)
phi (Ordinal n
idx, Right Int
nth) Sized n (Complex Double)
acc = Sized n (Complex Double)
acc Sized n (Complex Double)
-> (Sized n (Complex Double) -> Sized n (Complex Double))
-> Sized n (Complex Double)
forall a b. a -> (a -> b) -> b
& Index (Sized n (Complex Double))
-> Traversal'
(Sized n (Complex Double)) (IxValue (Sized n (Complex Double)))
forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix Ordinal n
Index (Sized n (Complex Double))
idx ((Complex Double -> Identity (Complex Double))
-> Sized n (Complex Double) -> Identity (Sized n (Complex Double)))
-> Complex Double
-> Sized n (Complex Double)
-> Sized n (Complex Double)
forall s t a b. ASetter s t a b -> b -> s -> t
.~ (Vector (Complex Double)
vec Vector (Complex Double) -> Int -> Complex Double
forall c t. Indexable c t => c -> Int -> t
LA.! Int
nth) Complex Double -> Complex Double -> Complex Double
forall a. Fractional a => a -> a -> a
P./ Complex Double
c
phi (Ordinal n
idx, Left OrderedPolynomial (Complex Double) ord n
g) Sized n (Complex Double)
acc = Sized n (Complex Double)
acc Sized n (Complex Double)
-> (Sized n (Complex Double) -> Sized n (Complex Double))
-> Sized n (Complex Double)
forall a b. a -> (a -> b) -> b
& Index (Sized n (Complex Double))
-> Traversal'
(Sized n (Complex Double)) (IxValue (Sized n (Complex Double)))
forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix Ordinal n
Index (Sized n (Complex Double))
idx ((Complex Double -> Identity (Complex Double))
-> Sized n (Complex Double) -> Identity (Sized n (Complex Double)))
-> Complex Double
-> Sized n (Complex Double)
-> Sized n (Complex Double)
forall s t a b. ASetter s t a b -> b -> s -> t
.~ (Coefficient (OrderedPolynomial (Complex Double) ord n)
-> Complex Double -> Complex Double)
-> Sized
(Arity (OrderedPolynomial (Complex Double) ord n)) (Complex Double)
-> OrderedPolynomial (Complex Double) ord n
-> Complex Double
forall poly m.
(IsPolynomial poly, Ring m) =>
(Coefficient poly -> m -> m) -> Sized (Arity poly) m -> poly -> m
substWith Coefficient (OrderedPolynomial (Complex Double) ord n)
-> Complex Double -> Complex Double
forall r. Multiplicative r => r -> r -> r
(*) Sized n (Complex Double)
Sized
(Arity (OrderedPolynomial (Complex Double) ord n)) (Complex Double)
acc OrderedPolynomial (Complex Double) ord n
g
in if Complex Double
c Complex Double -> Complex Double -> Bool
forall a. Eq a => a -> a -> Bool
== Complex Double
0
then Maybe (Sized n (Complex Double))
forall a. Maybe a
Nothing
else Sized n (Complex Double) -> Maybe (Sized n (Complex Double))
forall a. a -> Maybe a
Just (Sized n (Complex Double) -> Maybe (Sized n (Complex Double)))
-> Sized n (Complex Double) -> Maybe (Sized n (Complex Double))
forall a b. (a -> b) -> a -> b
$ ((Ordinal n, Either (OrderedPolynomial (Complex Double) ord n) Int)
-> Sized n (Complex Double) -> Sized n (Complex Double))
-> Sized n (Complex Double)
-> [(Ordinal n,
Either (OrderedPolynomial (Complex Double) ord n) Int)]
-> Sized n (Complex Double)
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr ({-# SCC "rewrite-answer" #-} (Ordinal n, Either (OrderedPolynomial (Complex Double) ord n) Int)
-> Sized n (Complex Double) -> Sized n (Complex Double)
phi) (SNat n -> Complex Double -> Sized n (Complex Double)
forall (f :: * -> *) (n :: Nat) a.
(CFreeMonoid f, Dom f a) =>
SNat n -> a -> Sized f n a
SV.replicate (OrderedPolynomial r ord n
-> SNat (Arity (OrderedPolynomial r ord n))
forall poly. IsPolynomial poly => poly -> SNat (Arity poly)
sArity' OrderedPolynomial r ord n
f0) ([Char] -> Complex Double
forall a. HasCallStack => [Char] -> a
error [Char]
"indec!")) [(Ordinal n,
Either (OrderedPolynomial (Complex Double) ord n) Int)]
inds
in [Maybe (Sized n (Complex Double))]
-> Maybe [Sized n (Complex Double)]
forall (t :: * -> *) (m :: * -> *) a.
(Traversable t, Monad m) =>
t (m a) -> m (t a)
sequence ([Maybe (Sized n (Complex Double))]
-> Maybe [Sized n (Complex Double)])
-> [Maybe (Sized n (Complex Double))]
-> Maybe [Sized n (Complex Double)]
forall a b. (a -> b) -> a -> b
$ (Vector (Complex Double) -> Maybe (Sized n (Complex Double)))
-> [Vector (Complex Double)] -> [Maybe (Sized n (Complex Double))]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
map Vector (Complex Double) -> Maybe (Sized n (Complex Double))
calc ([Vector (Complex Double)] -> [Maybe (Sized n (Complex Double))])
-> [Vector (Complex Double)] -> [Maybe (Sized n (Complex Double))]
forall a b. (a -> b) -> a -> b
$ Matrix (Complex Double) -> [Vector (Complex Double)]
forall t. Element t => Matrix t -> [Vector t]
LA.toColumns Matrix (Complex Double)
evecs
solve' ::
forall r n ord.
( Field r
, CoeffRing r
, KnownNat n
, 0 < n
, IsMonomialOrder n ord
, Convertible r Double
) =>
Double ->
Ideal (OrderedPolynomial r ord n) ->
[Sized n (Complex Double)]
solve' :: Double
-> Ideal (OrderedPolynomial r ord n) -> [Sized n (Complex Double)]
solve' Double
err Ideal (OrderedPolynomial r ord n)
ideal =
Ideal (OrderedPolynomial r ord n)
-> (forall ideal.
Reifies ideal (QIdeal (OrderedPolynomial r ord n)) =>
Proxy ideal -> [Sized n (Complex Double)])
-> [Sized n (Complex Double)]
forall poly a.
(IsOrderedPolynomial poly, Field (Coefficient poly)) =>
Ideal poly
-> (forall ideal. Reifies ideal (QIdeal poly) => Proxy ideal -> a)
-> a
reifyQuotient Ideal (OrderedPolynomial r ord n)
ideal ((forall ideal.
Reifies ideal (QIdeal (OrderedPolynomial r ord n)) =>
Proxy ideal -> [Sized n (Complex Double)])
-> [Sized n (Complex Double)])
-> (forall ideal.
Reifies ideal (QIdeal (OrderedPolynomial r ord n)) =>
Proxy ideal -> [Sized n (Complex Double)])
-> [Sized n (Complex Double)]
forall a b. (a -> b) -> a -> b
$ \Proxy ideal
ii ->
if Proxy ideal -> [OrderedPolynomial r ord n]
forall k (ideal :: k) poly.
Reifies ideal (QIdeal poly) =>
Proxy ideal -> [poly]
gBasis' Proxy ideal
ii [OrderedPolynomial r ord n] -> [OrderedPolynomial r ord n] -> Bool
forall a. Eq a => a -> a -> Bool
== [OrderedPolynomial r ord n
forall r. Unital r => r
one]
then []
else
let vs :: [[Complex Double]]
vs =
(OrderedPolynomial r ord n -> [Complex Double])
-> [OrderedPolynomial r ord n] -> [[Complex Double]]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
map ([Complex Double] -> [Complex Double]
forall a. Eq a => [a] -> [a]
nub ([Complex Double] -> [Complex Double])
-> (OrderedPolynomial r ord n -> [Complex Double])
-> OrderedPolynomial r ord n
-> [Complex Double]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector (Complex Double) -> [Complex Double]
forall a. Storable a => Vector a -> [a]
LA.toList (Vector (Complex Double) -> [Complex Double])
-> (OrderedPolynomial r ord n -> Vector (Complex Double))
-> OrderedPolynomial r ord n
-> [Complex Double]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Matrix (Complex Double) -> Vector (Complex Double)
forall t. Field t => Matrix t -> Vector (Complex Double)
LA.eigenvalues (Matrix (Complex Double) -> Vector (Complex Double))
-> (OrderedPolynomial r ord n -> Matrix (Complex Double))
-> OrderedPolynomial r ord n
-> Vector (Complex Double)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [[Complex Double]] -> Matrix (Complex Double)
forall (mat :: * -> *) a.
(Matrix mat, Elem mat a) =>
[[a]] -> mat a
AM.fromLists ([[Complex Double]] -> Matrix (Complex Double))
-> (OrderedPolynomial r ord n -> [[Complex Double]])
-> OrderedPolynomial r ord n
-> Matrix (Complex Double)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ([r] -> [Complex Double]) -> [[r]] -> [[Complex Double]]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
map ((r -> Complex Double) -> [r] -> [Complex Double]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
map r -> Complex Double
forall a. Convertible a Double => a -> Complex Double
toComplex) ([[r]] -> [[Complex Double]])
-> (OrderedPolynomial r ord n -> [[r]])
-> OrderedPolynomial r ord n
-> [[Complex Double]]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Quotient (OrderedPolynomial r ord n) ideal -> [[r]]
forall k t (n :: Nat) order (ideal :: k).
(DecidableZero t, Eq t, Field t, KnownNat n,
IsMonomialOrder n order,
Reifies ideal (QIdeal (OrderedPolynomial t order n))) =>
Quotient (OrderedPolynomial t order n) ideal -> [[t]]
matrixRep (Quotient (OrderedPolynomial r ord n) ideal -> [[r]])
-> (OrderedPolynomial r ord n
-> Quotient (OrderedPolynomial r ord n) ideal)
-> OrderedPolynomial r ord n
-> [[r]]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Proxy ideal
-> OrderedPolynomial r ord n
-> Quotient (OrderedPolynomial r ord n) ideal
forall k poly (ideal :: k).
(IsOrderedPolynomial poly, Field (Coefficient poly),
Reifies ideal (QIdeal poly)) =>
Proxy ideal -> poly -> Quotient poly ideal
modIdeal' Proxy ideal
ii) ([OrderedPolynomial r ord n] -> [[Complex Double]])
-> [OrderedPolynomial r ord n] -> [[Complex Double]]
forall a b. (a -> b) -> a -> b
$
Sized Vector n (OrderedPolynomial r ord n)
-> [OrderedPolynomial r ord n]
forall (f :: * -> *) (n :: Nat) a.
(CFoldable f, Dom f a) =>
Sized f n a -> [a]
SV.toList Sized Vector n (OrderedPolynomial r ord n)
forall k ord (n :: Nat).
(IsMonomialOrder n ord, CoeffRing k, KnownNat n) =>
Sized n (OrderedPolynomial k ord n)
allVars
mul :: a -> Complex Double -> Complex Double
mul a
p Complex Double
q = a -> Complex Double
forall a. Convertible a Double => a -> Complex Double
toComplex a
p Complex Double -> Complex Double -> Complex Double
forall r. Multiplicative r => r -> r -> r
* Complex Double
q
in [ Sized n (Complex Double)
xs
| [Complex Double]
xs0 <- [[Complex Double]] -> [[Complex Double]]
forall (t :: * -> *) (m :: * -> *) a.
(Traversable t, Monad m) =>
t (m a) -> m (t a)
sequence [[Complex Double]]
vs
, let xs :: Sized n (Complex Double)
xs = [Complex Double] -> Sized n (Complex Double)
forall (f :: * -> *) (n :: Nat) a.
(KnownNat n, CFreeMonoid f, Dom f a) =>
[a] -> Sized f n a
SV.unsafeFromList' [Complex Double]
xs0
, (OrderedPolynomial r ord n -> Bool)
-> [OrderedPolynomial r ord n] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all ((Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
err) (Double -> Bool)
-> (OrderedPolynomial r ord n -> Double)
-> OrderedPolynomial r ord n
-> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Complex Double -> Double
forall a. RealFloat a => Complex a -> a
magnitude (Complex Double -> Double)
-> (OrderedPolynomial r ord n -> Complex Double)
-> OrderedPolynomial r ord n
-> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Coefficient (OrderedPolynomial r ord n)
-> Complex Double -> Complex Double)
-> Sized (Arity (OrderedPolynomial r ord n)) (Complex Double)
-> OrderedPolynomial r ord n
-> Complex Double
forall poly m.
(IsPolynomial poly, Ring m) =>
(Coefficient poly -> m -> m) -> Sized (Arity poly) m -> poly -> m
substWith Coefficient (OrderedPolynomial r ord n)
-> Complex Double -> Complex Double
forall a.
Convertible a Double =>
a -> Complex Double -> Complex Double
mul Sized n (Complex Double)
Sized (Arity (OrderedPolynomial r ord n)) (Complex Double)
xs) ([OrderedPolynomial r ord n] -> Bool)
-> [OrderedPolynomial r ord n] -> Bool
forall a b. (a -> b) -> a -> b
$ Ideal (OrderedPolynomial r ord n) -> [OrderedPolynomial r ord n]
forall r. Ideal r -> [r]
generators Ideal (OrderedPolynomial r ord n)
ideal
]
subspMatrix ::
forall r n ord.
(Show r, Ord r, Field r, CoeffRing r, KnownNat n, IsMonomialOrder n ord) =>
Ordinal n ->
Ideal (OrderedPolynomial r ord n) ->
M.Matrix r
subspMatrix :: Ordinal n -> Ideal (OrderedPolynomial r ord n) -> Matrix r
subspMatrix Ordinal n
on Ideal (OrderedPolynomial r ord n)
ideal =
let poly :: OrderedPolynomial r ord n
poly = Ordinal n
-> Ideal (OrderedPolynomial r ord n) -> OrderedPolynomial r ord n
forall r ord (n :: Nat).
(Ord r, Field r, CoeffRing r, KnownNat n, IsMonomialOrder n ord) =>
Ordinal n
-> Ideal (OrderedPolynomial r ord n) -> OrderedPolynomial r ord n
univPoly Ordinal n
on Ideal (OrderedPolynomial r ord n)
ideal
v :: OrderedPolynomial r ord n
v = Ordinal (Arity (OrderedPolynomial r ord n))
-> OrderedPolynomial r ord n
forall poly. IsPolynomial poly => Ordinal (Arity poly) -> poly
var Ordinal n
Ordinal (Arity (OrderedPolynomial r ord n))
on :: OrderedPolynomial r ord n
dim :: Int
dim = OrderedPolynomial r ord n -> Int
forall poly. IsPolynomial poly => poly -> Int
totalDegree' OrderedPolynomial r ord n
poly
cfs :: [r]
cfs = [r -> r
forall r. Group r => r -> r
negate (r -> r) -> r -> r
forall a b. (a -> b) -> a -> b
$ OrderedMonomial
(MOrder (OrderedPolynomial r ord n))
(Arity (OrderedPolynomial r ord n))
-> OrderedPolynomial r ord n
-> Coefficient (OrderedPolynomial r ord n)
forall poly.
IsOrderedPolynomial poly =>
OrderedMonomial (MOrder poly) (Arity poly)
-> poly -> Coefficient poly
coeff (OrderedPolynomial r ord n
-> OrderedMonomial
(MOrder (OrderedPolynomial r ord n))
(Arity (OrderedPolynomial r ord n))
forall poly.
IsOrderedPolynomial poly =>
poly -> OrderedMonomial (MOrder poly) (Arity poly)
leadingMonomial (OrderedPolynomial r ord n
-> OrderedMonomial
(MOrder (OrderedPolynomial r ord n))
(Arity (OrderedPolynomial r ord n)))
-> OrderedPolynomial r ord n
-> OrderedMonomial
(MOrder (OrderedPolynomial r ord n))
(Arity (OrderedPolynomial r ord n))
forall a b. (a -> b) -> a -> b
$ OrderedPolynomial r ord n -> Natural -> OrderedPolynomial r ord n
forall r. Unital r => r -> Natural -> r
pow OrderedPolynomial r ord n
v (Natural
j :: Natural)) OrderedPolynomial r ord n
poly | Natural
j <- [Natural
0 .. Int -> Natural
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
dim Int -> Int -> Int
forall r. Group r => r -> r -> r
- Int
1)]]
leftTops :: Matrix r
leftTops = [[r]] -> Matrix r
forall a. [[a]] -> Matrix a
M.fromLists [Int -> r -> [r]
forall a. Int -> a -> [a]
replicate (Int
dim Int -> Int -> Int
forall r. Group r => r -> r -> r
- Int
1) r
forall m. Monoidal m => m
zero]
leftBots :: Matrix r
leftBots = (WrapAlgebra r -> r) -> Matrix (WrapAlgebra r) -> Matrix r
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap WrapAlgebra r -> r
forall a. WrapAlgebra a -> a
unwrapAlgebra (Int -> Matrix (WrapAlgebra r)
forall a. Num a => Int -> Matrix a
M.identity (Int
dim Int -> Int -> Int
forall r. Group r => r -> r -> r
- Int
1))
lfts :: Matrix r
lfts =
Matrix r
leftTops
Matrix r -> Matrix r -> Matrix r
forall a. Matrix a -> Matrix a -> Matrix a
M.<-> Matrix r
leftBots
rights :: Matrix r
rights = Vector r -> Matrix r
forall a. Vector a -> Matrix a
M.colVector ([r] -> Vector r
forall a. [a] -> Vector a
V.fromList [r]
cfs)
in Matrix r
lfts Matrix r -> Matrix r -> Matrix r
forall a. Matrix a -> Matrix a -> Matrix a
M.<|> Matrix r
rights
solveViaCompanion ::
forall r ord n.
(Show r, Ord r, Field r, CoeffRing r, KnownNat n, IsMonomialOrder n ord, Convertible r Double) =>
Double ->
Ideal (OrderedPolynomial r ord n) ->
[Sized n (Complex Double)]
solveViaCompanion :: Double
-> Ideal (OrderedPolynomial r ord n) -> [Sized n (Complex Double)]
solveViaCompanion Double
err Ideal (OrderedPolynomial r ord n)
ideal =
if Ideal (OrderedPolynomial r ord n) -> [OrderedPolynomial r ord n]
forall poly.
(Field (Coefficient poly), IsOrderedPolynomial poly) =>
Ideal poly -> [poly]
calcGroebnerBasis Ideal (OrderedPolynomial r ord n)
ideal [OrderedPolynomial r ord n] -> [OrderedPolynomial r ord n] -> Bool
forall a. Eq a => a -> a -> Bool
== [OrderedPolynomial r ord n
forall r. Unital r => r
one]
then []
else
let vs :: [[Complex Double]]
vs =
(Ordinal n -> [Complex Double])
-> [Ordinal n] -> [[Complex Double]]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
map ([Complex Double] -> [Complex Double]
forall a. Eq a => [a] -> [a]
nub ([Complex Double] -> [Complex Double])
-> (Ordinal n -> [Complex Double]) -> Ordinal n -> [Complex Double]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector (Complex Double) -> [Complex Double]
forall a. Storable a => Vector a -> [a]
LA.toList (Vector (Complex Double) -> [Complex Double])
-> (Ordinal n -> Vector (Complex Double))
-> Ordinal n
-> [Complex Double]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Matrix (Complex Double) -> Vector (Complex Double)
forall t. Field t => Matrix t -> Vector (Complex Double)
LA.eigenvalues (Matrix (Complex Double) -> Vector (Complex Double))
-> (Ordinal n -> Matrix (Complex Double))
-> Ordinal n
-> Vector (Complex Double)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [[Complex Double]] -> Matrix (Complex Double)
forall t. Element t => [[t]] -> Matrix t
LA.fromLists ([[Complex Double]] -> Matrix (Complex Double))
-> (Ordinal n -> [[Complex Double]])
-> Ordinal n
-> Matrix (Complex Double)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Matrix (Complex Double) -> [[Complex Double]]
forall a. Matrix a -> [[a]]
matToLists (Matrix (Complex Double) -> [[Complex Double]])
-> (Ordinal n -> Matrix (Complex Double))
-> Ordinal n
-> [[Complex Double]]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (r -> Complex Double) -> Matrix r -> Matrix (Complex Double)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap r -> Complex Double
forall a. Convertible a Double => a -> Complex Double
toComplex (Matrix r -> Matrix (Complex Double))
-> (Ordinal n -> Matrix r) -> Ordinal n -> Matrix (Complex Double)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Ordinal n -> Ideal (OrderedPolynomial r ord n) -> Matrix r)
-> Ideal (OrderedPolynomial r ord n) -> Ordinal n -> Matrix r
forall a b c. (a -> b -> c) -> b -> a -> c
flip Ordinal n -> Ideal (OrderedPolynomial r ord n) -> Matrix r
forall r (n :: Nat) ord.
(Show r, Ord r, Field r, CoeffRing r, KnownNat n,
IsMonomialOrder n ord) =>
Ordinal n -> Ideal (OrderedPolynomial r ord n) -> Matrix r
subspMatrix Ideal (OrderedPolynomial r ord n)
ideal) ([Ordinal n] -> [[Complex Double]])
-> [Ordinal n] -> [[Complex Double]]
forall a b. (a -> b) -> a -> b
$
SNat n -> [Ordinal n]
forall (n :: Nat). SNat n -> [Ordinal n]
enumOrdinal (SNat n
forall (n :: Nat). KnownNat n => SNat n
sNat :: SNat n)
mul :: a -> Complex Double -> Complex Double
mul a
p Complex Double
q = a -> Complex Double
forall a. Convertible a Double => a -> Complex Double
toComplex a
p Complex Double -> Complex Double -> Complex Double
forall r. Multiplicative r => r -> r -> r
* Complex Double
q
in [ Sized n (Complex Double)
xs
| [Complex Double]
xs0 <- [[Complex Double]] -> [[Complex Double]]
forall (t :: * -> *) (m :: * -> *) a.
(Traversable t, Monad m) =>
t (m a) -> m (t a)
sequence [[Complex Double]]
vs
, let xs :: Sized n (Complex Double)
xs = [Complex Double] -> Sized n (Complex Double)
forall (f :: * -> *) (n :: Nat) a.
(KnownNat n, CFreeMonoid f, Dom f a) =>
[a] -> Sized f n a
SV.unsafeFromList' [Complex Double]
xs0
, (OrderedPolynomial r ord n -> Bool)
-> [OrderedPolynomial r ord n] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all ((Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
err) (Double -> Bool)
-> (OrderedPolynomial r ord n -> Double)
-> OrderedPolynomial r ord n
-> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Complex Double -> Double
forall a. RealFloat a => Complex a -> a
magnitude (Complex Double -> Double)
-> (OrderedPolynomial r ord n -> Complex Double)
-> OrderedPolynomial r ord n
-> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Coefficient (OrderedPolynomial r ord n)
-> Complex Double -> Complex Double)
-> Sized (Arity (OrderedPolynomial r ord n)) (Complex Double)
-> OrderedPolynomial r ord n
-> Complex Double
forall poly m.
(IsPolynomial poly, Ring m) =>
(Coefficient poly -> m -> m) -> Sized (Arity poly) m -> poly -> m
substWith Coefficient (OrderedPolynomial r ord n)
-> Complex Double -> Complex Double
forall a.
Convertible a Double =>
a -> Complex Double -> Complex Double
mul Sized n (Complex Double)
Sized (Arity (OrderedPolynomial r ord n)) (Complex Double)
xs) ([OrderedPolynomial r ord n] -> Bool)
-> [OrderedPolynomial r ord n] -> Bool
forall a b. (a -> b) -> a -> b
$ Ideal (OrderedPolynomial r ord n) -> [OrderedPolynomial r ord n]
forall r. Ideal r -> [r]
generators Ideal (OrderedPolynomial r ord n)
ideal
]
matToLists :: M.Matrix a -> [[a]]
matToLists :: Matrix a -> [[a]]
matToLists Matrix a
mat = [Vector a -> [a]
forall a. Vector a -> [a]
V.toList (Vector a -> [a]) -> Vector a -> [a]
forall a b. (a -> b) -> a -> b
$ Int -> Matrix a -> Vector a
forall a. Int -> Matrix a -> Vector a
M.getRow Int
i Matrix a
mat | Int
i <- [Int
1 .. Matrix a -> Int
forall a. Matrix a -> Int
M.nrows Matrix a
mat]]
matrixRep ::
( DecidableZero t
, Eq t
, Field t
, KnownNat n
, IsMonomialOrder n order
, Reifies ideal (QIdeal (OrderedPolynomial t order n))
) =>
Quotient (OrderedPolynomial t order n) ideal ->
[[t]]
matrixRep :: Quotient (OrderedPolynomial t order n) ideal -> [[t]]
matrixRep Quotient (OrderedPolynomial t order n) ideal
f =
{-# SCC "matrixRep" #-}
case Maybe [Quotient (OrderedPolynomial t order n) ideal]
forall k poly (ideal :: k).
(IsOrderedPolynomial poly, Field (Coefficient poly),
Reifies ideal (QIdeal poly)) =>
Maybe [Quotient poly ideal]
standardMonomials of
Just [] -> []
Just [Quotient (OrderedPolynomial t order n) ideal]
bases ->
let anss :: [OrderedPolynomial t order n]
anss = (Quotient (OrderedPolynomial t order n) ideal
-> OrderedPolynomial t order n)
-> [Quotient (OrderedPolynomial t order n) ideal]
-> [OrderedPolynomial t order n]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
map (Quotient (OrderedPolynomial t order n) ideal
-> OrderedPolynomial t order n
forall k poly (ideal :: k). Quotient poly ideal -> poly
quotRepr (Quotient (OrderedPolynomial t order n) ideal
-> OrderedPolynomial t order n)
-> (Quotient (OrderedPolynomial t order n) ideal
-> Quotient (OrderedPolynomial t order n) ideal)
-> Quotient (OrderedPolynomial t order n) ideal
-> OrderedPolynomial t order n
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Quotient (OrderedPolynomial t order n) ideal
f Quotient (OrderedPolynomial t order n) ideal
-> Quotient (OrderedPolynomial t order n) ideal
-> Quotient (OrderedPolynomial t order n) ideal
forall r. Multiplicative r => r -> r -> r
*)) [Quotient (OrderedPolynomial t order n) ideal]
bases
in [[t]] -> [[t]]
forall a. [[a]] -> [[a]]
transpose ([[t]] -> [[t]]) -> [[t]] -> [[t]]
forall a b. (a -> b) -> a -> b
$ (OrderedPolynomial t order n -> [t])
-> [OrderedPolynomial t order n] -> [[t]]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
map (\OrderedPolynomial t order n
a -> (Quotient (OrderedPolynomial t order n) ideal -> t)
-> [Quotient (OrderedPolynomial t order n) ideal] -> [t]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
map ((OrderedMonomial order n -> OrderedPolynomial t order n -> t)
-> OrderedPolynomial t order n -> OrderedMonomial order n -> t
forall a b c. (a -> b -> c) -> b -> a -> c
flip OrderedMonomial order n -> OrderedPolynomial t order n -> t
forall poly.
IsOrderedPolynomial poly =>
OrderedMonomial (MOrder poly) (Arity poly)
-> poly -> Coefficient poly
coeff OrderedPolynomial t order n
a (OrderedMonomial order n -> t)
-> (Quotient (OrderedPolynomial t order n) ideal
-> OrderedMonomial order n)
-> Quotient (OrderedPolynomial t order n) ideal
-> t
forall b c a. (b -> c) -> (a -> b) -> a -> c
. OrderedPolynomial t order n -> OrderedMonomial order n
forall poly.
IsOrderedPolynomial poly =>
poly -> OrderedMonomial (MOrder poly) (Arity poly)
leadingMonomial (OrderedPolynomial t order n -> OrderedMonomial order n)
-> (Quotient (OrderedPolynomial t order n) ideal
-> OrderedPolynomial t order n)
-> Quotient (OrderedPolynomial t order n) ideal
-> OrderedMonomial order n
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Quotient (OrderedPolynomial t order n) ideal
-> OrderedPolynomial t order n
forall k poly (ideal :: k). Quotient poly ideal -> poly
quotRepr) [Quotient (OrderedPolynomial t order n) ideal]
bases) [OrderedPolynomial t order n]
anss
Maybe [Quotient (OrderedPolynomial t order n) ideal]
Nothing -> [Char] -> [[t]]
forall a. HasCallStack => [Char] -> a
error [Char]
"Not finite dimension"
toComplex :: Convertible a Double => a -> Complex Double
toComplex :: a -> Complex Double
toComplex a
a = a -> Double
forall a b. Convertible a b => a -> b
convert a
a Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ Double
0
reduction ::
(CoeffRing r, KnownNat n, IsMonomialOrder n ord, Field r) =>
Ordinal n ->
OrderedPolynomial r ord n ->
OrderedPolynomial r ord n
reduction :: Ordinal n -> OrderedPolynomial r ord n -> OrderedPolynomial r ord n
reduction Ordinal n
on OrderedPolynomial r ord n
f =
{-# SCC "reduction" #-}
let df :: OrderedPolynomial r ord n
df = {-# SCC "differentiate" #-} Ordinal (Arity (OrderedPolynomial r ord n))
-> OrderedPolynomial r ord n -> OrderedPolynomial r ord n
forall poly.
IsOrderedPolynomial poly =>
Ordinal (Arity poly) -> poly -> poly
diff Ordinal n
Ordinal (Arity (OrderedPolynomial r ord n))
on OrderedPolynomial r ord n
f
in (OrderedPolynomial r ord n, OrderedPolynomial r ord n)
-> OrderedPolynomial r ord n
forall a b. (a, b) -> b
snd ((OrderedPolynomial r ord n, OrderedPolynomial r ord n)
-> OrderedPolynomial r ord n)
-> (OrderedPolynomial r ord n, OrderedPolynomial r ord n)
-> OrderedPolynomial r ord n
forall a b. (a -> b) -> a -> b
$ [(OrderedPolynomial r ord n, OrderedPolynomial r ord n)]
-> (OrderedPolynomial r ord n, OrderedPolynomial r ord n)
forall a. [a] -> a
head ([(OrderedPolynomial r ord n, OrderedPolynomial r ord n)]
-> (OrderedPolynomial r ord n, OrderedPolynomial r ord n))
-> [(OrderedPolynomial r ord n, OrderedPolynomial r ord n)]
-> (OrderedPolynomial r ord n, OrderedPolynomial r ord n)
forall a b. (a -> b) -> a -> b
$ OrderedPolynomial r ord n
f OrderedPolynomial r ord n
-> [OrderedPolynomial r ord n]
-> [(OrderedPolynomial r ord n, OrderedPolynomial r ord n)]
forall poly.
(IsOrderedPolynomial poly, Field (Coefficient poly)) =>
poly -> [poly] -> [(poly, poly)]
`divPolynomial` Ideal (OrderedPolynomial r ord n) -> [OrderedPolynomial r ord n]
forall poly.
(Field (Coefficient poly), IsOrderedPolynomial poly) =>
Ideal poly -> [poly]
calcGroebnerBasis ([OrderedPolynomial r ord n] -> Ideal (OrderedPolynomial r ord n)
forall r. (DecidableZero r, Monoidal r) => [r] -> Ideal r
toIdeal [OrderedPolynomial r ord n
f, OrderedPolynomial r ord n
df])
univPoly ::
forall r ord n.
(Ord r, Field r, CoeffRing r, KnownNat n, IsMonomialOrder n ord) =>
Ordinal n ->
Ideal (OrderedPolynomial r ord n) ->
OrderedPolynomial r ord n
univPoly :: Ordinal n
-> Ideal (OrderedPolynomial r ord n) -> OrderedPolynomial r ord n
univPoly Ordinal n
nth Ideal (OrderedPolynomial r ord n)
ideal =
{-# SCC "univPoly" #-}
Ideal (OrderedPolynomial r ord n)
-> (forall ideal.
Reifies ideal (QIdeal (OrderedPolynomial r ord n)) =>
Proxy ideal -> OrderedPolynomial r ord n)
-> OrderedPolynomial r ord n
forall poly a.
(IsOrderedPolynomial poly, Field (Coefficient poly)) =>
Ideal poly
-> (forall ideal. Reifies ideal (QIdeal poly) => Proxy ideal -> a)
-> a
reifyQuotient Ideal (OrderedPolynomial r ord n)
ideal ((forall ideal.
Reifies ideal (QIdeal (OrderedPolynomial r ord n)) =>
Proxy ideal -> OrderedPolynomial r ord n)
-> OrderedPolynomial r ord n)
-> (forall ideal.
Reifies ideal (QIdeal (OrderedPolynomial r ord n)) =>
Proxy ideal -> OrderedPolynomial r ord n)
-> OrderedPolynomial r ord n
forall a b. (a -> b) -> a -> b
$ \Proxy ideal
pxy ->
if Proxy ideal -> [OrderedPolynomial r ord n]
forall k (ideal :: k) poly.
Reifies ideal (QIdeal poly) =>
Proxy ideal -> [poly]
gBasis' Proxy ideal
pxy [OrderedPolynomial r ord n] -> [OrderedPolynomial r ord n] -> Bool
forall a. Eq a => a -> a -> Bool
== [OrderedPolynomial r ord n
forall r. Unital r => r
one]
then OrderedPolynomial r ord n
forall r. Unital r => r
one
else
let x :: OrderedPolynomial r ord n
x = Ordinal (Arity (OrderedPolynomial r ord n))
-> OrderedPolynomial r ord n
forall poly. IsPolynomial poly => Ordinal (Arity poly) -> poly
var Ordinal n
Ordinal (Arity (OrderedPolynomial r ord n))
nth
monomDeg :: Natural
monomDeg =
Int -> Natural
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> Natural) -> Int -> Natural
forall a b. (a -> b) -> a -> b
$
[Quotient (OrderedPolynomial r ord n) ideal] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ([Quotient (OrderedPolynomial r ord n) ideal] -> Int)
-> [Quotient (OrderedPolynomial r ord n) ideal] -> Int
forall a b. (a -> b) -> a -> b
$
Maybe [Quotient (OrderedPolynomial r ord n) ideal]
-> [Quotient (OrderedPolynomial r ord n) ideal]
forall a. HasCallStack => Maybe a -> a
fromJust (Maybe [Quotient (OrderedPolynomial r ord n) ideal]
-> [Quotient (OrderedPolynomial r ord n) ideal])
-> Maybe [Quotient (OrderedPolynomial r ord n) ideal]
-> [Quotient (OrderedPolynomial r ord n) ideal]
forall a b. (a -> b) -> a -> b
$
Proxy ideal -> Maybe [Quotient (OrderedPolynomial r ord n) ideal]
forall k poly (ideal :: k).
(IsOrderedPolynomial poly, Field (Coefficient poly),
Reifies ideal (QIdeal poly)) =>
Proxy ideal -> Maybe [Quotient poly ideal]
standardMonomials' Proxy ideal
pxy
Vector (WrapAlgebra r)
p0 : [Vector (WrapAlgebra r)]
pows =
[ (r -> WrapAlgebra r) -> Vector r -> Vector (WrapAlgebra r)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap r -> WrapAlgebra r
forall a. a -> WrapAlgebra a
WrapAlgebra (Vector r -> Vector (WrapAlgebra r))
-> Vector r -> Vector (WrapAlgebra r)
forall a b. (a -> b) -> a -> b
$ Quotient (OrderedPolynomial r ord n) ideal
-> Vector (Coefficient (OrderedPolynomial r ord n))
forall k poly (ideal :: k).
(IsOrderedPolynomial poly, Reifies ideal (QIdeal poly)) =>
Quotient poly ideal -> Vector (Coefficient poly)
vectorRep (Quotient (OrderedPolynomial r ord n) ideal
-> Vector (Coefficient (OrderedPolynomial r ord n)))
-> Quotient (OrderedPolynomial r ord n) ideal
-> Vector (Coefficient (OrderedPolynomial r ord n))
forall a b. (a -> b) -> a -> b
$ Proxy ideal
-> OrderedPolynomial r ord n
-> Quotient (OrderedPolynomial r ord n) ideal
forall k poly (ideal :: k).
(IsOrderedPolynomial poly, Field (Coefficient poly),
Reifies ideal (QIdeal poly)) =>
Proxy ideal -> poly -> Quotient poly ideal
modIdeal' Proxy ideal
pxy (OrderedPolynomial r ord n -> Natural -> OrderedPolynomial r ord n
forall r. Unital r => r -> Natural -> r
pow OrderedPolynomial r ord n
x Natural
i)
| Natural
i <- [Natural
0 :: Natural .. Natural
monomDeg]
]
step :: Matrix (WrapAlgebra r)
-> [Vector (WrapAlgebra r)] -> OrderedPolynomial r ord n
step Matrix (WrapAlgebra r)
m ~(Vector (WrapAlgebra r)
p : [Vector (WrapAlgebra r)]
ps) =
{-# SCC "univPoly/step" #-}
case Matrix (WrapAlgebra r)
-> Vector (WrapAlgebra r) -> Maybe (Vector (WrapAlgebra r))
forall r.
(Ord r, Fractional r) =>
Matrix r -> Vector r -> Maybe (Vector r)
solveLinear Matrix (WrapAlgebra r)
m Vector (WrapAlgebra r)
p of
Maybe (Vector (WrapAlgebra r))
Nothing -> {-# SCC "recur" #-} Matrix (WrapAlgebra r)
-> [Vector (WrapAlgebra r)] -> OrderedPolynomial r ord n
step ({-# SCC "consCol" #-} Matrix (WrapAlgebra r)
m Matrix (WrapAlgebra r)
-> Matrix (WrapAlgebra r) -> Matrix (WrapAlgebra r)
forall a. Matrix a -> Matrix a -> Matrix a
M.<|> Vector (WrapAlgebra r) -> Matrix (WrapAlgebra r)
forall a. Vector a -> Matrix a
M.colVector Vector (WrapAlgebra r)
p) [Vector (WrapAlgebra r)]
ps
Just Vector (WrapAlgebra r)
ans ->
let cur :: Natural
cur = Int -> Natural
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> Natural) -> Int -> Natural
forall a b. (a -> b) -> a -> b
$ Vector (WrapAlgebra r) -> Int
forall a. Vector a -> Int
V.length Vector (WrapAlgebra r)
ans :: Natural
in {-# SCC "buildRelation" #-}
OrderedPolynomial r ord n -> Natural -> OrderedPolynomial r ord n
forall r. Unital r => r -> Natural -> r
pow OrderedPolynomial r ord n
x Natural
cur
OrderedPolynomial r ord n
-> OrderedPolynomial r ord n -> OrderedPolynomial r ord n
forall r. Group r => r -> r -> r
- [OrderedPolynomial r ord n] -> OrderedPolynomial r ord n
forall (f :: * -> *) m. (Foldable f, Monoidal m) => f m -> m
sum
( (r -> OrderedPolynomial r ord n -> OrderedPolynomial r ord n)
-> [r]
-> [OrderedPolynomial r ord n]
-> [OrderedPolynomial r ord n]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith
r -> OrderedPolynomial r ord n -> OrderedPolynomial r ord n
forall r m. Module (Scalar r) m => r -> m -> m
(.*.)
((WrapAlgebra r -> r) -> [WrapAlgebra r] -> [r]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap WrapAlgebra r -> r
forall a. WrapAlgebra a -> a
unwrapAlgebra ([WrapAlgebra r] -> [r]) -> [WrapAlgebra r] -> [r]
forall a b. (a -> b) -> a -> b
$ Vector (WrapAlgebra r) -> [WrapAlgebra r]
forall a. Vector a -> [a]
V.toList Vector (WrapAlgebra r)
ans)
[OrderedPolynomial r ord n -> Natural -> OrderedPolynomial r ord n
forall r. Unital r => r -> Natural -> r
pow OrderedPolynomial r ord n
x Natural
i | Natural
i <- [Natural
0 :: Natural .. Natural
cur Natural -> Natural -> Natural
forall a. Num a => a -> a -> a
P.- Natural
1]]
)
in Matrix (WrapAlgebra r)
-> [Vector (WrapAlgebra r)] -> OrderedPolynomial r ord n
step (Vector (WrapAlgebra r) -> Matrix (WrapAlgebra r)
forall a. Vector a -> Matrix a
M.colVector Vector (WrapAlgebra r)
p0) [Vector (WrapAlgebra r)]
pows
solveLinear ::
(Ord r, P.Fractional r) =>
M.Matrix r ->
V.Vector r ->
Maybe (V.Vector r)
solveLinear :: Matrix r -> Vector r -> Maybe (Vector r)
solveLinear Matrix r
mat Vector r
vec =
{-# SCC "solveLinear" #-}
if ({-# SCC "uRank" #-} Matrix r -> Int
forall a. (Eq a, Num a) => Matrix a -> Int
uRank Matrix r
u) Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Matrix r -> Int
forall a. (Eq a, Num a) => Matrix a -> Int
uRank Matrix r
u' Bool -> Bool -> Bool
|| Matrix r -> r
forall a. Num a => Matrix a -> a
M.diagProd Matrix r
u r -> r -> Bool
forall a. Eq a => a -> a -> Bool
== r
0 Bool -> Bool -> Bool
|| Matrix r -> Int
forall a. (Eq a, Num a) => Matrix a -> Int
uRank Matrix r
u Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Matrix r -> Int
forall a. Matrix a -> Int
M.ncols Matrix r
mat
then Maybe (Vector r)
forall a. Maybe a
Nothing
else
let ans :: Vector r
ans = Int -> Matrix r -> Vector r
forall a. Int -> Matrix a -> Vector a
M.getCol Int
1 (Matrix r -> Vector r) -> Matrix r -> Vector r
forall a b. (a -> b) -> a -> b
$ Matrix r
p Matrix r -> Matrix r -> Matrix r
forall a. Num a => a -> a -> a
P.* Vector r -> Matrix r
forall a. Vector a -> Matrix a
M.colVector Vector r
vec
lsol :: Vector r
lsol = {-# SCC "solveL" #-} Vector r -> Vector r
solveL Vector r
ans
cfs :: Vector r
cfs = Int -> Matrix r -> Vector r
forall a. Int -> Matrix a -> Vector a
M.getCol Int
1 (Matrix r -> Vector r) -> Matrix r -> Vector r
forall a b. (a -> b) -> a -> b
$ Matrix r
q Matrix r -> Matrix r -> Matrix r
forall a. Num a => a -> a -> a
P.* Vector r -> Matrix r
forall a. Vector a -> Matrix a
M.colVector ({-# SCC "solveU" #-} Vector r -> Vector r
solveU Vector r
lsol)
in Vector r -> Maybe (Vector r)
forall a. a -> Maybe a
Just Vector r
cfs
where
Just (Matrix r
u, Matrix r
l, Matrix r
p, Matrix r
q, r
_, r
_) = Matrix r -> Maybe (Matrix r, Matrix r, Matrix r, Matrix r, r, r)
forall a.
(Ord a, Fractional a) =>
Matrix a -> Maybe (Matrix a, Matrix a, Matrix a, Matrix a, a, a)
M.luDecomp' Matrix r
mat
Just (Matrix r
u', Matrix r
_, Matrix r
_, Matrix r
_, r
_, r
_) = Matrix r -> Maybe (Matrix r, Matrix r, Matrix r, Matrix r, r, r)
forall a.
(Ord a, Fractional a) =>
Matrix a -> Maybe (Matrix a, Matrix a, Matrix a, Matrix a, a, a)
M.luDecomp' (Matrix r
mat Matrix r -> Matrix r -> Matrix r
forall a. Matrix a -> Matrix a -> Matrix a
M.<|> Vector r -> Matrix r
forall a. Vector a -> Matrix a
M.colVector Vector r
vec)
uRank :: Matrix a -> Int
uRank = (a -> Int -> Int) -> Int -> Vector a -> Int
forall a b. (a -> b -> b) -> b -> Vector a -> b
V.foldr (\a
a Int
acc -> if a
a a -> a -> Bool
forall a. Eq a => a -> a -> Bool
/= a
0 then Int
acc Int -> Int -> Int
forall r. Additive r => r -> r -> r
+ Int
1 else Int
acc) (Int
0 :: Int) (Vector a -> Int) -> (Matrix a -> Vector a) -> Matrix a -> Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Matrix a -> Vector a
forall a. Matrix a -> Vector a
M.getDiag
solveL :: Vector r -> Vector r
solveL Vector r
v = (forall s. ST s (MVector s r)) -> Vector r
forall a. (forall s. ST s (MVector s a)) -> Vector a
V.create ((forall s. ST s (MVector s r)) -> Vector r)
-> (forall s. ST s (MVector s r)) -> Vector r
forall a b. (a -> b) -> a -> b
$ do
let stop :: Int
stop = Int -> Int -> Int
forall a. Ord a => a -> a -> a
min (Matrix r -> Int
forall a. Matrix a -> Int
M.ncols Matrix r
l) (Matrix r -> Int
forall a. Matrix a -> Int
M.nrows Matrix r
l)
MVector s r
mv <- Int -> r -> ST s (MVector (PrimState (ST s)) r)
forall (m :: * -> *) a.
PrimMonad m =>
Int -> a -> m (MVector (PrimState m) a)
MV.replicate (Matrix r -> Int
forall a. Matrix a -> Int
M.ncols Matrix r
l) r
0
[Int] -> (Int -> ST s ()) -> ST s ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [Int
0 .. Int
stop Int -> Int -> Int
forall r. Group r => r -> r -> r
- Int
1] ((Int -> ST s ()) -> ST s ()) -> (Int -> ST s ()) -> ST s ()
forall a b. (a -> b) -> a -> b
$ \Int
i -> do
MVector (PrimState (ST s)) r -> Int -> r -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
MV.write MVector s r
MVector (PrimState (ST s)) r
mv Int
i (r -> ST s ()) -> r -> ST s ()
forall a b. (a -> b) -> a -> b
$ Vector r
v Vector r -> Int -> r
forall a. Vector a -> Int -> a
V.! Int
i
[Int] -> (Int -> ST s ()) -> ST s ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [Int
0, Int
1 .. Int -> Int -> Int
forall a. Ord a => a -> a -> a
min (Int
i Int -> Int -> Int
forall r. Group r => r -> r -> r
-Int
1) (Matrix r -> Int
forall a. Matrix a -> Int
M.ncols Matrix r
l Int -> Int -> Int
forall r. Group r => r -> r -> r
- Int
1)] ((Int -> ST s ()) -> ST s ()) -> (Int -> ST s ()) -> ST s ()
forall a b. (a -> b) -> a -> b
$ \Int
j -> do
r
a <- MVector (PrimState (ST s)) r -> Int -> ST s r
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> m a
MV.read MVector s r
MVector (PrimState (ST s)) r
mv Int
i
r
b <- MVector (PrimState (ST s)) r -> Int -> ST s r
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> m a
MV.read MVector s r
MVector (PrimState (ST s)) r
mv Int
j
MVector (PrimState (ST s)) r -> Int -> r -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
MV.write MVector s r
MVector (PrimState (ST s)) r
mv Int
i (r -> ST s ()) -> r -> ST s ()
forall a b. (a -> b) -> a -> b
$ r
a r -> r -> r
forall a. Num a => a -> a -> a
P.- (Matrix r
l Matrix r -> (Int, Int) -> r
forall a. Matrix a -> (Int, Int) -> a
M.! (Int
i Int -> Int -> Int
forall r. Additive r => r -> r -> r
+ Int
1, Int
j Int -> Int -> Int
forall r. Additive r => r -> r -> r
+ Int
1)) r -> r -> r
forall a. Num a => a -> a -> a
P.* r
b
MVector s r -> ST s (MVector s r)
forall (m :: * -> *) a. Monad m => a -> m a
return MVector s r
mv
solveU :: Vector r -> Vector r
solveU Vector r
v = (forall s. ST s (MVector s r)) -> Vector r
forall a. (forall s. ST s (MVector s a)) -> Vector a
V.create ((forall s. ST s (MVector s r)) -> Vector r)
-> (forall s. ST s (MVector s r)) -> Vector r
forall a b. (a -> b) -> a -> b
$ do
let stop :: Int
stop = Int -> Int -> Int
forall a. Ord a => a -> a -> a
min (Matrix r -> Int
forall a. Matrix a -> Int
M.ncols Matrix r
u) (Matrix r -> Int
forall a. Matrix a -> Int
M.nrows Matrix r
u)
MVector s r
mv <- Int -> r -> ST s (MVector (PrimState (ST s)) r)
forall (m :: * -> *) a.
PrimMonad m =>
Int -> a -> m (MVector (PrimState m) a)
MV.replicate (Matrix r -> Int
forall a. Matrix a -> Int
M.ncols Matrix r
u) r
0
[Int] -> (Int -> ST s ()) -> ST s ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [Int
stop Int -> Int -> Int
forall r. Group r => r -> r -> r
- Int
1, Int
stop Int -> Int -> Int
forall r. Group r => r -> r -> r
- Int
2 .. Int
0] ((Int -> ST s ()) -> ST s ()) -> (Int -> ST s ()) -> ST s ()
forall a b. (a -> b) -> a -> b
$ \Int
i -> do
MVector (PrimState (ST s)) r -> Int -> r -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
MV.write MVector s r
MVector (PrimState (ST s)) r
mv Int
i (r -> ST s ()) -> r -> ST s ()
forall a b. (a -> b) -> a -> b
$ Vector r
v Vector r -> Int -> r
forall a. Vector a -> Int -> a
V.! Int
i
[Int] -> (Int -> ST s ()) -> ST s ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [Int
i Int -> Int -> Int
forall r. Additive r => r -> r -> r
+ Int
1, Int
i Int -> Int -> Int
forall r. Additive r => r -> r -> r
+ Int
2 .. Matrix r -> Int
forall a. Matrix a -> Int
M.ncols Matrix r
u Int -> Int -> Int
forall r. Group r => r -> r -> r
-Int
1] ((Int -> ST s ()) -> ST s ()) -> (Int -> ST s ()) -> ST s ()
forall a b. (a -> b) -> a -> b
$ \Int
j -> do
r
a <- MVector (PrimState (ST s)) r -> Int -> ST s r
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> m a
MV.read MVector s r
MVector (PrimState (ST s)) r
mv Int
i
r
b <- MVector (PrimState (ST s)) r -> Int -> ST s r
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> m a
MV.read MVector s r
MVector (PrimState (ST s)) r
mv Int
j
MVector (PrimState (ST s)) r -> Int -> r -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
MV.write MVector s r
MVector (PrimState (ST s)) r
mv Int
i (r -> ST s ()) -> r -> ST s ()
forall a b. (a -> b) -> a -> b
$ r
a r -> r -> r
forall a. Num a => a -> a -> a
P.- (Matrix r
u Matrix r -> (Int, Int) -> r
forall a. Matrix a -> (Int, Int) -> a
M.! (Int
i Int -> Int -> Int
forall r. Additive r => r -> r -> r
+ Int
1, Int
j Int -> Int -> Int
forall r. Additive r => r -> r -> r
+ Int
1)) r -> r -> r
forall a. Num a => a -> a -> a
P.* r
b
r
a0 <- MVector (PrimState (ST s)) r -> Int -> ST s r
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> m a
MV.read MVector s r
MVector (PrimState (ST s)) r
mv Int
i
MVector (PrimState (ST s)) r -> Int -> r -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
MV.write MVector s r
MVector (PrimState (ST s)) r
mv Int
i (r -> ST s ()) -> r -> ST s ()
forall a b. (a -> b) -> a -> b
$ r
a0 r -> r -> r
forall a. Fractional a => a -> a -> a
P./ (Matrix r
u Matrix r -> (Int, Int) -> r
forall a. Matrix a -> (Int, Int) -> a
M.! (Int
i Int -> Int -> Int
forall r. Additive r => r -> r -> r
+ Int
1, Int
i Int -> Int -> Int
forall r. Additive r => r -> r -> r
+ Int
1))
MVector s r -> ST s (MVector s r)
forall (m :: * -> *) a. Monad m => a -> m a
return MVector s r
mv
radical ::
forall r ord n.
( Ord r
, CoeffRing r
, KnownNat n
, Field r
, IsMonomialOrder n ord
) =>
Ideal (OrderedPolynomial r ord n) ->
Ideal (OrderedPolynomial r ord n)
radical :: Ideal (OrderedPolynomial r ord n)
-> Ideal (OrderedPolynomial r ord n)
radical Ideal (OrderedPolynomial r ord n)
ideal =
{-# SCC "radical" #-}
let gens :: [OrderedPolynomial r ord n]
gens = {-# SCC "calcGens" #-} (Ordinal n -> OrderedPolynomial r ord n)
-> [Ordinal n] -> [OrderedPolynomial r ord n]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
map (\Ordinal n
on -> Ordinal n -> OrderedPolynomial r ord n -> OrderedPolynomial r ord n
forall r (n :: Nat) ord.
(CoeffRing r, KnownNat n, IsMonomialOrder n ord, Field r) =>
Ordinal n -> OrderedPolynomial r ord n -> OrderedPolynomial r ord n
reduction Ordinal n
on (OrderedPolynomial r ord n -> OrderedPolynomial r ord n)
-> OrderedPolynomial r ord n -> OrderedPolynomial r ord n
forall a b. (a -> b) -> a -> b
$ Ordinal n
-> Ideal (OrderedPolynomial r ord n) -> OrderedPolynomial r ord n
forall r ord (n :: Nat).
(Ord r, Field r, CoeffRing r, KnownNat n, IsMonomialOrder n ord) =>
Ordinal n
-> Ideal (OrderedPolynomial r ord n) -> OrderedPolynomial r ord n
univPoly Ordinal n
on Ideal (OrderedPolynomial r ord n)
ideal) ([Ordinal n] -> [OrderedPolynomial r ord n])
-> [Ordinal n] -> [OrderedPolynomial r ord n]
forall a b. (a -> b) -> a -> b
$ SNat n -> [Ordinal n]
forall (n :: Nat). SNat n -> [Ordinal n]
enumOrdinal (SNat n
forall (n :: Nat). KnownNat n => SNat n
sNat :: SNat n)
in [OrderedPolynomial r ord n] -> Ideal (OrderedPolynomial r ord n)
forall r. (DecidableZero r, Monoidal r) => [r] -> Ideal r
toIdeal ([OrderedPolynomial r ord n] -> Ideal (OrderedPolynomial r ord n))
-> [OrderedPolynomial r ord n] -> Ideal (OrderedPolynomial r ord n)
forall a b. (a -> b) -> a -> b
$ Ideal (OrderedPolynomial r ord n) -> [OrderedPolynomial r ord n]
forall poly.
(Field (Coefficient poly), IsOrderedPolynomial poly) =>
Ideal poly -> [poly]
calcGroebnerBasis (Ideal (OrderedPolynomial r ord n) -> [OrderedPolynomial r ord n])
-> Ideal (OrderedPolynomial r ord n) -> [OrderedPolynomial r ord n]
forall a b. (a -> b) -> a -> b
$ [OrderedPolynomial r ord n] -> Ideal (OrderedPolynomial r ord n)
forall r. (DecidableZero r, Monoidal r) => [r] -> Ideal r
toIdeal ([OrderedPolynomial r ord n] -> Ideal (OrderedPolynomial r ord n))
-> [OrderedPolynomial r ord n] -> Ideal (OrderedPolynomial r ord n)
forall a b. (a -> b) -> a -> b
$ Ideal (OrderedPolynomial r ord n) -> [OrderedPolynomial r ord n]
forall r. Ideal r -> [r]
generators Ideal (OrderedPolynomial r ord n)
ideal [OrderedPolynomial r ord n]
-> [OrderedPolynomial r ord n] -> [OrderedPolynomial r ord n]
forall w. Monoid w => w -> w -> w
++ [OrderedPolynomial r ord n]
gens
isRadical ::
forall r ord n.
( Ord r
, CoeffRing r
, KnownNat n
, 0 < n
, Field r
, IsMonomialOrder n ord
) =>
Ideal (OrderedPolynomial r ord n) ->
Bool
isRadical :: Ideal (OrderedPolynomial r ord n) -> Bool
isRadical Ideal (OrderedPolynomial r ord n)
ideal =
let gens :: [OrderedPolynomial r ord n]
gens =
(Ordinal n -> OrderedPolynomial r ord n)
-> [Ordinal n] -> [OrderedPolynomial r ord n]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
map (\Ordinal n
on -> Ordinal n -> OrderedPolynomial r ord n -> OrderedPolynomial r ord n
forall r (n :: Nat) ord.
(CoeffRing r, KnownNat n, IsMonomialOrder n ord, Field r) =>
Ordinal n -> OrderedPolynomial r ord n -> OrderedPolynomial r ord n
reduction Ordinal n
on (OrderedPolynomial r ord n -> OrderedPolynomial r ord n)
-> OrderedPolynomial r ord n -> OrderedPolynomial r ord n
forall a b. (a -> b) -> a -> b
$ Ordinal n
-> Ideal (OrderedPolynomial r ord n) -> OrderedPolynomial r ord n
forall r ord (n :: Nat).
(Ord r, Field r, CoeffRing r, KnownNat n, IsMonomialOrder n ord) =>
Ordinal n
-> Ideal (OrderedPolynomial r ord n) -> OrderedPolynomial r ord n
univPoly Ordinal n
on Ideal (OrderedPolynomial r ord n)
ideal) ([Ordinal n] -> [OrderedPolynomial r ord n])
-> [Ordinal n] -> [OrderedPolynomial r ord n]
forall a b. (a -> b) -> a -> b
$
SNat n -> [Ordinal n]
forall (n :: Nat). SNat n -> [Ordinal n]
enumOrdinal (SNat n
forall (n :: Nat). KnownNat n => SNat n
sNat :: SNat n)
in (OrderedPolynomial r ord n -> Bool)
-> [OrderedPolynomial r ord n] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (OrderedPolynomial r ord n
-> Ideal (OrderedPolynomial r ord n) -> Bool
forall poly.
(Field (Coefficient poly), IsOrderedPolynomial poly) =>
poly -> Ideal poly -> Bool
`isIdealMember` Ideal (OrderedPolynomial r ord n)
ideal) [OrderedPolynomial r ord n]
gens
fglm ::
( Ord r
, KnownNat n
, Field r
, IsMonomialOrder n ord
, 0 < n
) =>
Ideal (OrderedPolynomial r ord n) ->
([OrderedPolynomial r Lex n], [OrderedPolynomial r Lex n])
fglm :: Ideal (OrderedPolynomial r ord n)
-> ([OrderedPolynomial r Lex n], [OrderedPolynomial r Lex n])
fglm Ideal (OrderedPolynomial r ord n)
ideal = Ideal (OrderedPolynomial r ord n)
-> (forall ideal.
Reifies ideal (QIdeal (OrderedPolynomial r ord n)) =>
Proxy ideal
-> ([OrderedPolynomial r Lex n], [OrderedPolynomial r Lex n]))
-> ([OrderedPolynomial r Lex n], [OrderedPolynomial r Lex n])
forall poly a.
(IsOrderedPolynomial poly, Field (Coefficient poly)) =>
Ideal poly
-> (forall ideal. Reifies ideal (QIdeal poly) => Proxy ideal -> a)
-> a
reifyQuotient Ideal (OrderedPolynomial r ord n)
ideal ((forall ideal.
Reifies ideal (QIdeal (OrderedPolynomial r ord n)) =>
Proxy ideal
-> ([OrderedPolynomial r Lex n], [OrderedPolynomial r Lex n]))
-> ([OrderedPolynomial r Lex n], [OrderedPolynomial r Lex n]))
-> (forall ideal.
Reifies ideal (QIdeal (OrderedPolynomial r ord n)) =>
Proxy ideal
-> ([OrderedPolynomial r Lex n], [OrderedPolynomial r Lex n]))
-> ([OrderedPolynomial r Lex n], [OrderedPolynomial r Lex n])
forall a b. (a -> b) -> a -> b
$ \Proxy ideal
pxy ->
(OrderedPolynomial r ord n -> Vector r)
-> ([OrderedPolynomial r Lex n], [OrderedPolynomial r Lex n])
forall k ord (n :: Nat).
(Ord k, Field k, 0 < n, IsMonomialOrder n ord, CoeffRing k,
KnownNat n) =>
(OrderedPolynomial k ord n -> Vector k)
-> ([OrderedPolynomial k Lex n], [OrderedPolynomial k Lex n])
fglmMap (Quotient (OrderedPolynomial r ord n) ideal -> Vector r
forall k poly (ideal :: k).
(IsOrderedPolynomial poly, Reifies ideal (QIdeal poly)) =>
Quotient poly ideal -> Vector (Coefficient poly)
vectorRep (Quotient (OrderedPolynomial r ord n) ideal -> Vector r)
-> (OrderedPolynomial r ord n
-> Quotient (OrderedPolynomial r ord n) ideal)
-> OrderedPolynomial r ord n
-> Vector r
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Proxy ideal
-> OrderedPolynomial r ord n
-> Quotient (OrderedPolynomial r ord n) ideal
forall k poly (ideal :: k).
(IsOrderedPolynomial poly, Field (Coefficient poly),
Reifies ideal (QIdeal poly)) =>
Proxy ideal -> poly -> Quotient poly ideal
modIdeal' Proxy ideal
pxy)
fglmMap ::
forall k ord n.
( Ord k
, Field k
, 0 < n
, IsMonomialOrder n ord
, CoeffRing k
, KnownNat n
) =>
(OrderedPolynomial k ord n -> V.Vector k) ->
( [OrderedPolynomial k Lex n]
, [OrderedPolynomial k Lex n]
)
fglmMap :: (OrderedPolynomial k ord n -> Vector k)
-> ([OrderedPolynomial k Lex n], [OrderedPolynomial k Lex n])
fglmMap OrderedPolynomial k ord n -> Vector k
l = (forall s.
ST s ([OrderedPolynomial k Lex n], [OrderedPolynomial k Lex n]))
-> ([OrderedPolynomial k Lex n], [OrderedPolynomial k Lex n])
forall a. (forall s. ST s a) -> a
runST ((forall s.
ST s ([OrderedPolynomial k Lex n], [OrderedPolynomial k Lex n]))
-> ([OrderedPolynomial k Lex n], [OrderedPolynomial k Lex n]))
-> (forall s.
ST s ([OrderedPolynomial k Lex n], [OrderedPolynomial k Lex n]))
-> ([OrderedPolynomial k Lex n], [OrderedPolynomial k Lex n])
forall a b. (a -> b) -> a -> b
$ do
FGLMEnv s k ord n
env <- (OrderedPolynomial k ord n -> Vector k)
-> STRef s [OrderedPolynomial k Lex n]
-> STRef s [OrderedPolynomial k ord n]
-> STRef s (Maybe (OrderedPolynomial k Lex n))
-> STRef s (OrderedMonomial Lex n)
-> FGLMEnv s k ord n
forall k s r (ord :: k) (n :: Nat).
(OrderedPolynomial r ord n -> Vector r)
-> STRef s [OrderedPolynomial r Lex n]
-> STRef s [OrderedPolynomial r ord n]
-> STRef s (Maybe (OrderedPolynomial r Lex n))
-> STRef s (OrderedMonomial Lex n)
-> FGLMEnv s r ord n
FGLMEnv OrderedPolynomial k ord n -> Vector k
l (STRef s [OrderedPolynomial k Lex n]
-> STRef s [OrderedPolynomial k ord n]
-> STRef s (Maybe (OrderedPolynomial k Lex n))
-> STRef s (OrderedMonomial Lex n)
-> FGLMEnv s k ord n)
-> ST s (STRef s [OrderedPolynomial k Lex n])
-> ST
s
(STRef s [OrderedPolynomial k ord n]
-> STRef s (Maybe (OrderedPolynomial k Lex n))
-> STRef s (OrderedMonomial Lex n)
-> FGLMEnv s k ord n)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> [OrderedPolynomial k Lex n]
-> ST s (STRef s [OrderedPolynomial k Lex n])
forall a s. a -> ST s (STRef s a)
newSTRef [] ST
s
(STRef s [OrderedPolynomial k ord n]
-> STRef s (Maybe (OrderedPolynomial k Lex n))
-> STRef s (OrderedMonomial Lex n)
-> FGLMEnv s k ord n)
-> ST s (STRef s [OrderedPolynomial k ord n])
-> ST
s
(STRef s (Maybe (OrderedPolynomial k Lex n))
-> STRef s (OrderedMonomial Lex n) -> FGLMEnv s k ord n)
forall (f :: * -> *) a b. Applicative f => f (a -> b) -> f a -> f b
<*> [OrderedPolynomial k ord n]
-> ST s (STRef s [OrderedPolynomial k ord n])
forall a s. a -> ST s (STRef s a)
newSTRef [] ST
s
(STRef s (Maybe (OrderedPolynomial k Lex n))
-> STRef s (OrderedMonomial Lex n) -> FGLMEnv s k ord n)
-> ST s (STRef s (Maybe (OrderedPolynomial k Lex n)))
-> ST s (STRef s (OrderedMonomial Lex n) -> FGLMEnv s k ord n)
forall (f :: * -> *) a b. Applicative f => f (a -> b) -> f a -> f b
<*> Maybe (OrderedPolynomial k Lex n)
-> ST s (STRef s (Maybe (OrderedPolynomial k Lex n)))
forall a s. a -> ST s (STRef s a)
newSTRef Maybe (OrderedPolynomial k Lex n)
forall a. Maybe a
Nothing ST s (STRef s (OrderedMonomial Lex n) -> FGLMEnv s k ord n)
-> ST s (STRef s (OrderedMonomial Lex n))
-> ST s (FGLMEnv s k ord n)
forall (f :: * -> *) a b. Applicative f => f (a -> b) -> f a -> f b
<*> OrderedMonomial Lex n -> ST s (STRef s (OrderedMonomial Lex n))
forall a s. a -> ST s (STRef s a)
newSTRef OrderedMonomial Lex n
forall r. Unital r => r
one
(ReaderT
(FGLMEnv s k ord n)
(ST s)
([OrderedPolynomial k Lex n], [OrderedPolynomial k Lex n])
-> FGLMEnv s k ord n
-> ST s ([OrderedPolynomial k Lex n], [OrderedPolynomial k Lex n]))
-> FGLMEnv s k ord n
-> ReaderT
(FGLMEnv s k ord n)
(ST s)
([OrderedPolynomial k Lex n], [OrderedPolynomial k Lex n])
-> ST s ([OrderedPolynomial k Lex n], [OrderedPolynomial k Lex n])
forall a b c. (a -> b -> c) -> b -> a -> c
flip ReaderT
(FGLMEnv s k ord n)
(ST s)
([OrderedPolynomial k Lex n], [OrderedPolynomial k Lex n])
-> FGLMEnv s k ord n
-> ST s ([OrderedPolynomial k Lex n], [OrderedPolynomial k Lex n])
forall r (m :: * -> *) a. ReaderT r m a -> r -> m a
runReaderT FGLMEnv s k ord n
env (ReaderT
(FGLMEnv s k ord n)
(ST s)
([OrderedPolynomial k Lex n], [OrderedPolynomial k Lex n])
-> ST s ([OrderedPolynomial k Lex n], [OrderedPolynomial k Lex n]))
-> ReaderT
(FGLMEnv s k ord n)
(ST s)
([OrderedPolynomial k Lex n], [OrderedPolynomial k Lex n])
-> ST s ([OrderedPolynomial k Lex n], [OrderedPolynomial k Lex n])
forall a b. (a -> b) -> a -> b
$ do
Machine s k ord n ()
forall r (n :: Nat) o s.
(DecidableZero r, Ord r, KnownNat n, Field r,
IsMonomialOrder n o) =>
Machine s r o n ()
mainLoop
ReaderT (FGLMEnv s k ord n) (ST s) Bool
-> Machine s k ord n () -> Machine s k ord n ()
forall (m :: * -> *) a. Monad m => m Bool -> m a -> m ()
whileM_ ReaderT (FGLMEnv s k ord n) (ST s) Bool
forall k s r (o :: k) (n :: Nat).
(0 < n, Ord r, KnownNat n, Field r) =>
Machine s r o n Bool
toContinue (Machine s k ord n () -> Machine s k ord n ())
-> Machine s k ord n () -> Machine s k ord n ()
forall a b. (a -> b) -> a -> b
$ Machine s k ord n ()
forall k s r (ord :: k) (n :: Nat).
(CoeffRing r, KnownNat n) =>
Machine s r ord n ()
nextMonomial Machine s k ord n ()
-> Machine s k ord n () -> Machine s k ord n ()
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> Machine s k ord n ()
forall r (n :: Nat) o s.
(DecidableZero r, Ord r, KnownNat n, Field r,
IsMonomialOrder n o) =>
Machine s r o n ()
mainLoop
(,) ([OrderedPolynomial k Lex n]
-> [OrderedPolynomial k Lex n]
-> ([OrderedPolynomial k Lex n], [OrderedPolynomial k Lex n]))
-> ReaderT (FGLMEnv s k ord n) (ST s) [OrderedPolynomial k Lex n]
-> ReaderT
(FGLMEnv s k ord n)
(ST s)
([OrderedPolynomial k Lex n]
-> ([OrderedPolynomial k Lex n], [OrderedPolynomial k Lex n]))
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Getting
(STRef s [OrderedPolynomial k Lex n])
(FGLMEnv s k ord n)
(STRef s [OrderedPolynomial k Lex n])
-> ReaderT (FGLMEnv s k ord n) (ST s) [OrderedPolynomial k Lex n]
forall k s b r (ord :: k) (n :: Nat).
Getting (STRef s b) (FGLMEnv s r ord n) (STRef s b)
-> Machine s r ord n b
look Getting
(STRef s [OrderedPolynomial k Lex n])
(FGLMEnv s k ord n)
(STRef s [OrderedPolynomial k Lex n])
forall k s r (ord :: k) (n :: Nat).
Lens' (FGLMEnv s r ord n) (STRef s [OrderedPolynomial r Lex n])
gLex ReaderT
(FGLMEnv s k ord n)
(ST s)
([OrderedPolynomial k Lex n]
-> ([OrderedPolynomial k Lex n], [OrderedPolynomial k Lex n]))
-> ReaderT (FGLMEnv s k ord n) (ST s) [OrderedPolynomial k Lex n]
-> ReaderT
(FGLMEnv s k ord n)
(ST s)
([OrderedPolynomial k Lex n], [OrderedPolynomial k Lex n])
forall (f :: * -> *) a b. Applicative f => f (a -> b) -> f a -> f b
<*> ((OrderedPolynomial k ord n -> OrderedPolynomial k Lex n)
-> [OrderedPolynomial k ord n] -> [OrderedPolynomial k Lex n]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
map (Lex -> OrderedPolynomial k ord n -> OrderedPolynomial k Lex n
forall k (n :: Nat) o o'.
(CoeffRing k, Eq (Monomial n), IsMonomialOrder n o,
IsMonomialOrder n o', KnownNat n) =>
o' -> OrderedPolynomial k o n -> OrderedPolynomial k o' n
changeOrder Lex
Lex) ([OrderedPolynomial k ord n] -> [OrderedPolynomial k Lex n])
-> ReaderT (FGLMEnv s k ord n) (ST s) [OrderedPolynomial k ord n]
-> ReaderT (FGLMEnv s k ord n) (ST s) [OrderedPolynomial k Lex n]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Getting
(STRef s [OrderedPolynomial k ord n])
(FGLMEnv s k ord n)
(STRef s [OrderedPolynomial k ord n])
-> ReaderT (FGLMEnv s k ord n) (ST s) [OrderedPolynomial k ord n]
forall k s b r (ord :: k) (n :: Nat).
Getting (STRef s b) (FGLMEnv s r ord n) (STRef s b)
-> Machine s r ord n b
look Getting
(STRef s [OrderedPolynomial k ord n])
(FGLMEnv s k ord n)
(STRef s [OrderedPolynomial k ord n])
forall k s r (ord :: k) (n :: Nat).
Lens' (FGLMEnv s r ord n) (STRef s [OrderedPolynomial r ord n])
bLex)
mainLoop ::
(DecidableZero r, Ord r, KnownNat n, Field r, IsMonomialOrder n o) =>
Machine s r o n ()
mainLoop :: Machine s r o n ()
mainLoop = do
OrderedMonomial Lex n
m <- Getting
(STRef s (OrderedMonomial Lex n))
(FGLMEnv s r o n)
(STRef s (OrderedMonomial Lex n))
-> Machine s r o n (OrderedMonomial Lex n)
forall k s b r (ord :: k) (n :: Nat).
Getting (STRef s b) (FGLMEnv s r ord n) (STRef s b)
-> Machine s r ord n b
look Getting
(STRef s (OrderedMonomial Lex n))
(FGLMEnv s r o n)
(STRef s (OrderedMonomial Lex n))
forall k s r (ord :: k) (n :: Nat).
Lens' (FGLMEnv s r ord n) (STRef s (OrderedMonomial Lex n))
monomial
let f :: OrderedPolynomial r o n
f = (Coefficient (OrderedPolynomial r o n),
OrderedMonomial
(MOrder (OrderedPolynomial r o n))
(Arity (OrderedPolynomial r o n)))
-> OrderedPolynomial r o n
forall poly.
IsOrderedPolynomial poly =>
(Coefficient poly, OrderedMonomial (MOrder poly) (Arity poly))
-> poly
toPolynomial (Coefficient (OrderedPolynomial r o n)
forall r. Unital r => r
one, Proxy o -> OrderedMonomial Lex n -> OrderedMonomial o n
forall k1 k2 (o' :: k1) (ord :: k2) (n :: Nat).
Proxy o' -> OrderedMonomial ord n -> OrderedMonomial o' n
changeMonomialOrderProxy Proxy o
forall k (t :: k). Proxy t
Proxy OrderedMonomial Lex n
m)
Vector r
lx <- OrderedPolynomial r o n
-> ReaderT (FGLMEnv s r o n) (ST s) (Vector r)
forall k s r (ord :: k) (n :: Nat) (f :: * -> *).
MonadReader (FGLMEnv s r ord n) f =>
OrderedPolynomial r ord n -> f (Vector r)
image OrderedPolynomial r o n
f
[Vector r]
bs <- (OrderedPolynomial r o n
-> ReaderT (FGLMEnv s r o n) (ST s) (Vector r))
-> [OrderedPolynomial r o n]
-> ReaderT (FGLMEnv s r o n) (ST s) [Vector r]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
mapM OrderedPolynomial r o n
-> ReaderT (FGLMEnv s r o n) (ST s) (Vector r)
forall k s r (ord :: k) (n :: Nat) (f :: * -> *).
MonadReader (FGLMEnv s r ord n) f =>
OrderedPolynomial r ord n -> f (Vector r)
image ([OrderedPolynomial r o n]
-> ReaderT (FGLMEnv s r o n) (ST s) [Vector r])
-> ReaderT (FGLMEnv s r o n) (ST s) [OrderedPolynomial r o n]
-> ReaderT (FGLMEnv s r o n) (ST s) [Vector r]
forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
=<< Getting
(STRef s [OrderedPolynomial r o n])
(FGLMEnv s r o n)
(STRef s [OrderedPolynomial r o n])
-> ReaderT (FGLMEnv s r o n) (ST s) [OrderedPolynomial r o n]
forall k s b r (ord :: k) (n :: Nat).
Getting (STRef s b) (FGLMEnv s r ord n) (STRef s b)
-> Machine s r ord n b
look Getting
(STRef s [OrderedPolynomial r o n])
(FGLMEnv s r o n)
(STRef s [OrderedPolynomial r o n])
forall k s r (ord :: k) (n :: Nat).
Lens' (FGLMEnv s r ord n) (STRef s [OrderedPolynomial r ord n])
bLex
let mat :: Matrix (WrapAlgebra r)
mat = (Matrix (WrapAlgebra r)
-> Matrix (WrapAlgebra r) -> Matrix (WrapAlgebra r))
-> Matrix (WrapAlgebra r)
-> [Matrix (WrapAlgebra r)]
-> Matrix (WrapAlgebra r)
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr Matrix (WrapAlgebra r)
-> Matrix (WrapAlgebra r) -> Matrix (WrapAlgebra r)
forall a. Matrix a -> Matrix a -> Matrix a
(M.<|>) (Int -> Int -> [WrapAlgebra r] -> Matrix (WrapAlgebra r)
forall a. Int -> Int -> [a] -> Matrix a
M.fromList Int
0 Int
0 []) ([Matrix (WrapAlgebra r)] -> Matrix (WrapAlgebra r))
-> [Matrix (WrapAlgebra r)] -> Matrix (WrapAlgebra r)
forall a b. (a -> b) -> a -> b
$ (Vector r -> Matrix (WrapAlgebra r))
-> [Vector r] -> [Matrix (WrapAlgebra r)]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
map (Vector (WrapAlgebra r) -> Matrix (WrapAlgebra r)
forall a. Vector a -> Matrix a
M.colVector (Vector (WrapAlgebra r) -> Matrix (WrapAlgebra r))
-> (Vector r -> Vector (WrapAlgebra r))
-> Vector r
-> Matrix (WrapAlgebra r)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (r -> WrapAlgebra r) -> Vector r -> Vector (WrapAlgebra r)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap r -> WrapAlgebra r
forall a. a -> WrapAlgebra a
WrapAlgebra) [Vector r]
bs
cond :: Maybe (Vector (WrapAlgebra r))
cond
| [Vector r] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Vector r]
bs =
if (r -> Bool) -> Vector r -> Bool
forall a. (a -> Bool) -> Vector a -> Bool
V.all (r -> r -> Bool
forall a. Eq a => a -> a -> Bool
== r
forall m. Monoidal m => m
zero) Vector r
lx
then Vector (WrapAlgebra r) -> Maybe (Vector (WrapAlgebra r))
forall a. a -> Maybe a
Just (Vector (WrapAlgebra r) -> Maybe (Vector (WrapAlgebra r)))
-> Vector (WrapAlgebra r) -> Maybe (Vector (WrapAlgebra r))
forall a b. (a -> b) -> a -> b
$ Int -> WrapAlgebra r -> Vector (WrapAlgebra r)
forall a. Int -> a -> Vector a
V.replicate ([Vector r] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Vector r]
bs) WrapAlgebra r
0
else Maybe (Vector (WrapAlgebra r))
forall a. Maybe a
Nothing
| Bool
otherwise = Matrix (WrapAlgebra r)
-> Vector (WrapAlgebra r) -> Maybe (Vector (WrapAlgebra r))
forall r.
(Ord r, Fractional r) =>
Matrix r -> Vector r -> Maybe (Vector r)
solveLinear Matrix (WrapAlgebra r)
mat ((r -> WrapAlgebra r) -> Vector r -> Vector (WrapAlgebra r)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap r -> WrapAlgebra r
forall a. a -> WrapAlgebra a
WrapAlgebra Vector r
lx)
case Maybe (Vector (WrapAlgebra r))
cond of
Maybe (Vector (WrapAlgebra r))
Nothing -> do
(STRef s (Maybe (OrderedPolynomial r Lex n))
-> Const
(STRef s (Maybe (OrderedPolynomial r Lex n)))
(STRef s (Maybe (OrderedPolynomial r Lex n))))
-> FGLMEnv s r o n
-> Const
(STRef s (Maybe (OrderedPolynomial r Lex n))) (FGLMEnv s r o n)
forall k s r (ord :: k) (n :: Nat).
Lens'
(FGLMEnv s r ord n) (STRef s (Maybe (OrderedPolynomial r Lex n)))
proced ((STRef s (Maybe (OrderedPolynomial r Lex n))
-> Const
(STRef s (Maybe (OrderedPolynomial r Lex n)))
(STRef s (Maybe (OrderedPolynomial r Lex n))))
-> FGLMEnv s r o n
-> Const
(STRef s (Maybe (OrderedPolynomial r Lex n))) (FGLMEnv s r o n))
-> Maybe (OrderedPolynomial r Lex n) -> Machine s r o n ()
forall (t :: (* -> *) -> * -> *) s s1 a.
(MonadTrans t, MonadReader s (t (ST s1))) =>
Getting (STRef s1 a) s (STRef s1 a) -> a -> t (ST s1) ()
.== Maybe (OrderedPolynomial r Lex n)
forall a. Maybe a
Nothing
Getting
(STRef s [OrderedPolynomial r o n])
(FGLMEnv s r o n)
(STRef s [OrderedPolynomial r o n])
forall k s r (ord :: k) (n :: Nat).
Lens' (FGLMEnv s r ord n) (STRef s [OrderedPolynomial r ord n])
bLex Getting
(STRef s [OrderedPolynomial r o n])
(FGLMEnv s r o n)
(STRef s [OrderedPolynomial r o n])
-> ([OrderedPolynomial r o n] -> [OrderedPolynomial r o n])
-> Machine s r o n ()
forall (t :: (* -> *) -> * -> *) s s1 a.
(MonadTrans t, MonadReader s (t (ST s1))) =>
Getting (STRef s1 a) s (STRef s1 a) -> (a -> a) -> t (ST s1) ()
@== (OrderedPolynomial r o n
f OrderedPolynomial r o n
-> [OrderedPolynomial r o n] -> [OrderedPolynomial r o n]
forall a. a -> [a] -> [a]
:)
Just Vector (WrapAlgebra r)
cs -> do
[OrderedPolynomial r o n]
bps <- Getting
(STRef s [OrderedPolynomial r o n])
(FGLMEnv s r o n)
(STRef s [OrderedPolynomial r o n])
-> ReaderT (FGLMEnv s r o n) (ST s) [OrderedPolynomial r o n]
forall k s b r (ord :: k) (n :: Nat).
Getting (STRef s b) (FGLMEnv s r ord n) (STRef s b)
-> Machine s r ord n b
look Getting
(STRef s [OrderedPolynomial r o n])
(FGLMEnv s r o n)
(STRef s [OrderedPolynomial r o n])
forall k s r (ord :: k) (n :: Nat).
Lens' (FGLMEnv s r ord n) (STRef s [OrderedPolynomial r ord n])
bLex
let g :: OrderedPolynomial r Lex n
g = Lex -> OrderedPolynomial r o n -> OrderedPolynomial r Lex n
forall k (n :: Nat) o o'.
(CoeffRing k, Eq (Monomial n), IsMonomialOrder n o,
IsMonomialOrder n o', KnownNat n) =>
o' -> OrderedPolynomial k o n -> OrderedPolynomial k o' n
changeOrder Lex
Lex (OrderedPolynomial r o n -> OrderedPolynomial r Lex n)
-> OrderedPolynomial r o n -> OrderedPolynomial r Lex n
forall a b. (a -> b) -> a -> b
$ OrderedPolynomial r o n
f OrderedPolynomial r o n
-> OrderedPolynomial r o n -> OrderedPolynomial r o n
forall r. Group r => r -> r -> r
- [OrderedPolynomial r o n] -> OrderedPolynomial r o n
forall (f :: * -> *) m. (Foldable f, Monoidal m) => f m -> m
sum ((r -> OrderedPolynomial r o n -> OrderedPolynomial r o n)
-> [r] -> [OrderedPolynomial r o n] -> [OrderedPolynomial r o n]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith r -> OrderedPolynomial r o n -> OrderedPolynomial r o n
forall r m. Module (Scalar r) m => r -> m -> m
(.*.) (Vector r -> [r]
forall a. Vector a -> [a]
V.toList (Vector r -> [r]) -> Vector r -> [r]
forall a b. (a -> b) -> a -> b
$ (WrapAlgebra r -> r) -> Vector (WrapAlgebra r) -> Vector r
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap WrapAlgebra r -> r
forall a. WrapAlgebra a -> a
unwrapAlgebra Vector (WrapAlgebra r)
cs) [OrderedPolynomial r o n]
bps)
(STRef s (Maybe (OrderedPolynomial r Lex n))
-> Const
(STRef s (Maybe (OrderedPolynomial r Lex n)))
(STRef s (Maybe (OrderedPolynomial r Lex n))))
-> FGLMEnv s r o n
-> Const
(STRef s (Maybe (OrderedPolynomial r Lex n))) (FGLMEnv s r o n)
forall k s r (ord :: k) (n :: Nat).
Lens'
(FGLMEnv s r ord n) (STRef s (Maybe (OrderedPolynomial r Lex n)))
proced ((STRef s (Maybe (OrderedPolynomial r Lex n))
-> Const
(STRef s (Maybe (OrderedPolynomial r Lex n)))
(STRef s (Maybe (OrderedPolynomial r Lex n))))
-> FGLMEnv s r o n
-> Const
(STRef s (Maybe (OrderedPolynomial r Lex n))) (FGLMEnv s r o n))
-> Maybe (OrderedPolynomial r Lex n) -> Machine s r o n ()
forall (t :: (* -> *) -> * -> *) s s1 a.
(MonadTrans t, MonadReader s (t (ST s1))) =>
Getting (STRef s1 a) s (STRef s1 a) -> a -> t (ST s1) ()
.== OrderedPolynomial r Lex n -> Maybe (OrderedPolynomial r Lex n)
forall a. a -> Maybe a
Just (Lex -> OrderedPolynomial r o n -> OrderedPolynomial r Lex n
forall k (n :: Nat) o o'.
(CoeffRing k, Eq (Monomial n), IsMonomialOrder n o,
IsMonomialOrder n o', KnownNat n) =>
o' -> OrderedPolynomial k o n -> OrderedPolynomial k o' n
changeOrder Lex
Lex OrderedPolynomial r o n
f)
(STRef s [OrderedPolynomial r Lex n]
-> Const
(STRef s [OrderedPolynomial r Lex n])
(STRef s [OrderedPolynomial r Lex n]))
-> FGLMEnv s r o n
-> Const (STRef s [OrderedPolynomial r Lex n]) (FGLMEnv s r o n)
forall k s r (ord :: k) (n :: Nat).
Lens' (FGLMEnv s r ord n) (STRef s [OrderedPolynomial r Lex n])
gLex ((STRef s [OrderedPolynomial r Lex n]
-> Const
(STRef s [OrderedPolynomial r Lex n])
(STRef s [OrderedPolynomial r Lex n]))
-> FGLMEnv s r o n
-> Const (STRef s [OrderedPolynomial r Lex n]) (FGLMEnv s r o n))
-> ([OrderedPolynomial r Lex n] -> [OrderedPolynomial r Lex n])
-> Machine s r o n ()
forall (t :: (* -> *) -> * -> *) s s1 a.
(MonadTrans t, MonadReader s (t (ST s1))) =>
Getting (STRef s1 a) s (STRef s1 a) -> (a -> a) -> t (ST s1) ()
@== (OrderedPolynomial r Lex n
g OrderedPolynomial r Lex n
-> [OrderedPolynomial r Lex n] -> [OrderedPolynomial r Lex n]
forall a. a -> [a] -> [a]
:)
toContinue ::
forall s r o n.
( 0 < n
, Ord r
, KnownNat n
, Field r
) =>
Machine s r o n Bool
toContinue :: Machine s r o n Bool
toContinue = do
Maybe (OrderedPolynomial r Lex n)
mans <- Getting
(STRef s (Maybe (OrderedPolynomial r Lex n)))
(FGLMEnv s r o n)
(STRef s (Maybe (OrderedPolynomial r Lex n)))
-> Machine s r o n (Maybe (OrderedPolynomial r Lex n))
forall k s b r (ord :: k) (n :: Nat).
Getting (STRef s b) (FGLMEnv s r ord n) (STRef s b)
-> Machine s r ord n b
look Getting
(STRef s (Maybe (OrderedPolynomial r Lex n)))
(FGLMEnv s r o n)
(STRef s (Maybe (OrderedPolynomial r Lex n)))
forall k s r (ord :: k) (n :: Nat).
Lens'
(FGLMEnv s r ord n) (STRef s (Maybe (OrderedPolynomial r Lex n)))
proced
case Maybe (OrderedPolynomial r Lex n)
mans of
Maybe (OrderedPolynomial r Lex n)
Nothing -> Bool -> Machine s r o n Bool
forall (m :: * -> *) a. Monad m => a -> m a
return Bool
True
Just OrderedPolynomial r Lex n
g -> do
let xLast :: OrderedPolynomial r Lex n
xLast = Sized Vector n (OrderedPolynomial r Lex n)
-> OrderedPolynomial r Lex n
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
P.maximum Sized Vector n (OrderedPolynomial r Lex n)
forall k ord (n :: Nat).
(IsMonomialOrder n ord, CoeffRing k, KnownNat n) =>
Sized n (OrderedPolynomial k ord n)
allVars OrderedPolynomial r Lex n
-> OrderedPolynomial r Lex n -> OrderedPolynomial r Lex n
forall a. a -> a -> a
`asTypeOf` OrderedPolynomial r Lex n
g
Bool -> Machine s r o n Bool
forall (m :: * -> *) a. Monad m => a -> m a
return (Bool -> Machine s r o n Bool) -> Bool -> Machine s r o n Bool
forall a b. (a -> b) -> a -> b
$ Bool -> Bool
not (Bool -> Bool) -> Bool -> Bool
forall a b. (a -> b) -> a -> b
$ OrderedPolynomial r Lex n
-> OrderedMonomial
(MOrder (OrderedPolynomial r Lex n))
(Arity (OrderedPolynomial r Lex n))
forall poly.
IsOrderedPolynomial poly =>
poly -> OrderedMonomial (MOrder poly) (Arity poly)
leadingMonomial OrderedPolynomial r Lex n
g OrderedMonomial Lex n -> OrderedMonomial Lex n -> Bool
forall k (n :: Nat) (ord :: k).
KnownNat n =>
OrderedMonomial ord n -> OrderedMonomial ord n -> Bool
`isPowerOf` OrderedPolynomial r Lex n
-> OrderedMonomial
(MOrder (OrderedPolynomial r Lex n))
(Arity (OrderedPolynomial r Lex n))
forall poly.
IsOrderedPolynomial poly =>
poly -> OrderedMonomial (MOrder poly) (Arity poly)
leadingMonomial OrderedPolynomial r Lex n
xLast
nextMonomial ::
forall s r ord n.
(CoeffRing r, KnownNat n) =>
Machine s r ord n ()
nextMonomial :: Machine s r ord n ()
nextMonomial = do
OrderedMonomial Lex n
m <- Getting
(STRef s (OrderedMonomial Lex n))
(FGLMEnv s r ord n)
(STRef s (OrderedMonomial Lex n))
-> Machine s r ord n (OrderedMonomial Lex n)
forall k s b r (ord :: k) (n :: Nat).
Getting (STRef s b) (FGLMEnv s r ord n) (STRef s b)
-> Machine s r ord n b
look Getting
(STRef s (OrderedMonomial Lex n))
(FGLMEnv s r ord n)
(STRef s (OrderedMonomial Lex n))
forall k s r (ord :: k) (n :: Nat).
Lens' (FGLMEnv s r ord n) (STRef s (OrderedMonomial Lex n))
monomial
[OrderedMonomial Lex n]
gs <- (OrderedPolynomial r Lex n -> OrderedMonomial Lex n)
-> [OrderedPolynomial r Lex n] -> [OrderedMonomial Lex n]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
map OrderedPolynomial r Lex n -> OrderedMonomial Lex n
forall poly.
IsOrderedPolynomial poly =>
poly -> OrderedMonomial (MOrder poly) (Arity poly)
leadingMonomial ([OrderedPolynomial r Lex n] -> [OrderedMonomial Lex n])
-> ReaderT (FGLMEnv s r ord n) (ST s) [OrderedPolynomial r Lex n]
-> ReaderT (FGLMEnv s r ord n) (ST s) [OrderedMonomial Lex n]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Getting
(STRef s [OrderedPolynomial r Lex n])
(FGLMEnv s r ord n)
(STRef s [OrderedPolynomial r Lex n])
-> ReaderT (FGLMEnv s r ord n) (ST s) [OrderedPolynomial r Lex n]
forall k s b r (ord :: k) (n :: Nat).
Getting (STRef s b) (FGLMEnv s r ord n) (STRef s b)
-> Machine s r ord n b
look Getting
(STRef s [OrderedPolynomial r Lex n])
(FGLMEnv s r ord n)
(STRef s [OrderedPolynomial r Lex n])
forall k s r (ord :: k) (n :: Nat).
Lens' (FGLMEnv s r ord n) (STRef s [OrderedPolynomial r Lex n])
gLex
let next :: OrderedMonomial Lex n
next =
(OrderedMonomial Lex n, Int) -> OrderedMonomial Lex n
forall a b. (a, b) -> a
fst ((OrderedMonomial Lex n, Int) -> OrderedMonomial Lex n)
-> (OrderedMonomial Lex n, Int) -> OrderedMonomial Lex n
forall a b. (a -> b) -> a -> b
$
((OrderedMonomial Lex n, Int)
-> (OrderedMonomial Lex n, Int) -> Ordering)
-> [(OrderedMonomial Lex n, Int)] -> (OrderedMonomial Lex n, Int)
forall (t :: * -> *) a.
Foldable t =>
(a -> a -> Ordering) -> t a -> a
maximumBy
(((OrderedMonomial Lex n, Int) -> Int)
-> (OrderedMonomial Lex n, Int)
-> (OrderedMonomial Lex n, Int)
-> Ordering
forall a b. Ord a => (b -> a) -> b -> b -> Ordering
comparing (OrderedMonomial Lex n, Int) -> Int
forall a b. (a, b) -> b
snd)
[ (Monomial n -> OrderedMonomial Lex n
forall k (ordering :: k) (n :: Nat).
Monomial n -> OrderedMonomial ordering n
OrderedMonomial Monomial n
monom, Ordinal n -> Int
forall a. Enum a => a -> Int
fromEnum Ordinal n
od)
| Ordinal n
od <- SNat n -> [Ordinal n]
forall (n :: Nat). SNat n -> [Ordinal n]
enumOrdinal (SNat n
forall (n :: Nat). KnownNat n => SNat n
sNat :: SNat n)
, let monom :: Monomial n
monom = Monomial n -> Ordinal n -> Monomial n
forall (n :: Nat).
KnownNat n =>
Monomial n -> Ordinal n -> Monomial n
beta (OrderedMonomial Lex n -> Monomial n
forall k (ordering :: k) (n :: Nat).
OrderedMonomial ordering n -> Monomial n
getMonomial OrderedMonomial Lex n
m) Ordinal n
od
, (OrderedMonomial Lex n -> Bool) -> [OrderedMonomial Lex n] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (Bool -> Bool
not (Bool -> Bool)
-> (OrderedMonomial Lex n -> Bool) -> OrderedMonomial Lex n -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (OrderedMonomial Lex n -> OrderedMonomial Lex n -> Bool
forall k (n :: Nat) (ord :: k).
KnownNat n =>
OrderedMonomial ord n -> OrderedMonomial ord n -> Bool
`divs` Monomial n -> OrderedMonomial Lex n
forall k (ordering :: k) (n :: Nat).
Monomial n -> OrderedMonomial ordering n
OrderedMonomial Monomial n
monom)) [OrderedMonomial Lex n]
gs
]
Getting
(STRef s (OrderedMonomial Lex n))
(FGLMEnv s r ord n)
(STRef s (OrderedMonomial Lex n))
forall k s r (ord :: k) (n :: Nat).
Lens' (FGLMEnv s r ord n) (STRef s (OrderedMonomial Lex n))
monomial Getting
(STRef s (OrderedMonomial Lex n))
(FGLMEnv s r ord n)
(STRef s (OrderedMonomial Lex n))
-> OrderedMonomial Lex n -> Machine s r ord n ()
forall (t :: (* -> *) -> * -> *) s s1 a.
(MonadTrans t, MonadReader s (t (ST s1))) =>
Getting (STRef s1 a) s (STRef s1 a) -> a -> t (ST s1) ()
.== OrderedMonomial Lex n
next
beta :: KnownNat n => Monomial n -> Ordinal n -> Monomial n
beta :: Monomial n -> Ordinal n -> Monomial n
beta Monomial n
xs o :: Ordinal n
o@(OLt SNat n1
k) =
let n :: SNat n
n = Monomial n -> SNat n
forall (n :: Nat) (f :: * -> *) a.
(KnownNat n, DomC f a) =>
Sized f n a -> SNat n
sizedLength Monomial n
xs
in (SNat (n1 + 1) -> Monomial n -> Sized Vector (n1 + 1) Int
forall (n :: Nat) (f :: * -> *) (m :: Nat) a.
(CFreeMonoid f, Dom f a, n <= m) =>
SNat n -> Sized f m a -> Sized f n a
SV.take (SNat n1 -> SNat (n1 + 1)
forall (n :: Nat). SNat n -> SNat (Succ n)
sSucc SNat n1
k) (Monomial n -> Sized Vector (n1 + 1) Int)
-> Monomial n -> Sized Vector (n1 + 1) Int
forall a b. (a -> b) -> a -> b
$ Monomial n
xs Monomial n -> (Monomial n -> Monomial n) -> Monomial n
forall a b. a -> (a -> b) -> b
& Index (Monomial n)
-> Traversal' (Monomial n) (IxValue (Monomial n))
forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix Ordinal n
Index (Monomial n)
o ((Int -> Identity Int) -> Monomial n -> Identity (Monomial n))
-> Int -> Monomial n -> Monomial n
forall a s t. Num a => ASetter s t a a -> a -> s -> t
+~ Int
1) Sized Vector (n1 + 1) Int
-> Sized Vector (n - (n1 + 1)) Int
-> Sized Vector ((n1 + 1) + (n - (n1 + 1))) Int
forall (f :: * -> *) (n :: Nat) (m :: Nat) a.
(CFreeMonoid f, Dom f a) =>
Sized f n a -> Sized f m a -> Sized f (n + m) a
SV.++ SNat (n - (n1 + 1)) -> Int -> Sized Vector (n - (n1 + 1)) Int
forall (f :: * -> *) (n :: Nat) a.
(CFreeMonoid f, Dom f a) =>
SNat n -> a -> Sized f n a
SV.replicate (SNat n
n SNat n -> SNat (n1 + 1) -> SNat (n - (n1 + 1))
forall (n :: Nat) (m :: Nat). SNat n -> SNat m -> SNat (n - m)
%- SNat n1 -> SNat (n1 + 1)
forall (n :: Nat). SNat n -> SNat (Succ n)
sSucc SNat n1
k) Int
0
beta Monomial n
_ Ordinal n
_ = [Char] -> Monomial n
forall a. HasCallStack => [Char] -> a
error [Char]
"beta: Bug in ghc!"