{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE CPP #-}
{-# LANGUAGE DeriveDataTypeable #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE FunctionalDependencies #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeFamilies #-}
{-# OPTIONS_GHC -Wall #-}
{-# OPTIONS_HADDOCK show-extensions #-}
-----------------------------------------------------------------------------
-- |
-- Module      :  ToySolver.Data.Polynomial.Base
-- Copyright   :  (c) Masahiro Sakai 2012-2013
-- License     :  BSD-style
--
-- Maintainer  :  masahiro.sakai@gmail.com
-- Stability   :  provisional
-- Portability :  non-portable
--
-- Polynomials
--
-- References:
--
-- * Monomial order <http://en.wikipedia.org/wiki/Monomial_order>
--
-- * Polynomial class for Ruby <http://www.math.kobe-u.ac.jp/~kodama/tips-RubyPoly.html>
--
-- * constructive-algebra package <http://hackage.haskell.org/package/constructive-algebra>
--
-----------------------------------------------------------------------------
module ToySolver.Data.Polynomial.Base
  (
  -- * Polynomial type
    Polynomial

  -- * Conversion
  , Var (..)
  , constant
  , terms
  , fromTerms
  , coeffMap
  , fromCoeffMap
  , fromTerm

  -- * Query
  , Degree (..)
  , Vars (..)
  , lt
  , lc
  , lm
  , coeff
  , lookupCoeff
  , isPrimitive

  -- * Operations
  , Factor (..)
  , SQFree (..)
  , ContPP (..)
  , deriv
  , integral
  , eval
  , subst
  , mapCoeff
  , toMonic
  , toUPolynomialOf
  , divModMP
  , reduce

  -- * Univariate polynomials
  , UPolynomial
  , X (..)
  , UTerm
  , UMonomial
  , div
  , mod
  , divMod
  , divides
  , gcd
  , lcm
  , exgcd
  , pdivMod
  , pdiv
  , pmod
  , gcd'
  , isRootOf
  , isSquareFree
  , nat
  , eisensteinsCriterion

  -- * Term
  , Term
  , tdeg
  , tscale
  , tmult
  , tdivides
  , tdiv
  , tderiv
  , tintegral

  -- * Monic monomial
  , Monomial
  , mone
  , mfromIndices
  , mfromIndicesMap
  , mindices
  , mindicesMap
  , mmult
  , mpow
  , mdivides
  , mdiv
  , mderiv
  , mintegral
  , mlcm
  , mgcd
  , mcoprime

  -- * Monomial order
  , MonomialOrder
  , lex
  , revlex
  , grlex
  , grevlex

  -- * Pretty Printing
  , PrintOptions (..)
  , prettyPrint
  , PrettyCoeff (..)
  , PrettyVar (..)
  ) where

import Prelude hiding (lex, div, mod, divMod, gcd, lcm)
import qualified Prelude
import Control.DeepSeq
import Control.Exception (assert)
import Control.Monad
import Data.Default.Class
import Data.Data
import qualified Data.FiniteField as FF
import Data.Function
import Data.Hashable
import Data.List
#if !MIN_VERSION_base(4,11,0)
import Data.Monoid
#endif
import Data.Numbers.Primes (primeFactors)
import Data.Ratio
import Data.String (IsString (..))
import Data.Map.Strict (Map)
import qualified Data.Map.Strict as Map
import Data.Set (Set)
import qualified Data.Set as Set
import Data.IntMap.Strict (IntMap)
import qualified Data.IntMap.Strict as IntMap
import Data.VectorSpace
import qualified Text.PrettyPrint.HughesPJClass as PP
import Text.PrettyPrint.HughesPJClass (Doc, PrettyLevel, Pretty (..), maybeParens)

infixl 7  `div`, `mod`

{--------------------------------------------------------------------
  Classes
--------------------------------------------------------------------}

class Vars a v => Var a v | a -> v where
  var :: v -> a

class Ord v => Vars a v | a -> v where
  vars :: a -> Set v

-- | total degree of a given polynomial
class Degree t where
  deg :: t -> Integer

{--------------------------------------------------------------------
  Polynomial type
--------------------------------------------------------------------}

-- | Polynomial over commutative ring r
newtype Polynomial r v = Polynomial{ Polynomial r v -> Map (Monomial v) r
coeffMap :: Map (Monomial v) r }
  deriving (Polynomial r v -> Polynomial r v -> Bool
(Polynomial r v -> Polynomial r v -> Bool)
-> (Polynomial r v -> Polynomial r v -> Bool)
-> Eq (Polynomial r v)
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
forall r v.
(Eq v, Eq r) =>
Polynomial r v -> Polynomial r v -> Bool
/= :: Polynomial r v -> Polynomial r v -> Bool
$c/= :: forall r v.
(Eq v, Eq r) =>
Polynomial r v -> Polynomial r v -> Bool
== :: Polynomial r v -> Polynomial r v -> Bool
$c== :: forall r v.
(Eq v, Eq r) =>
Polynomial r v -> Polynomial r v -> Bool
Eq, Eq (Polynomial r v)
Eq (Polynomial r v)
-> (Polynomial r v -> Polynomial r v -> Ordering)
-> (Polynomial r v -> Polynomial r v -> Bool)
-> (Polynomial r v -> Polynomial r v -> Bool)
-> (Polynomial r v -> Polynomial r v -> Bool)
-> (Polynomial r v -> Polynomial r v -> Bool)
-> (Polynomial r v -> Polynomial r v -> Polynomial r v)
-> (Polynomial r v -> Polynomial r v -> Polynomial r v)
-> Ord (Polynomial r v)
Polynomial r v -> Polynomial r v -> Bool
Polynomial r v -> Polynomial r v -> Ordering
Polynomial r v -> Polynomial r v -> Polynomial r v
forall a.
Eq a
-> (a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
forall r v. (Ord v, Ord r) => Eq (Polynomial r v)
forall r v.
(Ord v, Ord r) =>
Polynomial r v -> Polynomial r v -> Bool
forall r v.
(Ord v, Ord r) =>
Polynomial r v -> Polynomial r v -> Ordering
forall r v.
(Ord v, Ord r) =>
Polynomial r v -> Polynomial r v -> Polynomial r v
min :: Polynomial r v -> Polynomial r v -> Polynomial r v
$cmin :: forall r v.
(Ord v, Ord r) =>
Polynomial r v -> Polynomial r v -> Polynomial r v
max :: Polynomial r v -> Polynomial r v -> Polynomial r v
$cmax :: forall r v.
(Ord v, Ord r) =>
Polynomial r v -> Polynomial r v -> Polynomial r v
>= :: Polynomial r v -> Polynomial r v -> Bool
$c>= :: forall r v.
(Ord v, Ord r) =>
Polynomial r v -> Polynomial r v -> Bool
> :: Polynomial r v -> Polynomial r v -> Bool
$c> :: forall r v.
(Ord v, Ord r) =>
Polynomial r v -> Polynomial r v -> Bool
<= :: Polynomial r v -> Polynomial r v -> Bool
$c<= :: forall r v.
(Ord v, Ord r) =>
Polynomial r v -> Polynomial r v -> Bool
< :: Polynomial r v -> Polynomial r v -> Bool
$c< :: forall r v.
(Ord v, Ord r) =>
Polynomial r v -> Polynomial r v -> Bool
compare :: Polynomial r v -> Polynomial r v -> Ordering
$ccompare :: forall r v.
(Ord v, Ord r) =>
Polynomial r v -> Polynomial r v -> Ordering
$cp1Ord :: forall r v. (Ord v, Ord r) => Eq (Polynomial r v)
Ord, Typeable)

instance (Eq k, Num k, Ord v) => Num (Polynomial k v) where
  + :: Polynomial k v -> Polynomial k v -> Polynomial k v
(+)      = Polynomial k v -> Polynomial k v -> Polynomial k v
forall k v.
(Eq k, Num k, Ord v) =>
Polynomial k v -> Polynomial k v -> Polynomial k v
plus
  * :: Polynomial k v -> Polynomial k v -> Polynomial k v
(*)      = Polynomial k v -> Polynomial k v -> Polynomial k v
forall k v.
(Eq k, Num k, Ord v) =>
Polynomial k v -> Polynomial k v -> Polynomial k v
mult
  negate :: Polynomial k v -> Polynomial k v
negate   = Polynomial k v -> Polynomial k v
forall k v. Num k => Polynomial k v -> Polynomial k v
neg
  abs :: Polynomial k v -> Polynomial k v
abs Polynomial k v
x    = Polynomial k v
x -- OK?
  signum :: Polynomial k v -> Polynomial k v
signum Polynomial k v
_ = Polynomial k v
1 -- OK?
  fromInteger :: Integer -> Polynomial k v
fromInteger Integer
x = k -> Polynomial k v
forall k v. (Eq k, Num k) => k -> Polynomial k v
constant (Integer -> k
forall a. Num a => Integer -> a
fromInteger Integer
x)

instance (Eq k, Num k, Ord v, IsString v) => IsString (Polynomial k v) where
  fromString :: String -> Polynomial k v
fromString String
s = v -> Polynomial k v
forall a v. Var a v => v -> a
var (String -> v
forall a. IsString a => String -> a
fromString String
s)

instance (Eq k, Num k, Ord v) => AdditiveGroup (Polynomial k v) where
  ^+^ :: Polynomial k v -> Polynomial k v -> Polynomial k v
(^+^)   = Polynomial k v -> Polynomial k v -> Polynomial k v
forall k v.
(Eq k, Num k, Ord v) =>
Polynomial k v -> Polynomial k v -> Polynomial k v
plus
  zeroV :: Polynomial k v
zeroV   = Polynomial k v
forall k v. Polynomial k v
zero
  negateV :: Polynomial k v -> Polynomial k v
negateV = Polynomial k v -> Polynomial k v
forall k v. Num k => Polynomial k v -> Polynomial k v
neg

instance (Eq k, Num k, Ord v) => VectorSpace (Polynomial k v) where
  type Scalar (Polynomial k v) = k
  Scalar (Polynomial k v)
k *^ :: Scalar (Polynomial k v) -> Polynomial k v -> Polynomial k v
*^ Polynomial k v
p = k -> Polynomial k v -> Polynomial k v
forall k v. (Eq k, Num k) => k -> Polynomial k v -> Polynomial k v
scale k
Scalar (Polynomial k v)
k Polynomial k v
p

instance (Show v, Show k) => Show (Polynomial k v) where
  showsPrec :: Int -> Polynomial k v -> ShowS
showsPrec Int
d Polynomial k v
p  = Bool -> ShowS -> ShowS
showParen (Int
d Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
10) (ShowS -> ShowS) -> ShowS -> ShowS
forall a b. (a -> b) -> a -> b
$
    String -> ShowS
showString String
"fromTerms " ShowS -> ShowS -> ShowS
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Term k v] -> ShowS
forall a. Show a => a -> ShowS
shows (Polynomial k v -> [Term k v]
forall k v. Polynomial k v -> [Term k v]
terms Polynomial k v
p)

instance (NFData k, NFData v) => NFData (Polynomial k v) where
  rnf :: Polynomial k v -> ()
rnf (Polynomial Map (Monomial v) k
m) = Map (Monomial v) k -> ()
forall a. NFData a => a -> ()
rnf Map (Monomial v) k
m

instance (Hashable k, Hashable v) => Hashable (Polynomial k v) where
  hashWithSalt :: Int -> Polynomial k v -> Int
hashWithSalt = (Polynomial k v -> [(Monomial v, k)])
-> Int -> Polynomial k v -> Int
forall b a. Hashable b => (a -> b) -> Int -> a -> Int
hashUsing (Map (Monomial v) k -> [(Monomial v, k)]
forall k a. Map k a -> [(k, a)]
Map.toList (Map (Monomial v) k -> [(Monomial v, k)])
-> (Polynomial k v -> Map (Monomial v) k)
-> Polynomial k v
-> [(Monomial v, k)]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Polynomial k v -> Map (Monomial v) k
forall r v. Polynomial r v -> Map (Monomial v) r
coeffMap)

instance (Eq k, Num k, Ord v) => Var (Polynomial k v) v where
  var :: v -> Polynomial k v
var v
x = Term k v -> Polynomial k v
forall k v. (Eq k, Num k) => Term k v -> Polynomial k v
fromTerm (k
1, v -> Monomial v
forall a v. Var a v => v -> a
var v
x)

instance (Ord v) => Vars (Polynomial k v) v where
  vars :: Polynomial k v -> Set v
vars Polynomial k v
p = [Set v] -> Set v
forall (f :: * -> *) a. (Foldable f, Ord a) => f (Set a) -> Set a
Set.unions ([Set v] -> Set v) -> [Set v] -> Set v
forall a b. (a -> b) -> a -> b
$ [Monomial v -> Set v
forall a v. Vars a v => a -> Set v
vars Monomial v
mm | (k
_, Monomial v
mm) <- Polynomial k v -> [(k, Monomial v)]
forall k v. Polynomial k v -> [Term k v]
terms Polynomial k v
p]

instance Degree (Polynomial k v) where
  deg :: Polynomial k v -> Integer
deg Polynomial k v
p
    | Polynomial k v -> Bool
forall k v. Polynomial k v -> Bool
isZero Polynomial k v
p  = -Integer
1
    | Bool
otherwise = [Integer] -> Integer
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [Monomial v -> Integer
forall t. Degree t => t -> Integer
deg Monomial v
mm | (k
_,Monomial v
mm) <- Polynomial k v -> [(k, Monomial v)]
forall k v. Polynomial k v -> [Term k v]
terms Polynomial k v
p]

normalize :: (Eq k, Num k) => Polynomial k v -> Polynomial k v
normalize :: Polynomial k v -> Polynomial k v
normalize (Polynomial Map (Monomial v) k
m) = Map (Monomial v) k -> Polynomial k v
forall r v. Map (Monomial v) r -> Polynomial r v
Polynomial ((k -> Bool) -> Map (Monomial v) k -> Map (Monomial v) k
forall a k. (a -> Bool) -> Map k a -> Map k a
Map.filter (k
0k -> k -> Bool
forall a. Eq a => a -> a -> Bool
/=) Map (Monomial v) k
m)

asConstant :: Num k => Polynomial k v -> Maybe k
asConstant :: Polynomial k v -> Maybe k
asConstant Polynomial k v
p =
  case Polynomial k v -> [Term k v]
forall k v. Polynomial k v -> [Term k v]
terms Polynomial k v
p of
    [] -> k -> Maybe k
forall a. a -> Maybe a
Just k
0
    [(k
c,Monomial v
xs)] | Map v Integer -> Bool
forall k a. Map k a -> Bool
Map.null (Monomial v -> Map v Integer
forall v. Monomial v -> Map v Integer
mindicesMap Monomial v
xs) -> k -> Maybe k
forall a. a -> Maybe a
Just k
c
    [Term k v]
_ -> Maybe k
forall a. Maybe a
Nothing

scale :: (Eq k, Num k) => k -> Polynomial k v -> Polynomial k v
scale :: k -> Polynomial k v -> Polynomial k v
scale k
0 Polynomial k v
_ = Polynomial k v
forall k v. Polynomial k v
zero
scale k
1 Polynomial k v
p = Polynomial k v
p
scale k
a (Polynomial Map (Monomial v) k
m) = Map (Monomial v) k -> Polynomial k v
forall r v. Map (Monomial v) r -> Polynomial r v
Polynomial ((k -> Maybe k) -> Map (Monomial v) k -> Map (Monomial v) k
forall a b k. (a -> Maybe b) -> Map k a -> Map k b
Map.mapMaybe k -> Maybe k
f Map (Monomial v) k
m)
  where
    f :: k -> Maybe k
f k
b = if k
c k -> k -> Bool
forall a. Eq a => a -> a -> Bool
== k
0 then Maybe k
forall a. Maybe a
Nothing else k -> Maybe k
forall a. a -> Maybe a
Just k
c
      where c :: k
c = k
a k -> k -> k
forall a. Num a => a -> a -> a
* k
b

zero :: Polynomial k v
zero :: Polynomial k v
zero = Map (Monomial v) k -> Polynomial k v
forall r v. Map (Monomial v) r -> Polynomial r v
Polynomial (Map (Monomial v) k -> Polynomial k v)
-> Map (Monomial v) k -> Polynomial k v
forall a b. (a -> b) -> a -> b
$ Map (Monomial v) k
forall k a. Map k a
Map.empty

plus :: (Eq k, Num k, Ord v) => Polynomial k v -> Polynomial k v -> Polynomial k v
plus :: Polynomial k v -> Polynomial k v -> Polynomial k v
plus (Polynomial Map (Monomial v) k
m1) (Polynomial Map (Monomial v) k
m2) = Map (Monomial v) k -> Polynomial k v
forall r v. Map (Monomial v) r -> Polynomial r v
Polynomial (Map (Monomial v) k -> Polynomial k v)
-> Map (Monomial v) k -> Polynomial k v
forall a b. (a -> b) -> a -> b
$ (Monomial v -> k -> k -> Maybe k)
-> (Map (Monomial v) k -> Map (Monomial v) k)
-> (Map (Monomial v) k -> Map (Monomial v) k)
-> Map (Monomial v) k
-> Map (Monomial v) k
-> Map (Monomial v) k
forall k a b c.
Ord k =>
(k -> a -> b -> Maybe c)
-> (Map k a -> Map k c)
-> (Map k b -> Map k c)
-> Map k a
-> Map k b
-> Map k c
Map.mergeWithKey Monomial v -> k -> k -> Maybe k
forall a p. (Num a, Eq a) => p -> a -> a -> Maybe a
f Map (Monomial v) k -> Map (Monomial v) k
forall a. a -> a
id Map (Monomial v) k -> Map (Monomial v) k
forall a. a -> a
id Map (Monomial v) k
m1 Map (Monomial v) k
m2
  where
    f :: p -> a -> a -> Maybe a
f p
_ a
a a
b = if a
c a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 then Maybe a
forall a. Maybe a
Nothing else a -> Maybe a
forall a. a -> Maybe a
Just a
c
      where
        c :: a
c = a
a a -> a -> a
forall a. Num a => a -> a -> a
+ a
b

neg :: (Num k) => Polynomial k v -> Polynomial k v
neg :: Polynomial k v -> Polynomial k v
neg (Polynomial Map (Monomial v) k
m) = Map (Monomial v) k -> Polynomial k v
forall r v. Map (Monomial v) r -> Polynomial r v
Polynomial (Map (Monomial v) k -> Polynomial k v)
-> Map (Monomial v) k -> Polynomial k v
forall a b. (a -> b) -> a -> b
$ (k -> k) -> Map (Monomial v) k -> Map (Monomial v) k
forall a b k. (a -> b) -> Map k a -> Map k b
Map.map k -> k
forall a. Num a => a -> a
negate Map (Monomial v) k
m

mult :: (Eq k, Num k, Ord v) => Polynomial k v -> Polynomial k v -> Polynomial k v
mult :: Polynomial k v -> Polynomial k v -> Polynomial k v
mult Polynomial k v
a Polynomial k v
b
  | Just k
c <- Polynomial k v -> Maybe k
forall k v. Num k => Polynomial k v -> Maybe k
asConstant Polynomial k v
a = k -> Polynomial k v -> Polynomial k v
forall k v. (Eq k, Num k) => k -> Polynomial k v -> Polynomial k v
scale k
c Polynomial k v
b
  | Just k
c <- Polynomial k v -> Maybe k
forall k v. Num k => Polynomial k v -> Maybe k
asConstant Polynomial k v
b = k -> Polynomial k v -> Polynomial k v
forall k v. (Eq k, Num k) => k -> Polynomial k v -> Polynomial k v
scale k
c Polynomial k v
a
mult (Polynomial Map (Monomial v) k
m1) (Polynomial Map (Monomial v) k
m2) = Map (Monomial v) k -> Polynomial k v
forall k v. (Eq k, Num k) => Map (Monomial v) k -> Polynomial k v
fromCoeffMap (Map (Monomial v) k -> Polynomial k v)
-> Map (Monomial v) k -> Polynomial k v
forall a b. (a -> b) -> a -> b
$ (k -> k -> k) -> [(Monomial v, k)] -> Map (Monomial v) k
forall k a. Ord k => (a -> a -> a) -> [(k, a)] -> Map k a
Map.fromListWith k -> k -> k
forall a. Num a => a -> a -> a
(+)
      [ (Monomial v
xs1 Monomial v -> Monomial v -> Monomial v
forall v. Ord v => Monomial v -> Monomial v -> Monomial v
`mmult` Monomial v
xs2, k
c1k -> k -> k
forall a. Num a => a -> a -> a
*k
c2)
      | (Monomial v
xs1,k
c1) <- Map (Monomial v) k -> [(Monomial v, k)]
forall k a. Map k a -> [(k, a)]
Map.toList Map (Monomial v) k
m1, (Monomial v
xs2,k
c2) <- Map (Monomial v) k -> [(Monomial v, k)]
forall k a. Map k a -> [(k, a)]
Map.toList Map (Monomial v) k
m2
      ]

isZero :: Polynomial k v -> Bool
isZero :: Polynomial k v -> Bool
isZero (Polynomial Map (Monomial v) k
m) = Map (Monomial v) k -> Bool
forall k a. Map k a -> Bool
Map.null Map (Monomial v) k
m

-- | construct a polynomial from a constant
constant :: (Eq k, Num k) => k -> Polynomial k v
constant :: k -> Polynomial k v
constant k
c = Term k v -> Polynomial k v
forall k v. (Eq k, Num k) => Term k v -> Polynomial k v
fromTerm (k
c, Monomial v
forall v. Monomial v
mone)

-- | construct a polynomial from a list of terms
fromTerms :: (Eq k, Num k, Ord v) => [Term k v] -> Polynomial k v
fromTerms :: [Term k v] -> Polynomial k v
fromTerms = Map (Monomial v) k -> Polynomial k v
forall k v. (Eq k, Num k) => Map (Monomial v) k -> Polynomial k v
fromCoeffMap (Map (Monomial v) k -> Polynomial k v)
-> ([Term k v] -> Map (Monomial v) k)
-> [Term k v]
-> Polynomial k v
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (k -> k -> k) -> [(Monomial v, k)] -> Map (Monomial v) k
forall k a. Ord k => (a -> a -> a) -> [(k, a)] -> Map k a
Map.fromListWith k -> k -> k
forall a. Num a => a -> a -> a
(+) ([(Monomial v, k)] -> Map (Monomial v) k)
-> ([Term k v] -> [(Monomial v, k)])
-> [Term k v]
-> Map (Monomial v) k
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Term k v -> (Monomial v, k)) -> [Term k v] -> [(Monomial v, k)]
forall a b. (a -> b) -> [a] -> [b]
map (\(k
c,Monomial v
xs) -> (Monomial v
xs,k
c))

fromCoeffMap :: (Eq k, Num k) => Map (Monomial v) k -> Polynomial k v
fromCoeffMap :: Map (Monomial v) k -> Polynomial k v
fromCoeffMap Map (Monomial v) k
m = Polynomial k v -> Polynomial k v
forall k v. (Eq k, Num k) => Polynomial k v -> Polynomial k v
normalize (Polynomial k v -> Polynomial k v)
-> Polynomial k v -> Polynomial k v
forall a b. (a -> b) -> a -> b
$ Map (Monomial v) k -> Polynomial k v
forall r v. Map (Monomial v) r -> Polynomial r v
Polynomial Map (Monomial v) k
m

-- | construct a polynomial from a singlet term
fromTerm :: (Eq k, Num k) => Term k v -> Polynomial k v
fromTerm :: Term k v -> Polynomial k v
fromTerm (k
0,Monomial v
_) = Polynomial k v
forall k v. Polynomial k v
zero
fromTerm (k
c,Monomial v
xs) = Map (Monomial v) k -> Polynomial k v
forall r v. Map (Monomial v) r -> Polynomial r v
Polynomial (Map (Monomial v) k -> Polynomial k v)
-> Map (Monomial v) k -> Polynomial k v
forall a b. (a -> b) -> a -> b
$ Monomial v -> k -> Map (Monomial v) k
forall k a. k -> a -> Map k a
Map.singleton Monomial v
xs k
c

-- | list of terms
terms :: Polynomial k v -> [Term k v]
terms :: Polynomial k v -> [Term k v]
terms (Polynomial Map (Monomial v) k
m) = [(k
c,Monomial v
xs) | (Monomial v
xs,k
c) <- Map (Monomial v) k -> [(Monomial v, k)]
forall k a. Map k a -> [(k, a)]
Map.toList Map (Monomial v) k
m]

-- | leading term with respect to a given monomial order
lt :: (Num k) => MonomialOrder v -> Polynomial k v -> Term k v
lt :: MonomialOrder v -> Polynomial k v -> Term k v
lt MonomialOrder v
cmp Polynomial k v
p =
  case Polynomial k v -> [Term k v]
forall k v. Polynomial k v -> [Term k v]
terms Polynomial k v
p of
    [] -> (k
0, Monomial v
forall v. Monomial v
mone) -- should be error?
    [Term k v]
ms -> (Term k v -> Term k v -> Ordering) -> [Term k v] -> Term k v
forall (t :: * -> *) a.
Foldable t =>
(a -> a -> Ordering) -> t a -> a
maximumBy (MonomialOrder v
cmp MonomialOrder v
-> (Term k v -> Monomial v) -> Term k v -> Term k v -> Ordering
forall b c a. (b -> b -> c) -> (a -> b) -> a -> a -> c
`on` Term k v -> Monomial v
forall a b. (a, b) -> b
snd) [Term k v]
ms

-- | leading coefficient with respect to a given monomial order
lc :: (Num k) => MonomialOrder v -> Polynomial k v -> k
lc :: MonomialOrder v -> Polynomial k v -> k
lc MonomialOrder v
cmp = (k, Monomial v) -> k
forall a b. (a, b) -> a
fst ((k, Monomial v) -> k)
-> (Polynomial k v -> (k, Monomial v)) -> Polynomial k v -> k
forall b c a. (b -> c) -> (a -> b) -> a -> c
. MonomialOrder v -> Polynomial k v -> (k, Monomial v)
forall k v. Num k => MonomialOrder v -> Polynomial k v -> Term k v
lt MonomialOrder v
cmp

-- | leading monomial with respect to a given monomial order
lm :: (Num k) => MonomialOrder v -> Polynomial k v -> Monomial v
lm :: MonomialOrder v -> Polynomial k v -> Monomial v
lm MonomialOrder v
cmp = (k, Monomial v) -> Monomial v
forall a b. (a, b) -> b
snd ((k, Monomial v) -> Monomial v)
-> (Polynomial k v -> (k, Monomial v))
-> Polynomial k v
-> Monomial v
forall b c a. (b -> c) -> (a -> b) -> a -> c
. MonomialOrder v -> Polynomial k v -> (k, Monomial v)
forall k v. Num k => MonomialOrder v -> Polynomial k v -> Term k v
lt MonomialOrder v
cmp

-- | Look up a coefficient for a given monomial.
-- If no such term exists, it returns @0@.
coeff :: (Num k, Ord v) => Monomial v -> Polynomial k v -> k
coeff :: Monomial v -> Polynomial k v -> k
coeff Monomial v
xs (Polynomial Map (Monomial v) k
m) = k -> Monomial v -> Map (Monomial v) k -> k
forall k a. Ord k => a -> k -> Map k a -> a
Map.findWithDefault k
0 Monomial v
xs Map (Monomial v) k
m

-- | Look up a coefficient for a given monomial.
-- If no such term exists, it returns @Nothing@.
lookupCoeff :: Ord v => Monomial v -> Polynomial k v -> Maybe k
lookupCoeff :: Monomial v -> Polynomial k v -> Maybe k
lookupCoeff Monomial v
xs (Polynomial Map (Monomial v) k
m) = Monomial v -> Map (Monomial v) k -> Maybe k
forall k a. Ord k => k -> Map k a -> Maybe a
Map.lookup Monomial v
xs Map (Monomial v) k
m

contI :: (Integral r, Ord v) => Polynomial r v -> r
contI :: Polynomial r v -> r
contI Polynomial r v
0 = r
1
contI Polynomial r v
p = (r -> r -> r) -> [r] -> r
forall a. (a -> a -> a) -> [a] -> a
foldl1' r -> r -> r
forall a. Integral a => a -> a -> a
Prelude.gcd [r -> r
forall a. Num a => a -> a
abs r
c | (r
c,Monomial v
_) <- Polynomial r v -> [(r, Monomial v)]
forall k v. Polynomial k v -> [Term k v]
terms Polynomial r v
p]

ppI :: (Integral r, Ord v) => Polynomial r v -> Polynomial r v
ppI :: Polynomial r v -> Polynomial r v
ppI Polynomial r v
p = (r -> r) -> Polynomial r v -> Polynomial r v
forall k1 k v.
(Eq k1, Num k1) =>
(k -> k1) -> Polynomial k v -> Polynomial k1 v
mapCoeff r -> r
f Polynomial r v
p
  where
    c :: r
c = Polynomial r v -> r
forall r v. (Integral r, Ord v) => Polynomial r v -> r
contI Polynomial r v
p
    f :: r -> r
f r
x = Bool -> r -> r
forall a. (?callStack::CallStack) => Bool -> a -> a
assert (r
x r -> r -> r
forall a. Integral a => a -> a -> a
`Prelude.mod` r
c r -> r -> Bool
forall a. Eq a => a -> a -> Bool
== r
0) (r -> r) -> r -> r
forall a b. (a -> b) -> a -> b
$ r
x r -> r -> r
forall a. Integral a => a -> a -> a
`Prelude.div` r
c

class ContPP k where
  type PPCoeff k

  -- | Content of a polynomial
  cont :: (Ord v) => Polynomial k v -> k
  -- constructive-algebra-0.3.0 では cont 0 は error になる

  -- | Primitive part of a polynomial
  pp :: (Ord v) => Polynomial k v -> Polynomial (PPCoeff k) v

instance ContPP Integer where
  type PPCoeff Integer = Integer
  cont :: Polynomial Integer v -> Integer
cont = Polynomial Integer v -> Integer
forall r v. (Integral r, Ord v) => Polynomial r v -> r
contI
  pp :: Polynomial Integer v -> Polynomial (PPCoeff Integer) v
pp   = Polynomial Integer v -> Polynomial (PPCoeff Integer) v
forall r v. (Integral r, Ord v) => Polynomial r v -> Polynomial r v
ppI

instance Integral r => ContPP (Ratio r) where
  {-# SPECIALIZE instance ContPP (Ratio Integer) #-}

  type PPCoeff (Ratio r) = r

  cont :: Polynomial (Ratio r) v -> Ratio r
cont Polynomial (Ratio r) v
0 = Ratio r
1
  cont Polynomial (Ratio r) v
p = (r -> r -> r) -> [r] -> r
forall a. (a -> a -> a) -> [a] -> a
foldl1' r -> r -> r
forall a. Integral a => a -> a -> a
Prelude.gcd [r]
ns r -> r -> Ratio r
forall a. Integral a => a -> a -> Ratio a
% (r -> r -> r) -> r -> [r] -> r
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl' r -> r -> r
forall a. Integral a => a -> a -> a
Prelude.lcm r
1 [r]
ds
    where
      ns :: [r]
ns = [r -> r
forall a. Num a => a -> a
abs (Ratio r -> r
forall a. Ratio a -> a
numerator Ratio r
c) | (Ratio r
c,Monomial v
_) <- Polynomial (Ratio r) v -> [(Ratio r, Monomial v)]
forall k v. Polynomial k v -> [Term k v]
terms Polynomial (Ratio r) v
p]
      ds :: [r]
ds = [Ratio r -> r
forall a. Ratio a -> a
denominator Ratio r
c     | (Ratio r
c,Monomial v
_) <- Polynomial (Ratio r) v -> [(Ratio r, Monomial v)]
forall k v. Polynomial k v -> [Term k v]
terms Polynomial (Ratio r) v
p]

  pp :: Polynomial (Ratio r) v -> Polynomial (PPCoeff (Ratio r)) v
pp Polynomial (Ratio r) v
p = (Ratio r -> r) -> Polynomial (Ratio r) v -> Polynomial r v
forall k1 k v.
(Eq k1, Num k1) =>
(k -> k1) -> Polynomial k v -> Polynomial k1 v
mapCoeff (Ratio r -> r
forall a. Ratio a -> a
numerator (Ratio r -> r) -> (Ratio r -> Ratio r) -> Ratio r -> r
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Ratio r -> Ratio r -> Ratio r
forall a. Fractional a => a -> a -> a
/ Ratio r
c)) Polynomial (Ratio r) v
p
    where
      c :: Ratio r
c = Polynomial (Ratio r) v -> Ratio r
forall k v. (ContPP k, Ord v) => Polynomial k v -> k
cont Polynomial (Ratio r) v
p

-- | a polynomial is called primitive if the greatest common divisor of its coefficients is 1.
isPrimitive :: (Eq k, Num k, ContPP k, Ord v) => Polynomial k v -> Bool
isPrimitive :: Polynomial k v -> Bool
isPrimitive Polynomial k v
p = Polynomial k v -> Bool
forall k v. Polynomial k v -> Bool
isZero Polynomial k v
p Bool -> Bool -> Bool
|| Polynomial k v -> k
forall k v. (ContPP k, Ord v) => Polynomial k v -> k
cont Polynomial k v
p k -> k -> Bool
forall a. Eq a => a -> a -> Bool
== k
1

-- | Formal derivative of polynomials
deriv :: (Eq k, Num k, Ord v) => Polynomial k v -> v -> Polynomial k v
deriv :: Polynomial k v -> v -> Polynomial k v
deriv Polynomial k v
p v
x = [Polynomial k v] -> Polynomial k v
forall (f :: * -> *) v. (Foldable f, AdditiveGroup v) => f v -> v
sumV [Term k v -> Polynomial k v
forall k v. (Eq k, Num k) => Term k v -> Polynomial k v
fromTerm (Term k v -> v -> Term k v
forall k v. (Num k, Ord v) => Term k v -> v -> Term k v
tderiv Term k v
m v
x) | Term k v
m <- Polynomial k v -> [Term k v]
forall k v. Polynomial k v -> [Term k v]
terms Polynomial k v
p]

-- | Formal integral of polynomials
integral :: (Eq k, Fractional k, Ord v) => Polynomial k v -> v -> Polynomial k v
integral :: Polynomial k v -> v -> Polynomial k v
integral Polynomial k v
p v
x = [Polynomial k v] -> Polynomial k v
forall (f :: * -> *) v. (Foldable f, AdditiveGroup v) => f v -> v
sumV [Term k v -> Polynomial k v
forall k v. (Eq k, Num k) => Term k v -> Polynomial k v
fromTerm (Term k v -> v -> Term k v
forall k v. (Fractional k, Ord v) => Term k v -> v -> Term k v
tintegral Term k v
m v
x) | Term k v
m <- Polynomial k v -> [Term k v]
forall k v. Polynomial k v -> [Term k v]
terms Polynomial k v
p]

-- | Evaluation
eval :: (Num k) => (v -> k) -> Polynomial k v -> k
eval :: (v -> k) -> Polynomial k v -> k
eval v -> k
env Polynomial k v
p = [k] -> k
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [k
c k -> k -> k
forall a. Num a => a -> a -> a
* [k] -> k
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [(v -> k
env v
x) k -> Integer -> k
forall a b. (Num a, Integral b) => a -> b -> a
^ Integer
e | (v
x,Integer
e) <- Monomial v -> [(v, Integer)]
forall v. Monomial v -> [(v, Integer)]
mindices Monomial v
xs] | (k
c,Monomial v
xs) <- Polynomial k v -> [(k, Monomial v)]
forall k v. Polynomial k v -> [Term k v]
terms Polynomial k v
p]

-- | Substitution or bind
subst
  :: (Eq k, Num k, Ord v2)
  => Polynomial k v1 -> (v1 -> Polynomial k v2) -> Polynomial k v2
subst :: Polynomial k v1 -> (v1 -> Polynomial k v2) -> Polynomial k v2
subst Polynomial k v1
p v1 -> Polynomial k v2
s =
  [Polynomial k v2] -> Polynomial k v2
forall (f :: * -> *) v. (Foldable f, AdditiveGroup v) => f v -> v
sumV [k -> Polynomial k v2
forall k v. (Eq k, Num k) => k -> Polynomial k v
constant k
c Polynomial k v2 -> Polynomial k v2 -> Polynomial k v2
forall a. Num a => a -> a -> a
* [Polynomial k v2] -> Polynomial k v2
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [(v1 -> Polynomial k v2
s v1
x)Polynomial k v2 -> Integer -> Polynomial k v2
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
e | (v1
x,Integer
e) <- Monomial v1 -> [(v1, Integer)]
forall v. Monomial v -> [(v, Integer)]
mindices Monomial v1
xs] | (k
c, Monomial v1
xs) <- Polynomial k v1 -> [(k, Monomial v1)]
forall k v. Polynomial k v -> [Term k v]
terms Polynomial k v1
p]

-- | Transform polynomial coefficients.
mapCoeff :: (Eq k1, Num k1) => (k -> k1) -> Polynomial k v -> Polynomial k1 v
mapCoeff :: (k -> k1) -> Polynomial k v -> Polynomial k1 v
mapCoeff k -> k1
f (Polynomial Map (Monomial v) k
m) = Map (Monomial v) k1 -> Polynomial k1 v
forall r v. Map (Monomial v) r -> Polynomial r v
Polynomial (Map (Monomial v) k1 -> Polynomial k1 v)
-> Map (Monomial v) k1 -> Polynomial k1 v
forall a b. (a -> b) -> a -> b
$ (k -> Maybe k1) -> Map (Monomial v) k -> Map (Monomial v) k1
forall a b k. (a -> Maybe b) -> Map k a -> Map k b
Map.mapMaybe k -> Maybe k1
g Map (Monomial v) k
m
  where
    g :: k -> Maybe k1
g k
x = if k1
y k1 -> k1 -> Bool
forall a. Eq a => a -> a -> Bool
== k1
0 then Maybe k1
forall a. Maybe a
Nothing else k1 -> Maybe k1
forall a. a -> Maybe a
Just k1
y
      where
        y :: k1
y = k -> k1
f k
x

-- | Transform a polynomial into a monic polynomial with respect to the given monomial order.
toMonic :: (Eq r, Fractional r) => MonomialOrder v -> Polynomial r v -> Polynomial r v
toMonic :: MonomialOrder v -> Polynomial r v -> Polynomial r v
toMonic MonomialOrder v
cmp Polynomial r v
p
  | r
c r -> r -> Bool
forall a. Eq a => a -> a -> Bool
== r
0 Bool -> Bool -> Bool
|| r
c r -> r -> Bool
forall a. Eq a => a -> a -> Bool
== r
1 = Polynomial r v
p
  | Bool
otherwise = (r -> r) -> Polynomial r v -> Polynomial r v
forall k1 k v.
(Eq k1, Num k1) =>
(k -> k1) -> Polynomial k v -> Polynomial k1 v
mapCoeff (r -> r -> r
forall a. Fractional a => a -> a -> a
/r
c) Polynomial r v
p
  where
    c :: r
c = MonomialOrder v -> Polynomial r v -> r
forall k v. Num k => MonomialOrder v -> Polynomial k v -> k
lc MonomialOrder v
cmp Polynomial r v
p

-- | Convert /K[x,x1,x2,…]/ into /K[x1,x2,…][x]/.
toUPolynomialOf :: (Ord k, Num k, Ord v) => Polynomial k v -> v -> UPolynomial (Polynomial k v)
toUPolynomialOf :: Polynomial k v -> v -> UPolynomial (Polynomial k v)
toUPolynomialOf Polynomial k v
p v
v = [Term (Polynomial k v) X] -> UPolynomial (Polynomial k v)
forall k v. (Eq k, Num k, Ord v) => [Term k v] -> Polynomial k v
fromTerms ([Term (Polynomial k v) X] -> UPolynomial (Polynomial k v))
-> [Term (Polynomial k v) X] -> UPolynomial (Polynomial k v)
forall a b. (a -> b) -> a -> b
$ do
  (k
c,Monomial v
mm) <- Polynomial k v -> [(k, Monomial v)]
forall k v. Polynomial k v -> [Term k v]
terms Polynomial k v
p
  let m :: Map v Integer
m = Monomial v -> Map v Integer
forall v. Monomial v -> Map v Integer
mindicesMap Monomial v
mm
  Term (Polynomial k v) X -> [Term (Polynomial k v) X]
forall (m :: * -> *) a. Monad m => a -> m a
return ( [(k, Monomial v)] -> Polynomial k v
forall k v. (Eq k, Num k, Ord v) => [Term k v] -> Polynomial k v
fromTerms [(k
c, Map v Integer -> Monomial v
forall v. Ord v => Map v Integer -> Monomial v
mfromIndicesMap (v -> Map v Integer -> Map v Integer
forall k a. Ord k => k -> Map k a -> Map k a
Map.delete v
v Map v Integer
m))]
         , X -> Monomial X
forall a v. Var a v => v -> a
var X
X Monomial X -> Integer -> Monomial X
forall v. Ord v => Monomial v -> Integer -> Monomial v
`mpow` Integer -> v -> Map v Integer -> Integer
forall k a. Ord k => a -> k -> Map k a -> a
Map.findWithDefault Integer
0 v
v Map v Integer
m
         )

-- | Multivariate division algorithm
--
-- @divModMP cmp f [g1,g2,..]@ returns a pair @([q1,q2,…],r)@ such that
--
--   * @f = g1*q1 + g2*q2 + .. + r@ and
--
--   * @g1, g2, ..@ do not divide @r@.
divModMP
  :: forall k v. (Eq k, Fractional k, Ord v)
  => MonomialOrder v -> Polynomial k v -> [Polynomial k v] -> ([Polynomial k v], Polynomial k v)
divModMP :: MonomialOrder v
-> Polynomial k v
-> [Polynomial k v]
-> ([Polynomial k v], Polynomial k v)
divModMP MonomialOrder v
cmp Polynomial k v
p [Polynomial k v]
fs = IntMap (Polynomial k v)
-> [Term k v] -> ([Polynomial k v], Polynomial k v)
go IntMap (Polynomial k v)
forall a. IntMap a
IntMap.empty (Polynomial k v -> [Term k v]
terms' Polynomial k v
p)
  where
    terms' :: Polynomial k v -> [Term k v]
    terms' :: Polynomial k v -> [Term k v]
terms' Polynomial k v
g = (Term k v -> Term k v -> Ordering) -> [Term k v] -> [Term k v]
forall a. (a -> a -> Ordering) -> [a] -> [a]
sortBy (MonomialOrder v -> MonomialOrder v
forall a b c. (a -> b -> c) -> b -> a -> c
flip MonomialOrder v
cmp MonomialOrder v
-> (Term k v -> Monomial v) -> Term k v -> Term k v -> Ordering
forall b c a. (b -> b -> c) -> (a -> b) -> a -> a -> c
`on` Term k v -> Monomial v
forall a b. (a, b) -> b
snd) (Polynomial k v -> [Term k v]
forall k v. Polynomial k v -> [Term k v]
terms Polynomial k v
g)

    merge :: [Term k v] -> [Term k v] -> [Term k v]
    merge :: [Term k v] -> [Term k v] -> [Term k v]
merge [] [Term k v]
ys = [Term k v]
ys
    merge [Term k v]
xs [] = [Term k v]
xs
    merge xxs :: [Term k v]
xxs@(Term k v
x:[Term k v]
xs) yys :: [Term k v]
yys@(Term k v
y:[Term k v]
ys) =
      case MonomialOrder v
cmp (Term k v -> Monomial v
forall a b. (a, b) -> b
snd Term k v
x) (Term k v -> Monomial v
forall a b. (a, b) -> b
snd Term k v
y) of
        Ordering
GT -> Term k v
x Term k v -> [Term k v] -> [Term k v]
forall a. a -> [a] -> [a]
: [Term k v] -> [Term k v] -> [Term k v]
merge [Term k v]
xs [Term k v]
yys
        Ordering
LT -> Term k v
y Term k v -> [Term k v] -> [Term k v]
forall a. a -> [a] -> [a]
: [Term k v] -> [Term k v] -> [Term k v]
merge [Term k v]
xxs [Term k v]
ys
        Ordering
EQ ->
          case Term k v -> k
forall a b. (a, b) -> a
fst Term k v
x k -> k -> k
forall a. Num a => a -> a -> a
+ Term k v -> k
forall a b. (a, b) -> a
fst Term k v
y of
            k
0 -> [Term k v] -> [Term k v] -> [Term k v]
merge [Term k v]
xs [Term k v]
ys
            k
c -> (k
c, Term k v -> Monomial v
forall a b. (a, b) -> b
snd Term k v
x) Term k v -> [Term k v] -> [Term k v]
forall a. a -> [a] -> [a]
: [Term k v] -> [Term k v] -> [Term k v]
merge [Term k v]
xs [Term k v]
ys

    ls :: [(Int, (Term k v, [Term k v]))]
ls = [Int]
-> [(Term k v, [Term k v])] -> [(Int, (Term k v, [Term k v]))]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0..] [(MonomialOrder v -> Polynomial k v -> Term k v
forall k v. Num k => MonomialOrder v -> Polynomial k v -> Term k v
lt MonomialOrder v
cmp Polynomial k v
f, Polynomial k v -> [Term k v]
terms' Polynomial k v
f) | Polynomial k v
f <- [Polynomial k v]
fs]

    go :: IntMap (Polynomial k v) -> [Term k v] -> ([Polynomial k v], Polynomial k v)
    go :: IntMap (Polynomial k v)
-> [Term k v] -> ([Polynomial k v], Polynomial k v)
go IntMap (Polynomial k v)
qs [Term k v]
g =
      case [(Int, Polynomial k v, [Term k v])]
xs of
        [] -> ([Polynomial k v -> Int -> IntMap (Polynomial k v) -> Polynomial k v
forall a. a -> Int -> IntMap a -> a
IntMap.findWithDefault Polynomial k v
0 Int
i IntMap (Polynomial k v)
qs | Int
i <- [Int
0 .. [Polynomial k v] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Polynomial k v]
fs Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]], [Term k v] -> Polynomial k v
forall k v. (Eq k, Num k, Ord v) => [Term k v] -> Polynomial k v
fromTerms [Term k v]
g)
        (Int
i, Polynomial k v
b, [Term k v]
g') : [(Int, Polynomial k v, [Term k v])]
_ -> IntMap (Polynomial k v)
-> [Term k v] -> ([Polynomial k v], Polynomial k v)
go ((Polynomial k v -> Polynomial k v -> Polynomial k v)
-> Int
-> Polynomial k v
-> IntMap (Polynomial k v)
-> IntMap (Polynomial k v)
forall a. (a -> a -> a) -> Int -> a -> IntMap a -> IntMap a
IntMap.insertWith Polynomial k v -> Polynomial k v -> Polynomial k v
forall a. Num a => a -> a -> a
(+) Int
i Polynomial k v
b IntMap (Polynomial k v)
qs) [Term k v]
g'
      where
        xs :: [(Int, Polynomial k v, [Term k v])]
xs = do
          (Int
i,(Term k v
a,[Term k v]
f)) <- [(Int, (Term k v, [Term k v]))]
ls
          Term k v
h <- [Term k v]
g
          Bool -> [()]
forall (f :: * -> *). Alternative f => Bool -> f ()
guard (Bool -> [()]) -> Bool -> [()]
forall a b. (a -> b) -> a -> b
$ Term k v
a Term k v -> Term k v -> Bool
forall v k. Ord v => Term k v -> Term k v -> Bool
`tdivides` Term k v
h
          let b :: Term k v
b = Term k v -> Term k v -> Term k v
forall k v.
(Fractional k, Ord v) =>
Term k v -> Term k v -> Term k v
tdiv Term k v
h Term k v
a
          (Int, Polynomial k v, [Term k v])
-> [(Int, Polynomial k v, [Term k v])]
forall (m :: * -> *) a. Monad m => a -> m a
return (Int
i, Term k v -> Polynomial k v
forall k v. (Eq k, Num k) => Term k v -> Polynomial k v
fromTerm Term k v
b, [Term k v] -> [Term k v] -> [Term k v]
merge [Term k v]
g [(k -> Term k v -> Term k v
forall k v. Num k => k -> Term k v -> Term k v
tscale (-k
1) Term k v
b Term k v -> Term k v -> Term k v
forall k v. (Num k, Ord v) => Term k v -> Term k v -> Term k v
`tmult` Term k v
m) | Term k v
m <- [Term k v]
f])

-- | Multivariate division algorithm
--
-- @reduce cmp f gs = snd (divModMP cmp f gs)@
reduce
  :: (Eq k, Fractional k, Ord v)
  => MonomialOrder v -> Polynomial k v -> [Polynomial k v] -> Polynomial k v
reduce :: MonomialOrder v
-> Polynomial k v -> [Polynomial k v] -> Polynomial k v
reduce MonomialOrder v
cmp Polynomial k v
p [Polynomial k v]
fs = Polynomial k v -> Polynomial k v
go Polynomial k v
p
  where
    ls :: [(Term k v, Polynomial k v)]
ls = [(MonomialOrder v -> Polynomial k v -> Term k v
forall k v. Num k => MonomialOrder v -> Polynomial k v -> Term k v
lt MonomialOrder v
cmp Polynomial k v
f, Polynomial k v
f) | Polynomial k v
f <- [Polynomial k v]
fs]
    go :: Polynomial k v -> Polynomial k v
go Polynomial k v
g = if [Polynomial k v] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Polynomial k v]
xs then Polynomial k v
g else Polynomial k v -> Polynomial k v
go ([Polynomial k v] -> Polynomial k v
forall a. [a] -> a
head [Polynomial k v]
xs)
      where
        ms :: [Term k v]
ms = (Term k v -> Term k v -> Ordering) -> [Term k v] -> [Term k v]
forall a. (a -> a -> Ordering) -> [a] -> [a]
sortBy (MonomialOrder v -> MonomialOrder v
forall a b c. (a -> b -> c) -> b -> a -> c
flip MonomialOrder v
cmp MonomialOrder v
-> (Term k v -> Monomial v) -> Term k v -> Term k v -> Ordering
forall b c a. (b -> b -> c) -> (a -> b) -> a -> a -> c
`on` Term k v -> Monomial v
forall a b. (a, b) -> b
snd) (Polynomial k v -> [Term k v]
forall k v. Polynomial k v -> [Term k v]
terms Polynomial k v
g)
        xs :: [Polynomial k v]
xs = do
          (Term k v
a,Polynomial k v
f) <- [(Term k v, Polynomial k v)]
ls
          Term k v
h <- [Term k v]
ms
          Bool -> [()]
forall (f :: * -> *). Alternative f => Bool -> f ()
guard (Bool -> [()]) -> Bool -> [()]
forall a b. (a -> b) -> a -> b
$ Term k v
a Term k v -> Term k v -> Bool
forall v k. Ord v => Term k v -> Term k v -> Bool
`tdivides` Term k v
h
          Polynomial k v -> [Polynomial k v]
forall (m :: * -> *) a. Monad m => a -> m a
return (Polynomial k v
g Polynomial k v -> Polynomial k v -> Polynomial k v
forall a. Num a => a -> a -> a
- Term k v -> Polynomial k v
forall k v. (Eq k, Num k) => Term k v -> Polynomial k v
fromTerm (Term k v -> Term k v -> Term k v
forall k v.
(Fractional k, Ord v) =>
Term k v -> Term k v -> Term k v
tdiv Term k v
h Term k v
a) Polynomial k v -> Polynomial k v -> Polynomial k v
forall a. Num a => a -> a -> a
* Polynomial k v
f)

-- | Prime factorization of polynomials
class Factor a where
  -- | factor a polynomial @p@ into @p1 ^ n1 * p2 ^ n2 * ..@ and
  -- return a list @[(p1,n1), (p2,n2), ..]@.
  factor :: a -> [(a, Integer)]

-- | Square-free factorization of polynomials
class SQFree a where
  -- | factor a polynomial @p@ into @p1 ^ n1 * p2 ^ n2 * ..@ and
  -- return a list @[(p1,n1), (p2,n2), ..]@.
  sqfree :: a -> [(a, Integer)]

-- | Eisenstein's irreducibility criterion.
--
-- Given a polynomial with integer coefficients @P[x] = an x^n + .. + a1*x + a0@,
-- it is irreducible over rational numbers if there exists a prime number @p@
-- satisfying the following conditions:
--
-- * @p@ divides ai for @i /= n@,
--
-- * @p@ does not divides @an@, and
--
-- * @p^2@ does not divides @a0@.
--
eisensteinsCriterion
  :: (ContPP k, PPCoeff k ~ Integer)
  => UPolynomial k -> Bool
eisensteinsCriterion :: UPolynomial k -> Bool
eisensteinsCriterion UPolynomial k
p
  | UPolynomial k -> Integer
forall t. Degree t => t -> Integer
deg UPolynomial k
p Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
1 = Bool
True
  | Bool
otherwise  = UPolynomial Integer -> Bool
eisensteinsCriterion' (UPolynomial k -> Polynomial (PPCoeff k) X
forall k v.
(ContPP k, Ord v) =>
Polynomial k v -> Polynomial (PPCoeff k) v
pp UPolynomial k
p)

