{-# 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 #-}

{- | Algorithms for zero-dimensional ideals.

   Since 0.4.0.0
-}
module Algebra.Algorithms.ZeroDim
  ( -- * Root finding for zero-dimensional ideal
    solveM,
    solve',
    solveViaCompanion,
    solveLinear,

    -- * Radical computation
    radical,
    isRadical,

    -- * Converting monomial ordering to Lex using FGLM algorithm
    fglm,
    fglmMap,

    -- ** Internal helper function
    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

{- | Finds complex approximate  roots of given zero-dimensional ideal,
   using randomized altorithm.

   See also @'solve''@ and @'solveViaCompanion'@.
-}
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 f is@ finds complex approximate roots of the given zero-dimensional @n@-variate polynomial system @is@,
   using the given relatively prime polynomial @f@.
-}
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'' err is@ finds numeric approximate root of the
   given zero-dimensional polynomial system @is@,
   with error <@err@.

   See also @'solveViaCompanion'@ and @'solveM'@.
-}
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
            ]

{- | Given a zero-dimensional ideal \(I\),
   \('subspMatrix' i I\) computes a multiplication matrix
   \(m_{x_i} \mathbin{\upharpoonright} V_i \) by \(x_i\)
   restricted to the subspace \(V_i = \mathop{\mathrm{span}}(\left\{\ x_i^n \ \middle|\ 1 \leq n \right\})\) of \(k[\mathbf{X}]\).
-}
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' err is@ finds numeric approximate root of the
   given zero-dimensional polynomial system @is@,
   with error <@err@.

   See also @'solve''@ and @'solveM'@.
-}
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

-- | Calculates n-th reduction of f: @f `div` < f, ∂_{x_n} f >@.
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])

-- | Calculate the monic generator of \( k[X_0, ..., X_n] \cap k[X_i]\).
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

-- | Solves linear system. If the given matrix is degenerate, this returns @Nothing@.
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

-- | Calculate the radical of the given zero-dimensional ideal.
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

-- | Test if the given zero-dimensional ideal is radical or not.
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

{- | Calculate the Groebner basis w.r.t. lex ordering of the zero-dimensional ideal using FGLM algorithm.
   If the given ideal is not zero-dimensional this function may diverge.
-}
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)

-- | Compute the kernel and image of the given linear map using generalized FGLM algorithm.
fglmMap ::
  forall k ord n.
  ( Ord k
  , Field k
  , 0 < n
  , IsMonomialOrder n ord
  , CoeffRing k
  , KnownNat n
  ) =>
  -- | Linear map from polynomial ring.
  (OrderedPolynomial k ord n -> V.Vector k) ->
  -- | The tuple of:
  --
  --     * lex-Groebner basis of the kernel of the given linear map.
  --
  --     * The vector basis of the image of the linear map.
  ( [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!"