eisensteinsCriterion' :: UPolynomial Integer -> Bool
eisensteinsCriterion' :: UPolynomial Integer -> Bool
eisensteinsCriterion' UPolynomial Integer
p = [Bool] -> Bool
forall (t :: * -> *). Foldable t => t Bool -> Bool
or [Integer -> Bool
criterion Integer
prime | Integer
prime <- ([Integer] -> Integer) -> [[Integer]] -> [Integer]
forall a b. (a -> b) -> [a] -> [b]
map [Integer] -> Integer
forall a. [a] -> a
head ([[Integer]] -> [Integer]) -> [[Integer]] -> [Integer]
forall a b. (a -> b) -> a -> b
$ [Integer] -> [[Integer]]
forall a. Eq a => [a] -> [[a]]
group ([Integer] -> [[Integer]]) -> [Integer] -> [[Integer]]
forall a b. (a -> b) -> a -> b
$ Integer -> [Integer]
forall int. Integral int => int -> [int]
primeFactors Integer
c]
  where
    Just ((Monomial X
_,Integer
an), Map (Monomial X) Integer
ts) = Map (Monomial X) Integer
-> Maybe ((Monomial X, Integer), Map (Monomial X) Integer)
forall k a. Map k a -> Maybe ((k, a), Map k a)
Map.maxViewWithKey (UPolynomial Integer -> Map (Monomial X) Integer
forall r v. Polynomial r v -> Map (Monomial v) r
coeffMap UPolynomial Integer
p)
    a0 :: Integer
a0 = Monomial X -> UPolynomial Integer -> Integer
forall k v. (Num k, Ord v) => Monomial v -> Polynomial k v -> k
coeff Monomial X
forall v. Monomial v
mone UPolynomial Integer
p
    c :: Integer
c = (Integer -> Integer -> Integer) -> Integer -> [Integer] -> Integer
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl' Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
Prelude.gcd Integer
0 [Integer
ai | (Monomial X
_,Integer
ai) <- Map (Monomial X) Integer -> [(Monomial X, Integer)]
forall k a. Map k a -> [(k, a)]
Map.toList Map (Monomial X) Integer
ts]
    criterion :: Integer -> Bool
criterion Integer
prime = Integer
an Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`Prelude.mod` Integer
prime Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
/= Integer
0 Bool -> Bool -> Bool
&& Integer
a0 Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`Prelude.mod` (Integer
primeInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
prime) Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
/= Integer
0

{--------------------------------------------------------------------
  Pretty printing
--------------------------------------------------------------------}

-- | Options for pretty printing polynomials
--
-- The default value can be obtained by 'def'.
data PrintOptions k v
  = PrintOptions
  { PrintOptions k v -> PrettyLevel -> Rational -> v -> Doc
pOptPrintVar        :: PrettyLevel -> Rational -> v -> Doc
  , PrintOptions k v -> PrettyLevel -> Rational -> k -> Doc
pOptPrintCoeff      :: PrettyLevel -> Rational -> k -> Doc
  , PrintOptions k v -> k -> Bool
pOptIsNegativeCoeff :: k -> Bool
  , PrintOptions k v -> MonomialOrder v
pOptMonomialOrder   :: MonomialOrder v
  }

instance (PrettyCoeff k, PrettyVar v, Ord v) => Default (PrintOptions k v) where
  def :: PrintOptions k v
def =
    PrintOptions :: forall k v.
(PrettyLevel -> Rational -> v -> Doc)
-> (PrettyLevel -> Rational -> k -> Doc)
-> (k -> Bool)
-> MonomialOrder v
-> PrintOptions k v
PrintOptions
    { pOptPrintVar :: PrettyLevel -> Rational -> v -> Doc
pOptPrintVar        = PrettyLevel -> Rational -> v -> Doc
forall a. PrettyVar a => PrettyLevel -> Rational -> a -> Doc
pPrintVar
    , pOptPrintCoeff :: PrettyLevel -> Rational -> k -> Doc
pOptPrintCoeff      = PrettyLevel -> Rational -> k -> Doc
forall a. PrettyCoeff a => PrettyLevel -> Rational -> a -> Doc
pPrintCoeff
    , pOptIsNegativeCoeff :: k -> Bool
pOptIsNegativeCoeff = k -> Bool
forall a. PrettyCoeff a => a -> Bool
isNegativeCoeff
    , pOptMonomialOrder :: MonomialOrder v
pOptMonomialOrder   = MonomialOrder v
forall v. Ord v => MonomialOrder v
grlex
    }

instance (Ord k, Num k, Ord v, PrettyCoeff k, PrettyVar v) => Pretty (Polynomial k v) where
  pPrintPrec :: PrettyLevel -> Rational -> Polynomial k v -> Doc
pPrintPrec = PrintOptions k v
-> PrettyLevel -> Rational -> Polynomial k v -> Doc
forall k v.
(Ord k, Num k) =>
PrintOptions k v
-> PrettyLevel -> Rational -> Polynomial k v -> Doc
prettyPrint PrintOptions k v
forall a. Default a => a
def

prettyPrint
  :: (Ord k, Num k)
  => PrintOptions k v
  -> PrettyLevel -> Rational -> Polynomial k v -> Doc
prettyPrint :: PrintOptions k v
-> PrettyLevel -> Rational -> Polynomial k v -> Doc
prettyPrint PrintOptions k v
opt PrettyLevel
lv Rational
prec Polynomial k v
p =
    case ((k, Monomial v) -> (k, Monomial v) -> Ordering)
-> [(k, Monomial v)] -> [(k, Monomial v)]
forall a. (a -> a -> Ordering) -> [a] -> [a]
sortBy ((Monomial v -> Monomial v -> Ordering)
-> Monomial v -> Monomial v -> Ordering
forall a b c. (a -> b -> c) -> b -> a -> c
flip (PrintOptions k v -> Monomial v -> Monomial v -> Ordering
forall k v. PrintOptions k v -> MonomialOrder v
pOptMonomialOrder PrintOptions k v
opt) (Monomial v -> Monomial v -> Ordering)
-> ((k, Monomial v) -> Monomial v)
-> (k, Monomial v)
-> (k, Monomial v)
-> Ordering
forall b c a. (b -> b -> c) -> (a -> b) -> a -> a -> c
`on` (k, Monomial v) -> Monomial v
forall a b. (a, b) -> b
snd) ([(k, Monomial v)] -> [(k, Monomial v)])
-> [(k, Monomial v)] -> [(k, Monomial v)]
forall a b. (a -> b) -> a -> b
$ Polynomial k v -> [(k, Monomial v)]
forall k v. Polynomial k v -> [Term k v]
terms Polynomial k v
p of
      [] -> Int -> Doc
PP.int Int
0
      [(k, Monomial v)
t] -> Rational -> (k, Monomial v) -> Doc
pLeadingTerm Rational
prec (k, Monomial v)
t
      (k, Monomial v)
t:[(k, Monomial v)]
ts ->
        Bool -> Doc -> Doc
maybeParens (Rational
prec Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
> Rational
addPrec) (Doc -> Doc) -> Doc -> Doc
forall a b. (a -> b) -> a -> b
$
          [Doc] -> Doc
PP.hcat (Rational -> (k, Monomial v) -> Doc
pLeadingTerm Rational
addPrec (k, Monomial v)
t Doc -> [Doc] -> [Doc]
forall a. a -> [a] -> [a]
: ((k, Monomial v) -> Doc) -> [(k, Monomial v)] -> [Doc]
forall a b. (a -> b) -> [a] -> [b]
map (k, Monomial v) -> Doc
pTrailingTerm [(k, Monomial v)]
ts)
    where
      pLeadingTerm :: Rational -> (k, Monomial v) -> Doc
pLeadingTerm Rational
prec' (k
c,Monomial v
xs) =
        if PrintOptions k v -> k -> Bool
forall k v. PrintOptions k v -> k -> Bool
pOptIsNegativeCoeff PrintOptions k v
opt k
c
        then Bool -> Doc -> Doc
maybeParens (Rational
prec' Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
> Rational
addPrec) (Doc -> Doc) -> Doc -> Doc
forall a b. (a -> b) -> a -> b
$
               Char -> Doc
PP.char Char
'-' Doc -> Doc -> Doc
forall a. Semigroup a => a -> a -> a
<> PrintOptions k v
-> PrettyLevel -> Rational -> (k, Monomial v) -> Doc
forall k v.
(Ord k, Num k) =>
PrintOptions k v -> PrettyLevel -> Rational -> Term k v -> Doc
prettyPrintTerm PrintOptions k v
opt PrettyLevel
lv (Rational
addPrecRational -> Rational -> Rational
forall a. Num a => a -> a -> a
+Rational
1) (-k
c,Monomial v
xs)
        else PrintOptions k v
-> PrettyLevel -> Rational -> (k, Monomial v) -> Doc
forall k v.
(Ord k, Num k) =>
PrintOptions k v -> PrettyLevel -> Rational -> Term k v -> Doc
prettyPrintTerm PrintOptions k v
opt PrettyLevel
lv Rational
prec' (k
c,Monomial v
xs)

      pTrailingTerm :: (k, Monomial v) -> Doc
pTrailingTerm (k
c,Monomial v
xs) =
        if PrintOptions k v -> k -> Bool
forall k v. PrintOptions k v -> k -> Bool
pOptIsNegativeCoeff PrintOptions k v
opt k
c
        then Doc
PP.space Doc -> Doc -> Doc
forall a. Semigroup a => a -> a -> a
<> Char -> Doc
PP.char Char
'-' Doc -> Doc -> Doc
forall a. Semigroup a => a -> a -> a
<> Doc
PP.space Doc -> Doc -> Doc
forall a. Semigroup a => a -> a -> a
<> PrintOptions k v
-> PrettyLevel -> Rational -> (k, Monomial v) -> Doc
forall k v.
(Ord k, Num k) =>
PrintOptions k v -> PrettyLevel -> Rational -> Term k v -> Doc
prettyPrintTerm PrintOptions k v
opt PrettyLevel
lv (Rational
addPrecRational -> Rational -> Rational
forall a. Num a => a -> a -> a
+Rational
1) (-k
c,Monomial v
xs)
        else Doc
PP.space Doc -> Doc -> Doc
forall a. Semigroup a => a -> a -> a
<> Char -> Doc
PP.char Char
'+' Doc -> Doc -> Doc
forall a. Semigroup a => a -> a -> a
<> Doc
PP.space Doc -> Doc -> Doc
forall a. Semigroup a => a -> a -> a
<> PrintOptions k v
-> PrettyLevel -> Rational -> (k, Monomial v) -> Doc
forall k v.
(Ord k, Num k) =>
PrintOptions k v -> PrettyLevel -> Rational -> Term k v -> Doc
prettyPrintTerm PrintOptions k v
opt PrettyLevel
lv (Rational
addPrecRational -> Rational -> Rational
forall a. Num a => a -> a -> a
+Rational
1) (k
c,Monomial v
xs)

prettyPrintTerm
  :: (Ord k, Num k)
  => PrintOptions k v
  -> PrettyLevel -> Rational -> Term k v -> Doc
prettyPrintTerm :: PrintOptions k v -> PrettyLevel -> Rational -> Term k v -> Doc
prettyPrintTerm PrintOptions k v
opt PrettyLevel
lv Rational
prec (k
c,Monomial v
xs)
  | Int
len Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0  = PrintOptions k v -> PrettyLevel -> Rational -> k -> Doc
forall k v. PrintOptions k v -> PrettyLevel -> Rational -> k -> Doc
pOptPrintCoeff PrintOptions k v
opt PrettyLevel
lv (Rational
appPrecRational -> Rational -> Rational
forall a. Num a => a -> a -> a
+Rational
1) k
c
    -- intentionally specify (appPrec+1) to parenthesize any composite expression
  | Int
len Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1 Bool -> Bool -> Bool
&& k
c k -> k -> Bool
forall a. Eq a => a -> a -> Bool
== k
1 = Rational -> (v, Integer) -> Doc
pPow Rational
prec ((v, Integer) -> Doc) -> (v, Integer) -> Doc
forall a b. (a -> b) -> a -> b
$ [(v, Integer)] -> (v, Integer)
forall a. [a] -> a
head (Monomial v -> [(v, Integer)]
forall v. Monomial v -> [(v, Integer)]
mindices Monomial v
xs)
  | Bool
otherwise =
      Bool -> Doc -> Doc
maybeParens (Rational
prec Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
> Rational
mulPrec) (Doc -> Doc) -> Doc -> Doc
forall a b. (a -> b) -> a -> b
$
        [Doc] -> Doc
PP.hcat ([Doc] -> Doc) -> [Doc] -> Doc
forall a b. (a -> b) -> a -> b
$ Doc -> [Doc] -> [Doc]
forall a. a -> [a] -> [a]
intersperse (Char -> Doc
PP.char Char
'*') [Doc]
fs
    where
      len :: Int
len = Map v Integer -> Int
forall k a. Map k a -> Int
Map.size (Map v Integer -> Int) -> Map v Integer -> Int
forall a b. (a -> b) -> a -> b
$ Monomial v -> Map v Integer
forall v. Monomial v -> Map v Integer
mindicesMap Monomial v
xs
      fs :: [Doc]
fs  = [PrintOptions k v -> PrettyLevel -> Rational -> k -> Doc
forall k v. PrintOptions k v -> PrettyLevel -> Rational -> k -> Doc
pOptPrintCoeff PrintOptions k v
opt PrettyLevel
lv (Rational
appPrecRational -> Rational -> Rational
forall a. Num a => a -> a -> a
+Rational
1) k
c | k
c k -> k -> Bool
forall a. Eq a => a -> a -> Bool
/= k
1] [Doc] -> [Doc] -> [Doc]
forall a. [a] -> [a] -> [a]
++ [Rational -> (v, Integer) -> Doc
pPow (Rational
mulPrecRational -> Rational -> Rational
forall a. Num a => a -> a -> a
+Rational
1) (v, Integer)
p | (v, Integer)
p <- Monomial v -> [(v, Integer)]
forall v. Monomial v -> [(v, Integer)]
mindices Monomial v
xs]
      -- intentionally specify (appPrec+1) to parenthesize any composite expression

      pPow :: Rational -> (v, Integer) -> Doc
pPow Rational
prec' (v
x,Integer
1) = PrintOptions k v -> PrettyLevel -> Rational -> v -> Doc
forall k v. PrintOptions k v -> PrettyLevel -> Rational -> v -> Doc
pOptPrintVar PrintOptions k v
opt PrettyLevel
lv Rational
prec' v
x
      pPow Rational
prec' (v
x,Integer
n) =
        Bool -> Doc -> Doc
maybeParens (Rational
prec' Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
> Rational
expPrec) (Doc -> Doc) -> Doc -> Doc
forall a b. (a -> b) -> a -> b
$
          PrintOptions k v -> PrettyLevel -> Rational -> v -> Doc
forall k v. PrintOptions k v -> PrettyLevel -> Rational -> v -> Doc
pOptPrintVar PrintOptions k v
opt PrettyLevel
lv (Rational
expPrecRational -> Rational -> Rational
forall a. Num a => a -> a -> a
+Rational
1) v
x Doc -> Doc -> Doc
forall a. Semigroup a => a -> a -> a
<> Char -> Doc
PP.char Char
'^' Doc -> Doc -> Doc
forall a. Semigroup a => a -> a -> a
<> Integer -> Doc
PP.integer Integer
n

class PrettyCoeff a where
  pPrintCoeff :: PrettyLevel -> Rational -> a -> Doc
  isNegativeCoeff :: a -> Bool
  isNegativeCoeff a
_ = Bool
False

instance PrettyCoeff Integer where
  pPrintCoeff :: PrettyLevel -> Rational -> Integer -> Doc
pPrintCoeff = PrettyLevel -> Rational -> Integer -> Doc
forall a. Pretty a => PrettyLevel -> Rational -> a -> Doc
pPrintPrec
  isNegativeCoeff :: Integer -> Bool
isNegativeCoeff = (Integer
0Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
>)

instance (PrettyCoeff a, Integral a) => PrettyCoeff (Ratio a) where
  pPrintCoeff :: PrettyLevel -> Rational -> Ratio a -> Doc
pPrintCoeff PrettyLevel
lv Rational
p Ratio a
r
    | Ratio a -> a
forall a. Ratio a -> a
denominator Ratio a
r a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
1 = PrettyLevel -> Rational -> a -> Doc
forall a. PrettyCoeff a => PrettyLevel -> Rational -> a -> Doc
pPrintCoeff PrettyLevel
lv Rational
p (Ratio a -> a
forall a. Ratio a -> a
numerator Ratio a
r)
    | Bool
otherwise =
        Bool -> Doc -> Doc
maybeParens (Rational
p Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
> Rational
ratPrec) (Doc -> Doc) -> Doc -> Doc
forall a b. (a -> b) -> a -> b
$
          PrettyLevel -> Rational -> a -> Doc
forall a. PrettyCoeff a => PrettyLevel -> Rational -> a -> Doc
pPrintCoeff PrettyLevel
lv (Rational
ratPrecRational -> Rational -> Rational
forall a. Num a => a -> a -> a
+Rational
1) (Ratio a -> a
forall a. Ratio a -> a
numerator Ratio a
r) Doc -> Doc -> Doc
forall a. Semigroup a => a -> a -> a
<>
          Char -> Doc
PP.char Char
'/' Doc -> Doc -> Doc
forall a. Semigroup a => a -> a -> a
<>
          PrettyLevel -> Rational -> a -> Doc
forall a. PrettyCoeff a => PrettyLevel -> Rational -> a -> Doc
pPrintCoeff PrettyLevel
lv (Rational
ratPrecRational -> Rational -> Rational
forall a. Num a => a -> a -> a
+Rational
1) (Ratio a -> a
forall a. Ratio a -> a
denominator Ratio a
r)
  isNegativeCoeff :: Ratio a -> Bool
isNegativeCoeff Ratio a
x = a -> Bool
forall a. PrettyCoeff a => a -> Bool
isNegativeCoeff (Ratio a -> a
forall a. Ratio a -> a
numerator Ratio a
x)

instance PrettyCoeff (FF.PrimeField a) where
  pPrintCoeff :: PrettyLevel -> Rational -> PrimeField a -> Doc
pPrintCoeff PrettyLevel
lv Rational
p PrimeField a
a = PrettyLevel -> Rational -> Integer -> Doc
forall a. PrettyCoeff a => PrettyLevel -> Rational -> a -> Doc
pPrintCoeff PrettyLevel
lv Rational
p (PrimeField a -> Integer
forall (p :: Nat). PrimeField p -> Integer
FF.toInteger PrimeField a
a)
  isNegativeCoeff :: PrimeField a -> Bool
isNegativeCoeff PrimeField a
_  = Bool
False

instance (Num c, Ord c, PrettyCoeff c, Ord v, PrettyVar v) => PrettyCoeff (Polynomial c v) where
  pPrintCoeff :: PrettyLevel -> Rational -> Polynomial c v -> Doc
pPrintCoeff = PrettyLevel -> Rational -> Polynomial c v -> Doc
forall a. Pretty a => PrettyLevel -> Rational -> a -> Doc
pPrintPrec

class PrettyVar a where
  pPrintVar :: PrettyLevel -> Rational -> a -> Doc

instance PrettyVar Int where
  pPrintVar :: PrettyLevel -> Rational -> Int -> Doc
pPrintVar PrettyLevel
_ Rational
_ Int
n = Char -> Doc
PP.char Char
'x' Doc -> Doc -> Doc
forall a. Semigroup a => a -> a -> a
<> Int -> Doc
PP.int Int
n

instance PrettyVar X where
  pPrintVar :: PrettyLevel -> Rational -> X -> Doc
pPrintVar PrettyLevel
_ Rational
_ X
X = Char -> Doc
PP.char Char
'x'

addPrec, mulPrec, ratPrec, expPrec, appPrec :: Rational
addPrec :: Rational
addPrec = Rational
6 -- Precedence of '+'
mulPrec :: Rational
mulPrec = Rational
7 -- Precedence of '*'
ratPrec :: Rational
ratPrec = Rational
7 -- Precedence of '/'
expPrec :: Rational
expPrec = Rational
8 -- Precedence of '^'
appPrec :: Rational
appPrec = Rational
10 -- Precedence of function application

{--------------------------------------------------------------------
  Univariate polynomials
--------------------------------------------------------------------}

-- | Univariate polynomials over commutative ring r
type UPolynomial r = Polynomial r X

-- | Variable "x"
data X = X
  deriving (X -> X -> Bool
(X -> X -> Bool) -> (X -> X -> Bool) -> Eq X
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: X -> X -> Bool
$c/= :: X -> X -> Bool
== :: X -> X -> Bool
$c== :: X -> X -> Bool
Eq, Eq X
Eq X
-> (X -> X -> Ordering)
-> (X -> X -> Bool)
-> (X -> X -> Bool)
-> (X -> X -> Bool)
-> (X -> X -> Bool)
-> (X -> X -> X)
-> (X -> X -> X)
-> Ord X
X -> X -> Bool
X -> X -> Ordering
X -> X -> X
forall a.
Eq a
-> (a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
min :: X -> X -> X
$cmin :: X -> X -> X
max :: X -> X -> X
$cmax :: X -> X -> X
>= :: X -> X -> Bool
$c>= :: X -> X -> Bool
> :: X -> X -> Bool
$c> :: X -> X -> Bool
<= :: X -> X -> Bool
$c<= :: X -> X -> Bool
< :: X -> X -> Bool
$c< :: X -> X -> Bool
compare :: X -> X -> Ordering
$ccompare :: X -> X -> Ordering
$cp1Ord :: Eq X
Ord, X
X -> X -> Bounded X
forall a. a -> a -> Bounded a
maxBound :: X
$cmaxBound :: X
minBound :: X
$cminBound :: X
Bounded, Int -> X
X -> Int
X -> [X]
X -> X
X -> X -> [X]
X -> X -> X -> [X]
(X -> X)
-> (X -> X)
-> (Int -> X)
-> (X -> Int)
-> (X -> [X])
-> (X -> X -> [X])
-> (X -> X -> [X])
-> (X -> X -> X -> [X])
-> Enum X
forall a.
(a -> a)
-> (a -> a)
-> (Int -> a)
-> (a -> Int)
-> (a -> [a])
-> (a -> a -> [a])
-> (a -> a -> [a])
-> (a -> a -> a -> [a])
-> Enum a
enumFromThenTo :: X -> X -> X -> [X]
$cenumFromThenTo :: X -> X -> X -> [X]
enumFromTo :: X -> X -> [X]
$cenumFromTo :: X -> X -> [X]
enumFromThen :: X -> X -> [X]
$cenumFromThen :: X -> X -> [X]
enumFrom :: X -> [X]
$cenumFrom :: X -> [X]
fromEnum :: X -> Int
$cfromEnum :: X -> Int
toEnum :: Int -> X
$ctoEnum :: Int -> X
pred :: X -> X
$cpred :: X -> X
succ :: X -> X
$csucc :: X -> X
Enum, Int -> X -> ShowS
[X] -> ShowS
X -> String
(Int -> X -> ShowS) -> (X -> String) -> ([X] -> ShowS) -> Show X
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [X] -> ShowS
$cshowList :: [X] -> ShowS
show :: X -> String
$cshow :: X -> String
showsPrec :: Int -> X -> ShowS
$cshowsPrec :: Int -> X -> ShowS
Show, ReadPrec [X]
ReadPrec X
Int -> ReadS X
ReadS [X]
(Int -> ReadS X)
-> ReadS [X] -> ReadPrec X -> ReadPrec [X] -> Read X
forall a.
(Int -> ReadS a)
-> ReadS [a] -> ReadPrec a -> ReadPrec [a] -> Read a
readListPrec :: ReadPrec [X]
$creadListPrec :: ReadPrec [X]
readPrec :: ReadPrec X
$creadPrec :: ReadPrec X
readList :: ReadS [X]
$creadList :: ReadS [X]
readsPrec :: Int -> ReadS X
$creadsPrec :: Int -> ReadS X
Read, Typeable, Typeable X
DataType
Constr
Typeable X
-> (forall (c :: * -> *).
    (forall d b. Data d => c (d -> b) -> d -> c b)
    -> (forall g. g -> c g) -> X -> c X)
-> (forall (c :: * -> *).
    (forall b r. Data b => c (b -> r) -> c r)
    -> (forall r. r -> c r) -> Constr -> c X)
-> (X -> Constr)
-> (X -> DataType)
-> (forall (t :: * -> *) (c :: * -> *).
    Typeable t =>
    (forall d. Data d => c (t d)) -> Maybe (c X))
-> (forall (t :: * -> * -> *) (c :: * -> *).
    Typeable t =>
    (forall d e. (Data d, Data e) => c (t d e)) -> Maybe (c X))
-> ((forall b. Data b => b -> b) -> X -> X)
-> (forall r r'.
    (r -> r' -> r) -> r -> (forall d. Data d => d -> r') -> X -> r)
-> (forall r r'.
    (r' -> r -> r) -> r -> (forall d. Data d => d -> r') -> X -> r)
-> (forall u. (forall d. Data d => d -> u) -> X -> [u])
-> (forall u. Int -> (forall d. Data d => d -> u) -> X -> u)
-> (forall (m :: * -> *).
    Monad m =>
    (forall d. Data d => d -> m d) -> X -> m X)
-> (forall (m :: * -> *).
    MonadPlus m =>
    (forall d. Data d => d -> m d) -> X -> m X)
-> (forall (m :: * -> *).
    MonadPlus m =>
    (forall d. Data d => d -> m d) -> X -> m X)
-> Data X
X -> DataType
X -> Constr
(forall b. Data b => b -> b) -> X -> X
(forall d b. Data d => c (d -> b) -> d -> c b)
-> (forall g. g -> c g) -> X -> c X
(forall b r. Data b => c (b -> r) -> c r)
-> (forall r. r -> c r) -> Constr -> c X
forall a.
Typeable a
-> (forall (c :: * -> *).
    (forall d b. Data d => c (d -> b) -> d -> c b)
    -> (forall g. g -> c g) -> a -> c a)
-> (forall (c :: * -> *).
    (forall b r. Data b => c (b -> r) -> c r)
    -> (forall r. r -> c r) -> Constr -> c a)
-> (a -> Constr)
-> (a -> DataType)
-> (forall (t :: * -> *) (c :: * -> *).
    Typeable t =>
    (forall d. Data d => c (t d)) -> Maybe (c a))
-> (forall (t :: * -> * -> *) (c :: * -> *).
    Typeable t =>
    (forall d e. (Data d, Data e) => c (t d e)) -> Maybe (c a))
-> ((forall b. Data b => b -> b) -> a -> a)
-> (forall r r'.
    (r -> r' -> r) -> r -> (forall d. Data d => d -> r') -> a -> r)
-> (forall r r'.
    (r' -> r -> r) -> r -> (forall d. Data d => d -> r') -> a -> r)
-> (forall u. (forall d. Data d => d -> u) -> a -> [u])
-> (forall u. Int -> (forall d. Data d => d -> u) -> a -> u)
-> (forall (m :: * -> *).
    Monad m =>
    (forall d. Data d => d -> m d) -> a -> m a)
-> (forall (m :: * -> *).
    MonadPlus m =>
    (forall d. Data d => d -> m d) -> a -> m a)
-> (forall (m :: * -> *).
    MonadPlus m =>
    (forall d. Data d => d -> m d) -> a -> m a)
-> Data a
forall u. Int -> (forall d. Data d => d -> u) -> X -> u
forall u. (forall d. Data d => d -> u) -> X -> [u]
forall r r'.
(r -> r' -> r) -> r -> (forall d. Data d => d -> r') -> X -> r
forall r r'.
(r' -> r -> r) -> r -> (forall d. Data d => d -> r') -> X -> r
forall (m :: * -> *).
Monad m =>
(forall d. Data d => d -> m d) -> X -> m X
forall (m :: * -> *).
MonadPlus m =>
(forall d. Data d => d -> m d) -> X -> m X
forall (c :: * -> *).
(forall b r. Data b => c (b -> r) -> c r)
-> (forall r. r -> c r) -> Constr -> c X
forall (c :: * -> *).
(forall d b. Data d => c (d -> b) -> d -> c b)
-> (forall g. g -> c g) -> X -> c X
forall (t :: * -> *) (c :: * -> *).
Typeable t =>
(forall d. Data d => c (t d)) -> Maybe (c X)
forall (t :: * -> * -> *) (c :: * -> *).
Typeable t =>
(forall d e. (Data d, Data e) => c (t d e)) -> Maybe (c X)
$cX :: Constr
$tX :: DataType
gmapMo :: (forall d. Data d => d -> m d) -> X -> m X
$cgmapMo :: forall (m :: * -> *).
MonadPlus m =>
(forall d. Data d => d -> m d) -> X -> m X
gmapMp :: (forall d. Data d => d -> m d) -> X -> m X
$cgmapMp :: forall (m :: * -> *).
MonadPlus m =>
(forall d. Data d => d -> m d) -> X -> m X
gmapM :: (forall d. Data d => d -> m d) -> X -> m X
$cgmapM :: forall (m :: * -> *).
Monad m =>
(forall d. Data d => d -> m d) -> X -> m X
gmapQi :: Int -> (forall d. Data d => d -> u) -> X -> u
$cgmapQi :: forall u. Int -> (forall d. Data d => d -> u) -> X -> u
gmapQ :: (forall d. Data d => d -> u) -> X -> [u]
$cgmapQ :: forall u. (forall d. Data d => d -> u) -> X -> [u]
gmapQr :: (r' -> r -> r) -> r -> (forall d. Data d => d -> r') -> X -> r
$cgmapQr :: forall r r'.
(r' -> r -> r) -> r -> (forall d. Data d => d -> r') -> X -> r
gmapQl :: (r -> r' -> r) -> r -> (forall d. Data d => d -> r') -> X -> r
$cgmapQl :: forall r r'.
(r -> r' -> r) -> r -> (forall d. Data d => d -> r') -> X -> r
gmapT :: (forall b. Data b => b -> b) -> X -> X
$cgmapT :: (forall b. Data b => b -> b) -> X -> X
dataCast2 :: (forall d e. (Data d, Data e) => c (t d e)) -> Maybe (c X)
$cdataCast2 :: forall (t :: * -> * -> *) (c :: * -> *).
Typeable t =>
(forall d e. (Data d, Data e) => c (t d e)) -> Maybe (c X)
dataCast1 :: (forall d. Data d => c (t d)) -> Maybe (c X)
$cdataCast1 :: forall (t :: * -> *) (c :: * -> *).
Typeable t =>
(forall d. Data d => c (t d)) -> Maybe (c X)
dataTypeOf :: X -> DataType
$cdataTypeOf :: X -> DataType
toConstr :: X -> Constr
$ctoConstr :: X -> Constr
gunfold :: (forall b r. Data b => c (b -> r) -> c r)
-> (forall r. r -> c r) -> Constr -> c X
$cgunfold :: forall (c :: * -> *).
(forall b r. Data b => c (b -> r) -> c r)
-> (forall r. r -> c r) -> Constr -> c X
gfoldl :: (forall d b. Data d => c (d -> b) -> d -> c b)
-> (forall g. g -> c g) -> X -> c X
$cgfoldl :: forall (c :: * -> *).
(forall d b. Data d => c (d -> b) -> d -> c b)
-> (forall g. g -> c g) -> X -> c X
$cp1Data :: Typeable X
Data)

instance NFData X where
   rnf :: X -> ()
rnf X
a = X
a X -> () -> ()
`seq` ()

instance Hashable X where
  hashWithSalt :: Int -> X -> Int
hashWithSalt = (X -> Int) -> Int -> X -> Int
forall b a. Hashable b => (a -> b) -> Int -> a -> Int
hashUsing X -> Int
forall a. Enum a => a -> Int
fromEnum

-- | Natural ordering /x^0 < x^1 < x^2 ../ is the unique monomial order for
-- univariate polynomials.
nat :: MonomialOrder X
nat :: MonomialOrder X
nat = Integer -> Integer -> Ordering
forall a. Ord a => a -> a -> Ordering
compare (Integer -> Integer -> Ordering)
-> (Monomial X -> Integer) -> MonomialOrder X
forall b c a. (b -> b -> c) -> (a -> b) -> a -> a -> c
`on` Monomial X -> Integer
forall t. Degree t => t -> Integer
deg

-- | division of univariate polynomials
div :: (Eq k, Fractional k) => UPolynomial k -> UPolynomial k -> UPolynomial k
div :: UPolynomial k -> UPolynomial k -> UPolynomial k
div UPolynomial k
f1 UPolynomial k
f2 = (UPolynomial k, UPolynomial k) -> UPolynomial k
forall a b. (a, b) -> a
fst (UPolynomial k -> UPolynomial k -> (UPolynomial k, UPolynomial k)
forall k.
(Eq k, Fractional k) =>
UPolynomial k -> UPolynomial k -> (UPolynomial k, UPolynomial k)
divMod UPolynomial k
f1 UPolynomial k
f2)

-- | division of univariate polynomials
mod :: (Eq k, Fractional k) => UPolynomial k -> UPolynomial k -> UPolynomial k
mod :: UPolynomial k -> UPolynomial k -> UPolynomial k
mod UPolynomial k
f1 UPolynomial k
f2 = (UPolynomial k, UPolynomial k) -> UPolynomial k
forall a b. (a, b) -> b
snd (UPolynomial k -> UPolynomial k -> (UPolynomial k, UPolynomial k)
forall k.
(Eq k, Fractional k) =>
UPolynomial k -> UPolynomial k -> (UPolynomial k, UPolynomial k)
divMod UPolynomial k
f1 UPolynomial k
f2)

-- | division of univariate polynomials
divMod :: (Eq k, Fractional k) => UPolynomial k -> UPolynomial k -> (UPolynomial k, UPolynomial k)
divMod :: UPolynomial k -> UPolynomial k -> (UPolynomial k, UPolynomial k)
divMod UPolynomial k
f UPolynomial k
g
  | UPolynomial k -> Bool
forall k v. Polynomial k v -> Bool
isZero UPolynomial k
g  = String -> (UPolynomial k, UPolynomial k)
forall a. (?callStack::CallStack) => String -> a
error String
"ToySolver.Data.Polynomial.divMod: division by zero"
  | Bool
otherwise = UPolynomial k -> UPolynomial k -> (UPolynomial k, UPolynomial k)
go UPolynomial k
0 UPolynomial k
f
  where
    lt_g :: Term k X
lt_g = MonomialOrder X -> UPolynomial k -> Term k X
forall k v. Num k => MonomialOrder v -> Polynomial k v -> Term k v
lt MonomialOrder X
nat UPolynomial k
g
    go :: UPolynomial k -> UPolynomial k -> (UPolynomial k, UPolynomial k)
go !UPolynomial k
q !UPolynomial k
r
      | UPolynomial k -> Integer
forall t. Degree t => t -> Integer
deg UPolynomial k
r Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< UPolynomial k -> Integer
forall t. Degree t => t -> Integer
deg UPolynomial k
g = (UPolynomial k
q,UPolynomial k
r)
      | Bool
otherwise     = UPolynomial k -> UPolynomial k -> (UPolynomial k, UPolynomial k)
go (UPolynomial k
q UPolynomial k -> UPolynomial k -> UPolynomial k
forall a. Num a => a -> a -> a
+ UPolynomial k
t) (UPolynomial k
r UPolynomial k -> UPolynomial k -> UPolynomial k
forall a. Num a => a -> a -> a
- UPolynomial k
t UPolynomial k -> UPolynomial k -> UPolynomial k
forall a. Num a => a -> a -> a
* UPolynomial k
g)
        where
          lt_r :: Term k X
lt_r = MonomialOrder X -> UPolynomial k -> Term k X
forall k v. Num k => MonomialOrder v -> Polynomial k v -> Term k v
lt MonomialOrder X
nat UPolynomial k
r
          t :: UPolynomial k
t    = Term k X -> UPolynomial k
forall k v. (Eq k, Num k) => Term k v -> Polynomial k v
fromTerm (Term k X -> UPolynomial k) -> Term k X -> UPolynomial k
forall a b. (a -> b) -> a -> b
$ Term k X
lt_r Term k X -> Term k X -> Term k X
forall k v.
(Fractional k, Ord v) =>
Term k v -> Term k v -> Term k v
`tdiv` Term k X
lt_g

divides :: (Eq k, Fractional k) => UPolynomial k -> UPolynomial k -> Bool
divides :: UPolynomial k -> UPolynomial k -> Bool
divides UPolynomial k
f1 UPolynomial k
f2 = UPolynomial k
f2 UPolynomial k -> UPolynomial k -> UPolynomial k
forall k.
(Eq k, Fractional k) =>
UPolynomial k -> UPolynomial k -> UPolynomial k
`mod` UPolynomial k
f1 UPolynomial k -> UPolynomial k -> Bool
forall a. Eq a => a -> a -> Bool
== UPolynomial k
0

-- | GCD of univariate polynomials
gcd :: (Eq k, Fractional k) => UPolynomial k -> UPolynomial k -> UPolynomial k
gcd :: UPolynomial k -> UPolynomial k -> UPolynomial k
gcd UPolynomial k
f1 UPolynomial k
0  = MonomialOrder X -> UPolynomial k -> UPolynomial k
forall r v.
(Eq r, Fractional r) =>
MonomialOrder v -> Polynomial r v -> Polynomial r v
toMonic MonomialOrder X
nat UPolynomial k
f1
gcd UPolynomial k
f1 UPolynomial k
f2 = UPolynomial k -> UPolynomial k -> UPolynomial k
forall k.
(Eq k, Fractional k) =>
UPolynomial k -> UPolynomial k -> UPolynomial k
gcd UPolynomial k
f2 (UPolynomial k
f1 UPolynomial k -> UPolynomial k -> UPolynomial k
forall k.
(Eq k, Fractional k) =>
UPolynomial k -> UPolynomial k -> UPolynomial k
`mod` UPolynomial k
f2)

-- | LCM of univariate polynomials
lcm :: (Eq k, Fractional k) => UPolynomial k -> UPolynomial k -> UPolynomial k
lcm :: UPolynomial k -> UPolynomial k -> UPolynomial k
lcm UPolynomial k
_ UPolynomial k
0 = UPolynomial k
0
lcm UPolynomial k
0 UPolynomial k
_ = UPolynomial k
0
lcm UPolynomial k
f1 UPolynomial k
f2 = MonomialOrder X -> UPolynomial k -> UPolynomial k
forall r v.
(Eq r, Fractional r) =>
MonomialOrder v -> Polynomial r v -> Polynomial r v
toMonic MonomialOrder X
nat (UPolynomial k -> UPolynomial k) -> UPolynomial k -> UPolynomial k
forall a b. (a -> b) -> a -> b
$ (UPolynomial k
f1 UPolynomial k -> UPolynomial k -> UPolynomial k
forall k.
(Eq k, Fractional k) =>
UPolynomial k -> UPolynomial k -> UPolynomial k
`mod` (UPolynomial k -> UPolynomial k -> UPolynomial k
forall k.
(Eq k, Fractional k) =>
UPolynomial k -> UPolynomial k -> UPolynomial k
gcd UPolynomial k
f1 UPolynomial k
f2)) UPolynomial k -> UPolynomial k -> UPolynomial k
forall a. Num a => a -> a -> a
* UPolynomial k
f2

-- | Extended GCD algorithm
exgcd
  :: (Eq k, Fractional k)
  => UPolynomial k
  -> UPolynomial k
  -> (UPolynomial k, UPolynomial k, UPolynomial k)
exgcd :: UPolynomial k
-> UPolynomial k -> (UPolynomial k, UPolynomial k, UPolynomial k)
exgcd UPolynomial k
f1 UPolynomial k
f2 = (UPolynomial k, UPolynomial k, UPolynomial k)
-> (UPolynomial k, UPolynomial k, UPolynomial k)
forall k1 v v.
(Eq k1, Fractional k1) =>
(Polynomial k1 X, Polynomial k1 v, Polynomial k1 v)
-> (Polynomial k1 X, Polynomial k1 v, Polynomial k1 v)
f ((UPolynomial k, UPolynomial k, UPolynomial k)
 -> (UPolynomial k, UPolynomial k, UPolynomial k))
-> (UPolynomial k, UPolynomial k, UPolynomial k)
-> (UPolynomial k, UPolynomial k, UPolynomial k)
forall a b. (a -> b) -> a -> b
$ UPolynomial k
-> UPolynomial k
-> UPolynomial k
-> UPolynomial k
-> UPolynomial k
-> UPolynomial k
-> (UPolynomial k, UPolynomial k, UPolynomial k)
forall k.
(Eq k, Fractional k) =>
UPolynomial k
-> UPolynomial k
-> UPolynomial k
-> UPolynomial k
-> UPolynomial k
-> UPolynomial k
-> (UPolynomial k, UPolynomial k, UPolynomial k)
go UPolynomial k
f1 UPolynomial k
f2 UPolynomial k
1 UPolynomial k
0 UPolynomial k
0 UPolynomial k
1
  where
    go :: UPolynomial k
-> UPolynomial k
-> UPolynomial k
-> UPolynomial k
-> UPolynomial k
-> UPolynomial k
-> (UPolynomial k, UPolynomial k, UPolynomial k)
go !UPolynomial k
r0 !UPolynomial k
r1 !UPolynomial k
s0 !UPolynomial k
s1 !UPolynomial k
t0 !UPolynomial k
t1
      | UPolynomial k
r1 UPolynomial k -> UPolynomial k -> Bool
forall a. Eq a => a -> a -> Bool
== UPolynomial k
0   = (UPolynomial k
r0, UPolynomial k
s0, UPolynomial k
t0)
      | Bool
otherwise = UPolynomial k
-> UPolynomial k
-> UPolynomial k
-> UPolynomial k
-> UPolynomial k
-> UPolynomial k
-> (UPolynomial k, UPolynomial k, UPolynomial k)
go UPolynomial k
r1 UPolynomial k
r2 UPolynomial k
s1 UPolynomial k
s2 UPolynomial k
t1 UPolynomial k
t2
      where
        (UPolynomial k
q, UPolynomial k
r2) = UPolynomial k
r0 UPolynomial k -> UPolynomial k -> (UPolynomial k, UPolynomial k)
forall k.
(Eq k, Fractional k) =>
UPolynomial k -> UPolynomial k -> (UPolynomial k, UPolynomial k)
`divMod` UPolynomial k
r1
        s2 :: UPolynomial k
s2 = UPolynomial k
s0 UPolynomial k -> UPolynomial k -> UPolynomial k
forall a. Num a => a -> a -> a
- UPolynomial k
qUPolynomial k -> UPolynomial k -> UPolynomial k
forall a. Num a => a -> a -> a
*UPolynomial k
s1
        t2 :: UPolynomial k
t2 = UPolynomial k
t0 UPolynomial k -> UPolynomial k -> UPolynomial k
forall a. Num a => a -> a -> a
- UPolynomial k
qUPolynomial k -> UPolynomial k -> UPolynomial k
forall a. Num a => a -> a -> a
*UPolynomial k
t1
    f :: (Polynomial k1 X, Polynomial k1 v, Polynomial k1 v)
-> (Polynomial k1 X, Polynomial k1 v, Polynomial k1 v)
f (Polynomial k1 X
g,Polynomial k1 v
u,Polynomial k1 v
v)
      | k1
lc_g k1 -> k1 -> Bool
forall a. Eq a => a -> a -> Bool
== k1
0 = (Polynomial k1 X
g, Polynomial k1 v
u, Polynomial k1 v
v)
      | Bool
otherwise = ((k1 -> k1) -> Polynomial k1 X -> Polynomial k1 X
forall k1 k v.
(Eq k1, Num k1) =>
(k -> k1) -> Polynomial k v -> Polynomial k1 v
mapCoeff (k1 -> k1 -> k1
forall a. Fractional a => a -> a -> a
/k1
lc_g) Polynomial k1 X
g, (k1 -> k1) -> Polynomial k1 v -> Polynomial k1 v
forall k1 k v.
(Eq k1, Num k1) =>
(k -> k1) -> Polynomial k v -> Polynomial k1 v
mapCoeff (k1 -> k1 -> k1
forall a. Fractional a => a -> a -> a
/k1
lc_g) Polynomial k1 v
u, (k1 -> k1) -> Polynomial k1 v -> Polynomial k1 v
forall k1 k v.
(Eq k1, Num k1) =>
(k -> k1) -> Polynomial k v -> Polynomial k1 v
mapCoeff (k1 -> k1 -> k1
forall a. Fractional a => a -> a -> a
/k1
lc_g) Polynomial k1 v
v)
      where
        lc_g :: k1
lc_g = MonomialOrder X -> Polynomial k1 X -> k1
forall k v. Num k => MonomialOrder v -> Polynomial k v -> k
lc MonomialOrder X
nat Polynomial k1 X
g

-- | pseudo division
pdivMod :: (Eq r, Num r) => UPolynomial r -> UPolynomial r -> (r, UPolynomial r, UPolynomial r)
pdivMod :: UPolynomial r -> UPolynomial r -> (r, UPolynomial r, UPolynomial r)
pdivMod UPolynomial r
_ UPolynomial r
0 = String -> (r, UPolynomial r, UPolynomial r)
forall a. (?callStack::CallStack) => String -> a
error String
"ToySolver.Data.Polynomial.pdivMod: division by 0"
pdivMod UPolynomial r
f UPolynomial r
g
  | UPolynomial r -> Integer
forall t. Degree t => t -> Integer
deg UPolynomial r
f Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< UPolynomial r -> Integer
forall t. Degree t => t -> Integer
deg UPolynomial r
g = (r
1, UPolynomial r
0, UPolynomial r
f)
  | Bool
otherwise     = Integer
-> UPolynomial r
-> UPolynomial r
-> (r, UPolynomial r, UPolynomial r)
go (UPolynomial r -> Integer
forall t. Degree t => t -> Integer
deg UPolynomial r
f Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- UPolynomial r -> Integer
forall t. Degree t => t -> Integer
deg UPolynomial r
g Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
1) UPolynomial r
f UPolynomial r
0
  where
    (r
lc_g, Monomial X
lm_g) = MonomialOrder X -> UPolynomial r -> (r, Monomial X)
forall k v. Num k => MonomialOrder v -> Polynomial k v -> Term k v
lt MonomialOrder X
nat UPolynomial r
g
    b :: r
b = r
lc_g r -> Integer -> r
forall a b. (Num a, Integral b) => a -> b -> a
^ (UPolynomial r -> Integer
forall t. Degree t => t -> Integer
deg UPolynomial r
f Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
deg_g Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
1)
    deg_g :: Integer
deg_g = UPolynomial r -> Integer
forall t. Degree t => t -> Integer
deg UPolynomial r
g
    go :: Integer
-> UPolynomial r
-> UPolynomial r
-> (r, UPolynomial r, UPolynomial r)
go !Integer
n !UPolynomial r
f1 !UPolynomial r
q
      | Integer
deg_g Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> UPolynomial r -> Integer
forall t. Degree t => t -> Integer
deg UPolynomial r
f1 = (r
b, UPolynomial r
q, r -> UPolynomial r -> UPolynomial r
forall k v. (Eq k, Num k) => k -> Polynomial k v -> Polynomial k v
scale (r
lc_g r -> Integer -> r
forall a b. (Num a, Integral b) => a -> b -> a
^ Integer
n) UPolynomial r
f1)
      | Bool
otherwise      = Integer
-> UPolynomial r
-> UPolynomial r
-> (r, UPolynomial r, UPolynomial r)
go (Integer
n Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
1) (r -> UPolynomial r -> UPolynomial r
forall k v. (Eq k, Num k) => k -> Polynomial k v -> Polynomial k v
scale r
lc_g UPolynomial r
f1 UPolynomial r -> UPolynomial r -> UPolynomial r
forall a. Num a => a -> a -> a
- UPolynomial r
s UPolynomial r -> UPolynomial r -> UPolynomial r
forall a. Num a => a -> a -> a
* UPolynomial r
g) (UPolynomial r
q UPolynomial r -> UPolynomial r -> UPolynomial r
forall a. Num a => a -> a -> a
+ r -> UPolynomial r -> UPolynomial r
forall k v. (Eq k, Num k) => k -> Polynomial k v -> Polynomial k v
scale (r
lc_g r -> Integer -> r
forall a b. (Num a, Integral b) => a -> b -> a
^ (Integer
nInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
1)) UPolynomial r
s)
          where
            (r
lc_f1, Monomial X
lm_f1) = MonomialOrder X -> UPolynomial r -> (r, Monomial X)
forall k v. Num k => MonomialOrder v -> Polynomial k v -> Term k v
lt MonomialOrder X
nat UPolynomial r
f1
            s :: UPolynomial r
s = (r, Monomial X) -> UPolynomial r
forall k v. (Eq k, Num k) => Term k v -> Polynomial k v
fromTerm (r
lc_f1, Monomial X
lm_f1 Monomial X -> Monomial X -> Monomial X
forall v. Ord v => Monomial v -> Monomial v -> Monomial v
`mdiv` Monomial X
lm_g)

-- | pseudo quotient
pdiv :: (Eq r, Num r) => UPolynomial r -> UPolynomial r -> UPolynomial r
pdiv :: UPolynomial r -> UPolynomial r -> UPolynomial r
pdiv UPolynomial r
f UPolynomial r
g =
  case UPolynomial r
f UPolynomial r -> UPolynomial r -> (r, UPolynomial r, UPolynomial r)
forall r.
(Eq r, Num r) =>
UPolynomial r -> UPolynomial r -> (r, UPolynomial r, UPolynomial r)
`pdivMod` UPolynomial r
g of
    (r
_, UPolynomial r
q, UPolynomial r
_) -> UPolynomial r
q

-- | pseudo reminder
pmod :: (Eq r, Num r) => UPolynomial r -> UPolynomial r -> UPolynomial r
pmod :: UPolynomial r -> UPolynomial r -> UPolynomial r
pmod UPolynomial r
_ UPolynomial r
0 = String -> UPolynomial r
forall a. (?callStack::CallStack) => String -> a
error String
"ToySolver.Data.Polynomial.pmod: division by 0"
pmod UPolynomial r
f UPolynomial r
g
  | UPolynomial r -> Integer
forall t. Degree t => t -> Integer
deg UPolynomial r
f Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< UPolynomial r -> Integer
forall t. Degree t => t -> Integer
deg UPolynomial r
g = UPolynomial r
f
  | Bool
otherwise     = Integer -> UPolynomial r -> UPolynomial r
go (UPolynomial r -> Integer
forall t. Degree t => t -> Integer
deg UPolynomial r
f Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- UPolynomial r -> Integer
forall t. Degree t => t -> Integer
deg UPolynomial r
g Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
1) UPolynomial r
f
  where
    (r
lc_g, Monomial X
lm_g) = MonomialOrder X -> UPolynomial r -> (r, Monomial X)
forall k v. Num k => MonomialOrder v -> Polynomial k v -> Term k v
lt MonomialOrder X
nat UPolynomial r
g
    deg_g :: Integer
deg_g = UPolynomial r -> Integer
forall t. Degree t => t -> Integer
deg UPolynomial r
g
    go :: Integer -> UPolynomial r -> UPolynomial r
go !Integer
n !UPolynomial r
f1
      | Integer
deg_g Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> UPolynomial r -> Integer
forall t. Degree t => t -> Integer
deg UPolynomial r
f1 = r -> UPolynomial r -> UPolynomial r
forall k v. (Eq k, Num k) => k -> Polynomial k v -> Polynomial k v
scale (r
lc_g r -> Integer -> r
forall a b. (Num a, Integral b) => a -> b -> a
^ Integer
n) UPolynomial r
f1
      | Bool
otherwise      = Integer -> UPolynomial r -> UPolynomial r
go (Integer
n Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
1) (r -> UPolynomial r -> UPolynomial r
forall k v. (Eq k, Num k) => k -> Polynomial k v -> Polynomial k v
scale r
lc_g UPolynomial r
f1 UPolynomial r -> UPolynomial r -> UPolynomial r
forall a. Num a => a -> a -> a
- UPolynomial r
s UPolynomial r -> UPolynomial r -> UPolynomial r
forall a. Num a => a -> a -> a
* UPolynomial r
g)
          where
            (r
lc_f1, Monomial X
lm_f1) = MonomialOrder X -> UPolynomial r -> (r, Monomial X)
forall k v. Num k => MonomialOrder v -> Polynomial k v -> Term k v
lt MonomialOrder X
nat UPolynomial r
f1
            s :: UPolynomial r
s = (r, Monomial X) -> UPolynomial r
forall k v. (Eq k, Num k) => Term k v -> Polynomial k v
fromTerm (r
lc_f1, Monomial X
lm_f1 Monomial X -> Monomial X -> Monomial X
forall v. Ord v => Monomial v -> Monomial v -> Monomial v
`mdiv` Monomial X
lm_g)

-- | GCD of univariate polynomials
gcd' :: (Integral r) => UPolynomial r -> UPolynomial r -> UPolynomial r
gcd' :: UPolynomial r -> UPolynomial r -> UPolynomial r
gcd' UPolynomial r
f1 UPolynomial r
0  = UPolynomial r -> UPolynomial r
forall r v. (Integral r, Ord v) => Polynomial r v -> Polynomial r v
ppI UPolynomial r
f1
gcd' UPolynomial r
f1 UPolynomial r
f2 = UPolynomial r -> UPolynomial r -> UPolynomial r
forall r.
Integral r =>
UPolynomial r -> UPolynomial r -> UPolynomial r
gcd' UPolynomial r
f2 (UPolynomial r
f1 UPolynomial r -> UPolynomial r -> UPolynomial r
forall r.
(Eq r, Num r) =>
UPolynomial r -> UPolynomial r -> UPolynomial r
`pmod` UPolynomial r
f2)

-- | Is the number a root of the polynomial?
isRootOf :: (Eq k, Num k) => k -> UPolynomial k -> Bool
isRootOf :: k -> UPolynomial k -> Bool
isRootOf k
x UPolynomial k
p = (X -> k) -> UPolynomial k -> k
forall k v. Num k => (v -> k) -> Polynomial k v -> k
eval (\X
_ -> k
x) UPolynomial k
p k -> k -> Bool
forall a. Eq a => a -> a -> Bool
== k
0

-- | Is the polynomial square free?
isSquareFree :: (Eq k, Fractional k) => UPolynomial k -> Bool
isSquareFree :: UPolynomial k -> Bool
isSquareFree UPolynomial k
p = UPolynomial k -> UPolynomial k -> UPolynomial k
forall k.
(Eq k, Fractional k) =>
UPolynomial k -> UPolynomial k -> UPolynomial k
gcd UPolynomial k
p (UPolynomial k -> X -> UPolynomial k
forall k v.
(Eq k, Num k, Ord v) =>
Polynomial k v -> v -> Polynomial k v
deriv UPolynomial k
p X
X) UPolynomial k -> UPolynomial k -> Bool
forall a. Eq a => a -> a -> Bool
== UPolynomial k
1

{--------------------------------------------------------------------
  Term
--------------------------------------------------------------------}

type Term k v = (k, Monomial v)

type UTerm k = Term k X

tdeg :: Term k v -> Integer
tdeg :: Term k v -> Integer
tdeg (k
_,Monomial v
xs) = Monomial v -> Integer
forall t. Degree t => t -> Integer
deg Monomial v
xs

tscale :: (Num k) => k -> Term k v -> Term k v
tscale :: k -> Term k v -> Term k v
tscale k
a (k
c, Monomial v
xs) = (k
ak -> k -> k
forall a. Num a => a -> a -> a
*k
c, Monomial v
xs)

tmult :: (Num k, Ord v) => Term k v -> Term k v -> Term k v
tmult :: Term k v -> Term k v -> Term k v
tmult (k
c1,Monomial v
xs1) (k
c2,Monomial v
xs2) = (k
c1k -> k -> k
forall a. Num a => a -> a -> a
*k
c2, Monomial v
xs1 Monomial v -> Monomial v -> Monomial v
forall v. Ord v => Monomial v -> Monomial v -> Monomial v
`mmult` Monomial v
xs2)

tdivides :: (Ord v) => Term k v -> Term k v -> Bool
tdivides :: Term k v -> Term k v -> Bool
tdivides (k
_,Monomial v
xs1) (k
_,Monomial v
xs2) = Monomial v
xs1 Monomial v -> Monomial v -> Bool
forall v. Ord v => Monomial v -> Monomial v -> Bool
`mdivides` Monomial v
xs2

tdiv :: (Fractional k, Ord v) => Term k v -> Term k v -> Term k v
tdiv :: Term k v -> Term k v -> Term k v
tdiv (k
c1,Monomial v
xs1) (k
c2,Monomial v
xs2) = (k
c1 k -> k -> k
forall a. Fractional a => a -> a -> a
/ k
c2, Monomial v
xs1 Monomial v -> Monomial v -> Monomial v
forall v. Ord v => Monomial v -> Monomial v -> Monomial v
`mdiv` Monomial v
xs2)

tderiv :: (Num k, Ord v) => Term k v -> v -> Term k v
tderiv :: Term k v -> v -> Term k v
tderiv (k
c,Monomial v
xs) v
x =
  case Monomial v -> v -> (Integer, Monomial v)
forall v. Ord v => Monomial v -> v -> (Integer, Monomial v)
mderiv Monomial v
xs v
x of
    (Integer
s,Monomial v
ys) -> (k
c k -> k -> k
forall a. Num a => a -> a -> a
* Integer -> k
forall a b. (Integral a, Num b) => a -> b
fromIntegral Integer
s, Monomial v
ys)

tintegral :: (Fractional k, Ord v) => Term k v -> v -> Term k v
tintegral :: Term k v -> v -> Term k v
tintegral (k
c,Monomial v
xs) v
x =
  case Monomial v -> v -> (Rational, Monomial v)
forall v. Ord v => Monomial v -> v -> (Rational, Monomial v)
mintegral Monomial v
xs v
x of
    (Rational
s,Monomial v
ys) -> (k
c k -> k -> k
forall a. Num a => a -> a -> a
* Rational -> k
forall a. Fractional a => Rational -> a
fromRational Rational
s, Monomial v
ys)

{--------------------------------------------------------------------
  Monic Monomial
--------------------------------------------------------------------}

-- 本当は変数の型に応じて type family で表現を変えたい

-- | Monic monomials
newtype Monomial v = Monomial{ Monomial v -> Map v Integer
mindicesMap :: Map v Integer }
  deriving (Monomial v -> Monomial v -> Bool
(Monomial v -> Monomial v -> Bool)
-> (Monomial v -> Monomial v -> Bool) -> Eq (Monomial v)
forall v. Eq v => Monomial v -> Monomial v -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: Monomial v -> Monomial v -> Bool
$c/= :: forall v. Eq v => Monomial v -> Monomial v -> Bool
== :: Monomial v -> Monomial v -> Bool
$c== :: forall v. Eq v => Monomial v -> Monomial v -> Bool
Eq, Eq (Monomial v)
Eq (Monomial v)
-> (Monomial v -> Monomial v -> Ordering)
-> (Monomial v -> Monomial v -> Bool)
-> (Monomial v -> Monomial v -> Bool)
-> (Monomial v -> Monomial v -> Bool)
-> (Monomial v -> Monomial v -> Bool)
-> (Monomial v -> Monomial v -> Monomial v)
-> (Monomial v -> Monomial v -> Monomial v)
-> Ord (Monomial v)
Monomial v -> Monomial v -> Bool
Monomial v -> Monomial v -> Ordering
Monomial v -> Monomial v -> Monomial v
forall a.
Eq a
-> (a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
forall v. Ord v => Eq (Monomial v)
forall v. Ord v => Monomial v -> Monomial v -> Bool
forall v. Ord v => MonomialOrder v
forall v. Ord v => Monomial v -> Monomial v -> Monomial v
min :: Monomial v -> Monomial v -> Monomial v
$cmin :: forall v. Ord v => Monomial v -> Monomial v -> Monomial v
max :: Monomial v -> Monomial v -> Monomial v
$cmax :: forall v. Ord v => Monomial v -> Monomial v -> Monomial v
>= :: Monomial v -> Monomial v -> Bool
$c>= :: forall v. Ord v => Monomial v -> Monomial v -> Bool
> :: Monomial v -> Monomial v -> Bool
$c> :: forall v. Ord v => Monomial v -> Monomial v -> Bool
<= :: Monomial v -> Monomial v -> Bool
$c<= :: forall v. Ord v => Monomial v -> Monomial v -> Bool
< :: Monomial v -> Monomial v -> Bool
$c< :: forall v. Ord v => Monomial v -> Monomial v -> Bool
compare :: Monomial v -> Monomial v -> Ordering
$ccompare :: forall v. Ord v => MonomialOrder v
$cp1Ord :: forall v. Ord v => Eq (Monomial v)
Ord, Typeable)

type UMonomial = Monomial X

instance (Show v) => Show (Monomial v) where
  showsPrec :: Int -> Monomial v -> ShowS
showsPrec Int
d Monomial v
m  = Bool -> ShowS -> ShowS
showParen (Int
d Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
10) (ShowS -> ShowS) -> ShowS -> ShowS
forall a b. (a -> b) -> a -> b
$
    String -> ShowS
showString String
"mfromIndices " ShowS -> ShowS -> ShowS
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [(v, Integer)] -> ShowS
forall a. Show a => a -> ShowS
shows (Monomial v -> [(v, Integer)]
forall v. Monomial v -> [(v, Integer)]
mindices Monomial v
m)

instance (Ord v, IsString v) => IsString (Monomial v) where
  fromString :: String -> Monomial v
fromString String
s = v -> Monomial v
forall a v. Var a v => v -> a
var (String -> v
forall a. IsString a => String -> a
fromString String
s)

instance (NFData v) => NFData (Monomial v) where
  rnf :: Monomial v -> ()
rnf (Monomial Map v Integer
m) = Map v Integer -> ()
forall a. NFData a => a -> ()
rnf Map v Integer
m

instance Degree (Monomial v) where
  deg :: Monomial v -> Integer
deg (Monomial Map v Integer
m) = [Integer] -> Integer
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ([Integer] -> Integer) -> [Integer] -> Integer
forall a b. (a -> b) -> a -> b
$ Map v Integer -> [Integer]
forall k a. Map k a -> [a]
Map.elems Map v Integer
m

instance Ord v => Var (Monomial v) v where
  var :: v -> Monomial v
var v
x = Map v Integer -> Monomial v
forall v. Map v Integer -> Monomial v
Monomial (Map v Integer -> Monomial v) -> Map v Integer -> Monomial v
forall a b. (a -> b) -> a -> b
$ v -> Integer -> Map v Integer
forall k a. k -> a -> Map k a
Map.singleton v
x Integer
1

instance Ord v => Vars (Monomial v) v where
  vars :: Monomial v -> Set v
vars Monomial v
mm = Map v Integer -> Set v
forall k a. Map k a -> Set k
Map.keysSet (Monomial v -> Map v Integer
forall v. Monomial v -> Map v Integer
mindicesMap Monomial v
mm)

instance Hashable v => Hashable (Monomial v) where
  hashWithSalt :: Int -> Monomial v -> Int
hashWithSalt = (Monomial v -> [(v, Integer)]) -> Int -> Monomial v -> Int
forall b a. Hashable b => (a -> b) -> Int -> a -> Int
hashUsing (Map v Integer -> [(v, Integer)]
forall k a. Map k a -> [(k, a)]
Map.toList (Map v Integer -> [(v, Integer)])
-> (Monomial v -> Map v Integer) -> Monomial v -> [(v, Integer)]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Monomial v -> Map v Integer
forall v. Monomial v -> Map v Integer
mindicesMap)

mone :: Monomial v
mone :: Monomial v
mone = Map v Integer -> Monomial v
forall v. Map v Integer -> Monomial v
Monomial (Map v Integer -> Monomial v) -> Map v Integer -> Monomial v
forall a b. (a -> b) -> a -> b
$ Map v Integer
forall k a. Map k a
Map.empty

mfromIndices :: Ord v => [(v, Integer)] -> Monomial v
mfromIndices :: [(v, Integer)] -> Monomial v
mfromIndices [(v, Integer)]
xs
  | ((v, Integer) -> Bool) -> [(v, Integer)] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any (\(v
_,Integer
e) -> Integer
0Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
>Integer
e) [(v, Integer)]
xs = String -> Monomial v
forall a. (?callStack::CallStack) => String -> a
error String
"ToySolver.Data.Polynomial.mfromIndices: negative exponent"
  | Bool
otherwise = Map v Integer -> Monomial v
forall v. Map v Integer -> Monomial v
Monomial (Map v Integer -> Monomial v) -> Map v Integer -> Monomial v
forall a b. (a -> b) -> a -> b
$ (Integer -> Integer -> Integer) -> [(v, Integer)] -> Map v Integer
forall k a. Ord k => (a -> a -> a) -> [(k, a)] -> Map k a
Map.fromListWith Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
(+) [(v
x,Integer
e) | (v
x,Integer
e) <- [(v, Integer)]
xs, Integer
e Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> Integer
0]

mfromIndicesMap :: Ord v => Map v Integer -> Monomial v
mfromIndicesMap :: Map v Integer -> Monomial v
mfromIndicesMap Map v Integer
m
  | ((v, Integer) -> Bool) -> [(v, Integer)] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any (\(v
_,Integer
e) -> Integer
0Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
>Integer
e) (Map v Integer -> [(v, Integer)]
forall k a. Map k a -> [(k, a)]
Map.toList Map v Integer
m) = String -> Monomial v
forall a. (?callStack::CallStack) => String -> a
error String
"ToySolver.Data.Polynomial.mfromIndicesMap: negative exponent"
  | Bool
otherwise = Map v Integer -> Monomial v
forall v. Map v Integer -> Monomial v
mfromIndicesMap' Map v Integer
m

mfromIndicesMap' :: Map v Integer -> Monomial v
mfromIndicesMap' :: Map v Integer -> Monomial v
mfromIndicesMap' Map v Integer
m = Map v Integer -> Monomial v
forall v. Map v Integer -> Monomial v
Monomial (Map v Integer -> Monomial v) -> Map v Integer -> Monomial v
forall a b. (a -> b) -> a -> b
$ (Integer -> Bool) -> Map v Integer -> Map v Integer
forall a k. (a -> Bool) -> Map k a -> Map k a
Map.filter (Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
>Integer
0) Map v Integer
m

mindices :: Monomial v -> [(v, Integer)]
mindices :: Monomial v -> [(v, Integer)]
mindices = Map v Integer -> [(v, Integer)]
forall k a. Map k a -> [(k, a)]
Map.toAscList (Map v Integer -> [(v, Integer)])
-> (Monomial v -> Map v Integer) -> Monomial v -> [(v, Integer)]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Monomial v -> Map v Integer
forall v. Monomial v -> Map v Integer
mindicesMap

mmult :: Ord v => Monomial v -> Monomial v -> Monomial v
mmult :: Monomial v -> Monomial v -> Monomial v
mmult (Monomial Map v Integer
xs1) (Monomial Map v Integer
xs2) = Map v Integer -> Monomial v
forall v. Map v Integer -> Monomial v
mfromIndicesMap' (Map v Integer -> Monomial v) -> Map v Integer -> Monomial v
forall a b. (a -> b) -> a -> b
$ (Integer -> Integer -> Integer)
-> Map v Integer -> Map v Integer -> Map v Integer
forall k a. Ord k => (a -> a -> a) -> Map k a -> Map k a -> Map k a
Map.unionWith Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
(+) Map v Integer
xs1 Map v Integer
xs2

mpow :: Ord v => Monomial v -> Integer -> Monomial v
mpow :: Monomial v -> Integer -> Monomial v
mpow Monomial v
_ Integer
0 = Monomial v
forall v. Monomial v
mone
mpow Monomial v
m Integer
1 = Monomial v
m
mpow (Monomial Map v Integer
xs) Integer
e
  | Integer
0 Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> Integer
e     = String -> Monomial v
forall a. (?callStack::CallStack) => String -> a
error String
"ToySolver.Data.Polynomial.mpow: negative exponent"
  | Bool
otherwise = Map v Integer -> Monomial v
forall v. Map v Integer -> Monomial v
Monomial (Map v Integer -> Monomial v) -> Map v Integer -> Monomial v
forall a b. (a -> b) -> a -> b
$ (Integer -> Integer) -> Map v Integer -> Map v Integer
forall a b k. (a -> b) -> Map k a -> Map k b
Map.map (Integer
eInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*) Map v Integer
xs

mdivides :: Ord v => Monomial v -> Monomial v -> Bool
mdivides :: Monomial v -> Monomial v -> Bool
mdivides (Monomial Map v Integer
xs1) (Monomial Map v Integer
xs2) = (Integer -> Integer -> Bool)
-> Map v Integer -> Map v Integer -> Bool
forall k a b.
Ord k =>
(a -> b -> Bool) -> Map k a -> Map k b -> Bool
Map.isSubmapOfBy Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
(<=) Map v Integer
xs1 Map v Integer
xs2

mdiv :: Ord v => Monomial v -> Monomial v -> Monomial v
mdiv :: Monomial v -> Monomial v -> Monomial v
mdiv (Monomial Map v Integer
xs1) (Monomial Map v Integer
xs2) = Map v Integer -> Monomial v
forall v. Map v Integer -> Monomial v
Monomial (Map v Integer -> Monomial v) -> Map v Integer -> Monomial v
forall a b. (a -> b) -> a -> b
$ (Integer -> Integer -> Maybe Integer)
-> Map v Integer -> Map v Integer -> Map v Integer
forall k a b.
Ord k =>
(a -> b -> Maybe a) -> Map k a -> Map k b -> Map k a
Map.differenceWith Integer -> Integer -> Maybe Integer
forall a. (Ord a, Num a) => a -> a -> Maybe a
f Map v Integer
xs1 Map v Integer
xs2
  where
    f :: a -> a -> Maybe a
f a
m a
n
      | a
m a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a
n    = Maybe a
forall a. Maybe a
Nothing
      | Bool
otherwise = a -> Maybe a
forall a. a -> Maybe a
Just (a
m a -> a -> a
forall a. Num a => a -> a -> a
- a
n)

mderiv :: Ord v => Monomial v -> v -> (Integer, Monomial v)
mderiv :: Monomial v -> v -> (Integer, Monomial v)
mderiv (Monomial Map v Integer
xs) v
x
  | Integer
nInteger -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
==Integer
0      = (Integer
0, Monomial v
forall v. Monomial v
mone)
  | Bool
otherwise = (Integer
n, Map v Integer -> Monomial v
forall v. Map v Integer -> Monomial v
Monomial (Map v Integer -> Monomial v) -> Map v Integer -> Monomial v
forall a b. (a -> b) -> a -> b
$ (Integer -> Maybe Integer) -> v -> Map v Integer -> Map v Integer
forall k a. Ord k => (a -> Maybe a) -> k -> Map k a -> Map k a
Map.update Integer -> Maybe Integer
forall a. (Ord a, Num a) => a -> Maybe a
f v
x Map v Integer
xs)
  where
    n :: Integer
n = Integer -> v -> Map v Integer -> Integer
forall k a. Ord k => a -> k -> Map k a -> a
Map.findWithDefault Integer
0 v
x Map v Integer
xs
    f :: a -> Maybe a
f a
m
      | a
m a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a
1    = Maybe a
forall a. Maybe a
Nothing
      | Bool
otherwise = a -> Maybe a
forall a. a -> Maybe a
Just (a -> Maybe a) -> a -> Maybe a
forall a b. (a -> b) -> a -> b
$! a
m a -> a -> a
forall a. Num a => a -> a -> a
- a
1

mintegral :: Ord v => Monomial v -> v -> (Rational, Monomial v)
mintegral :: Monomial v -> v -> (Rational, Monomial v)
mintegral (Monomial Map v Integer
xs) v
x =
  (Integer
1 Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
% Integer -> Integer
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Integer
n Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
1), Map v Integer -> Monomial v
forall v. Map v Integer -> Monomial v
Monomial (Map v Integer -> Monomial v) -> Map v Integer -> Monomial v
forall a b. (a -> b) -> a -> b
$ v -> Integer -> Map v Integer -> Map v Integer
forall k a. Ord k => k -> a -> Map k a -> Map k a
Map.insert v
x (Integer
nInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
1) Map v Integer
xs)
  where
    n :: Integer
n = Integer -> v -> Map v Integer -> Integer
forall k a. Ord k => a -> k -> Map k a -> a
Map.findWithDefault Integer
0 v
x Map v Integer
xs

mlcm :: Ord v => Monomial v -> Monomial v -> Monomial v
mlcm :: Monomial v -> Monomial v -> Monomial v
mlcm (Monomial Map v Integer
m1) (Monomial Map v Integer
m2) = Map v Integer -> Monomial v
forall v. Map v Integer -> Monomial v
Monomial (Map v Integer -> Monomial v) -> Map v Integer -> Monomial v
forall a b. (a -> b) -> a -> b
$ (Integer -> Integer -> Integer)
-> Map v Integer -> Map v Integer -> Map v Integer
forall k a. Ord k => (a -> a -> a) -> Map k a -> Map k a -> Map k a
Map.unionWith Integer -> Integer -> Integer
forall a. Ord a => a -> a -> a
max Map v Integer
m1 Map v Integer
m2

mgcd :: Ord v => Monomial v -> Monomial v -> Monomial v
mgcd :: Monomial v -> Monomial v -> Monomial v
mgcd (Monomial Map v Integer
m1) (Monomial Map v Integer
m2) = Map v Integer -> Monomial v
forall v. Map v Integer -> Monomial v
Monomial (Map v Integer -> Monomial v) -> Map v Integer -> Monomial v
forall a b. (a -> b) -> a -> b
$ (Integer -> Integer -> Integer)
-> Map v Integer -> Map v Integer -> Map v Integer
forall k a b c.
Ord k =>
(a -> b -> c) -> Map k a -> Map k b -> Map k c
Map.intersectionWith Integer -> Integer -> Integer
forall a. Ord a => a -> a -> a
min Map v Integer
m1 Map v Integer
m2

mcoprime :: Ord v => Monomial v -> Monomial v -> Bool
mcoprime :: Monomial v -> Monomial v -> Bool
mcoprime Monomial v
m1 Monomial v
m2 = Monomial v -> Monomial v -> Monomial v
forall v. Ord v => Monomial v -> Monomial v -> Monomial v
mgcd Monomial v
m1 Monomial v
m2 Monomial v -> Monomial v -> Bool
forall a. Eq a => a -> a -> Bool
== Monomial v
forall v. Monomial v
mone

{--------------------------------------------------------------------
  Monomial Order
--------------------------------------------------------------------}

type MonomialOrder v = Monomial v -> Monomial v -> Ordering

-- | Lexicographic order
lex :: Ord v => MonomialOrder v
lex :: MonomialOrder v
lex = [(v, Integer)] -> [(v, Integer)] -> Ordering
forall a a. (Ord a, Ord a) => [(a, a)] -> [(a, a)] -> Ordering
go ([(v, Integer)] -> [(v, Integer)] -> Ordering)
-> (Monomial v -> [(v, Integer)]) -> MonomialOrder v
forall b c a. (b -> b -> c) -> (a -> b) -> a -> a -> c
`on` Monomial v -> [(v, Integer)]
forall v. Monomial v -> [(v, Integer)]
mindices
  where
    go :: [(a, a)] -> [(a, a)] -> Ordering
go [] [] = Ordering
EQ
    go [] [(a, a)]
_  = Ordering
LT -- = compare 0 n2
    go [(a, a)]
_ []  = Ordering
GT -- = compare n1 0
    go ((a
x1,a
n1):[(a, a)]
xs1) ((a
x2,a
n2):[(a, a)]
xs2) =
      case a -> a -> Ordering
forall a. Ord a => a -> a -> Ordering
compare a
x1 a
x2 of
        Ordering
LT -> Ordering
GT -- = compare n1 0
        Ordering
GT -> Ordering
LT -- = compare 0 n2
        Ordering
EQ -> a -> a -> Ordering
forall a. Ord a => a -> a -> Ordering
compare a
n1 a
n2 Ordering -> Ordering -> Ordering
forall a. Monoid a => a -> a -> a
`mappend` [(a, a)] -> [(a, a)] -> Ordering
go [(a, a)]
xs1 [(a, a)]
xs2

-- | Reverse lexicographic order.
--
-- Note that revlex is /NOT/ a monomial order.
revlex :: Ord v => Monomial v -> Monomial v -> Ordering
revlex :: Monomial v -> Monomial v -> Ordering
revlex = [(v, Integer)] -> [(v, Integer)] -> Ordering
forall a a. (Ord a, Ord a) => [(a, a)] -> [(a, a)] -> Ordering
go ([(v, Integer)] -> [(v, Integer)] -> Ordering)
-> (Monomial v -> [(v, Integer)])
-> Monomial v
-> Monomial v
-> Ordering
forall b c a. (b -> b -> c) -> (a -> b) -> a -> a -> c
`on` (Map v Integer -> [(v, Integer)]
forall k a. Map k a -> [(k, a)]
Map.toDescList (Map v Integer -> [(v, Integer)])
-> (Monomial v -> Map v Integer) -> Monomial v -> [(v, Integer)]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Monomial v -> Map v Integer
forall v. Monomial v -> Map v Integer
mindicesMap)
  where
    go :: [(a, a)] -> [(a, a)] -> Ordering
go [] [] = Ordering
EQ
    go [] [(a, a)]
_  = Ordering
GT -- = cmp 0 n2
    go [(a, a)]
_ []  = Ordering
LT -- = cmp n1 0
    go ((a
x1,a
n1):[(a, a)]
xs1) ((a
x2,a
n2):[(a, a)]
xs2) =
      case a -> a -> Ordering
forall a. Ord a => a -> a -> Ordering
compare a
x1 a
x2 of
        Ordering
LT -> Ordering
GT -- = cmp 0 n2
        Ordering
GT -> Ordering
LT -- = cmp n1 0
        Ordering
EQ -> a -> a -> Ordering
forall a. Ord a => a -> a -> Ordering
cmp a
n1 a
n2 Ordering -> Ordering -> Ordering
forall a. Monoid a => a -> a -> a
`mappend` [(a, a)] -> [(a, a)] -> Ordering
go [(a, a)]
xs1 [(a, a)]
xs2
    cmp :: a -> a -> Ordering
cmp a
n1 a
n2 = a -> a -> Ordering
forall a. Ord a => a -> a -> Ordering
compare a
n2 a
n1

-- | Graded lexicographic order
grlex :: Ord v => MonomialOrder v
grlex :: MonomialOrder v
grlex = (Integer -> Integer -> Ordering
forall a. Ord a => a -> a -> Ordering
compare (Integer -> Integer -> Ordering)
-> (Monomial v -> Integer) -> MonomialOrder v
forall b c a. (b -> b -> c) -> (a -> b) -> a -> a -> c
`on` Monomial v -> Integer
forall t. Degree t => t -> Integer
deg) MonomialOrder v -> MonomialOrder v -> MonomialOrder v
forall a. Monoid a => a -> a -> a
`mappend` MonomialOrder v
forall v. Ord v => MonomialOrder v
lex

-- | Graded reverse lexicographic order
grevlex :: Ord v => MonomialOrder v
grevlex :: MonomialOrder v
grevlex = (Integer -> Integer -> Ordering
forall a. Ord a => a -> a -> Ordering
compare (Integer -> Integer -> Ordering)
-> (Monomial v -> Integer) -> MonomialOrder v
forall b c a. (b -> b -> c) -> (a -> b) -> a -> a -> c
`on` Monomial v -> Integer
forall t. Degree t => t -> Integer
deg) MonomialOrder v -> MonomialOrder v -> MonomialOrder v
forall a. Monoid a => a -> a -> a
`mappend` MonomialOrder v
forall v. Ord v => MonomialOrder v
revlex