packages feed

polynomial-algebra (empty) → 0.1

raw patch · 28 files changed

+4351/−0 lines, 28 filesdep +arraydep +basedep +compact-word-vectorssetup-changed

Dependencies added: array, base, compact-word-vectors, containers

Files

+ LICENSE view
@@ -0,0 +1,29 @@+Copyright (c) 2018-2019, Balazs Komuves+All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++- Redistributions of source code must retain the above copyright notice,+this list of conditions and the following disclaimer.+ +- Redistributions in binary form must reproduce the above copyright notice,+this list of conditions and the following disclaimer in the documentation+and/or other materials provided with the distribution.+ +- Neither names of the copyright holders nor the names of the contributors+may be used to endorse or promote products derived from this software without+specific prior written permission. ++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER +OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,+EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,+PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR+PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF+LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING+NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS+SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.+
+ README.md view
@@ -0,0 +1,22 @@+polynomial-algebra Haskell library+==================================++This is a Haskell library to compute with multivariate polynomials.++Polynomials are implemented as free modules (with a coefficient ring)+over the monoid of monomials. The free module implementation is basically+a map from monomials to coefficients, with the invariant that zero +coefficients should be never present. Different implementations of monomials+are available with different speed and usability tradeoffs:++* generic monomial over a variable set given by inhabitants of a type+* monomials over x1, x2 ... xn (two different in-memory representations)+* monomials over an infinite number of variables x1, x2, ...+* univariate monomial (basically, an integer exponent)+* exterior monomial (for exterior algebra)++Type level parameters are used for the variable names (used for pretty-printing)+and number of variables where possible.++A type class interface allows the user to work uniformly over different+implementations.
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ polynomial-algebra.cabal view
@@ -0,0 +1,87 @@+cabal-version:        2.4+Name:                 polynomial-algebra+Version:              0.1+Synopsis:             Multivariate polynomial rings+Description:          Multivariate and univariate polynomial rings, with several different representations+License:              BSD-3-Clause+License-file:         LICENSE+Author:               Balazs Komuves+Copyright:            (c) 2018-2019 Balazs Komuves+Maintainer:           bkomuves (plus) hackage (at) gmail (dot) com+Homepage:             https://github.com/bkomuves/polynomial-algebra+Stability:            Experimental+Category:             Math+Tested-With:          GHC == 8.6.5+Build-Type:           Simple++extra-source-files:   README.md++source-repository head+  type:                 git+  location:             https://github.com/bkomuves/polynomial-algebra++--------------------------------------------------------------------------------++Library++  Build-Depends:        base >= 4 && < 5, +                        array >= 0.5, containers >= 0.6, +                        compact-word-vectors >= 0.2.0.2++  Exposed-Modules:      Math.Algebra.Polynomial.Class+                        Math.Algebra.Polynomial.FreeModule+                        Math.Algebra.Polynomial.Monomial.Generic+                        Math.Algebra.Polynomial.Monomial.Indexed+                        Math.Algebra.Polynomial.Monomial.Infinite+                        Math.Algebra.Polynomial.Monomial.Compact+                        Math.Algebra.Polynomial.Monomial.Univariate+                        -- Math.Algebra.Polynomial.Univariate.Dense+                        -- Math.Algebra.Polynomial.Univariate.Sparse+                        Math.Algebra.Polynomial.Monomial.Tensor+                        Math.Algebra.Polynomial.Monomial.Exterior.Indexed+                        Math.Algebra.Polynomial.Univariate+                        Math.Algebra.Polynomial.Univariate.Pochhammer+                        Math.Algebra.Polynomial.Univariate.Lagrange+                        Math.Algebra.Polynomial.Univariate.Cyclotomic+                        Math.Algebra.Polynomial.Univariate.Chebysev+                        Math.Algebra.Polynomial.Univariate.Hermite+                        Math.Algebra.Polynomial.Univariate.Legendre+                        Math.Algebra.Polynomial.Univariate.Bernoulli+                        Math.Algebra.Polynomial.Multivariate.Generic+                        Math.Algebra.Polynomial.Multivariate.Compact+                        Math.Algebra.Polynomial.Multivariate.Indexed+                        Math.Algebra.Polynomial.Multivariate.Infinite+                        Math.Algebra.Polynomial.Exterior.Indexed+                        -- Math.Algebra.Polynomial.Exterior.Generic+                        Math.Algebra.Polynomial.Pretty+                        Math.Algebra.Polynomial.Misc++  Default-Extensions:   CPP, BangPatterns +  Other-Extensions:     ScopedTypeVariables, TypeSynonymInstances, FlexibleContexts, +                        GeneralizedNewtypeDeriving, TypeFamilies, DataKinds++  Default-Language:     Haskell2010++  Hs-Source-Dirs:       src++  ghc-options:          -fwarn-tabs -fno-warn-unused-matches -fno-warn-name-shadowing -fno-warn-unused-imports+    +--------------------------------------------------------------------------------+    +-- test-suite test+-- +--   default-language:     Haskell2010+--   type:                 exitcode-stdio-1.0+--   hs-source-dirs:       test+--   main-is:              testSuite.hs+--   +--   build-depends:        base >= 4 && < 5, containers >= 0.6, array >= 0.5,+--                         tasty >= 0.11, tasty-hunit >= 0.9, +--                         tasty-quickcheck >= 0.9, QuickCheck >= 2.6,+--                         polynomial-algebra >= 0.1+--                        +--   other-modules:        Tests.Common+--                         Tests.Conversion++--------------------------------------------------------------------------------+
+ src/Math/Algebra/Polynomial/Class.hs view
@@ -0,0 +1,190 @@++{-# LANGUAGE +      FlexibleContexts, TypeFamilies, TypeSynonymInstances, FlexibleInstances, +      GeneralizedNewtypeDeriving, ConstraintKinds+  #-}++module Math.Algebra.Polynomial.Class where++--------------------------------------------------------------------------------++import Data.List ( foldl' , foldl1' , maximum , null )+import Data.Typeable+import Data.Proxy++import Math.Algebra.Polynomial.Misc+import Math.Algebra.Polynomial.Pretty+import Math.Algebra.Polynomial.FreeModule ( FreeModule(..) ) ++--------------------------------------------------------------------------------+-- * Indices++-- | The index of a variable. This will be used as the variable type +-- when the set of variables is a continguous set like @{x_1, x_2, ... , x_N}@ +newtype Index +  = Index Int +  deriving (Eq,Ord,Show,Enum)++fromIndex :: Index -> Int+fromIndex (Index j) = j++instance Pretty Index where+  pretty (Index j) = "x_" ++ show j++--------------------------------------------------------------------------------+-- * Rings++-- | The class of coefficient rings. +--+-- Since base rings like integers or rational behave differently than say+-- another polynomial ring as a coefficient ring, we have to be explicit+-- about some things (mainly for pretty-printing purposes+--+-- TODO: clean this up!+class (Eq c, Ord c, Num c, IsSigned c, Show c, Pretty c, Typeable c) => Ring c where+  isZeroR   :: c -> Bool+  signumR   :: c -> Maybe Sign+  absR      :: c -> c+  isSignedR :: Proxy c -> Bool+  isAtomicR :: Proxy c -> Bool++  isZeroR   = (==0)+  signumR   = signOf+  absR      = abs+  isSignedR = const True+  isAtomicR = const True++instance Ring Int+instance Ring Integer+instance Ring Rational++-- | The class of coefficient fields (this is just a constraint synonym for now)+type Field c = (Ring c, Fractional c)++--------------------------------------------------------------------------------++-- | The class of types whose inhabitants can serve as variables+-- (this is just a constraint synonym for now)+type Variable v = (Ord v, Show v, Pretty v, Typeable v)++--------------------------------------------------------------------------------+-- * Monomials++-- | The class of (multivariate) monomials+-- +-- The @Maybe@-s are there to allow truncated and exterior polynomial rings+class (Pretty m) => Monomial m where+  -- | the type of variables+  type VarM m :: *                          +  +  -- checking the invariant+  normalizeM  :: m -> m                         -- ^ enforce the invariant +  isNormalM   :: m -> Bool                      -- ^ check the invariant+  -- construction and deconstruction+  fromListM   :: [(VarM m,Int)] -> m            -- ^ construction from @(variable,exponent)@ pairs+  toListM     :: m -> [(VarM m,Int)]            -- ^ extracting @(variable,exponent)@ pairs+  -- simple monomials+  emptyM      :: m                              -- ^ the empty monomial (corresponding to the polynomial 1)+  isEmptyM    :: m -> Bool                      -- ^ checks whether it is empty+  variableM   :: VarM m        -> m             -- ^ a single variable+  singletonM  :: VarM m -> Int -> m             -- ^ a single variable raised to a power+  -- algebra+  mulM        :: m -> m -> m                    -- ^ multiplication of monomials+  productM    :: [m] -> m                       -- ^ product of several monomials+  powM        :: m -> Int -> m                  -- ^ raising to a power+  divM        :: m -> m -> Maybe m              -- ^ division of monomials+  -- calculus+  diffM       :: Num c => VarM m -> Int -> m -> Maybe (m,c)       -- ^ differentiation+  -- degrees+  maxDegM     :: m -> Int                       -- ^ maximum degree+  totalDegM   :: m -> Int                       -- ^ total degree+  -- substitution and evaluation+  evalM       :: Num c => (VarM m -> c) -> m -> c                 -- ^ ring substitution (evaluation)+  varSubsM    :: (VarM m -> VarM m) -> m -> m                     -- ^ simple variable substitution+  termSubsM   :: Num c => (VarM m -> Maybe c) -> (m,c) -> (m,c)   -- ^ term substitution++{-+  -- some (inefficient) default implementations+  normalizeM     = id+  isNormalM      = const True+  productM       = foldl' mulM emptyM+  mulM a b       = productM [a,b]+  emptyM         = fromListM []+  variableM  v   = fromListM [(v,1)]+  singletonM v e = fromListM [(v,e)]+  maxDegM        = maximum      . map snd . toListM+  totalDegM      = foldl' (+) 0 . map snd . toListM+  isEmptyM       = null . toListM+-}++proxyVarM :: Monomial m => m -> Proxy (VarM m)+proxyVarM _ = Proxy++--------------------------------------------------------------------------------+-- * Polynomial rings++-- | The class of almost polynomial rings+class (Pretty p, Ring (CoeffP p), FreeModule p, CoeffP p ~ CoeffF p, MonomP p ~ BaseF p) => AlmostPolynomial p where++  -- | Type of coefficients+  type CoeffP p :: *+  -- | Type of monomials+  type MonomP p :: *+  -- | Type of variables+  type VarP   p :: *++  -- conversion+  fromListP     :: [(MonomP p, CoeffP p)] -> p       -- ^ construction from @(variable,exponent)@ pairs+  toListP       :: p -> [(MonomP p, CoeffP p)]       -- ^ extracting @(variable,exponent)@ pairs++  -- zero, one+  zeroP         :: p+  isZeroP       :: p -> Bool+  oneP          :: p++  -- construction+  variableP     :: VarP p        -> p                -- ^ a single variable+  singletonP    :: VarP p -> Int -> p                -- ^ a single variable raised to a power+  monomP        :: MonomP p -> p+  monomP'       :: MonomP p -> CoeffP p -> p+  scalarP       :: CoeffP p -> p++  -- algebra+  addP          :: p -> p -> p+  subP          :: p -> p -> p+  negP          :: p -> p +  sumP          :: [p] -> p++  mulP          :: p -> p -> p+  productP      :: [p] -> p++  coeffOfP      :: MonomP p -> p -> CoeffP p+  mulByMonomP   :: MonomP p -> p -> p+  scaleP        :: CoeffP p -> p -> p ++  -- default implementations+  sumP     ps = case ps of { [] -> zeroP ; _ -> foldl1' addP ps }+  productP ps = case ps of { [] -> oneP  ; _ -> foldl1' mulP ps }++--------------------------------------------------------------------------------++-- | The class of polynomial rings+class (AlmostPolynomial p, Num p, Monomial (MonomP p), VarM (MonomP p) ~ VarP p) => Polynomial p where++  evalP         :: Num d => (CoeffP p -> d) -> (VarP p -> d) -> p -> d+  varSubsP      :: (VarP p -> VarP p) -> p -> p+  coeffSubsP    :: (VarP p -> Maybe (CoeffP p)) -> p -> p+  subsP         :: (VarP p -> p) -> p -> p++--------------------------------------------------------------------------------++proxyCoeffP :: AlmostPolynomial p => p -> Proxy (CoeffP p)+proxyCoeffP _ = Proxy++proxyMonomP :: AlmostPolynomial p => p -> Proxy (MonomP p)+proxyMonomP _ = Proxy++proxyVarP :: AlmostPolynomial p => p -> Proxy (VarP p)+proxyVarP _ = Proxy++--------------------------------------------------------------------------------
+ src/Math/Algebra/Polynomial/Exterior/Indexed.hs view
@@ -0,0 +1,166 @@++-- | Exterior algebra where the variable set +-- looks like @{x_1, x_2, ... , x_N}@ +--+-- See <https://en.wikipedia.org/wiki/Exterior_algebra>++{-# LANGUAGE +      BangPatterns, TypeFamilies, DataKinds, KindSignatures, ScopedTypeVariables,+      FlexibleContexts, StandaloneDeriving+  #-}+module Math.Algebra.Polynomial.Exterior.Indexed+  (+    ExtAlg(..) , unExtAlg , polyVar , nOfExtAlg+  , ZExtAlg , QExtAlg+  , embed+  , Ext(..)+  )+  where++--------------------------------------------------------------------------------++import Data.Maybe+import Data.List+import Data.Array.Unboxed ++import Data.Typeable+import GHC.TypeLits+import Data.Proxy++import Data.Foldable as F ++import qualified Math.Algebra.Polynomial.FreeModule as ZMod+import Math.Algebra.Polynomial.FreeModule ( FreeMod , FreeModule(..) ) -- , ZMod , QMod )++import Math.Algebra.Polynomial.Monomial.Exterior.Indexed ++import Math.Algebra.Polynomial.Class+import Math.Algebra.Polynomial.Pretty+import Math.Algebra.Polynomial.Misc++--------------------------------------------------------------------------------+-- * Exterior algebra++-- | An exterior polynomial in with a given coefficient ring. +--+-- It is also indexed by the (shared) /name/ of the variables and the /number of/+-- variable. For example @ExtAlgn Rational "x" 3@ the type of polynomials in the+-- variables @x1, x2, x3@ with rational coefficients.+newtype ExtAlg (coeff :: *) (var :: Symbol) (n :: Nat) +  = ExtAlg (FreeMod coeff (Ext var n))+  deriving (Eq,Ord,Show,Typeable)++unExtAlg :: ExtAlg c v n -> FreeMod c (Ext v n) +unExtAlg (ExtAlg p) = p++-- | Name of the variables+polyVar :: KnownSymbol var => ExtAlg c var n -> String+polyVar = symbolVal . varProxy where+  varProxy :: ExtAlg c var n -> Proxy var+  varProxy _ = Proxy++-- | Number of variables+nOfExtAlg :: KnownNat n => ExtAlg c var n -> Int+nOfExtAlg = fromInteger . natVal . natProxy where+  natProxy :: ExtAlg c var n -> Proxy n+  natProxy _ = Proxy++instance FreeModule (ExtAlg c v n) where+  type BaseF  (ExtAlg c v n) = Ext v n+  type CoeffF (ExtAlg c v n) = c+  toFreeModule   = unExtAlg+  fromFreeModule = ExtAlg++--------------------------------------------------------------------------------++type ZExtAlg = ExtAlg Integer+type QExtAlg = ExtAlg Rational++--------------------------------------------------------------------------------++instance (Ring c, KnownSymbol v, KnownNat n) => AlmostPolynomial (ExtAlg c v n) where+  type CoeffP (ExtAlg c v n) = c+  type MonomP (ExtAlg c v n) = Ext v n+  type VarP   (ExtAlg c v n) = Index++  zeroP         = ExtAlg ZMod.zero+  isZeroP       = ZMod.isZero . unExtAlg+  oneP          = ExtAlg (ZMod.generator emptyExt)++  fromListP     = ExtAlg . ZMod.fromList+  toListP       = ZMod.toList . unExtAlg++  variableP     = ExtAlg . ZMod.generator . variableExt+  singletonP    = error "ExtAlg/singletonP: not implemented (because it is meaningless)"+  monomP        = \m     -> ExtAlg $ ZMod.generator m+  monomP'       = \m c   -> ExtAlg $ ZMod.singleton m c+  scalarP       = \s     -> ExtAlg $ ZMod.singleton emptyExt s++  addP          = \p1 p2 -> ExtAlg $ ZMod.add (unExtAlg p1) (unExtAlg p2)+  subP          = \p1 p2 -> ExtAlg $ ZMod.sub (unExtAlg p1) (unExtAlg p2)+  negP          = ExtAlg . ZMod.neg . unExtAlg+  mulP          = \p1 p2 -> ExtAlg $ ZMod.mulWith'' mulExtCoeff (unExtAlg p1) (unExtAlg p2)++  coeffOfP      = \m p   -> ZMod.coeffOf m (unExtAlg p)+  productP      = \ps    -> ExtAlg $ ZMod.productWith'' emptyExt mulExtCoeff $ map unExtAlg ps+  mulByMonomP   = \m p   -> ExtAlg $ ZMod.mapMaybeBaseCoeff (mulExtCoeff m) (unExtAlg p)    -- not injective!!!+  scaleP        = \s p   -> ExtAlg $ ZMod.scale s (unExtAlg p) ++{-+  evalP      = error "ExtAlg/evalP: not implemented"+  varSubsP   = error "ExtAlg/varSubsP: not implemented"+  coeffSubsP = error "ExtAlg/coeffSubsP: not implemented"+  subsP      = error "ExtAlg/subsP: not implemented"+-}+++instance (Ring c, KnownSymbol v, KnownNat n) => Num (ExtAlg c v n) where+  fromInteger = scalarP . fromInteger+  (+)    = addP+  (-)    = subP+  negate = negP+  (*)    = mulP+  abs    = id+  signum = \_ -> scalarP 1++instance (Ring c, KnownSymbol v, KnownNat n, Pretty (Ext v n)) => Pretty (ExtAlg c v n) where+  pretty poly@(ExtAlg fm) = if isSignedR (proxyCoeffP poly)+    then prettyFreeMod'  True   pretty fm+    else prettyFreeMod'' pretty pretty fm++-- hackety hack hack...+instance IsSigned (ExtAlg c v n) where+  signOf = const (Just Plus)++-- So that we can use it again as a coefficient ring+instance (Ring c, KnownSymbol v, KnownNat n) => Ring (ExtAlg c v n) where+  isZeroR   = ZMod.isZero . unExtAlg+  isAtomicR = const False+  isSignedR = const False+  absR      = id+  signumR   = const (Just Plus)++--------------------------------------------------------------------------------++-- | You can always increase the number of variables; +-- you can also decrease, if don't use the last few ones.+--+-- This will throw an error if you try to eliminate variables which are in fact used.+-- To do that, you can instead substitute 1 or 0 into them.+--+embed :: (Ring c, KnownNat n, KnownNat m) => ExtAlg c v n -> ExtAlg c v m+embed old@(ExtAlg old_fm) = new where+  n = nOfExtAlg old+  m = nOfExtAlg new+  new = ExtAlg $ case compare m n of +    LT -> ZMod.unsafeMapBase project old_fm+    EQ -> ZMod.unsafeMapBase keep    old_fm+    GT -> ZMod.unsafeMapBase extend  old_fm+  extend  (Ext int) = Ext int+  keep    (Ext int) = Ext int+  project (Ext int) = let new = Ext int +                      in  if isNormalExt new +                            then new+                            else error "Exterior/Indexed/embed: the projected variables are actually used!"++--------------------------------------------------------------------------------
+ src/Math/Algebra/Polynomial/FreeModule.hs view
@@ -0,0 +1,439 @@++-- | Free modules over some generator set.  +--+-- This module should be imported qualified.++{-# LANGUAGE +      BangPatterns, FlexibleContexts, FlexibleInstances, TypeSynonymInstances, TypeFamilies,+      ConstraintKinds, KindSignatures+  #-}+module Math.Algebra.Polynomial.FreeModule where++--------------------------------------------------------------------------------++import Prelude   hiding ( sum , product )+import Data.List hiding ( sum , product )++import Data.Monoid+import Data.Ratio+import Data.Maybe++import Data.Typeable+import Data.Proxy++import Control.Monad ( foldM )++import qualified Data.Map.Strict as Map+import Data.Map.Strict (Map)++import Data.Set (Set)++--------------------------------------------------------------------------------+-- * Partial monoids++class PartialMonoid a where+  pmUnit :: a+  pmAdd  :: a -> a -> Maybe a+  pmSum  :: [a] -> Maybe a++  pmSum xs  = case xs of { [] -> Just pmUnit ; (y:ys) -> foldM pmAdd y ys }+  pmAdd x y = pmSum [x,y]++{- undecidable instance+instance Monoid a => PartialMonoid a where+  pmUnit    = mempty+  pmAdd x y = Just $ mappend x y+  pmSum xs  = Just $ mconcat xs+-}++--------------------------------------------------------------------------------+-- * A type class++-- | The reason for this type class is to make using newtype wrappers more convenient+class (Ord (BaseF a)) => FreeModule a where+  type BaseF  a :: *+  type CoeffF a :: *+  toFreeModule   :: a -> FreeMod (CoeffF a) (BaseF a)+  fromFreeModule :: FreeMod (CoeffF a) (BaseF a) -> a+  +instance Ord b => FreeModule (FreeMod c b) where+  type BaseF  (FreeMod c b) = b+  type CoeffF (FreeMod c b) = c+  toFreeModule   = id+  fromFreeModule = id++--------------------------------------------------------------------------------+-- * Free modules++-- | Free module over a coefficient ring with the given base. Internally a map+-- storing the coefficients. We maintain the invariant that the coefficients+-- are never zero.+newtype FreeMod coeff base = FreeMod { unFreeMod :: Map base coeff } deriving (Eq,Ord,Show)++-- | Free module with integer coefficients+type ZMod base = FreeMod Integer base++-- | Free module with rational coefficients+type QMod base = FreeMod Rational base++--------------------------------------------------------------------------------+-- * Support++-- | Number of terms+size :: FreeMod c b -> Int+size (FreeMod table) = Map.size table++-- | The support as a list+supportList :: Ord b => FreeMod c b -> [b] +supportList (FreeMod table) = Map.keys table++-- | The support as a set+supportSet :: Ord b => FreeMod c b -> Set b +supportSet (FreeMod table) = Map.keysSet table++--------------------------------------------------------------------------------++instance (Monoid b, Ord b, Eq c, Num c) => Num (FreeMod c b) where+  (+)    = add+  (-)    = sub+  negate = neg+  (*)    = mul+  fromInteger = konst . fromInteger+  abs      = id       -- error "FreeMod/abs"+  signum _ = konst 1  -- error "FreeMod/signum"++--------------------------------------------------------------------------------+-- * Sanity checking++-- | Should be the identity function+normalize :: (Ord b, Eq c, Num c) => FreeMod c b -> FreeMod c b+normalize = FreeMod . Map.filter (/=0) . unFreeMod++-- | Safe equality testing (should be identical to @==@)+safeEq :: (Ord b, Eq b, Eq c, Num c) => FreeMod c b -> FreeMod c b -> Bool+safeEq x y = normalize x == normalize y++--------------------------------------------------------------------------------+-- * Constructing and deconstructing++-- | The additive unit+zero :: FreeMod c b+zero = FreeMod $ Map.empty++-- | Testing for equality with zero +-- (WARNING: this assumes that the invariant of never having zero coefficients actually holds!)+isZero :: Ord b => FreeMod c b -> Bool+isZero (FreeMod mp) = Map.null mp++-- | A module generator+generator :: Num c => b -> FreeMod c b +generator x = FreeMod $ Map.singleton x 1++-- | A single generator with a coefficient+singleton :: (Ord b, Num c, Eq c) => b -> c -> FreeMod c b+singleton b c = FreeMod $ if c/=0+  then Map.singleton b c+  else Map.empty++-- | Conversion from list. +-- This should handle repeated generators correctly (adding their coefficients).+fromList :: (Eq c, Num c, Ord b) => [(b,c)] -> FreeMod c b+fromList = FreeMod . Map.filter cond . Map.fromListWith (+) where+  cond x = (x /= 0)++fromMap :: (Eq c, Num c, Ord b) => Map b c -> FreeMod c b+fromMap = FreeMod . Map.filter (/=0) ++-- | Returns the sum of the given generator elements+fromGeneratorSet :: (Ord b, Num c) => Set b -> FreeMod c b+fromGeneratorSet = FreeMod . Map.fromSet (const 1)++-- | Returns the sum of the given generator elements +fromGeneratorList :: (Ord b, Eq c, Num c) => [b] -> FreeMod c b+fromGeneratorList bs = FreeMod $ foldl' f Map.empty bs where+  f !old !b = Map.alter g b old+  g !mb     = case mb of+    Nothing   -> Just 1+    Just k    -> let k' = k+1                                  -- when for example working in a finite field+                 in  if k' /= 0 then Just k' else Nothing      -- repeated adding of 1 can result in zero...++-- | Conversion to list +toList :: FreeMod c b -> [(b,c)]+toList = Map.toList . unFreeMod++-- | Extract the coefficient of a generator+coeffOf :: (Ord b, Num c) => b -> FreeMod c b -> c+coeffOf b (FreeMod x) = case Map.lookup b x of+  Just c  -> c+  Nothing -> 0++-- | Finds the term with the largest generator (in the natural ordering of the generators)+findMaxTerm :: (Ord b) => FreeMod c b -> Maybe (b,c)+findMaxTerm (FreeMod m) = if Map.null m+  then Nothing+  else Just (Map.findMax m)++-- | Finds the term with the smallest generator (in the natural ordering of the generators)+findMinTerm :: (Ord b) => FreeMod c b -> Maybe (b,c)+findMinTerm (FreeMod m) = if Map.null m+  then Nothing+  else Just (Map.findMin m)++--------------------------------------------------------------------------------+-- * Basic operations++-- | Negation+neg :: Num c => FreeMod c b -> FreeMod c b +neg (FreeMod m) = FreeMod (Map.map negate m)++-- | Additions+add :: (Ord b, Eq c, Num c) => FreeMod c b -> FreeMod c b -> FreeMod c b+add (FreeMod m1) (FreeMod m2) = FreeMod (Map.mergeWithKey f id id m1 m2) where+  f _ x y = case x+y of { 0 -> Nothing ; z -> Just z }++-- | Subtraction+sub :: (Ord b, Eq c, Num c) => FreeMod c b -> FreeMod c b -> FreeMod c b+sub (FreeMod m1) (FreeMod m2) = FreeMod (Map.mergeWithKey f id (Map.map negate) m1 m2) where+  f _ x y = case x-y of { 0 -> Nothing ; z -> Just z }++-- | Scaling by a number+scale :: (Ord b, Eq c, Num c) => c -> FreeMod c b -> FreeMod c b+scale 0 _           = zero+scale x (FreeMod m) = FreeMod (Map.mapMaybe f m) where+  f y = case x*y of { 0 -> Nothing ; z -> Just z }++-- | Dividing by a number (assuming that the coefficient ring is integral, and each coefficient+-- is divisible by the given number)+divideByConst :: (Ord b, Eq c, Integral c, Show c) => c -> FreeMod c b -> FreeMod c b+divideByConst d (FreeMod m) = FreeMod (Map.mapMaybe f m) where+  f a = case divMod a d of+    (b,0) -> case b of { 0 -> Nothing ; z -> Just z }+    _     -> error $ "FreeMod/divideByConst: not divisible by " ++ show d++-- | Addition after scaling: @A + c*B@. +addScale :: (Ord b, Eq c, Num c) => FreeMod c b -> c -> FreeMod c b -> FreeMod c b+addScale (FreeMod m1) cf (FreeMod m2) = +  if cf == 0 +    then FreeMod m1 +    else FreeMod (Map.mergeWithKey f id (Map.mapMaybe g) m1 m2) +  where+    g     y = case     cf*y of { 0 -> Nothing ; z -> Just z }+    f _ x y = case x + cf*y of { 0 -> Nothing ; z -> Just z }++-- | Subtraction after scaling: @A - c*B@. This is a handy optimization for conversion algorithms.+subScale :: (Ord b, Eq c, Num c) => FreeMod c b -> c -> FreeMod c b -> FreeMod c b+subScale (FreeMod m1) cf (FreeMod m2) = +  if cf == 0 +    then FreeMod m1 +    else FreeMod (Map.mergeWithKey f id (Map.mapMaybe g) m1 m2) +  where+    g     y = case   - cf*y of { 0 -> Nothing ; z -> Just z }+    f _ x y = case x - cf*y of { 0 -> Nothing ; z -> Just z }++--------------------------------------------------------------------------------++-- | Summation+sum :: (Ord b, Eq c, Num c) => [FreeMod c b] -> FreeMod c b+sum []  = zero+sum zms = (foldl1' add) zms++-- | Linear combination+linComb :: (Ord b, Eq c, Num c) => [(c, FreeMod c b)] -> FreeMod c b+linComb = sumWith where++   sumWith :: (Ord b, Eq c, Num c) => [(c, FreeMod c b)] -> FreeMod c b+   sumWith []  = zero+   sumWith zms = sum [ scale c zm | (c,zm) <- zms ]++-- | Expand each generator into a term in another module and then sum the results+flatMap :: (Ord b1, Ord b2, Eq c, Num c) => (b1 -> FreeMod c b2) -> FreeMod c b1 -> FreeMod c b2+flatMap f = sum . map g . Map.assocs . unFreeMod where+  g (x,c) = scale c (f x)++flatMap' :: (Ord b1, Ord b2, Eq c2, Num c2) => (c1 -> c2) -> (b1 -> FreeMod c2 b2) -> FreeMod c1 b1 -> FreeMod c2 b2+flatMap' embed f = sum . map g . Map.assocs . unFreeMod where+  g (x,c) = scale (embed c) (f x)++-- | The histogram of a multiset of generators is naturally an element of the given Z-module.+{-# SPECIALIZE histogram :: Ord b => [b] -> ZMod b #-} +histogram :: (Ord b, Num c) => [b] -> FreeMod c b+histogram xs = FreeMod $ foldl' f Map.empty xs where+  f old x = Map.insertWith (+) x 1 old+  +--------------------------------------------------------------------------------+-- * Rings (at least some simple ones, where the basis form a partial monoid)++-- | The multiplicative unit+one :: (Monoid b, Num c) => FreeMod c b+one = FreeMod (Map.singleton mempty 1)++-- | A constant+konst :: (Monoid b, Eq c, Num c) => c -> FreeMod c b+konst c = FreeMod $ if c/=0 +  then Map.singleton mempty c+  else Map.empty++-- | Multiplying two ring elements+mul :: (Ord b, Monoid b, Eq c, Num c) => FreeMod c b -> FreeMod c b -> FreeMod c b+mul = mulWith (<>)++-- | Multiplying two ring elements, when the base forms a partial monoid+mul' :: (Ord b, PartialMonoid b, Eq c, Num c) => FreeMod c b -> FreeMod c b -> FreeMod c b+mul' = mulWith' (pmAdd)++-- | Product+product :: (Ord b, Monoid b, Eq c, Num c) => [FreeMod c b] -> FreeMod c b+product []  = generator mempty+product xs  = foldl1' mul xs++-- | Product, when the base forms a partial monoid+product' :: (Ord b, PartialMonoid b, Eq c, Num c) => [FreeMod c b] -> FreeMod c b+product' []  = generator pmUnit+product' xs  = foldl1' mul' xs++-- | Multiplying two ring elements, using the given monoid operation on base+mulWith :: (Ord b, Eq c, Num c) => (b -> b -> b) -> FreeMod c b -> FreeMod c b -> FreeMod c b+mulWith op xx yy = normalize $ sum [ (f x c) | (x,c) <- toList xx ] where+  -- fromListWith can result in zero coefficients... +  -- and if the sum is one term only, then it won't rmeove them :(+  f !x !c = FreeMod $ Map.fromListWith (+) [ (op x y, cd) | (y,d) <- ylist , let cd = c*d , cd /= 0 ]+  ylist = toList yy++-- | Multiplication using a \"truncated\" operation on the base+mulWith' :: (Ord b, Eq c, Num c) => (b -> b -> Maybe b) -> FreeMod c b -> FreeMod c b -> FreeMod c b+mulWith' op xx yy = normalize $ sum [ (f x c) | (x,c) <- toList xx ] where+  f !x !c = FreeMod $ Map.fromListWith (+) [ (z, cd) | (y,d) <- ylist , Just z <- [op x y] , let cd = c*d , cd /= 0 ]+  ylist = toList yy++mulWith'' :: (Ord b, Eq c, Num c) => (b -> b -> Maybe (b,c)) -> FreeMod c b -> FreeMod c b -> FreeMod c b+mulWith'' op xx yy = normalize $ sum [ (f x c) | (x,c) <- toList xx ] where+  f !x !c = FreeMod $ Map.fromListWith (+) [ (z, cde) | (y,d) <- ylist , Just (z,e) <- [op x y] , let cde = c*d*e , cde /= 0 ]+  ylist = toList yy++-- | Product, using the given Monoid empty and operation. +--+-- Implementation note: we only use the user-supported +-- empty value in case of an empty product.+productWith :: (Ord b, Eq c, Num c) => b -> (b -> b -> b) -> [FreeMod c b] -> FreeMod c b+productWith empty op []  = generator empty+productWith empty op xs  = foldl1' (mulWith op) xs++productWith' :: (Ord b, Eq c, Num c) => b -> (b -> b -> Maybe b) -> [FreeMod c b] -> FreeMod c b+productWith' empty op []  = generator empty+productWith' empty op xs  = foldl1' (mulWith' op) xs++productWith'' :: (Ord b, Eq c, Num c) => b -> (b -> b -> Maybe (b,c)) -> [FreeMod c b] -> FreeMod c b+productWith'' empty op []  = generator empty+productWith'' empty op xs  = foldl1' (mulWith'' op) xs++-- | Multiplies by a monomial+mulByMonom :: (Eq c, Num c, Ord b, Monoid b) => b -> FreeMod c b -> FreeMod c b+mulByMonom monom = mapBase (mappend monom) ++-- | Multiplies by a monomial (NOTE: we assume that this is an injective operation!!!)+unsafeMulByMonom :: (Ord b, Monoid b) => b -> FreeMod c b -> FreeMod c b+unsafeMulByMonom monom = unsafeMapBase (mappend monom)++-- | Multiplies by a partial monomial +mulByMonom' :: (Eq c, Num c, Ord b, PartialMonoid b) => b -> FreeMod c b -> FreeMod c b+mulByMonom' monom = mapMaybeBase (pmAdd monom) ++unsafeMulByMonom' :: (Ord b, PartialMonoid b) => b -> FreeMod c b -> FreeMod c b+unsafeMulByMonom' monom = unsafeMapMaybeBase (pmAdd monom) ++--------------------------------------------------------------------------------+-- * Integer \/ Rational conversions++freeModCoeffProxy :: FreeMod c b -> Proxy c+freeModCoeffProxy _ = Proxy++{-+typeRepZ, typeRepQ :: TypeRep+typeRepZ = typeOf (0 :: Integer ) +typeRepQ = typeOf (0 :: Rational)+-}++-- | This is an optimized coefficient ring change function. It detects runtime whether the output +-- coefficient ring is also the integers, and does nothing in that case. +fromZMod :: (Num c, Typeable c, Eq c, Num c, Ord b, Typeable b) => ZMod b -> FreeMod c b+fromZMod = unsafeCoeffChange fromInteger++-- | This is an optimized coefficient ring change function. It detects runtime whether the output +-- coefficient ring is also the rational, and does nothing in that case. +fromQMod :: (Fractional c, Typeable c, Eq c, Num c, Ord b, Typeable b) => QMod b -> FreeMod c b+fromQMod = unsafeCoeffChange fromRational++-- | This is an optimized coefficient ring change function. It detects runtime whether the output +-- coefficient ring is the same as the input, and does nothing in that case. +--+-- For this to be valid, it is required that the supported function is identity in+-- the case @c1 ~ c2@ !!!+unsafeCoeffChange :: (Typeable c1, Typeable c2, Eq c2, Num c2, Ord b, Typeable b) => (c1 -> c2) -> FreeMod c1 b -> FreeMod c2 b+unsafeCoeffChange f input = case cast input of+  Just out -> {- trace "*** no cast" $ -} out+  Nothing  -> {- trace "!!! cast"    $ -} mapCoeff f input++-- | Given a polynomial with formally rational coefficients, but whose coeffiecients+-- are actually integers, we return the corresponding polynomial with integer coefficients+unsafeZModFromQMod :: Ord b => QMod b -> ZMod b+unsafeZModFromQMod = mapCoeff f where+  f :: Rational -> Integer+  f q = if denominator q == 1 then numerator q else error "unsafeZModFromQMod: coefficient is not integral"++--------------------------------------------------------------------------------+-- * Misc++-- | Changing the base set+mapBase :: (Ord a, Ord b, Eq c, Num c) => (a -> b) -> FreeMod c a -> FreeMod c b+mapBase f +  = normalize                            -- it can happen that we merge a (-1) and (+1) for example ...+  . onFreeMod (Map.mapKeysWith (+) f)++-- | Changing the base set (the user must guarantee that the map is injective!!)+unsafeMapBase :: (Ord a, Ord b) => (a -> b) -> FreeMod c a -> FreeMod c b+unsafeMapBase f = onFreeMod (Map.mapKeys f)++-- | Changing the coefficient ring+mapCoeff :: (Ord b, Eq c2, Num c2) => (c1 -> c2) -> FreeMod c1 b -> FreeMod c2 b+mapCoeff f = onFreeMod' (Map.mapMaybe mbf) where+  mbf x = case f x of { 0 -> Nothing ; y -> Just y }++mapCoeffWithKey :: (Ord b, Eq c2, Num c2) => (b -> c1 -> c2) -> FreeMod c1 b -> FreeMod c2 b+mapCoeffWithKey f = onFreeMod' (Map.mapMaybeWithKey mbf) where+  mbf k x = case f k x of { 0 -> Nothing ; y -> Just y }++-- | Extract a subset of terms+filterBase :: (Ord b) => (b -> Bool) -> FreeMod c b -> FreeMod c b+filterBase f = onFreeMod (Map.filterWithKey g) where g k _ = f k ++-- | Map and extract a subset of terms +mapMaybeBase :: (Ord a, Ord b, Eq c, Num c) => (a -> Maybe b) -> FreeMod c a -> FreeMod c b+mapMaybeBase f +  = normalize      -- it can happen that we merge a (-1) and (+1) for example ...+  . onFreeMod (Map.fromListWith (+) . mapMaybe g . Map.toList)+  where+    g (k,x) = case f k of { Just k' -> Just (k',x) ; Nothing -> Nothing }++-- | Like 'mapMaybeBase', but the user must guarantee that the @Just@ part of the map is injective!+unsafeMapMaybeBase :: (Ord a, Ord b) => (a -> Maybe b) -> FreeMod c a -> FreeMod c b+unsafeMapMaybeBase f = onFreeMod (Map.fromList . mapMaybe g . Map.toList)+  where+    g (k,x) = case f k of { Just k' -> Just (k',x) ; Nothing -> Nothing }++mapMaybeBaseCoeff :: (Ord a, Ord b, Eq c, Num c) => (a -> Maybe (b,c)) -> FreeMod c a -> FreeMod c b+mapMaybeBaseCoeff f +  = normalize      -- it can happen that we merge a (-1) and (+1) for example ...+  . onFreeMod (Map.fromListWith (+) . mapMaybe g . Map.toList)+  where+    g (k,x) = case f k of +        Just (k',y) -> let z = x*y in if z/=0 then Just (k',z) else Nothing+        Nothing     -> Nothing ++-- | NOTE: This is UNSAFE! The user must guarantee that the map respects the invariants!+onFreeMod :: (Ord a, Ord b) => (Map a c -> Map b c) -> FreeMod c a -> FreeMod c b+onFreeMod f = FreeMod . f . unFreeMod++onFreeMod' :: (Ord a, Ord b) => (Map a c -> Map b d) -> FreeMod c a -> FreeMod d b+onFreeMod' f = FreeMod . f . unFreeMod++--------------------------------------------------------------------------------
+ src/Math/Algebra/Polynomial/Misc.hs view
@@ -0,0 +1,358 @@++-- | Some auxilary functions used internally++{-# LANGUAGE CPP, BangPatterns, TypeSynonymInstances, FlexibleInstances, DeriveFunctor #-}+module Math.Algebra.Polynomial.Misc where++--------------------------------------------------------------------------------++import Data.List+import Data.Monoid+import Data.Ratio+import Data.Array++-- Semigroup became a superclass of Monoid+#if MIN_VERSION_base(4,11,0)     +import Data.Foldable+import Data.Semigroup+#endif++import Control.Monad++import qualified Data.Map.Strict as Map ; import Data.Map (Map)++--------------------------------------------------------------------------------++{-+-- * Pairs++data Pair a +  = Pair a a +  deriving (Eq,Ord,Show,Functor)+-}++--------------------------------------------------------------------------------++equating :: Eq b => (a -> b) -> a -> a -> Bool  +equating f x y = (f x == f y)++--------------------------------------------------------------------------------+-- * Lists++unique :: Ord a => [a] -> [a]+unique = map head . group . sort++-- | Synonym for histogram+count :: Ord b => [b] -> Map b Integer+count = histogram++histogram :: Ord b => [b] -> Map b Integer+histogram xs = foldl' f Map.empty xs where+  f old x = Map.insertWith (+) x 1 old++{-# SPECIALIZE sum' :: [Int]     -> Int     #-}+{-# SPECIALIZE sum' :: [Integer] -> Integer #-}+sum' :: Num a => [a] -> a+sum' = foldl' (+) 0++longZipWith :: (a -> c) -> (b -> c) -> (a -> b -> c) -> [a] -> [b] -> [c]+longZipWith f g h = go where+  go (x:xs) (y:ys) = h x y : go xs ys+  go xs     []     = map f xs+  go []     ys     = map g ys++longReplaceListElem :: a -> Int -> a -> [a] -> [a]+longReplaceListElem x0 j y xs = go j xs  where+  go  0 (x:xs) = y  : xs+  go !i (x:xs) = x  : go (i-1) xs+  go  0 []     = y  : []+  go !i []     = x0 : go (i-1) []++--------------------------------------------------------------------------------+-- * Maps+  +deleteLookup :: Ord a => a -> Map a b -> (Maybe b, Map a b)+deleteLookup k table = (Map.lookup k table, Map.delete k table)  ++unsafeDeleteLookup :: Ord a => a -> Map a b -> (b, Map a b)+unsafeDeleteLookup k table = (fromJust (Map.lookup k table), Map.delete k table) where+  fromJust mb = case mb of+    Just y  -> y+    Nothing -> error "unsafeDeleteLookup: key not found"++-- | Example usage: @insertMap (:[]) (:) ...@+insertMap :: Ord k => (b -> a) -> (b -> a -> a) -> k -> b -> Map k a -> Map k a+insertMap f g k y = Map.alter h k where+  h mb = case mb of+    Nothing -> Just (f y)+    Just x  -> Just (g y x)    ++-- | Example usage: @buildMap (:[]) (:) ...@+buildMap :: Ord k => (b -> a) -> (b -> a -> a) -> [(k,b)] -> Map k a+buildMap f g xs = foldl' worker Map.empty xs where+  worker !old (k,y) = Map.alter h k old where+    h mb = case mb of+      Nothing -> Just (f y)+      Just x  -> Just (g y x)    ++--------------------------------------------------------------------------------+-- * Signs++data Sign+  = Plus                            -- hmm, this way @Plus < Minus@, not sure about that+  | Minus+  deriving (Eq,Ord,Show,Read)++oppositeSign :: Sign -> Sign+oppositeSign s = case s of+  Plus  -> Minus+  Minus -> Plus++mulSign :: Sign -> Sign -> Sign+mulSign s1 s2 = case s1 of+  Plus  -> s2+  Minus -> oppositeSign s2++productOfSigns :: [Sign] -> Sign+productOfSigns = go Plus where+  go !acc []     = acc+  go !acc (x:xs) = case x of+    Plus  -> go acc xs+    Minus -> go (oppositeSign acc) xs++-- | Negate the second argument if the first is odd+negateIfOdd :: (Integral a, Num b) => a -> b -> b+negateIfOdd k y = if even k then y else negate y++--------------------------------------------------------------------------------++-- Semigroup became a superclass of Monoid+#if MIN_VERSION_base(4,11,0)        ++instance Semigroup Sign where+  (<>)    = mulSign+  sconcat = foldl1 mulSign++instance Monoid Sign where+  mempty  = Plus+  mconcat = productOfSigns++#else++instance Monoid Sign where+  mempty  = Plus+  mappend = mulSign+  mconcat = productOfSigns++#endif++--------------------------------------------------------------------------------++{-# SPECIALIZE signValue :: Sign -> Int     #-}+{-# SPECIALIZE signValue :: Sign -> Integer #-}++-- | @+1@ or @-1@+signValue :: Num a => Sign -> a+signValue s = case s of+  Plus  ->  1+  Minus -> -1++{-# SPECIALIZE signed :: Sign -> Int     -> Int     #-}+{-# SPECIALIZE signed :: Sign -> Integer -> Integer #-}++-- | Negate the second argument if the first is 'Minus'+signed :: Num a => Sign -> a -> a+signed s y = case s of+  Plus  -> y+  Minus -> negate y++class IsSigned a where+  signOf :: a -> Maybe Sign++signOfNum :: (Ord a, Num a) => a -> Maybe Sign +signOfNum x = case compare x 0 of+  LT -> Just Minus+  GT -> Just Plus+  EQ -> Nothing++instance IsSigned Int      where signOf = signOfNum+instance IsSigned Integer  where signOf = signOfNum+instance IsSigned Rational where signOf = signOfNum++--------------------------------------------------------------------------------+-- * Numbers++fromRat :: Rational -> Integer+fromRat r = case denominator r of+  1 -> numerator r+  _ -> error "fromRat: not an integer"    ++safeDiv :: Integer -> Integer -> Integer+safeDiv a b = case divMod a b of+  (q,0) -> q+  (q,r) -> error $ "saveDiv: " ++ show a ++ " = " ++ show b ++ " * " ++ show q ++ " + " ++ show r++--------------------------------------------------------------------------------+-- * Basic number theory++-- | A000142.+factorial :: Integral a => a -> Integer+factorial n+  | n <  0    = error "factorial: input should be nonnegative"+  | n == 0    = 1+  | otherwise = product [1..fromIntegral n]++-- | A007318. Note: This is zero for @n<0@ or @k<0@; see also 'signedBinomial' below.+binomial :: Integral a => a -> a -> Integer+binomial n k +  | k > n = 0+  | k < 0 = 0+  | k > (n `div` 2) = binomial n (n-k)+  | otherwise = (product [n'-k'+1 .. n']) `div` (product [1..k'])+  where +    k' = fromIntegral k+    n' = fromIntegral n++moebiusMu :: Num c => Int -> c+moebiusMu n +  | any (>1) expos       =  0+  | even (length primes) =  1+  | otherwise            = -1+  where+    factors = groupIntegerFactors $ integerFactorsTrialDivision (fromIntegral n)+    (primes,expos) = unzip factors++divisors :: Int -> [Int]+divisors n = [ f tup | tup <- tuples' expos ] where+  grps = groupIntegerFactors $ integerFactorsTrialDivision $ fromIntegral n+  (primes,expos) = unzip grps+  int_ps = map fromInteger primes :: [Int]+  f es = foldl' (*) 1 $ zipWith (^) int_ps es++-- | Square-free divisors together with their Mobius mu value+squareFreeDivisors :: Int -> [(Int,Sign)]+squareFreeDivisors n = map f (sublists int_ps) where+  grps = groupIntegerFactors $ integerFactorsTrialDivision $ fromIntegral n+  primes = map fst grps+  int_ps = map fromInteger primes :: [Int]+  f ps = ( foldl' (*) 1 ps , if even (length ps) then Plus else Minus)++-- | List of primes, using tree merge with wheel. Code by Will Ness.+primes :: [Integer]+primes = 2:3:5:7: gaps 11 wheel (fold3t $ roll 11 wheel primes') where                                                             ++  primes' = 11: gaps 13 (tail wheel) (fold3t $ roll 11 wheel primes')+  fold3t ((x:xs): ~(ys:zs:t)) +    = x : union xs (union ys zs) `union` fold3t (pairs t)            +  pairs ((x:xs):ys:t) = (x : union xs ys) : pairs t +  wheel = 2:4:2:4:6:2:6:4:2:4:6:6:2:6:4:2:6:4:6:8:4:2:4:2:  +          4:8:6:4:6:2:4:6:2:6:6:4:2:4:6:2:6:4:2:4:2:10:2:10:wheel +  gaps k ws@(w:t) cs@ ~(c:u) +    | k==c  = gaps (k+w) t u              +    | True  = k : gaps (k+w) t cs  +  roll k ws@(w:t) ps@ ~(p:u) +    | k==p  = scanl (\c d->c+p*d) (p*p) ws : roll (k+w) t u              +    | True  = roll (k+w) t ps   ++  minus xxs@(x:xs) yys@(y:ys) = case compare x y of +    LT -> x : minus xs  yys+    EQ ->     minus xs  ys +    GT ->     minus xxs ys+  minus xs [] = xs+  minus [] _  = []+  +  union xxs@(x:xs) yys@(y:ys) = case compare x y of +    LT -> x : union xs  yys+    EQ -> x : union xs  ys +    GT -> y : union xxs ys+  union xs [] = xs+  union [] ys =ys++--------------------------------------------------------------------------------+-- Prime factorization++-- | Groups integer factors. Example: from [2,2,2,3,3,5] we produce [(2,3),(3,2),(5,1)]  +groupIntegerFactors :: [Integer] -> [(Integer,Int)]+groupIntegerFactors = map f . group . sort where+  f xs = (head xs, length xs)++-- | The naive trial division algorithm.+integerFactorsTrialDivision :: Integer -> [Integer]+integerFactorsTrialDivision n +  | n<1 = error "integerFactorsTrialDivision: n should be at least 1"+  | otherwise = go primes n +  where+    go _  1 = []+    go rs k = sub ps k where+      sub [] k = [k]+      sub qqs@(q:qs) k = case mod k q of+        0 -> q : go qqs (div k q)+        _ -> sub qs k+      ps = takeWhile (\p -> p*p <= k) rs  ++--------------------------------------------------------------------------------+-- * Basic combinatorics++tuples' :: [Int] -> [[Int]]+tuples' [] = [[]]+tuples' (s:ss) = [ x:xs | x <- [0..s] , xs <- tuples' ss ] ++-- | All sublists of a list.+sublists :: [a] -> [[a]]+sublists [] = [[]]+sublists (x:xs) = sublists xs ++ map (x:) (sublists xs)  ++--------------------------------------------------------------------------------+-- * Integer-indexed cache++intCache :: ((Int -> a) -> (Int -> a)) -> (Int -> a)+intCache compute = cached where+  cached n = lkpITable n table+  table    = mkITable (map (compute cached) [0..])+  +newtype ITable a = ITable [Array Int a] ++mkITable :: [a] -> ITable a+mkITable = ITable . go 1 where+  go !siz list = arr : go (2*siz) rest where+    (this,rest) = splitAt siz list+    arr = listArray (0,siz-1) this++lkpITable :: Int -> ITable a -> a        +lkpITable idx (ITable list) = go 1 idx list where+  go !siz !idx (this:rest) = if idx < siz+    then (this ! idx)+    else go (2*siz) (idx-siz) rest++--------------------------------------------------------------------------------+-- * Stirling numbers++-- | Rows of (signed) Stirling numbers of the first kind. OEIS:A008275.+-- Coefficients of the polinomial @(x-1)*(x-2)*...*(x-n+1)@.+-- This function uses the recursion formula.+signedStirling1stArray :: Integral a => a -> Array Int Integer+signedStirling1stArray n+  | n <  1    = error "stirling1stArray: n should be at least 1"+  | n == 1    = listArray (1,1 ) [1]+  | otherwise = listArray (1,n') [ lkp (k-1) - fromIntegral (n-1) * lkp k | k<-[1..n'] ] +  where+    prev = signedStirling1stArray (n-1)+    n' = fromIntegral n :: Int+    lkp j | j <  1    = 0+          | j >= n'   = 0+          | otherwise = prev ! j ++-- | Stirling numbers of the second kind. OEIS:A008277.+-- This function uses an explicit formula.+-- +-- Argument order: @stirling2nd n k@+--+stirling2nd :: Integral a => a -> a -> Integer+stirling2nd n k +  | k==0 && n==0 = 1+  | k < 1        = 0+  | k > n        = 0+  | otherwise = sum xs `div` factorial k where+      xs = [ negateIfOdd (k-i) $ binomial k i * (fromIntegral i)^n | i<-[0..k] ]++--------------------------------------------------------------------------------
+ src/Math/Algebra/Polynomial/Monomial/Compact.hs view
@@ -0,0 +1,252 @@++-- | Multivariate compact monomials where the variable set +-- looks like @{x_1, x_2, ... , x_N}@. +--+-- This is very similar to the \"Indexed\" version, but should have much more+-- compact in-memory representation (which is useful in case of large or many +-- polynomials, and should be in theory also faster, because of cache friendlyness)+--++{-# LANGUAGE CPP, BangPatterns, TypeFamilies, DataKinds, KindSignatures, ScopedTypeVariables #-}+module Math.Algebra.Polynomial.Monomial.Compact where++--------------------------------------------------------------------------------++import Data.List+import Data.Word++import Data.Array.Unboxed  -- used only by compactFromList++#if MIN_VERSION_base(4,11,0)        +import Data.Semigroup+import Data.Monoid+#else+import Data.Monoid+#endif++import Data.Typeable+import GHC.TypeLits+import Data.Proxy++import Data.Foldable as F ++import qualified Data.Vector.Compact.WordVec as V++import Math.Algebra.Polynomial.Class+import Math.Algebra.Polynomial.Pretty+import Math.Algebra.Polynomial.Misc++import Math.Algebra.Polynomial.Monomial.Indexed ( XS , xsFromExponents , xsToExponents )++--------------------------------------------------------------------------------+-- * Monomials++-- | Monomials of the variables @x1,x2,...,xn@. The internal representation is a+-- compact vector of the exponents.+--+-- The type is indexed by the /name/ of the variables, and then the /number/ of variables.+--+-- Note that we assume here that the internal vector has length @n@.+newtype Compact (var :: Symbol) (n :: Nat) +  = Compact V.WordVec +  deriving (Eq,Show,Typeable)++--------------------------------------------------------------------------------++-- note: this must be a monomial ordering!+instance Ord (Compact var n) where +  compare (Compact a) (Compact b) = compare a b++instance KnownNat n => Semigroup (Compact var n) where+  (<>) = mulCompact++instance KnownNat n => Monoid (Compact var n) where+  mempty  = emptyCompact+  mappend = mulCompact++instance KnownSymbol var => Pretty (Compact var n) where +  pretty monom =   +    case [ showXPow i e | (i,e) <- zip [1..] es , e /= 0 ] of +      [] -> "(1)"+      xs -> intercalate "*" xs+    where+      es = compactToWordExpoList monom+      v  = compactVar monom+      showXPow !i !e = case e of+        0 -> "1"+        1 -> v ++ show i+        _ -> v ++ show i ++ "^" ++ show e++-- | Name of the variables+compactVar :: KnownSymbol var => Compact var n -> String+compactVar = symbolVal . varProxy where+  varProxy :: Compact var n -> Proxy var+  varProxy _ = Proxy++-- | Number of variables+nOfCompact :: KnownNat n => Compact var n -> Int+nOfCompact = fromInteger . natVal . natProxy where+  natProxy :: Compact var n -> Proxy n+  natProxy _ = Proxy++--------------------------------------------------------------------------------+-- * Conversion++-- | from @(variable,exponent)@ pairs+compactFromList :: KnownNat n => [(Index,Int)] -> Compact v n+compactFromList list = xs where+  xs  = Compact $ V.fromList {- n -} (elems arr)+  arr = accumArray (+) 0 (1,n) list' :: UArray Int Word+  n   = nOfCompact xs+  list' = map f list :: [(Int,Word)]+  f (Index j , e)+    | j < 1      = error "compactFromList: index out of bounds (too small)"+    | j > n      = error "compactFromList: index out of bounds (too big)"+    | e < 0      = error "compactFromList: negative exponent"+    | otherwise  = (j,fromIntegral e)+  +-- | to @(variable,exponent)@ pairs+compactToList :: Compact v n -> [(Index,Int)]+compactToList (Compact vec) = filter cond $ zipWith f [1..] (V.toList vec) where+  f j e = (Index j, fromIntegral e)+  cond (_,e) = e > 0++-- | from @Word@ exponent list+compactFromWordExpoList :: KnownNat n => [Word] -> Compact var n+compactFromWordExpoList ws = cpt where+  n   = nOfCompact cpt+  cpt = Compact vec+  vec = V.fromList {- n -} (take n (ws ++ repeat 0))++-- | to @Word@ exponent list+compactToWordExpoList :: Compact var n -> [Word]+compactToWordExpoList (Compact vec) = V.toList vec++-- | from @Int@ exponent list+compactFromExponents :: KnownNat n => [Int] -> Compact v n+compactFromExponents = compactFromWordExpoList . map fromIntegral++-- | to @Int@ exponent list+compactToExponents :: KnownNat n => Compact v n -> [Int]+compactToExponents = map fromIntegral . compactToWordExpoList++-- | from 'XS' exponent list+compactFromXS :: KnownNat n => XS v n -> Compact v n +compactFromXS = compactFromExponents . xsToExponents++-- | to 'XS' exponent list+compactToXS :: KnownNat n => Compact v n -> XS v n+compactToXS = xsFromExponents . compactToExponents++--------------------------------------------------------------------------------+-- * empty (all zero exponents)++emptyCompact :: KnownNat n => Compact v n+emptyCompact = xs where +  xs = Compact $ V.fromList' (V.Shape n 4) (replicate n (0::Word))+  n  = nOfCompact xs++isEmptyCompact :: Compact v n -> Bool+isEmptyCompact monom@(Compact vec) = (V.maximum vec == 0)+  -- all (==0) (compactToWordExpoList monom)++--------------------------------------------------------------------------------+-- * normalization++isNormalCompact :: KnownNat n => Compact v n -> Bool+isNormalCompact cpt@(Compact vec) = nOfCompact cpt == V.vecLen vec++--------------------------------------------------------------------------------+-- * creation++variableCompact :: KnownNat n => Index -> Compact v n +variableCompact idx = singletonCompact idx 1++singletonCompact :: KnownNat n => Index -> Int -> Compact v n +singletonCompact (Index j) e0+  | j < 1     = error "singletonCompact: index out of bounds (too small)"+  | j > n     = error "singletonCompact: index out of bounds (too big)"+  | e < 0     = error "singletonCompact: negative exponent"+  | otherwise = cpt+  where+    e    = fromIntegral e0 :: Word+    list = replicate (j-1) 0 ++ e : replicate (n-j) 0+    n    = nOfCompact cpt+    cpt  = Compact $ V.fromList' (V.Shape n (V.bitsNeededFor e)) list++--------------------------------------------------------------------------------+-- * products++mulCompact :: KnownNat n => Compact v n -> Compact v n -> Compact v n+mulCompact (Compact vec1) (Compact vec2) = Compact $ V.add vec1 vec2++productCompact :: (KnownNat n, Foldable f) => f (Compact v n) -> Compact v n+productCompact = F.foldl' mulCompact emptyCompact ++powCompact :: KnownNat n => Compact v n -> Int -> Compact v n+powCompact (Compact vec) e +  | e < 0     = error "powCompact: negative exponent"+  | e == 0    = emptyCompact+  | otherwise = Compact $ V.scale (fromIntegral e) vec+  +divCompact :: KnownNat n => Compact v n -> Compact v n -> Maybe (Compact v n)+divCompact (Compact vec1) (Compact vec2) = Compact <$> V.subtract vec1 vec2++--------------------------------------------------------------------------------+-- * degree++maxDegCompact :: Compact v n -> Int+maxDegCompact (Compact vec) = fromIntegral (V.maximum vec)++totalDegCompact :: Compact v n -> Int+totalDegCompact (Compact vec) = fromIntegral (V.sum vec)++--------------------------------------------------------------------------------+-- * differentiation++diffCompact :: Num c => Index -> Int -> Compact v n -> Maybe (Compact v n, c)+diffCompact = error "diffCompact: not implemented yet"++{-+diffCompact :: Num c => Index -> Int -> Compact v n -> Maybe (Compact v n, c)+diffCompact _         0 cpt          = Just (cpt,1)+diffCompact (Index j) k (Compact ba) =+  if k8 > m8+    then Nothing+    else Just (Compact ba' , fromInteger c) +  where+    k8   = fromIntegral k :: Word8+    m8   = indexByteArray ba (j-1) :: Word8+    m    = fromIntegral m8 :: Int+    ba'  = byteArrayFromList $ change $ byteArrayToList ba+    c    = product [ fromIntegral (m - i) | i<-[0..k-1] ] :: Integer+    change = go 1 where+      go i (x:xs) = if i == j then (x-k8) : xs else x : go (i+1) xs+      go i []     = [] +-}++--------------------------------------------------------------------------------++instance (KnownNat n, KnownSymbol v) => Monomial (Compact v n) where+  type VarM (Compact v n) = Index+  normalizeM = id+  isNormalM  = isNormalCompact+  fromListM  = compactFromList+  toListM    = compactToList+  emptyM     = emptyCompact+  isEmptyM   = isEmptyCompact+  variableM  = variableCompact+  singletonM = singletonCompact+  mulM       = mulCompact+  divM       = divCompact+  productM   = productCompact+  powM       = powCompact +  maxDegM    = maxDegCompact              +  totalDegM  = totalDegCompact+  diffM      = diffCompact+ +  evalM      = error "Compact/evalM: not yet implemented"+  varSubsM   = error "Compact/varSubsM: not yet implemented"+  termSubsM  = error "Compact/termSubsM: not yet implemented"++--------------------------------------------------------------------------------
+ src/Math/Algebra/Polynomial/Monomial/Exterior/Indexed.hs view
@@ -0,0 +1,305 @@++-- | Exterior monomials where the variable set +-- looks like @{x_1, x_2, ... , x_N}@ +--+-- The internal representation is a bit vector++{-# LANGUAGE +      CPP, BangPatterns, TypeFamilies, DataKinds, KindSignatures, ScopedTypeVariables,+      FlexibleContexts+  #-}+module Math.Algebra.Polynomial.Monomial.Exterior.Indexed where++--------------------------------------------------------------------------------++import Data.Maybe+import Data.List+import Data.Array.Unboxed +import Data.Ord+import Data.Bits++import Data.Typeable+import GHC.TypeLits+import Data.Proxy++import Data.Foldable as F ++import qualified Data.Set as Set ; import Data.Set (Set)   -- IntSet and+import qualified Data.Map as Map ; import Data.Map (Map)   -- IntMap would be more efficient++import Math.Algebra.Polynomial.Class+import Math.Algebra.Polynomial.Pretty+import Math.Algebra.Polynomial.Misc++import Math.Algebra.Polynomial.FreeModule ( PartialMonoid(..) )++--------------------------------------------------------------------------------+-- * Exterior monomials++-- | Exterior monomials of the variables @x1,x2,...,xn@. The internal representation is +-- a bit vector encoded as an @Integer@+newtype Ext (var :: Symbol) (n :: Nat) = Ext Integer deriving (Eq,Ord,Show)++-- | Signed exterior monomial+data SgnExt (var :: Symbol) (n :: Nat) = SgnExt !Sign !(Ext var n) deriving (Eq,Ord,Show)++unExt :: Ext v n -> Integer+unExt (Ext k) = k++instance KnownNat n => PartialMonoid (SgnExt var n) where+  pmUnit = emptySgnExt+  pmAdd  = mulSgnExt   ++instance KnownSymbol var => Pretty (Ext var n) where +  pretty monom =   +    case [ showX i | i <- extToList monom ] of +      [] -> "(1)"+      xs -> intercalate "/\\" xs+    where+      v = extVar monom+      showX !(Index i) = v ++ show i++instance KnownSymbol var => Pretty (SgnExt var n) where +  pretty (SgnExt sgn ext) = case sgn of+    Plus  -> '+' : pretty ext+    Minus -> '-' : pretty ext++-- | Name of the variables+extVar :: KnownSymbol var => Ext var n -> String+extVar = symbolVal . varProxy where+  varProxy :: Ext var n -> Proxy var+  varProxy _ = Proxy++-- | Number of variables+nOfExt :: KnownNat n => Ext var n -> Int+nOfExt = fromInteger . natVal . natProxy where+  natProxy :: Ext var n -> Proxy n+  natProxy _ = Proxy++nOfMbExt :: KnownNat n => Maybe (Ext var n) -> Int+nOfMbExt = fromInteger . natVal . natProxy where+  natProxy :: Maybe (Ext var n) -> Proxy n+  natProxy _ = Proxy++nOfSgnExt :: KnownNat n => SgnExt var n -> Int+nOfSgnExt = fromInteger . natVal . natProxy where+  natProxy :: SgnExt var n -> Proxy n+  natProxy _ = Proxy++nOfMbSgnExt :: KnownNat n => Maybe (SgnExt var n) -> Int+nOfMbSgnExt = fromInteger . natVal . natProxy where+  natProxy :: Maybe (SgnExt var n) -> Proxy n+  natProxy _ = Proxy++--------------------------------------------------------------------------------+-- * emptyness++emptyExt :: KnownNat n => Ext v n+emptyExt = Ext 0++emptySgnExt :: KnownNat n => SgnExt v n+emptySgnExt = SgnExt Plus (Ext 0)++isEmptyExt :: Ext v n -> Bool+isEmptyExt (Ext a) = (a==0)++--------------------------------------------------------------------------------+-- * normalization++isNormalExt :: KnownNat n => Ext v n -> Bool+isNormalExt ext@(Ext int) = shiftR int n == 0 where n = nOfExt ext++--------------------------------------------------------------------------------+-- * conversion++extFromList :: KnownNat n => [Index] -> Maybe (SgnExt v n)+extFromList list = result where+  result+    | Map.null multiset                     = Just emptySgnExt+    | fst (Map.findMin multiset) < Index 1  = error "extFromList: index out of bounds (too small)"+    | fst (Map.findMax multiset) > Index n  = Nothing  -- should this be an error ?+    | any (>1) (Map.elems multiset)         = Nothing+    | otherwise                             = Just (SgnExt sgn $ Ext int)+  n = nOfMbSgnExt result+  multiset = Map.fromListWith (+) [ (idx,1) | idx <- list ] +  perm = sortingPermutationAsc list+  sgn = if odd (numberOfInversionsMerge perm) then Minus else Plus+  int = sum' [ shiftL 1 (j-1) | Index j <- Map.keys multiset ]++extToList :: Ext v n -> [Index]+extToList (Ext int) = go 0 int where+  go _   0 = []+  go !i !a = if testBit a 0 +    then Index (i+1) : go (i+1) (shiftR a 1)+    else               go (i+1) (shiftR a 1)++extFromSet :: KnownNat n => Set Index -> Maybe (Ext v n)+extFromSet set = result where+  n = nOfMbExt result+  result = case Set.lookupMax set of+    Nothing        -> Just emptyExt+    Just (Index k) -> if k > n+      then Nothing   -- should this be an error ?+      else if (Set.findMin set < Index 1)+        then error "extFromSet: index smaller than 1"+        else Just $ Ext $ sum' [ shiftL 1 (j-1) | Index j <- Set.toList set ]+  +extToSet :: Ext v n -> Set Index+extToSet = Set.fromList . extToList++--------------------------------------------------------------------------------+-- * creation++variableExt :: KnownNat n => Index -> Ext v n +variableExt (Index idx) = Ext $ shiftL 1 (idx-1)++--------------------------------------------------------------------------------+-- * multiplication++mulExt :: KnownNat n => Ext v n -> Ext v n -> Maybe (SgnExt v n)+mulExt x y = extFromList (extToList x ++ extToList y)++mulExtCoeff :: (Num c, KnownNat n) => Ext v n -> Ext v n -> Maybe (Ext v n, c)+mulExtCoeff x y = case mulExt x y of+  Nothing             -> Nothing+  Just (SgnExt sgn z) -> case sgn of+    Plus  -> Just (z,  1)+    Minus -> Just (z, -1)++mulSgnExt :: KnownNat n => SgnExt v n -> SgnExt v n -> Maybe (SgnExt v n)+mulSgnExt (SgnExt s x) (SgnExt t y) = case mulExt x y of+  Nothing           -> Nothing+  Just (SgnExt u z) -> Just $ SgnExt (s `mappend` t `mappend` u) z++productExt :: (KnownNat n, Foldable f) => f (Ext v n) -> Maybe (SgnExt v n)+productExt t = go Plus (F.toList t) where+  go !sgn list = case list of+    []         -> Just (SgnExt sgn emptyExt)+    [x]        -> Just (SgnExt sgn x)+    (x:y:rest) -> case mulExt x y of+      Nothing           -> Nothing+      Just (SgnExt s z) -> go (sgn `mappend` s) (z:rest) ++--------------------------------------------------------------------------------++[v1,v2,v3,v4,v5,v6] = [ variableExt (Index i) | i<-[1..6] ] :: [Ext "a" 7]+Just a = productExt [v1,v2,v3]+Just b = productExt [v4,v5,v6]++prop_graded_anticomm :: KnownNat n => Ext v n -> Ext v n -> Bool+prop_graded_anticomm x y +  | isNothing mb1 && isNothing mb2   = True+  | isJust    mb1 && isJust    mb2   = u == v && maybeFlip s == t+  | otherwise                        = False+  where+    mb1 = x `mulExt` y+    mb2 = y `mulExt` x+    Just (SgnExt s u) = mb1+    Just (SgnExt t v) = mb2+    d1 = degreeExt x+    d2 = degreeExt y+    maybeFlip = if odd (d1*d2) then oppositeSign else id    ++prop_graded_anticomm_sgn :: KnownNat n => SgnExt v n -> SgnExt v n -> Bool+prop_graded_anticomm_sgn x y +  | isNothing mb1 && isNothing mb2   = True+  | isJust    mb1 && isJust    mb2   = u == v && maybeFlip s == t+  | otherwise                        = False+  where+    mb1 = x `mulSgnExt` y+    mb2 = y `mulSgnExt` x+    Just (SgnExt s u) = mb1+    Just (SgnExt t v) = mb2+    d1 = degreeSgnExt x+    d2 = degreeSgnExt y+    maybeFlip = if odd (d1*d2) then oppositeSign else id    ++--------------------------------------------------------------------------------+-- * degree++degreeExt :: Ext v n -> Int+degreeExt (Ext k) = popCount k++degreeSgnExt :: SgnExt v n -> Int+degreeSgnExt (SgnExt sgn ext) = degreeExt ext++--------------------------------------------------------------------------------+-- * Permutations++--+newtype Permutation = Permutation (UArray Int Int) deriving (Eq,Ord)++toPermutationUnsafe :: [Int] -> Permutation+toPermutationUnsafe xs = Permutation perm where+  n    = length xs+  perm = listArray (1,n) xs++sortingPermutationAsc :: Ord a => [a] -> Permutation+sortingPermutationAsc xs = toPermutationUnsafe (map fst sorted) where+  sorted = sortBy (comparing snd) $ zip [1..] xs++{-+isEvenPermutation :: Permutation -> Bool+isEvenPermutation (Permutation perm) = res where++  (1,n) = bounds perm+  res = runST $ do+    tag <- newArray (1,n) False +    cycles <- unfoldM (step tag) 1 +    return $ even (sum cycles)+    +  step :: STUArray s Int Bool -> Int -> ST s (Int,Maybe Int)+  step tag k = do+    cyclen <- worker tag k k 0+    m <- next tag (k+1)+    return (cyclen,m)+    +  next :: STUArray s Int Bool -> Int -> ST s (Maybe Int)+  next tag k = if k > n+    then return Nothing+    else readArray tag k >>= \b -> if b +      then next tag (k+1)  +      else return (Just k)+      +  worker :: STUArray s Int Bool -> Int -> Int -> Int -> ST s Int+  worker tag k l cyclen = do+    writeArray tag l True+    let m = perm ! l+    if m == k +      then return cyclen+      else worker tag k m (1+cyclen)      +-}++isEvenPermutation :: Permutation -> Bool+isEvenPermutation = even . numberOfInversionsMerge++-- | Returns the number of inversions, using the merge-sort algorithm.+-- This should be @O(n*log(n))@+--+numberOfInversionsMerge :: Permutation -> Int+numberOfInversionsMerge (Permutation arr) = fst (sortCnt n $ elems arr) where+  (_,n) = bounds arr+                                        +  -- | First argument is length of the list.+  -- Returns also the inversion count.+  sortCnt :: Int -> [Int] -> (Int,[Int])+  sortCnt 0 _     = (0,[] )+  sortCnt 1 [x]   = (0,[x])+  sortCnt 2 [x,y] = if x>y then (1,[y,x]) else (0,[x,y])+  sortCnt n xs    = mergeCnt (sortCnt k us) (sortCnt l vs) where+    k = div n 2+    l = n - k +    (us,vs) = splitAt k xs++  mergeCnt :: (Int,[Int]) -> (Int,[Int]) -> (Int,[Int])+  mergeCnt (!c,us) (!d,vs) = (c+d+e,ws) where++    (e,ws) = go 0 us vs ++    go !k xs [] = ( k*length xs , xs )+    go _  [] ys = ( 0 , ys)+    go !k xxs@(x:xs) yys@(y:ys) = if x < y+      then let (a,zs) = go  k     xs yys in (a+k, x:zs)+      else let (a,zs) = go (k+1) xxs  ys in (a  , y:zs)++--------------------------------------------------------------------------------
+ src/Math/Algebra/Polynomial/Monomial/Generic.hs view
@@ -0,0 +1,192 @@++-- | Multivariate monomials where the set of variables is given by +-- the inhabitants of a type++{-# LANGUAGE CPP, BangPatterns, TypeFamilies, KindSignatures, StandaloneDeriving, FlexibleContexts #-}+module Math.Algebra.Polynomial.Monomial.Generic where++--------------------------------------------------------------------------------++import Data.Maybe+import Data.List+import Data.Foldable as F++import Data.Proxy+import Data.Typeable++#if MIN_VERSION_base(4,11,0)        +import Data.Semigroup+import Data.Monoid+import Data.List.NonEmpty ( NonEmpty )+#else+import Data.Monoid+#endif++import Data.Map.Strict (Map) ; import qualified Data.Map.Strict as Map ++import Math.Algebra.Polynomial.Class+import Math.Algebra.Polynomial.FreeModule+import Math.Algebra.Polynomial.Pretty+import Math.Algebra.Polynomial.Misc++--------------------------------------------------------------------------------+-- * Monomials++-- | A monomial over the set of variables represented by the inhabitants +-- of the type @v@.+-- The invariant we keep is that the exponents present in the @Map@ are always positive.+newtype Monom v +  = Monom (Map v Int)+  deriving (Eq,Ord,Show,Typeable)++unMonom :: Monom v -> Map v Int+unMonom (Monom table) = table++normalizeMonom :: Ord v => Monom v -> Monom v +normalizeMonom (Monom m) = Monom $ Map.filter (>0) m++isNormalMonom :: Monom v -> Bool+isNormalMonom (Monom m) = all (>0) (Map.elems m)++monomToList :: Ord v => Monom v -> [(v,Int)]+monomToList (Monom m) = Map.toList m++monomFromList :: Ord v => [(v,Int)] -> Monom v+monomFromList list = Monom $ Map.fromList $ filter f list where+  f (v,e) | e <  0    = error "monomFromList: negative exponent"+          | e == 0    = False+          | otherwise = True++-- | Note: we can collapse variables together!+mapMonom :: (Ord v, Ord w) => (v -> w) -> Monom v -> Monom w+mapMonom f (Monom old) = Monom $ Map.mapKeysWith (+) f old ++emptyMonom :: Monom v+emptyMonom = Monom Map.empty++isEmptyMonom :: Monom v -> Bool+isEmptyMonom (Monom m) = Map.null m++mulMonom :: Ord v => Monom v -> Monom v -> Monom v+mulMonom (Monom m1) (Monom m2) = Monom (Map.unionWith (+) m1 m2)++divMonom :: Ord v => Monom v -> Monom v -> Maybe (Monom v)+divMonom (Monom m1) (Monom m2) = +  if all (>=0) (Map.elems pre_monom) +    then Just $ Monom pre_monom+    else Nothing+  where+    minus_m2  = Map.map negate m2+    pre_monom = Map.filter (/=0) +              $ Map.unionWith (+) m1 minus_m2++{-# SPECIALIZE prodMonoms :: Ord v => [Monom v] ->  Monom v #-}+prodMonoms :: (Foldable f, Ord v) => f (Monom v) -> Monom v+prodMonoms list = Monom $ Map.unionsWith (+) $ map unMonom $ F.toList list++powMonom :: Ord v => Monom v -> Int -> Monom v+powMonom (Monom m) e +  | e < 0     = error "powMonom: expecting a non-negative exponent"+  | e == 0    = emptyMonom+  | otherwise = Monom (Map.map (*e) m)++varMonom :: Ord v => v -> Monom v+varMonom v = Monom (Map.singleton v 1)++singletonMonom :: Ord v => v -> Int -> Monom v+singletonMonom v e +  | e < 0     = error "singletonMonom: expecting a non-negative exponent"+  | e == 0    = emptyMonom+  | otherwise = Monom (Map.singleton v e)++maxDegMonom :: Monom v -> Int+maxDegMonom (Monom m) = maximum (Map.elems m)++totalDegMonom :: Monom v -> Int+totalDegMonom (Monom m) = sum' (Map.elems m)++evalMonom :: (Num c, Ord v) => (v -> c) -> Monom v -> c+evalMonom f (Monom m) = foldl' (*) 1 (map g $ Map.toList m) where+  g (v,e) = f v ^ e++termSubsMonom :: (Num c, Ord v) => (v -> Maybe c) -> (Monom v, c) -> (Monom v, c)+termSubsMonom f (Monom m , c0) = (Monom m' , c0*proj) where+  m'   = Map.fromList $ catMaybes mbs+  list = Map.toList m+  (proj, mbs)  = mapAccumL g 1 list +  g !s (!v,!e) = case f v of+    Nothing     -> ( s       , Just (v,e) )   +    Just c      -> ( s * c^e , Nothing    )++--------------------------------------------------------------------------------+-- * differentiation++diffMonom :: (Ord v, Num c) => v -> Int -> Monom v -> Maybe (Monom v, c)+diffMonom _ 0 mon           = Just (mon,1)+diffMonom v k (Monom table) =+  if k > m +    then Nothing+    else Just (Monom table' , fromInteger c) +  where+    m      = Map.findWithDefault 0 v table+    table' = Map.insert v (m-k) table+    c      = Data.List.product [ fromIntegral (m-i) | i<-[0..k-1] ] :: Integer++--------------------------------------------------------------------------------++-- Semigroup became a superclass of Monoid+#if MIN_VERSION_base(4,11,0)        ++{-# SPECIALIZE prodMonoms :: Ord v => NonEmpty (Monom v) ->  Monom v #-}++instance Ord v => Semigroup (Monom v) where+  (<>)    = mulMonom+  sconcat = prodMonoms++instance Ord v => Monoid (Monom v)  where+  mempty  = emptyMonom+  mconcat = prodMonoms++#else++instance Ord v => Monoid (Monom v) where+  mempty  = emptyMonom+  mappend = mulMonom+  mconcat = prodMonoms++#endif++--------------------------------------------------------------------------------++instance Pretty v => Pretty (Monom v) where+  pretty (Monom m) = worker (Map.toList m) where+    worker []    = "1"+    worker pairs = intercalate "*" (map f pairs) where+      f (b,0) = "1"                                      -- shouldn't normall happen+      f (b,1) = pretty b+      f (b,k) = pretty b ++ "^" ++ show k++--------------------------------------------------------------------------------++instance (Ord v, Pretty v) => Monomial (Monom v) where+  type VarM (Monom v) = v+  normalizeM  = normalizeMonom+  isNormalM   = isNormalMonom+  fromListM   = monomFromList+  toListM     = monomToList+  emptyM      = emptyMonom+  isEmptyM    = isEmptyMonom+  variableM   = varMonom+  singletonM  = singletonMonom+  mulM        = mulMonom+  divM        = divMonom+  productM    = prodMonoms+  powM        = powMonom+  maxDegM     = maxDegMonom              +  totalDegM   = totalDegMonom+  diffM       = diffMonom+  evalM       = evalMonom+  varSubsM    = mapMonom+  termSubsM   = termSubsMonom++--------------------------------------------------------------------------------
+ src/Math/Algebra/Polynomial/Monomial/Indexed.hs view
@@ -0,0 +1,270 @@++-- | Multivariate monomials where the variable set +-- looks like @{x_1, x_2, ... , x_N}@ ++{-# LANGUAGE +      CPP, BangPatterns, TypeFamilies, DataKinds, KindSignatures, ScopedTypeVariables,+      FlexibleContexts+  #-}+module Math.Algebra.Polynomial.Monomial.Indexed where++--------------------------------------------------------------------------------++import Data.Maybe+import Data.List+import Data.Array.Unboxed ++#if MIN_VERSION_base(4,11,0)        +import Data.Semigroup+import Data.Monoid+#else+import Data.Monoid+#endif++import Data.Typeable+import GHC.TypeLits+import Data.Proxy++import Data.Foldable as F ++import Math.Algebra.Polynomial.Class+import Math.Algebra.Polynomial.Pretty+import Math.Algebra.Polynomial.Misc++--------------------------------------------------------------------------------+-- * Monomials++-- | Monomials of the variables @x1,x2,...,xn@. The internal representation is the+-- dense array of exponents: @x1^e1*x2^e2*...*xn^en@ is represented by @[e1,e2,...,en]@.+--+-- The type is indexed by the /name/ of the variables, and then the /number/ of variables.+--+-- Note that we require here that the array has bounds @(1,n)@+newtype XS (var :: Symbol) (n :: Nat) = XS (UArray Int Int) deriving (Eq,Show,Typeable)++-- | Note: this must be a monomial ordering!+instance Ord (XS var n) where compare (XS a) (XS b) = compare a b++instance KnownNat n => Semigroup (XS var n) where+  (<>) = mulXS++instance KnownNat n => Monoid (XS var n) where+  mempty  = emptyXS+  mappend = mulXS++instance KnownSymbol var => Pretty (XS var n) where +  pretty monom@(XS arr) =   +    case [ showXPow i e | (i,e) <- zip [1..] es , e /= 0 ] of +      [] -> "(1)"+      xs -> intercalate "*" xs+    where+      es = elems arr+      v = xsVar monom+      showXPow !i !e = case e of+        0 -> "1"+        1 -> v ++ show i+        _ -> v ++ show i ++ "^" ++ show e++-- | Name of the variables+xsVar :: KnownSymbol var => XS var n -> String+xsVar = symbolVal . varProxy where+  varProxy :: XS var n -> Proxy var+  varProxy _ = Proxy++-- | Number of variables+nOfXS :: KnownNat n => XS var n -> Int+nOfXS = fromInteger . natVal . natProxy where+  natProxy :: XS var n -> Proxy n+  natProxy _ = Proxy++--------------------------------------------------------------------------------+-- * emptyness++emptyXS :: KnownNat n => XS v n+emptyXS = xs where +  xs = XS $ accumArray const 0 (1,n) []+  n  = nOfXS xs++isEmptyXS :: XS v n -> Bool+isEmptyXS (XS arr) = all (==0) (elems arr)++--------------------------------------------------------------------------------+-- * normalization++isNormalXS :: KnownNat n => XS v n -> Bool+isNormalXS xs@(XS arr) = bounds arr == (1,n) where n = nOfXS xs++--------------------------------------------------------------------------------+-- * conversion++-- | from @(variable,exponent)@ pairs+xsFromList :: KnownNat n => [(Index,Int)] -> XS v n+xsFromList list = xs where+  xs = XS $ accumArray (+) 0 (1,n) list'+  n  = nOfXS xs+  list' = map f list +  f (Index j , e)+    | j < 1      = error "xsFromList: index out of bounds (too small)"+    | j > n      = error "xsFromList: index out of bounds (too big)"+    | e < 0      = error "xsFromList: negative exponent"+    | otherwise  = (j,e)+  +-- | to @(variable,exponent)@ pairs+xsToList :: XS v n -> [(Index,Int)]+xsToList (XS arr) = [ (Index j, e) | (j,e) <- assocs arr , e > 0 ]++--------------------------------------------------------------------------------++-- | from exponent list+xsFromExponents :: KnownNat n => [Int] -> XS v n+xsFromExponents expos = xs where+  xs   = XS $ listArray (1,n) list+  n    = nOfXS xs+  list = take n (expos ++ repeat 0)++-- | to exponent list+xsToExponents :: KnownNat n => XS v n -> [Int]+xsToExponents (XS arr) = elems arr++--------------------------------------------------------------------------------+-- * creation++variableXS :: KnownNat n => Index -> XS v n +variableXS idx = singletonXS idx 1++singletonXS :: KnownNat n => Index -> Int -> XS v n +singletonXS (Index j) e +  | j < 1     = error "singletonXS: index out of bounds (too small)"+  | j > n     = error "singletonXS: index out of bounds (too big)"+  | e < 0     = error "singletonXS: negative exponent"+  | otherwise = xs+  where+    xs = XS $ accumArray (+) 0 (1,n) [(j,e)]+    n = nOfXS xs++--------------------------------------------------------------------------------+-- * multiplication++mulXS :: KnownNat n => XS v n -> XS v n -> XS v n+mulXS xs1@(XS es) (XS fs) = ys where+  ys = XS $ listArray (1,n) $ zipWith (+) (elems es) (elems fs) where+  n  = nOfXS xs1++productXS :: (KnownNat n, Foldable f) => f (XS v n) -> XS v n+productXS = F.foldl' mulXS emptyXS ++powXS :: XS v n -> Int -> XS v n+powXS (XS arr) e +  | e < 0     = error "powXS: negative exponent"+  | e == 0    = XS (amap (const 0) arr)+  | otherwise = XS (amap (*e)      arr)++divXS :: KnownNat n => XS v n -> XS v n -> Maybe (XS v n)+divXS xs1@(XS es) (XS fs) +  | all (>=0) gs  = Just (XS $ listArray (1,n) gs)+  | otherwise     = Nothing+  where+    gs = zipWith (-) (elems es) (elems fs) where+    n  = nOfXS xs1++--------------------------------------------------------------------------------+-- * degree++maxDegXS :: XS v n -> Int+maxDegXS (XS arr) = maximum (elems arr)++totalDegXS :: XS v n -> Int+totalDegXS (XS arr) = sum' (elems arr)++--------------------------------------------------------------------------------+-- * evaluation and substituion++evalXS :: Num c => (Index -> c) -> XS v n -> c+evalXS f xs@(XS arr) = foldl' (*) 1 (map g $ assocs arr) where+  g (!j,!e) = case e of+    0 -> 1+    1 -> f (Index j) +    _ -> f (Index j) ^ e ++varSubsXS :: KnownNat n => (Index -> Index) -> XS v n -> XS v n+varSubsXS f xs@(XS arr) = XS arr' where+  n    = nOfXS xs+  arr' = accumArray (+) 0 (1,n) list+  list = [ ( myFromIndex (f (Index j)) , e ) | (j,e) <- assocs arr ]+  myFromIndex (Index j)  +    | j >= 1 && j <= 1  = j+    | otherwise         = error "varSubsXS: variable index out of bounds"++termSubsXS :: (KnownNat n, Num c) => (Index -> Maybe c) -> (XS v n, c) -> (XS v n, c) +termSubsXS f (xs@(XS arr) , c0) = (XS arr', c0*proj) where+  n    = nOfXS xs+  arr' = accumArray (+) 0 (1,n) $ catMaybes mbs+  (proj,mbs)   = mapAccumL g 1 (assocs arr)+  g !s (!j,!e) = case f (Index j) of+    Nothing     -> (s       , Just (j,e) )+    Just c      -> (s * c^e , Nothing    )+ +--------------------------------------------------------------------------------+-- * differentiation++diffXS :: Num c => Index -> Int -> XS v n -> Maybe (XS v n, c)+diffXS _         0 xs          = Just (xs,1)+diffXS (Index j) k xs@(XS arr) =+  if k > m +    then Nothing+    else Just (XS arr' , fromInteger c) +  where+    m    = arr!j+    arr' = arr // [(j,m-k)]+    c    = product [ fromIntegral (m-i) | i<-[0..k-1] ] :: Integer++--------------------------------------------------------------------------------++instance (KnownNat n, KnownSymbol v) => Monomial (XS v n) where+  type VarM (XS v n) = Index+  normalizeM  = id+  isNormalM   = isNormalXS+  fromListM   = xsFromList+  toListM     = xsToList+  emptyM      = emptyXS+  isEmptyM    = isEmptyXS+  variableM   = variableXS+  singletonM  = singletonXS+  mulM        = mulXS+  divM        = divXS+  productM    = productXS+  powM        = powXS+  maxDegM     = maxDegXS              +  totalDegM   = totalDegXS+  diffM       = diffXS+  evalM       = evalXS+  varSubsM    = varSubsXS+  termSubsM   = termSubsXS++--------------------------------------------------------------------------------+-- * changing the number of variables++-- | You can always increase the number of variables; +-- you can also decrease, if don't use the last few ones.+--+-- This will throw an error if you try to eliminate variables which are in fact used.+-- To do that, you can instead substitute 1 into them.+--+embedXS :: (KnownNat n, KnownNat m) => XS v n -> XS v m +embedXS old = new where+  n = nOfXS old+  m = nOfXS new+  new = case compare m n of +    LT -> project old+    EQ -> keep    old+    GT -> extend  old+  extend  (XS arr) = XS $ listArray (1,m) (elems arr ++ replicate (m-n) 0)+  keep    (XS arr) = XS arr+  project (XS arr) = +    let old = elems arr+        (new,rest) = splitAt m old+    in  if all (==0) rest +          then XS $ listArray (1,m) new+          else error "Indexed/embed: the projected variables are actually used!"++--------------------------------------------------------------------------------
+ src/Math/Algebra/Polynomial/Monomial/Infinite.hs view
@@ -0,0 +1,216 @@++-- | Multivariate monomials where the variable set is the countable infinite +-- set @{x_1, x_2, x_3,... }@ ++{-# LANGUAGE +      CPP, BangPatterns, TypeFamilies, DataKinds, KindSignatures, ScopedTypeVariables,+      FlexibleContexts+  #-}+module Math.Algebra.Polynomial.Monomial.Infinite where++--------------------------------------------------------------------------------++import Data.Maybe+import Data.List+import Data.Array.Unboxed ++#if MIN_VERSION_base(4,11,0)        +import Data.Semigroup+import Data.Monoid+#else+import Data.Monoid+#endif++import Data.Typeable+import GHC.TypeLits+import Data.Proxy++import Data.Foldable as F ++import Data.Map.Strict ( Map )+import qualified Data.Map.Strict as Map++import Math.Algebra.Polynomial.Class+import Math.Algebra.Polynomial.Pretty+import Math.Algebra.Polynomial.Misc++--------------------------------------------------------------------------------+-- * Monomials++-- | Monomials of the variables @x1,x2,x3,...@. The internal representation is a+-- list of exponents: @x1^e1*x2^e2*x3^e3...@ is represented by @[e1,e2,e3,...]@.+--+-- We assume that only finitely many nonzero exponents appear.+--+-- The type is indexed by the /name/ of the variables.+--+newtype XInf (var :: Symbol) = XInf [Int] deriving (Eq,Show,Typeable)++unXInf :: XInf var -> [Int]+unXInf (XInf ns) = ns++-- The opposite order does not makes sense here...+instance Ord (XInf var) where compare (XInf a) (XInf b) = compare a b++instance Semigroup (XInf var) where+  (<>) = mulXInf++instance Monoid (XInf var) where+  mempty  = emptyXInf+  mappend = mulXInf++instance KnownSymbol var => Pretty (XInf var) where +  pretty monom@(XInf es) =   +    case [ showXPow i e | (i,e) <- zip [1..] es , e /= 0 ] of +      [] -> "(1)"+      xs -> intercalate "*" xs+    where+      v  = xInfVar monom+      showXPow !i !e = case e of+        0 -> "1"+        1 -> v ++ show i+        _ -> v ++ show i ++ "^" ++ show e++-- | Name of the variables+xInfVar :: KnownSymbol var => XInf var -> String+xInfVar = symbolVal . varProxy where+  varProxy :: XInf var -> Proxy var+  varProxy _ = Proxy++--------------------------------------------------------------------------------++emptyXInf :: XInf v+emptyXInf = XInf []++isEmptyXInf :: XInf v -> Bool+isEmptyXInf (XInf arr) = all (==0) arr++mulXInf :: XInf v -> XInf v -> XInf v+mulXInf (XInf es) (XInf fs) = XInf $ longZipWith id id (+) es fs++productXInf :: (Foldable f) => f (XInf v) -> XInf v+productXInf = F.foldl' mulXInf emptyXInf ++divXInf :: XInf v -> XInf v -> Maybe (XInf v)+divXInf xs1@(XInf es) (XInf fs) +  | all (>=0) gs  = Just (XInf gs)+  | otherwise     = Nothing+  where+    gs = longZipWith id negate (-) es fs where++--------------------------------------------------------------------------------++xInfFromList :: [(Index,Int)] -> XInf v+xInfFromList list = +  case Map.lookupMax table of+    Nothing    -> XInf []+    Just (n,_) -> XInf [ Map.findWithDefault 0 i table | i<-[Index 1 .. n] ]+  where+    table = Map.fromListWith (+) list+  +xInfToList :: XInf v -> [(Index,Int)]+xInfToList (XInf arr) +  = filter cond +  $ zip [ Index j | j<-[1..] ] arr +  where+    cond (_,e) = e > 0++xInfToMap :: XInf var -> Map Index Int+xInfToMap = Map.fromList . xInfToList++--------------------------------------------------------------------------------++normalizeXInf :: XInf v -> XInf v+normalizeXInf (XInf arr) = XInf $ reverse $ dropWhile (==0) $ reverse arr++isNormalXInf :: XInf v -> Bool+isNormalXInf (XInf arr) = null (takeWhile (==0) $ reverse arr) ++--------------------------------------------------------------------------------++variableXInf :: Index -> XInf v +variableXInf idx = singletonXInf idx 1++singletonXInf :: Index -> Int -> XInf v +singletonXInf (Index j) e +  | j < 1     = error "singletonXInf: index out of bounds (smaller than 1)"+  | e < 0     = error "singletonXInf: negative exponent"+  | otherwise = XInf $ replicate (j-1) 0 ++ [e]++--------------------------------------------------------------------------------++powXInf :: XInf v -> Int -> XInf v+powXInf (XInf arr) e +  | e < 0     = error "powXInf: negative exponent"+  | e == 0    = XInf []+  | otherwise = XInf (map (*e)      arr)++--------------------------------------------------------------------------------++maxDegXInf :: XInf v -> Int+maxDegXInf (XInf arr) = maximum arr++totalDegXInf :: XInf v -> Int+totalDegXInf (XInf arr) = sum' arr++--------------------------------------------------------------------------------++evalXInf :: Num c => (Index -> c) -> XInf v -> c+evalXInf f xinf = foldl' (*) 1 (map g $ xInfToList xinf) where+  g (!j,!e) = case e of+    0 ->  1+    1 ->  f j +    _ -> (f j) ^ e ++varSubsXInf :: (Index -> Index) -> XInf v -> XInf v+varSubsXInf f xinf = new where+  table = xInfToMap xinf+  new   = xInfFromList [ (f v , e) | (v,e) <- Map.assocs table ] +  -- NOTE: ^^^^^^^^ xInfFromList handles repeated variables!++termSubsXInf :: (Num c) => (Index -> Maybe c) -> (XInf v, c) -> (XInf v, c) +termSubsXInf f (xinf, c0) = (xInfFromList list, c1) where+  (list,c1) = foldl' g ([],c0) (xInfToList xinf)+  g (old,c) (v,e) = case f v of+      Just d  -> (old , c * d^e)+      Nothing -> ((v,e):old , c)+  +--------------------------------------------------------------------------------+-- * differentiation+ +diffXInf :: Num c => Index -> Int -> XInf v -> Maybe (XInf v, c)+diffXInf _         0 xinf      = Just (xinf,1)+diffXInf (Index j) k (XInf es) =+  if k > m +    then Nothing+    else Just (XInf es' , fromInteger c) +  where+    m    = (es ++ repeat 0) !! (j-1)+    es'  = longReplaceListElem 0 (j-1) (m-k) es+    c    = product [ fromIntegral (m-i) | i<-[0..k-1] ] :: Integer++--------------------------------------------------------------------------------++instance (KnownSymbol v) => Monomial (XInf v) where+  type VarM (XInf v) = Index+  normalizeM  = normalizeXInf+  isNormalM   = isNormalXInf+  fromListM   = xInfFromList+  toListM     = xInfToList+  emptyM      = emptyXInf+  isEmptyM    = isEmptyXInf+  variableM   = variableXInf+  singletonM  = singletonXInf+  mulM        = mulXInf+  divM        = divXInf+  productM    = productXInf+  powM        = powXInf+  diffM       = diffXInf+  maxDegM     = maxDegXInf              +  totalDegM   = totalDegXInf+  evalM       = evalXInf+  varSubsM    = varSubsXInf+  termSubsM   = termSubsXInf++--------------------------------------------------------------------------------+
+ src/Math/Algebra/Polynomial/Monomial/Tensor.hs view
@@ -0,0 +1,149 @@++-- | Tensor product (that is, pairs) of monomials++{-# LANGUAGE CPP, BangPatterns, TypeFamilies, UnicodeSyntax, KindSignatures, DataKinds #-}+module Math.Algebra.Polynomial.Monomial.Tensor where++--------------------------------------------------------------------------------++import Data.Typeable+import Data.Either++import Data.Proxy+import GHC.TypeLits++#if MIN_VERSION_base(4,11,0)        +import Data.Semigroup+import Data.Monoid+#else+import Data.Monoid+#endif++import Math.Algebra.Polynomial.Class+import Math.Algebra.Polynomial.Pretty++--------------------------------------------------------------------------------++-- | Elementary tensors (basically pairs). The phantom type parameter+-- @symbol@ is used to render an infix symbol when pretty-printing+data Tensor (symbol :: Symbol) (a :: *) (b :: *) = Tensor !a !b deriving (Eq,Ord,Show,Typeable)++instance (Semigroup a, Semigroup b) => Semigroup (Tensor sym a b) where+  (<>) (Tensor x1 y1) (Tensor x2 y2) = Tensor (x1<>x2) (y1<>y2)+  +instance (Monoid a, Monoid b) => Monoid (Tensor sym a b) where+  mempty = Tensor mempty mempty+  mappend (Tensor x1 y1) (Tensor x2 y2) = Tensor (x1 `mappend` x2) (y1 `mappend` y2)++instance (KnownSymbol sym, Pretty a, Pretty b) => Pretty (Tensor sym a b) where+  pretty t@(Tensor a b) = pretty a ++ tensorSymbol t ++ pretty b+  +tensorSymbol :: KnownSymbol sym => Tensor sym a b -> String+tensorSymbol = symbolVal . symProxy where+  symProxy :: Tensor sym a b -> Proxy sym+  symProxy _ = Proxy++--------------------------------------------------------------------------------++flip :: Tensor sym a b -> Tensor sym b a+flip (Tensor x y) = Tensor y x++--------------------------------------------------------------------------------+-- * Injections++injLeft :: Monoid b => a -> Tensor sym a b+injLeft x = Tensor x mempty++injRight :: Monoid a => b -> Tensor sym a b+injRight x = Tensor mempty x++--------------------------------------------------------------------------------+-- * Projections++projLeft :: Tensor sym a b -> a+projLeft (Tensor x _) = x++projRight :: Tensor sym a b -> b+projRight (Tensor _ y) = y++--------------------------------------------------------------------------------+-- * differentiation+ +diffTensor :: (Monomial a, Monomial b, Num c) => Either (VarM a) (VarM b) -> Int -> Tensor sym a b -> Maybe (Tensor sym a b, c)+diffTensor ei k (Tensor left right) = case ei of+  Left v  -> case diffM v k left of+    Just (left' ,c) -> Just (Tensor left' right , c)+    Nothing         -> Nothing+  Right v -> case diffM v k right of+    Just (right',c) -> Just (Tensor left  right', c)+    Nothing         -> Nothing++--------------------------------------------------------------------------------++instance (KnownSymbol sym, Monomial a, Monomial b) => Monomial (Tensor sym a b) where+  type VarM (Tensor sym a b) = Either (VarM a) (VarM b)+  +  -- checking the invariant+  normalizeM  (Tensor x y) = Tensor (normalizeM x) (normalizeM y)+  isNormalM   (Tensor x y) = isNormalM x && isNormalM y++  -- construction and deconstruction+  fromListM   list = Tensor (fromListM list1) (fromListM list2) where+                (list1,list2) = partitionEithers $ map distEither list                                         +  toListM     (Tensor x y) = map f (toListM x) ++ map g (toListM y) where+                f (v,e) = (Left  v, e)+                g (v,e) = (Right v, e)++  -- simple monomials+  emptyM      = Tensor emptyM emptyM+  isEmptyM    (Tensor x y) = isEmptyM x && isEmptyM y+  variableM   ei = case ei of +                       Left  v -> Tensor (variableM v) emptyM+                       Right v -> Tensor emptyM (variableM v)+  singletonM  ei k = case ei of +                       Left  v -> Tensor (singletonM v k) emptyM+                       Right v -> Tensor emptyM (singletonM v k)+  -- algebra+  mulM        (Tensor x1 y1) (Tensor x2 y2) = Tensor (mulM x1 x2) (mulM y1 y2)+  productM    tensors = Tensor (productM $ map projLeft tensors) (productM $ map projRight tensors)+  powM        (Tensor x y) k = Tensor (powM x k) (powM y k)++  divM        (Tensor x1 y1) (Tensor x2 y2) = case (divM x1 x2, divM y1 y2) of+                  (Just z1 , Just z2) -> Just (Tensor z1 z2)+                  (_       , _      ) -> Nothing++  -- calculus+  diffM = diffTensor++  -- degrees+  maxDegM     (Tensor x y) = max (maxDegM x) (maxDegM y)+  totalDegM   (Tensor x y) = totalDegM x + totalDegM y++  -- substitution and evaluation+  evalM       f (Tensor x y) = evalM (f . Left) x * evalM (f . Right) y+  varSubsM    f (Tensor x y) = Tensor x' y' where+                  x' = varSubsM (unsafeFromLeft  . f . Left ) x+                  y' = varSubsM (unsafeFromRight . f . Right) y+  termSubsM   f (Tensor x y, c) = (Tensor x' y', c*d*e) where+                  (x',d) = termSubsM (f . Left ) (x,1)+                  (y',e) = termSubsM (f . Right) (y,1)++--------------------------------------------------------------------------------+-- * Helpers++distEither :: (Either a b, c) -> Either (a,c) (b,c)+distEither (ei, z) = case ei of+  Left  x -> Left  (x,z)+  Right y -> Right (y,z)++unsafeFromLeft :: Either a b -> a+unsafeFromLeft ei = case ei of +  Left  x -> x+  Right _ -> error "unsafeFromLeft: Right"++unsafeFromRight :: Either a b -> b+unsafeFromRight ei = case ei of +  Left  _ -> error "unsafeFromRight: Left"+  Right y -> y++--------------------------------------------------------------------------------
+ src/Math/Algebra/Polynomial/Monomial/Univariate.hs view
@@ -0,0 +1,117 @@++-- | Univariate \"monomials\" (basically the natural numbers)++{-# LANGUAGE BangPatterns, DataKinds, KindSignatures, TypeFamilies #-}+module Math.Algebra.Polynomial.Monomial.Univariate where++--------------------------------------------------------------------------------++import Data.Array ( assocs ) +import Data.List++#if MIN_VERSION_base(4,11,0)        +import Data.Semigroup+import Data.Monoid+#else+import Data.Monoid+#endif++import Data.Typeable+import GHC.TypeLits+import Data.Proxy++import Math.Algebra.Polynomial.Class+import Math.Algebra.Polynomial.Pretty+import Math.Algebra.Polynomial.Misc++--------------------------------------------------------------------------------+-- * Univariate monomials++-- | A monomial in a univariate polynomial, indexed by its name, eg @U "x"@+newtype U (var :: Symbol) = U Int deriving (Eq,Ord,Show,Typeable)++-- | Name of the variable+uVar :: KnownSymbol var => U var -> String+uVar = symbolVal . uproxy where+  uproxy :: U var -> Proxy var+  uproxy _ = Proxy++instance KnownSymbol var => Pretty (U var) where+  pretty u@(U e) = case e of+    0 -> "1"+    1 -> uVar u+    _ -> uVar u ++ "^" ++ show e++--------------------------------------------------------------------------------++#if MIN_VERSION_base(4,11,0)        ++instance Semigroup (U var) where+  (<>) (U e) (U f) = U (e+f)++instance Monoid (U var) where+  mempty = U 0+  mappend (U e) (U f) = U (e+f)+  mconcat us = U $ sum' [ e | U e <- us ]++#else++instance Monoid (U var) where+  mempty  = U 0+  mappend (U e) (U f) = U (e+f)+  mconcat us = U $ sum' [ e | U e <- us ]++#endif++--------------------------------------------------------------------------------++instance KnownSymbol var => Monomial (U var) where+  -- | the type of variables+  type VarM (U var) = ()+  +  -- checking the invariant+  normalizeM  = id+  isNormalM   = const True++  -- construction and deconstruction+  fromListM   ves = U $ sum' (map snd ves)+  toListM     (U e) = [((),e)]++  -- simple monomials+  emptyM      = U 0+  isEmptyM    (U e) = (e==0)+  variableM   _   = U 1+  singletonM  _ e = U e++  -- algebra+  mulM         = mappend+  productM     = mconcat+  divM (U e) (U f) = if e >= f then Just (U (e-f)) else Nothing+  powM (U e) k = U (k*e)++  -- degrees+  maxDegM     (U e) = e+  totalDegM   (U e) = e++  -- calculus+  diffM _ = diffU++  -- substitution and evaluation+  evalM       f (U e) = (f ())^e+  varSubsM    _ = id+  termSubsM   f (U e, c) = case f () of  +                Nothing  -> (U e, c      )+                (Just x) -> (U 0, c * x^e)++--------------------------------------------------------------------------------+-- * differentiation++diffU :: Num c => Int -> U v -> Maybe (U v, c)+diffU k (U m) =+  if k > m +    then Nothing+    else Just (U (m-k) , fromInteger c) +  where+    c = product [ fromIntegral (m-i) | i<-[0..k-1] ] :: Integer++--------------------------------------------------------------------------------
+ src/Math/Algebra/Polynomial/Multivariate/Compact.hs view
@@ -0,0 +1,154 @@++-- | Multivariate compact polynomials where the variable set +-- looks like @{x_1, x_2, ... , x_N}@.+--+-- This is very similar to the \"Indexed\" version, but should have much more+-- compact in-memory representation (which is useful in case of large or many +-- polynomials; and should be in theory also faster, because of cache-friendlyness)+--+--++{-# LANGUAGE BangPatterns, TypeFamilies, DataKinds, KindSignatures, ScopedTypeVariables, FlexibleContexts #-}+module Math.Algebra.Polynomial.Multivariate.Compact +  ( Poly(..) , unPoly , polyVar , nOfPoly , renamePolyVar+  , ZPoly , QPoly , fromZPoly, fromQPoly+  , Compact+  )+  where++--------------------------------------------------------------------------------++import Data.List+import Data.Word++import Data.Typeable+import GHC.TypeLits+import Data.Proxy+import Unsafe.Coerce as Unsafe++import Data.Foldable as F ++import qualified Math.Algebra.Polynomial.FreeModule as ZMod+import Math.Algebra.Polynomial.FreeModule ( FreeMod , FreeModule(..) ) -- , ZMod , QMod )++import Math.Algebra.Polynomial.Monomial.Compact++import Math.Algebra.Polynomial.Class+import Math.Algebra.Polynomial.Pretty+import Math.Algebra.Polynomial.Misc++--------------------------------------------------------------------------------++-- | A multivariate polynomial in with a given coefficient ring. +--+-- It is also indexed by the (shared) /name/ of the variables and the /number of/+-- variable. For example @Polyn Rational "x" 3@ the type of polynomials in the+-- variables @x1, x2, x3@ with rational coefficients.+newtype Poly (coeff :: *) (var :: Symbol) (n :: Nat) = Poly (FreeMod coeff (Compact var n))+  deriving (Eq,Ord,Show,Typeable)++unPoly :: Poly c v n -> FreeMod c (Compact v n)+unPoly (Poly x) = x++-- | Name of the variables+polyVar :: KnownSymbol var => Poly c var n -> String+polyVar = symbolVal . varProxy where+  varProxy :: Poly c var n -> Proxy var+  varProxy _ = Proxy++-- | Number of variables+nOfPoly :: KnownNat n => Poly c var n -> Int+nOfPoly = fromInteger . natVal . natProxy where+  natProxy :: Poly c var n -> Proxy n+  natProxy _ = Proxy++instance FreeModule (Poly c v n) where+  type BaseF  (Poly c v n) = Compact v n+  type CoeffF (Poly c v n) = c+  toFreeModule   = unPoly+  fromFreeModule = Poly++-- | Rename the variables (zero cost)+renamePolyVar :: Poly c var1 n -> Poly c var2 n+renamePolyVar = Unsafe.unsafeCoerce++--------------------------------------------------------------------------------++type ZPoly = Poly Integer+type QPoly = Poly Rational++-- | Change the coefficient ring (from integers)+fromZPoly :: (Ring c, KnownSymbol v, KnownNat n) => Poly Integer v n -> Poly c v n+fromZPoly= Poly . ZMod.fromZMod . unPoly++-- | Change the coefficient field (from rationals)+fromQPoly :: (Field c, KnownSymbol v, KnownNat n) => Poly Rational v n -> Poly c v n+fromQPoly = Poly . ZMod.fromQMod . unPoly++--------------------------------------------------------------------------------++instance (Ring c, KnownSymbol v, KnownNat n) => AlmostPolynomial (Poly c v n) where+  type CoeffP (Poly c v n) = c+  type MonomP (Poly c v n) = Compact v n+  type VarP   (Poly c v n) = Index++  zeroP         = Poly ZMod.zero+  isZeroP       = ZMod.isZero . unPoly+  oneP          = Poly (ZMod.generator emptyCompact)++  fromListP     = Poly . ZMod.fromList+  toListP       = ZMod.toList . unPoly++  variableP     = Poly . ZMod.generator . variableCompact+  singletonP    = \v e -> Poly (ZMod.generator (singletonCompact v e))+  monomP        = \m     -> Poly $ ZMod.generator m+  monomP'       = \m c   -> Poly $ ZMod.singleton m c+  scalarP       = \s     -> Poly $ ZMod.singleton emptyCompact s++  addP          = \p1 p2 -> Poly $ ZMod.add (unPoly p1) (unPoly p2)+  subP          = \p1 p2 -> Poly $ ZMod.sub (unPoly p1) (unPoly p2)+  negP          = Poly . ZMod.neg . unPoly+  mulP          = \p1 p2 -> Poly $ ZMod.mulWith     mulCompact (unPoly p1) (unPoly p2)++  coeffOfP      = \m p   -> ZMod.coeffOf m (unPoly p)+  productP      = \ps    -> Poly $ ZMod.productWith emptyCompact mulCompact $ map unPoly ps+  mulByMonomP   = \m p   -> Poly $ ZMod.mulByMonom  m (unPoly p)+  scaleP        = \s p   -> Poly $ ZMod.scale s (unPoly p) ++instance (Ring c, KnownSymbol v, KnownNat n) => Polynomial (Poly c v n) where+  evalP         = \g f p -> let { !z = evalM f ; h (!m,!c) = g c * z m } in sum' $ map h $ ZMod.toList $ unPoly p+  --varSubsP      = \f p   -> Poly $ ZMod.mapBase (varSubsCompact f) (unPoly p)+  --coeffSubsP    = \f p   -> Poly $ ZMod.fromList $ map (termSubsCompact f) $ ZMod.toList $ unPoly p +  --subsP         = \f p   -> Poly $ ZMod.flatMap (evalCompact (unPoly . f)) (unPoly p)+  varSubsP   = error "Compact/varSubsP: not yet implemented"+  coeffSubsP = error "Compact/coeffSubsP: not yet implemented"+  subsP      = error "Compact/subsP: not yet implemented"+  ++instance (Ring c, KnownSymbol v, KnownNat n) => Num (Poly c v n) where+  fromInteger = scalarP . fromInteger+  (+)    = addP+  (-)    = subP+  negate = negP+  (*)    = mulP+  abs    = id+  signum = \_ -> scalarP 1++instance (Ring c, KnownSymbol v, KnownNat n, Pretty (Compact v n)) => Pretty (Poly c v n) where+  pretty poly@(Poly fm) = if isSignedR (proxyCoeffP poly)+    then prettyFreeMod'  True   pretty fm+    else prettyFreeMod'' pretty pretty fm++-- hackety hack hack...+instance IsSigned (Poly c v n) where+  signOf = const (Just Plus)++-- So that we can use it again as a coefficient ring+instance (Ring c, KnownSymbol v, KnownNat n) => Ring (Poly c v n) where+  isZeroR   = ZMod.isZero . unPoly+  isAtomicR = const False+  isSignedR = const False+  absR      = id+  signumR   = const (Just Plus)++--------------------------------------------------------------------------------
+ src/Math/Algebra/Polynomial/Multivariate/Generic.hs view
@@ -0,0 +1,119 @@++-- | Multivariate polynomial algebra where the set of variables is given by +-- the inhabitants of a type++{-# LANGUAGE BangPatterns, TypeFamilies, KindSignatures, StandaloneDeriving, FlexibleContexts #-}+module Math.Algebra.Polynomial.Multivariate.Generic where++--------------------------------------------------------------------------------++import Data.Maybe+import Data.List+import Data.Foldable as F++import Data.Proxy+import Data.Typeable++import qualified Data.Map.Strict as Map ; import Data.Map.Strict (Map)++import qualified Math.Algebra.Polynomial.FreeModule as ZMod+import Math.Algebra.Polynomial.FreeModule ( FreeMod , FreeModule(..) ) -- , ZMod , QMod )++import Math.Algebra.Polynomial.Monomial.Generic ++import Math.Algebra.Polynomial.Class+-- import Math.Algebra.Polynomial.FreeModule+import Math.Algebra.Polynomial.Pretty+import Math.Algebra.Polynomial.Misc++--------------------------------------------------------------------------------++newtype Poly (coeff :: *) (var :: *)+  = Poly (FreeMod coeff (Monom var) )+  deriving (Eq,Ord,Show,Typeable)++unPoly :: Poly c v -> FreeMod c (Monom v) +unPoly (Poly p) = p++instance Ord v => FreeModule (Poly c v) where+  type BaseF  (Poly c v) = Monom v+  type CoeffF (Poly c v) = c+  toFreeModule   = unPoly+  fromFreeModule = Poly++--------------------------------------------------------------------------------++type ZPoly = Poly Integer+type QPoly = Poly Rational++-- | Change the coefficient ring (from integers)+fromZPoly :: (Ring c, Variable v) => Poly Integer v -> Poly c v +fromZPoly= Poly . ZMod.fromZMod . unPoly++-- | Change the coefficient ring (from rationals)+fromQPoly :: (Field c, Variable v) => Poly Rational v -> Poly c v +fromQPoly = Poly . ZMod.fromQMod . unPoly++--------------------------------------------------------------------------------++instance (Ring c, Ord v, Pretty v) => AlmostPolynomial (Poly c v) where+  type CoeffP (Poly c v) = c+  type MonomP (Poly c v) = Monom v+  type VarP   (Poly c v) = v++  fromListP     = Poly . ZMod.fromList+  toListP       = ZMod.toList . unPoly++  zeroP         = Poly ZMod.zero+  isZeroP       = ZMod.isZero . unPoly+  oneP          = Poly (ZMod.generator emptyMonom)++  variableP     = Poly . ZMod.generator . varMonom+  singletonP    = \v e -> Poly (ZMod.generator (singletonMonom v e))+  monomP        = \m     -> Poly $ ZMod.generator m+  monomP'       = \m c   -> Poly $ ZMod.singleton m c+  scalarP       = \s     -> Poly $ ZMod.singleton emptyMonom s++  addP          = \p1 p2 -> Poly $ ZMod.add (unPoly p1) (unPoly p2)+  subP          = \p1 p2 -> Poly $ ZMod.sub (unPoly p1) (unPoly p2)+  negP          = Poly . ZMod.neg . unPoly+  mulP          = \p1 p2 -> Poly $ ZMod.mulWith     mulMonom (unPoly p1) (unPoly p2)+  productP      = \ps    -> Poly $ ZMod.productWith emptyMonom mulMonom $ map unPoly ps++  coeffOfP      = \m p   -> ZMod.coeffOf m (unPoly p)+  mulByMonomP   = \m p   -> Poly $ ZMod.unsafeMulByMonom m (unPoly p)+  scaleP        = \s p   -> Poly $ ZMod.scale s (unPoly p) ++instance (Ring c, Ord v, Pretty v) => Polynomial (Poly c v) where+  evalP         = \g f p -> let { !z = evalM f ; h (!m,!c) = g c * z m } in sum' $ map h $ ZMod.toList $ unPoly p+  varSubsP      = \f p   -> Poly $ ZMod.mapBase (mapMonom f) (unPoly p)+  coeffSubsP    = \f p   -> Poly $ ZMod.fromList $ map (termSubsMonom f) $ ZMod.toList $ unPoly p +  subsP         = \f p   -> Poly $ ZMod.flatMap (evalMonom (unPoly . f)) (unPoly p)++instance (Ring c, Ord v, Pretty v) => Num (Poly c v) where+  fromInteger = scalarP . fromInteger+  (+)    = addP+  (-)    = subP+  negate = negP+  (*)    = mulP+  abs    = id+  signum = \_ -> scalarP 1++instance (Ring c, Ord v, Pretty v, Pretty (Monom v)) => Pretty (Poly c v) where+  pretty poly@(Poly fm) = if isSignedR (proxyCoeffP poly)+    then prettyFreeMod'  True   pretty fm+    else prettyFreeMod'' pretty pretty fm++-- hackety hack hack...+instance IsSigned (Poly c v) where+  signOf = const (Just Plus)++-- So that we can use it again as a coefficient ring+instance (Ring c, Variable v) => Ring (Poly c v) where+  isZeroR   = ZMod.isZero . unPoly+  isAtomicR = const False+  isSignedR = const False+  absR      = id+  signumR   = const (Just Plus)++--------------------------------------------------------------------------------
+ src/Math/Algebra/Polynomial/Multivariate/Indexed.hs view
@@ -0,0 +1,178 @@++-- | Multivariate polynomials where the variable set +-- looks like @{x_1, x_2, ... , x_N}@ ++{-# LANGUAGE +      BangPatterns, TypeFamilies, DataKinds, KindSignatures, ScopedTypeVariables,+      FlexibleContexts, StandaloneDeriving+  #-}+module Math.Algebra.Polynomial.Multivariate.Indexed+  (+    Poly(..) , unPoly , polyVar , nOfPoly , renamePolyVar+  , ZPoly , QPoly , fromZPoly, fromQPoly+  , embed+  , XS(..)+  )+  where++--------------------------------------------------------------------------------++import Data.Maybe+import Data.List+import Data.Array.Unboxed ++import Data.Typeable+import GHC.TypeLits+import Data.Proxy+import Unsafe.Coerce as Unsafe++import Data.Foldable as F ++import qualified Math.Algebra.Polynomial.FreeModule as ZMod+import Math.Algebra.Polynomial.FreeModule ( FreeMod , FreeModule(..) ) -- , ZMod , QMod )++import Math.Algebra.Polynomial.Monomial.Indexed ++import Math.Algebra.Polynomial.Class+import Math.Algebra.Polynomial.Pretty+import Math.Algebra.Polynomial.Misc++--------------------------------------------------------------------------------+-- * Polynomials++-- | A multivariate polynomial in with a given coefficient ring. +--+-- It is also indexed by the (shared) /name/ of the variables and the /number of/+-- variable. For example @Polyn Rational "x" 3@ the type of polynomials in the+-- variables @x1, x2, x3@ with rational coefficients.+newtype Poly (coeff :: *) (var :: Symbol) (n :: Nat) +  = Poly (FreeMod coeff (XS var n))+  deriving (Eq,Ord,Show,Typeable)++-- deriving instance (Ord coeff) => Ord (Poly coeff var n)++unPoly :: Poly c v n -> FreeMod c (XS v n) +unPoly (Poly p) = p++-- | Name of the variables+polyVar :: KnownSymbol var => Poly c var n -> String+polyVar = symbolVal . varProxy where+  varProxy :: Poly c var n -> Proxy var+  varProxy _ = Proxy++-- | Number of variables+nOfPoly :: KnownNat n => Poly c var n -> Int+nOfPoly = fromInteger . natVal . natProxy where+  natProxy :: Poly c var n -> Proxy n+  natProxy _ = Proxy++instance FreeModule (Poly c v n) where+  type BaseF  (Poly c v n) = XS v n+  type CoeffF (Poly c v n) = c+  toFreeModule   = unPoly+  fromFreeModule = Poly++renamePolyVar :: Poly c var1 n -> Poly c var2 n+renamePolyVar = Unsafe.unsafeCoerce++--------------------------------------------------------------------------------++type ZPoly = Poly Integer+type QPoly = Poly Rational++-- | Change the coefficient ring (from integers)+fromZPoly :: (Ring c, KnownSymbol v, KnownNat n) => Poly Integer v n -> Poly c v n+fromZPoly= Poly . ZMod.fromZMod . unPoly++-- | Change the coefficient field (from rationals)+fromQPoly :: (Field c, KnownSymbol v, KnownNat n) => Poly Rational v n -> Poly c v n+fromQPoly = Poly . ZMod.fromQMod . unPoly++--------------------------------------------------------------------------------++instance (Ring c, KnownSymbol v, KnownNat n) => AlmostPolynomial (Poly c v n) where+  type CoeffP (Poly c v n) = c+  type MonomP (Poly c v n) = XS v n+  type VarP   (Poly c v n) = Index++  zeroP         = Poly ZMod.zero+  isZeroP       = ZMod.isZero . unPoly+  oneP          = Poly (ZMod.generator emptyXS)++  fromListP     = Poly . ZMod.fromList+  toListP       = ZMod.toList . unPoly++  variableP     = Poly . ZMod.generator . variableXS+  singletonP    = \v e -> Poly (ZMod.generator (singletonXS v e))+  monomP        = \m     -> Poly $ ZMod.generator m+  monomP'       = \m c   -> Poly $ ZMod.singleton m c+  scalarP       = \s     -> Poly $ ZMod.singleton emptyXS s++  addP          = \p1 p2 -> Poly $ ZMod.add (unPoly p1) (unPoly p2)+  subP          = \p1 p2 -> Poly $ ZMod.sub (unPoly p1) (unPoly p2)+  negP          = Poly . ZMod.neg . unPoly+  mulP          = \p1 p2 -> Poly $ ZMod.mulWith mulXS (unPoly p1) (unPoly p2)++  coeffOfP      = \m p   -> ZMod.coeffOf m (unPoly p)+  productP      = \ps    -> Poly $ ZMod.productWith emptyXS mulXS $ map unPoly ps+  mulByMonomP   = \m p   -> Poly $ ZMod.unsafeMulByMonom m (unPoly p)+  scaleP        = \s p   -> Poly $ ZMod.scale s (unPoly p) ++instance (Ring c, KnownSymbol v, KnownNat n) => Polynomial (Poly c v n) where+  evalP         = \g f p -> let { !z = evalM f ; h (!m,!c) = g c * z m } in sum' $ map h $ ZMod.toList $ unPoly p+  varSubsP      = \f p   -> Poly $ ZMod.mapBase (varSubsXS f) (unPoly p)+  coeffSubsP    = \f p   -> Poly $ ZMod.fromList $ map (termSubsXS f) $ ZMod.toList $ unPoly p +  subsP         = \f p   -> Poly $ ZMod.flatMap (evalXS (unPoly . f)) (unPoly p)++instance (Ring c, KnownSymbol v, KnownNat n) => Num (Poly c v n) where+  fromInteger = scalarP . fromInteger+  (+)    = addP+  (-)    = subP+  negate = negP+  (*)    = mulP+  abs    = id+  signum = \_ -> scalarP 1++instance (Ring c, KnownSymbol v, KnownNat n, Pretty (XS v n)) => Pretty (Poly c v n) where+  pretty poly@(Poly fm) = if isSignedR (proxyCoeffP poly)+    then prettyFreeMod'  True   pretty fm+    else prettyFreeMod'' pretty pretty fm++-- hackety hack hack...+instance IsSigned (Poly c v n) where+  signOf = const (Just Plus)++-- So that we can use it again as a coefficient ring+instance (Ring c, KnownSymbol v, KnownNat n) => Ring (Poly c v n) where+  isZeroR   = ZMod.isZero . unPoly+  isAtomicR = const False+  isSignedR = const False+  absR      = id+  signumR   = const (Just Plus)++--------------------------------------------------------------------------------++-- | You can always increase the number of variables; +-- you can also decrease, if don't use the last few ones.+--+-- This will throw an error if you try to eliminate variables which are in fact used.+-- To do that, you can instead substitute 0 or 1 into them.+--+embed :: (Ring c, KnownNat n, KnownNat m) => Poly c v n -> Poly c v m+embed old@(Poly old_fm) = new where+  n = nOfPoly old+  m = nOfPoly new+  new = Poly $ case compare m n of +    LT -> ZMod.unsafeMapBase project old_fm+    EQ -> ZMod.unsafeMapBase keep    old_fm+    GT -> ZMod.unsafeMapBase extend  old_fm+  extend  (XS arr) = XS $ listArray (1,m) (elems arr ++ replicate (m-n) 0)+  keep    (XS arr) = XS arr+  project (XS arr) = +    let old = elems arr+        (new,rest) = splitAt m old+    in  if all (==0) rest +          then XS $ listArray (1,m) new+          else error "Indexed/embed: the projected variables are actually used!"++--------------------------------------------------------------------------------
+ src/Math/Algebra/Polynomial/Multivariate/Infinite.hs view
@@ -0,0 +1,165 @@++-- | Multivariate polynomials where the variable set is the countably infinite+-- set @{x_1, x_2, x_3, ...}@ ++{-# LANGUAGE +      BangPatterns, TypeFamilies, DataKinds, KindSignatures, ScopedTypeVariables,+      FlexibleContexts+  #-}+module Math.Algebra.Polynomial.Multivariate.Infinite+  (+    Poly(..) , unPoly , polyVar , renamePolyVar+  , ZPoly , QPoly , fromZPoly, fromQPoly+  , truncate+  , XInf(..)+  )+  where++--------------------------------------------------------------------------------++import Prelude hiding ( truncate )++import Data.Maybe+import Data.List+import Data.Array.Unboxed ++import Data.Typeable+import GHC.TypeLits+import Data.Proxy+import Unsafe.Coerce as Unsafe++import Data.Foldable as F ++import qualified Math.Algebra.Polynomial.FreeModule as ZMod+import Math.Algebra.Polynomial.FreeModule ( FreeMod , FreeModule(..) ) -- , ZMod , QMod )++import Math.Algebra.Polynomial.Monomial.Infinite++import Math.Algebra.Polynomial.Class+import Math.Algebra.Polynomial.Pretty+import Math.Algebra.Polynomial.Misc++import qualified Math.Algebra.Polynomial.Monomial.Indexed     as Fin+import qualified Math.Algebra.Polynomial.Multivariate.Indexed as Fin++--------------------------------------------------------------------------------+-- * Polynomials++-- | A multivariate polynomial in with a given coefficient ring. +--+-- It is also indexed by the (shared) /name/ of the variables and the /number of/+-- variable. For example @Polyn Rational "x" 3@ the type of polynomials in the+-- variables @x1, x2, x3@ with rational coefficients.+newtype Poly (coeff :: *) (var :: Symbol) +  = Poly (FreeMod coeff (XInf var))+  deriving (Eq,Ord,Show,Typeable)++unPoly :: Poly c v -> FreeMod c (XInf v) +unPoly (Poly p) = p++-- | Name of the variables+polyVar :: KnownSymbol var => Poly c var -> String+polyVar = symbolVal . varProxy where+  varProxy :: Poly c var -> Proxy var+  varProxy _ = Proxy++instance FreeModule (Poly c v) where+  type BaseF  (Poly c v) = XInf v +  type CoeffF (Poly c v) = c+  toFreeModule   = unPoly+  fromFreeModule = Poly++renamePolyVar :: Poly c var1 -> Poly c var2 +renamePolyVar = Unsafe.unsafeCoerce++--------------------------------------------------------------------------------++type ZPoly = Poly Integer+type QPoly = Poly Rational++-- | Change the coefficient ring (from integers)+fromZPoly :: (Ring c, KnownSymbol v) => Poly Integer v -> Poly c v +fromZPoly= Poly . ZMod.fromZMod . unPoly++-- | Change the coefficient field (from rationals)+fromQPoly :: (Field c, KnownSymbol v) => Poly Rational v -> Poly c v +fromQPoly = Poly . ZMod.fromQMod . unPoly++--------------------------------------------------------------------------------++instance (Ring c, KnownSymbol v) => AlmostPolynomial (Poly c v) where+  type CoeffP (Poly c v) = c+  type MonomP (Poly c v) = XInf v+  type VarP   (Poly c v) = Index++  zeroP         = Poly ZMod.zero+  isZeroP       = ZMod.isZero . unPoly+  oneP          = Poly (ZMod.generator emptyM)++  fromListP     = Poly . ZMod.fromList+  toListP       = ZMod.toList . unPoly++  variableP     = Poly . ZMod.generator . variableXInf+  singletonP    = \v e -> Poly (ZMod.generator (singletonXInf v e))+  monomP        = \m     -> Poly $ ZMod.generator m+  monomP'       = \m c   -> Poly $ ZMod.singleton m c+  scalarP       = \s     -> Poly $ ZMod.singleton emptyXInf s++  addP          = \p1 p2 -> Poly $ ZMod.add (unPoly p1) (unPoly p2)+  subP          = \p1 p2 -> Poly $ ZMod.sub (unPoly p1) (unPoly p2)+  negP          = Poly . ZMod.neg . unPoly+  mulP          = \p1 p2 -> Poly $ ZMod.mulWith     mulXInf (unPoly p1) (unPoly p2)++  coeffOfP      = \m p   -> ZMod.coeffOf m (unPoly p)+  productP      = \ps    -> Poly $ ZMod.productWith emptyXInf mulXInf $ map unPoly ps+  mulByMonomP   = \m p   -> Poly $ ZMod.mulByMonom  m (unPoly p)+  scaleP        = \s p   -> Poly $ ZMod.scale s (unPoly p) ++instance (Ring c, KnownSymbol v) => Polynomial (Poly c v) where+  evalP         = \g f p -> let { !z = evalM f ; h (!m,!c) = g c * z m } in sum' $ map h $ ZMod.toList $ unPoly p+  varSubsP      = \f p   -> Poly $ ZMod.mapBase (varSubsM f) (unPoly p)+  coeffSubsP    = \f p   -> Poly $ ZMod.fromList $ map (termSubsM f) $ ZMod.toList $ unPoly p +  subsP         = \f p   -> Poly $ ZMod.flatMap (evalM (unPoly . f)) (unPoly p)++instance (Ring c, KnownSymbol v) => Num (Poly c v) where+  fromInteger = scalarP . fromInteger+  (+)    = addP+  (-)    = subP+  negate = negP+  (*)    = mulP+  abs    = id+  signum = \_ -> scalarP 1++instance (Ring c, KnownSymbol v, Pretty (XInf v)) => Pretty (Poly c v) where+  pretty poly@(Poly fm) = if isSignedR (proxyCoeffP poly)+    then prettyFreeMod'  True   pretty fm+    else prettyFreeMod'' pretty pretty fm++-- hackety hack hack...+instance IsSigned (Poly c v) where+  signOf = const (Just Plus)++-- So that we can use it again as a coefficient ring+instance (Ring c, KnownSymbol v) => Ring (Poly c v) where+  isZeroR   = ZMod.isZero . unPoly+  isAtomicR = const False+  isSignedR = const False+  absR      = id+  signumR   = const (Just Plus)++--------------------------------------------------------------------------------++-- | We can always truncate to a given number of variables, simply +-- by substituting zero to the rest+truncate :: (Eq c, Num c, KnownNat n) => Poly c v -> Fin.Poly c v n+truncate input@(Poly inpoly) = output where+  n = Fin.nOfPoly output+  output = Fin.Poly $ ZMod.mapMaybeBase f inpoly+  f (XInf es) =+      let (fs,rest) = splitAt n es+      in  if all (==0) rest +            then Just (Fin.xsFromExponents fs) +            else Nothing ++--------------------------------------------------------------------------------+
+ src/Math/Algebra/Polynomial/Pretty.hs view
@@ -0,0 +1,171 @@++{-# LANGUAGE FlexibleInstances #-}++-- | Pretty-printing.+--+-- Tip: you can try putting something like this into your @.ghci@ file to+-- make life more convenient:+--+-- > :m    +Math.Algebra.Polynomial.Pretty  +-- > :seti -interactive-print=prettyPrint+--+ +module Math.Algebra.Polynomial.Pretty where++--------------------------------------------------------------------------------++import Data.List+import Data.Ratio++import Math.Algebra.Polynomial.FreeModule ( FreeMod, ZMod, QMod )+import qualified Math.Algebra.Polynomial.FreeModule as ZMod++import Math.Algebra.Polynomial.Misc++--------------------------------------------------------------------------------++class Pretty a where+  pretty :: a -> String++  prettyInParens :: a -> String+  prettyInParens = pretty++prettyPrint :: Pretty a => a -> IO ()+prettyPrint = putStrLn . pretty++--------------------------------------------------------------------------------++-- instance Pretty a => Pretty (ZMod a) where+--   pretty = prettyZMod pretty++instance (Num c, Eq c, Pretty c, IsSigned c, Pretty b) => Pretty (FreeMod c b) where+  pretty = prettyFreeMod' True pretty+  prettyInParens x = "(" ++ pretty x ++ ")"++--------------------------------------------------------------------------------++instance Pretty Int where+  pretty = show++instance Pretty Integer where+  pretty = show++instance (Eq a, Num a, Pretty a) => Pretty (Ratio a) where+  pretty q = case denominator q of+    1 -> prettyInParens (numerator q)+    _ -> prettyInParens (numerator q) ++ "/" ++ prettyInParens (denominator q)++--------------------------------------------------------------------------------+-- * Pretty printing elements of free modules++-- | Example: @showVarPower "x" 5 == "x^5"@+showVarPower :: String -> Int -> String+showVarPower name expo = case expo of+  0 -> "1"+  1 -> name+  _ -> name ++ "^" ++ show expo++--------------------------------------------------------------------------------++-- | no multiplication sign (ok for mathematica and humans)+prettyZMod_ :: (b -> String) -> ZMod b -> String+prettyZMod_ = prettyFreeMod' False+  +-- | multiplication sign (ok for maple etc)+prettyZMod :: (b -> String) -> ZMod b -> String+prettyZMod = prettyFreeMod' True++--------------------------------------------------------------------------------++prettyFreeMod' +  :: (Num c, Eq c, IsSigned c, Pretty c) +  => Bool                -- ^ use star for multiplication (@False@ means just concatenation)+  -> (b -> String)       -- ^ show base+  -> FreeMod c b +  -> String+prettyFreeMod' star showBase what = final where+  final = if take 3 stuff == " + " then drop 3 stuff else drop 1 stuff+  stuff = concatMap f (ZMod.toList what) +  f (g,  1) = plus  ++ showBase' g+  f (g, -1) = minus ++ showBase' g+  f (g, c)  = case showBase' g of+    "1" -> sgn c ++ (prettyInParens $ abs c)+    b   -> sgn c ++ (prettyInParens $ abs c) ++ starStr ++ b+  -- cond (_,c) = (c/=0)+  starStr = if star then "*" else " "+  showBase' g = case showBase g of+    "" -> "1"  -- "(1)"+    s  -> s+  sgn c = case signOf c of+    Just Minus -> minus+    _          -> plus+  plus  = " + "+  minus = " - "++prettyFreeMod'' +  :: (c -> String)    -- ^ show coefficient+  -> (b -> String)    -- ^ show base+  -> FreeMod c b +  -> String+prettyFreeMod'' showCoeff showBase what = result where+  result = intercalate " + " (map f $ ZMod.toList what) +  f (g, c) = case showBase g of+    ""  -> showCoeff c+    "1" -> showCoeff c+    b   -> showCoeff c ++ starStr ++ b+  starStr = "*"++--------------------------------------------------------------------------------++{-+-- * Utility++-- | Put into parentheses+paren :: String -> String+paren s = '(' : s ++ ")"++--------------------------------------------------------------------------------++{-++-- | Exponential form of a partition+expFormString :: Partition -> String+expFormString p = "(" ++ intercalate "," (map f ies) ++ ")" where+  ies = toExponentialForm p+  f (i,e) = show i ++ "^" ++ show e+-}++extendStringL :: Int -> String -> String+extendStringL k s = s ++ replicate (k - length s) ' '++extendStringR :: Int -> String -> String+extendStringR k s = replicate (k - length s) ' ' ++ s++--------------------------------------------------------------------------------+-- * Mathematica-formatted output++class Mathematica a where+  mathematica :: a -> String++instance Mathematica Int where+  mathematica = show++instance Mathematica Integer where+  mathematica = show++instance Mathematica String where+  mathematica = show++{-+instance Mathematica Partition where+  mathematica (Partition ps) = "{" ++ intercalate "," (map show ps) ++ "}"+-}++data Indexed a = Indexed String a++instance Mathematica a => Mathematica (Indexed a) where+  mathematica (Indexed x sub) = "Subscript[" ++ x ++ "," ++ mathematica sub ++ "]"++--------------------------------------------------------------------------------++-}
+ src/Math/Algebra/Polynomial/Univariate.hs view
@@ -0,0 +1,154 @@++-- | Univariate polynomials++{-# LANGUAGE BangPatterns, DataKinds, KindSignatures, GeneralizedNewtypeDeriving, TypeFamilies #-}+module Math.Algebra.Polynomial.Univariate+  ( -- * Univariate polynomials+    Univariate(..) ,  U(..) , unUni , uniVar , renameUniVar+  , ZUni , QUni , fromZUni , fromQUni+  , differentiateUni , integrateUni , integrateUni'+  )+  where++--------------------------------------------------------------------------------++import Data.Array ( Array , (!) , listArray , assocs ) +import Data.List++import GHC.TypeLits+import Data.Proxy+import Unsafe.Coerce as Unsafe++import Math.Algebra.Polynomial.Class+import Math.Algebra.Polynomial.Misc+import Math.Algebra.Polynomial.Pretty++import qualified Math.Algebra.Polynomial.FreeModule as ZMod+import Math.Algebra.Polynomial.FreeModule ( FreeMod , FreeModule(..) , ZMod , QMod )++import Math.Algebra.Polynomial.Monomial.Univariate++--------------------------------------------------------------------------------+-- * Univariate polynomials++-- | A univariate polynomial with the given coefficient ring. Note: this +-- is also indexed by the /name/ of the variable.+newtype Univariate (coeff :: *) (var :: Symbol) = Uni (FreeMod coeff (U var))+  deriving (Eq,Ord,Show)++unUni :: Univariate c v -> FreeMod c (U v)+unUni (Uni a) = a++instance FreeModule (Univariate c v) where+  type BaseF  (Univariate c v) = U v +  type CoeffF (Univariate c v) = c+  toFreeModule   = unUni+  fromFreeModule = Uni++-- | Name of the variable+uniVar :: KnownSymbol var => Univariate c var -> String+uniVar = symbolVal . varProxy where+  varProxy :: Univariate c var -> Proxy var+  varProxy _ = Proxy++-- | Rename the variable (zero cost)+renameUniVar :: Univariate c var1 -> Univariate c var2+renameUniVar = Unsafe.unsafeCoerce++--------------------------------------------------------------------------------++-- | An univariate polynomial integer coefficients+type ZUni var = Univariate Integer var++-- | An univariate polynomial with rational coefficients+type QUni var = Univariate Rational var++-- | Change the coefficient ring+fromZUni :: (Ring c, KnownSymbol v) => Univariate Integer v -> Univariate c v +fromZUni = Uni . ZMod.fromZMod . unUni++-- | Change the coefficient ring+fromQUni :: (Field c, KnownSymbol v) => Univariate Rational v -> Univariate c v +fromQUni = Uni . ZMod.fromQMod . unUni++--------------------------------------------------------------------------------++-- | Differentiation+differentiateUni :: (Ring c, KnownSymbol var) => Univariate c var -> Univariate c var+differentiateUni = Uni . ZMod.mapMaybeBaseCoeff f . unUni where+  f (U k) = case k of+    0 -> Nothing+    _ -> Just ( U (k-1) , fromIntegral k )++-- | Integration+integrateUni :: (Field c, KnownSymbol var) => Univariate c var -> Univariate c var+integrateUni = Uni . ZMod.mapMaybeBaseCoeff f . unUni where+  f (U k) = Just ( U (k+1) , 1 / fromIntegral (k+1) )++integrateUni' :: (Field c, KnownSymbol var) => c -> Univariate c var -> Univariate c var+integrateUni' c0 p = integrateUni p + scalarP c0++--------------------------------------------------------------------------------++instance (Ring coeff, KnownSymbol var) => AlmostPolynomial (Univariate coeff var) where+                                          +  type CoeffP (Univariate coeff var) = coeff+  type MonomP (Univariate coeff var) = U var+  type VarP   (Univariate coeff var) = ()++  fromListP     = Uni . ZMod.fromList+  toListP       = ZMod.toList . unUni++  zeroP         = Uni ZMod.zero+  isZeroP       = ZMod.isZero . unUni+  oneP          = Uni (ZMod.generator emptyM)++  variableP     = Uni . ZMod.generator . variableM+  singletonP    = \v e -> Uni (ZMod.generator (singletonM v e))+  monomP        = \m     -> Uni $ ZMod.generator m+  monomP'       = \m c   -> Uni $ ZMod.singleton m c+  scalarP       = \s     -> Uni $ ZMod.singleton emptyM s++  addP          = \p1 p2 -> Uni $ ZMod.add (unUni p1) (unUni p2)+  subP          = \p1 p2 -> Uni $ ZMod.sub (unUni p1) (unUni p2)+  negP          = Uni . ZMod.neg . unUni+  mulP          = \p1 p2 -> Uni $ ZMod.mulWith mulM (unUni p1) (unUni p2)+  productP      = \ps    -> Uni $ ZMod.productWith emptyM mulM $ map unUni ps++  coeffOfP      = \m p   -> ZMod.coeffOf m (unUni p)+  mulByMonomP   = \m p   -> Uni $ ZMod.mulByMonom  m (unUni p)+  scaleP        = \s p   -> Uni $ ZMod.scale s (unUni p) ++instance (Ring coeff, KnownSymbol var) => Polynomial (Univariate coeff var) where+  evalP         = \g f p -> let { !z = evalM f ; h (!m,!c) = g c * z m } in sum' $ map h $ ZMod.toList $ unUni p+  varSubsP      = \f p   -> Uni $ ZMod.mapBase (varSubsM f) (unUni p)+  coeffSubsP    = \f p   -> Uni $ ZMod.fromList $ map (termSubsM f) $ ZMod.toList $ unUni p +  subsP         = \f p   -> Uni $ ZMod.flatMap (evalM (unUni . f)) (unUni p)++instance (Ring c, KnownSymbol v) => Num (Univariate c v) where+  fromInteger = scalarP . fromInteger+  (+)    = addP+  (-)    = subP+  negate = negP+  (*)    = mulP+  abs    = id+  signum = \_ -> scalarP 1++instance (Ring c, KnownSymbol v) => Pretty (Univariate c v) where+  pretty poly@(Uni fm) = if isSignedR (proxyCoeffP poly)+    then prettyFreeMod'  True   pretty fm+    else prettyFreeMod'' pretty pretty fm++-- hackety hack hack...+instance IsSigned (Univariate c v) where+  signOf = const (Just Plus)++-- So that we can use it again as a coefficient ring+instance (Ring c, KnownSymbol v) => Ring (Univariate c v) where+  isZeroR   = ZMod.isZero . unUni+  isAtomicR = const False+  isSignedR = const False+  absR      = id+  signumR   = const (Just Plus)++--------------------------------------------------------------------------------
+ src/Math/Algebra/Polynomial/Univariate/Bernoulli.hs view
@@ -0,0 +1,135 @@++-- | Bernoulli and Euler polynomials+--+-- See <https://en.wikipedia.org/wiki/Bernoulli_polynomials>+-- ++{-# LANGUAGE DataKinds, TypeSynonymInstances, FlexibleContexts, FlexibleInstances, BangPatterns, ScopedTypeVariables #-}+module Math.Algebra.Polynomial.Univariate.Bernoulli+  ( bernoulliB, eulerE+  , rationalBernoulliB, rationalEulerE+  , bernoulliNumber, signedEulerNumber, unsignedEulerNumber+  , eulerianPolynomial+  ) +  where++--------------------------------------------------------------------------------++import Data.List+import Data.Ratio++import Data.Semigroup+import Data.Monoid++import GHC.TypeLits++import qualified Math.Algebra.Polynomial.FreeModule as ZMod+import Math.Algebra.Polynomial.FreeModule ( FreeMod , FreeModule(..) , ZMod , QMod )++import Math.Algebra.Polynomial.Univariate++import Math.Algebra.Polynomial.Class+import Math.Algebra.Polynomial.Pretty+import Math.Algebra.Polynomial.Misc++--------------------------------------------------------------------------------+-- * Bernoulli polynomials++-- | Bernoulli polynomials+bernoulliB :: (Field c, KnownSymbol v) => Int -> Univariate c v+bernoulliB = fromQUni . renameUniVar . rationalBernoulliB++-- | Euler polynomials (not to be confused with the related Eulerian polynomials!)+eulerE :: (Field c, KnownSymbol v) => Int -> Univariate c v+eulerE = fromQUni . renameUniVar . rationalEulerE++--------------------------------------------------------------------------------++rationalBernoulliB :: Int -> Univariate Rational "x"+rationalBernoulliB n = Uni $ ZMod.fromList+  [ (U k , fromInteger (binomial n k) * bernoulliNumber (n-k) ) | k<-[0..n] ]++rationalEulerE :: Int -> Univariate Rational "x"+rationalEulerE n = sumP+  [ scaleP coeff $ xMinusHalfPowN (n-k)+  | k <- [0,2..n]+  , let coeff = fromInteger (binomial n k * eulerNumber k) / 2^k +  ]+  where+    eulerNumber k = if even k +      then signedEulerNumber (div k 2) +      else 0++xMinusHalfPowN :: Int -> Univariate Rational "x"+xMinusHalfPowN = intCache compute where+  x_minus_half = Uni $ ZMod.fromList [ (U 1 , 1) , (U 0 , -1/2) ]+  compute recur n = case n of +    0 -> 1+    n -> x_minus_half * recur (n-1)++--------------------------------------------------------------------------------+-- * Bernoulli numbers++-- | Bernoulli numbers. @bernoulli 1 == -1%2@ and @bernoulli k == 0@ for+-- k>2 and /odd/. This function uses the formula involving Stirling numbers+-- of the second kind. Numerators: A027641, denominators: A027642.+bernoulliNumber :: Integral a => a -> Rational+bernoulliNumber n +  | n <  0    = error "bernoulli: n should be nonnegative"+  | n == 0    = 1+  | n == 1    = -1/2+  | otherwise = sum [ f k | k<-[1..n] ] +  where+    f k = toRational (negateIfOdd (n+k) $ factorial k * stirling2nd n k) +        / toRational (k+1)++--------------------------------------------------------------------------------+-- * Euler numbers++{-+x, oneminusx, xoneminusx, halfoneminusx :: Univariate Rational "x"+x = variableP ()+oneminusx = 1 - x+xoneminusx = x * (1 - x)+halfoneminusx = scaleP (1/2) oneminusx+nx :: Int -> Univariate Rational "x"+nx n = monomP' (U 1) (fromIntegral n)++scaledEulerianPolynomial :: Int -> Univariate Rational "x"+scaledEulerianPolynomial = intCache compute where+  compute recur n = case n of +    0 -> 1+    n -> (nx n + halfoneminusx) * recur (n-1) + xoneminusx * differentiateUni (recur (n-1))+-}++x, oneminusx, two_xoneminusx :: Univariate Integer "x"+x = variableP ()+oneminusx = 1 - x+two_xoneminusx = 2 * x * (1 - x)+two_nx :: Int -> Univariate Integer "x"+two_nx n = monomP' (U 1) (2 * fromIntegral n)++-- | Eulerian polynomials (row polynomials of A060187)+eulerianPolynomial :: Int -> Univariate Integer "x"+eulerianPolynomial = intCache compute where+  compute recur n = case n of +    0 -> 1+    n -> (two_nx n + oneminusx) * recur (n-1) + two_xoneminusx * differentiateUni (recur (n-1))++-- | Signed Euler numbers (unsigned version: A000364)+--+-- See <https://en.wikipedia.org/wiki/Euler_number>+--+-- NOTE: we skip the zeros (every other index)+signedEulerNumber :: Int -> Integer+signedEulerNumber = intCache compute where+  compute _ n = div (evalP id (\_ -> -1) $ eulerianPolynomial (2*n)) (4^n)++-- | unsigned Euler numbers (A000364)+--+-- NOTE: we skip the zeros (every other index)+unsignedEulerNumber :: Int -> Integer+unsignedEulerNumber n = negateIfOdd n $ signedEulerNumber n++--------------------------------------------------------------------------------+
+ src/Math/Algebra/Polynomial/Univariate/Chebysev.hs view
@@ -0,0 +1,63 @@++-- | Chebysev polynomials+-- +-- See <https://en.wikipedia.org/wiki/Chebyshev_polynomials>++{-# LANGUAGE DataKinds, TypeSynonymInstances, FlexibleContexts, FlexibleInstances, BangPatterns, ScopedTypeVariables #-}+module Math.Algebra.Polynomial.Univariate.Chebysev+  ( chebysevT+  , chebysevU+  , integralChebysevT +  , integralChebysevU+  ) +  where++--------------------------------------------------------------------------------++import Data.List+import Data.Ratio++import Data.Semigroup+import Data.Monoid++import GHC.TypeLits++import qualified Math.Algebra.Polynomial.FreeModule as ZMod+import Math.Algebra.Polynomial.FreeModule ( FreeMod , FreeModule(..) , ZMod , QMod )++import Math.Algebra.Polynomial.Univariate++import Math.Algebra.Polynomial.Class+import Math.Algebra.Polynomial.Pretty+import Math.Algebra.Polynomial.Misc++--------------------------------------------------------------------------------++chebysevT :: (Ring c, KnownSymbol v) => Int -> Univariate c v+chebysevT = fromZUni . renameUniVar . integralChebysevT++chebysevU :: (Ring c, KnownSymbol v) => Int -> Univariate c v+chebysevU = fromZUni . renameUniVar . integralChebysevU++--------------------------------------------------------------------------------++x, twox :: Univariate Integer "x"+x    = variableP ()+twox = monomP'   (U 1) 2++integralChebysevT :: Int -> Univariate Integer "x"+integralChebysevT = intCache compute where+  compute recur n = case n of +    0 -> 1+    1 -> x+    n -> twox * recur (n-1) - recur (n-2)++integralChebysevU :: Int -> Univariate Integer "x"+integralChebysevU = intCache compute where+  compute recur n = case n of +    0 -> 1+    1 -> twox+    n -> twox * recur (n-1) - recur (n-2)++--------------------------------------------------------------------------------+
+ src/Math/Algebra/Polynomial/Univariate/Cyclotomic.hs view
@@ -0,0 +1,155 @@++-- | Cyclotomic polynomials+-- +-- See <https://en.wikipedia.org/wiki/Cyclotomic_polynomial>+--++{-# LANGUAGE DataKinds, TypeSynonymInstances, FlexibleContexts, FlexibleInstances, BangPatterns, ScopedTypeVariables #-}+module Math.Algebra.Polynomial.Univariate.Cyclotomic +  ( cyclotomic+  , cyclotomicMoebius+  , cyclotomicNaive+  ) +  where++--------------------------------------------------------------------------------++import Data.Ratio++import Data.Semigroup+import Data.Monoid++import qualified Math.Algebra.Polynomial.FreeModule as ZMod+import Math.Algebra.Polynomial.FreeModule ( FreeMod , FreeModule(..) , ZMod , QMod )++import Math.Algebra.Polynomial.Univariate++import Math.Algebra.Polynomial.Class+import Math.Algebra.Polynomial.Pretty+import Math.Algebra.Polynomial.Misc++--------------------------------------------------------------------------------++-- | test whether the product of cyclotomic polynomials @Phi_d$ for the divisors @d@ of @n@ is @x^n-1@+test_product :: Int -> Bool+test_product n = lhs == rhs where+  lhs = product [ cyclotomic d | d <- divisors n ] +  rhs = Uni $ ZMod.fromList [ (U n , 1) , (U 0 , -1) ]++-- | Synonym to 'cyclotomicMoebius'+cyclotomic :: Int -> Univariate Integer "x" +cyclotomic = cyclotomicMoebius ++--------------------------------------------------------------------------------+-- * Implementation via Moebius inversion++-- | Cyclotomic polynomials via Moebius inversion+cyclotomicMoebius :: Int -> Univariate Integer "x"+cyclotomicMoebius = Uni . ZMod.unsafeZModFromQMod . unUni . cyclotomicMoebiusQ++cyclotomicMoebiusQ :: Int -> Univariate Rational "x" +cyclotomicMoebiusQ n = unsafeDivide (Uni numer) (Uni denom) where+  sqfd   = squareFreeDivisors n +  numer  = ZMod.product [ term d | (d,Plus ) <- sqfd ]+  denom  = ZMod.product [ term d | (d,Minus) <- sqfd ]+  term d = ZMod.fromList [ (U (div n d) , 1) , (U 0 , -1) ]++polynomialLongDivision :: forall p. (Polynomial p, Fractional (CoeffP p)) => p -> p -> (p,p)+polynomialLongDivision p0 q = go zeroP p0 where++  (bq,cq) = case findMaxTerm q of+    Just bc -> bc+    Nothing -> error "polynomialLongDivision: division by zero"++  go !acc !p = case findMaxTerm p of+    Nothing      -> (acc,zeroP)+    Just (bp,cp) -> case divM bp bq of+      Nothing      -> (acc,p)+      Just br      -> let cr = (cp/cq)+                          u  = scaleP cr (mulByMonomP br q)+                          p' = p - u+                          acc' = (acc + monomP' br cr) +                      in  go acc' p' ++divideMaybe :: (Polynomial p, Fractional (CoeffP p)) => p -> p -> Maybe p+divideMaybe p q = case polynomialLongDivision p q of+  (s,r) -> if isZeroP r then Just s else Nothing++unsafeDivide :: (Polynomial p, Fractional (CoeffP p)) => p -> p -> p+unsafeDivide p q = case divideMaybe p q of+  Just s  -> s+  Nothing -> error "Polynomial/unsafeDivide: not divisible"++findMaxTerm :: FreeModule f => f -> Maybe (BaseF f, CoeffF f)+findMaxTerm = ZMod.findMaxTerm . toFreeModule++--------------------------------------------------------------------------------+-- * Naive implementation++newtype E2PiI = E2PiI Rational deriving (Eq,Ord,Show)++instance Pretty E2PiI where+  pretty (E2PiI y) = "e^{2*pi*i*" ++ show p ++ "/" ++ show q ++ "}" where+    p = numerator   y+    q = denominator y++mod1 :: Rational -> Rational+mod1 x = x - fromInteger (floor x)++mulE2PiI :: E2PiI -> E2PiI -> E2PiI+mulE2PiI (E2PiI x) (E2PiI y) = E2PiI (mod1 $ x+y)++instance Semigroup E2PiI where+  (<>) = mulE2PiI++instance Monoid E2PiI where+  mempty  = E2PiI 0+  mappend = mulE2PiI++half :: Rational+half = 1/2++reduce :: ZMod E2PiI -> Either (ZMod E2PiI) Integer+reduce input = +  case ZMod.findMaxTerm input of+    Nothing           -> Right 0+    Just (E2PiI y, c) +      | y == 0        -> Right   c+      | otherwise     -> case minimalZeroSumCircle y of+          Just circle     -> reduce $ ZMod.sub input (ZMod.scale c circle)+          Nothing         -> Left input++minimalZeroSumCircle :: Rational -> Maybe (ZMod E2PiI)+minimalZeroSumCircle y +  | y < half     = Nothing+  | r == 1       = Just $ ZMod.fromList [ (E2PiI (i % q) , 1) | i<-[0..q-1] ]+  | otherwise    = Nothing+  where+    p = numerator   y+    q = denominator y+    r = q - p++--------------------------------------------------------------------------------++type X   = U "x"+type Pre = ZMod E2PiI++instance IsSigned Pre where+  signOf _ = Just Plus++preCyclotomic :: Int -> FreeMod Pre X+preCyclotomic n = ZMod.product terms where+  terms = [ ZMod.sub x (ZMod.konst (e k)) | k <- [1..n] , gcd k n == 1 ] where+  x   = ZMod.generator (U 1)+  qn  = fromIntegral n :: Rational+  e k = ZMod.generator $ E2PiI (mod1 $ fromIntegral k / qn) :: ZMod E2PiI++-- | Naive algorithm (using the direct definition of cyclotomic polynomials, and reducing+-- sums of roots of unity till they become integers)+cyclotomicNaive :: Int -> Univariate Integer "x"+cyclotomicNaive = Uni . ZMod.mapCoeff fun . preCyclotomic where+  fun term = case reduce term of+    Right c -> c+    Left {} -> error "cyclotomicNaive: shouldn't happen" ++--------------------------------------------------------------------------------
+ src/Math/Algebra/Polynomial/Univariate/Hermite.hs view
@@ -0,0 +1,68 @@++-- | Hermit polynomials+--+-- See <https://en.wikipedia.org/wiki/Hermite_polynomials>+-- ++{-# LANGUAGE DataKinds, TypeSynonymInstances, FlexibleContexts, FlexibleInstances, BangPatterns, ScopedTypeVariables #-}+module Math.Algebra.Polynomial.Univariate.Hermite+  ( hermiteH+  , hermiteHe+  , integralHermiteH+  , integralHermiteHe+  ) +  where++--------------------------------------------------------------------------------++import Data.List+import Data.Ratio++import Data.Semigroup+import Data.Monoid++import GHC.TypeLits++import qualified Math.Algebra.Polynomial.FreeModule as ZMod+import Math.Algebra.Polynomial.FreeModule ( FreeMod , FreeModule(..) , ZMod , QMod )++import Math.Algebra.Polynomial.Univariate++import Math.Algebra.Polynomial.Class+import Math.Algebra.Polynomial.Pretty+import Math.Algebra.Polynomial.Misc++--------------------------------------------------------------------------------++-- | Probabilists\' Hermite polynomials+hermiteHe :: (Ring c, KnownSymbol v) => Int -> Univariate c v+hermiteHe = fromZUni . renameUniVar . integralHermiteHe++-- | Physicists\' Hermite polynomials+hermiteH :: (Ring c, KnownSymbol v) => Int -> Univariate c v+hermiteH = fromZUni . renameUniVar . integralHermiteH++--------------------------------------------------------------------------------++x, twox :: Univariate Integer "x"+x    = variableP ()+twox = monomP'   (U 1) 2++-- | Probabilists\' Hermite polynomials+integralHermiteHe :: Int -> Univariate Integer "x"+integralHermiteHe = intCache compute where+  compute recur n = case n of +    0 -> 1+    1 -> x+    n -> x * recur (n-1) - differentiateUni (recur (n-1))++-- | Physicists\' Hermite polynomials+integralHermiteH :: Int -> Univariate Integer "x"+integralHermiteH = intCache compute where+  compute recur n = case n of +    0 -> 1+    1 -> twox+    n -> twox * recur (n-1) - differentiateUni (recur (n-1))++--------------------------------------------------------------------------------+
+ src/Math/Algebra/Polynomial/Univariate/Lagrange.hs view
@@ -0,0 +1,54 @@++-- | Lagrange interpolation+--+-- See <https://en.wikipedia.org/wiki/Lagrange_polynomial>++{-# LANGUAGE BangPatterns, DataKinds, KindSignatures, GeneralizedNewtypeDeriving, TypeFamilies #-}+module Math.Algebra.Polynomial.Univariate.Lagrange where++--------------------------------------------------------------------------------++import GHC.TypeLits++{-+import Data.Array ( Array , (!) , listArray , assocs ) +import Data.List++import Data.Proxy++import Math.Algebra.Polynomial.Misc+import Math.Algebra.Polynomial.Pretty+-}++import Math.Algebra.Polynomial.Class++import qualified Math.Algebra.Polynomial.FreeModule as ZMod+import Math.Algebra.Polynomial.FreeModule ( FreeMod , ZMod , QMod )++import Math.Algebra.Polynomial.Univariate++--------------------------------------------------------------------------------+-- * Lagrange interpolation++-- | The Lagrange interpolation for a list of @(x_i,y_i)@ pairs: the resulting polynomial P+-- is the one with minimal degree for which @P(x_i) = y_i@+lagrangeInterp :: KnownSymbol var => [(Rational,Rational)] -> Univariate Rational var +lagrangeInterp xys = final where+  final = sumP [ scaleP (ys!!j) (lagrangePoly xs j) | j<-[0..m-1] ] where+  m = length xys+  (xs,ys) = unzip xys++-- | The minimal degree polynomial which evaluates to zero at the given inputs, except a single one+-- on which it evaluates to 1+lagrangePoly :: KnownSymbol var => [Rational] -> Int -> Univariate Rational var+lagrangePoly xs j = Uni $ ZMod.scale (1/denom) numer where+  numer  = ZMod.product [ term i    | i<-[0..m-1] , i /= j ]+  denom  = product      [ x j - x i | i<-[0..m-1] , i /= j ]+  m      = length xs+  x i    = xs !! i+  term i = ZMod.fromList +    [ (U 1 ,     1 )+    , (U 0 , - x i )+    ]++--------------------------------------------------------------------------------
+ src/Math/Algebra/Polynomial/Univariate/Legendre.hs view
@@ -0,0 +1,55 @@++-- | Legendre polynomials+--+-- See <https://en.wikipedia.org/wiki/Legendre_polynomials>+-- ++{-# LANGUAGE DataKinds, TypeSynonymInstances, FlexibleContexts, FlexibleInstances, BangPatterns, ScopedTypeVariables #-}+module Math.Algebra.Polynomial.Univariate.Legendre+  ( legendreP+  , rationalLegendreP+  ) +  where++--------------------------------------------------------------------------------++import Data.List+import Data.Ratio++import Data.Semigroup+import Data.Monoid++import GHC.TypeLits++import qualified Math.Algebra.Polynomial.FreeModule as ZMod+import Math.Algebra.Polynomial.FreeModule ( FreeMod , FreeModule(..) , ZMod , QMod )++import Math.Algebra.Polynomial.Univariate++import Math.Algebra.Polynomial.Class+import Math.Algebra.Polynomial.Pretty+import Math.Algebra.Polynomial.Misc++--------------------------------------------------------------------------------++-- | Legendre polynomials+legendreP :: (Field c, KnownSymbol v) => Int -> Univariate c v+legendreP = fromQUni . renameUniVar . rationalLegendreP++--------------------------------------------------------------------------------++x :: Univariate Rational "x"+x = variableP ()++rationalLegendreP :: Int -> Univariate Rational "x"+rationalLegendreP = intCache compute where+  fi = (fromIntegral :: Int -> Rational)+  compute recur n = case n of +    0 -> 1+    1 -> x+    n -> scaleP (1 / fi n) +       $ scaleP (fi $ 2*n-1) (x * recur (n-1)) +       - scaleP (fi $   n-1) (    recur (n-2))++--------------------------------------------------------------------------------+
+ src/Math/Algebra/Polynomial/Univariate/Pochhammer.hs view
@@ -0,0 +1,86 @@++-- | Rising and falling factorials+--+-- See <https://en.wikipedia.org/wiki/Falling_and_rising_factorials>++{-# LANGUAGE BangPatterns, DataKinds, KindSignatures #-}+module Math.Algebra.Polynomial.Univariate.Pochhammer where++--------------------------------------------------------------------------------++import Data.Array ( Array , (!) , listArray , assocs ) +import Data.List++import Data.Typeable+import GHC.TypeLits+import Data.Proxy++import Math.Algebra.Polynomial.Pretty+import Math.Algebra.Polynomial.Misc++import qualified Math.Algebra.Polynomial.FreeModule as ZMod+import Math.Algebra.Polynomial.FreeModule ( FreeMod , ZMod , QMod )++import Math.Algebra.Polynomial.Monomial.Univariate ( U(..) )+import Math.Algebra.Polynomial.Univariate++--------------------------------------------------------------------------------+-- * Rising and Falling factorials++risingFactorial :: Int -> Univariate Integer "x"+risingFactorial n = expandRisingFactorial (RF n)++fallingFactorial :: Int -> Univariate Integer "x"+fallingFactorial n = expandFallingFactorial (FF n)++--------------------------------------------------------------------------------+-- * Polynomials using rising or falling factorials as bases++-- | Univariate polynomials using /rising factorials/ as a basis function+newtype RisingPoly  coeff = RisingPoly  { fromRisingPoly  :: FreeMod coeff RisingF }  deriving (Eq,Show)++-- | Univariate polynomials using /falling factorials/ as a basis function+newtype FallingPoly coeff = FallingPoly { fromFallingPoly :: FreeMod coeff FallingF } deriving (Eq,Show)++expandRisingPoly :: (KnownSymbol var, Typeable c, Eq c, Num c) => FreeMod c RisingF -> Univariate c var+expandRisingPoly = Uni . ZMod.flatMap (unUni . expandRisingFactorial)++expandFallingPoly :: (KnownSymbol var, Typeable c, Eq c, Num c) => FreeMod c FallingF -> Univariate c var+expandFallingPoly = Uni . ZMod.flatMap (unUni . expandFallingFactorial)++--------------------------------------------------------------------------------+-- * Rising and falling factorial types++-- | Rising factorial @x^(k) = x(x+1)(x+2)...(x+k-1)@+newtype RisingF = RF Int deriving (Eq,Ord,Show)++-- | Falling factorial @x_(k) = x(x-1)(x-2)...(x-k+1)@+newtype FallingF = FF Int deriving (Eq,Ord,Show)++instance Pretty RisingF where+  pretty (RF k) = case k of+    0 -> "1"+    1 -> "x"+    _ -> "x^(" ++ show k ++ ")"++instance Pretty FallingF where+  pretty (FF k) = case k of+    0 -> "1"+    1 -> "x"+    _ -> "x_(" ++ show k ++ ")"++expandRisingFactorial :: (KnownSymbol var, Typeable c, Eq c, Num c) => RisingF -> Univariate c var+expandRisingFactorial = Uni . ZMod.fromZMod . unUni . expandRisingFactorialZ++expandFallingFactorial ::(KnownSymbol var, Typeable c, Eq c, Num c) => FallingF -> Univariate c var+expandFallingFactorial = Uni . ZMod.fromZMod . unUni . expandFallingFactorialZ++expandRisingFactorialZ :: RisingF -> Univariate Integer var+expandRisingFactorialZ (RF k) = Uni $ ZMod.fromList+  [ (U p, abs c) | (p,c) <- assocs (signedStirling1stArray k) ]++expandFallingFactorialZ :: FallingF -> Univariate Integer var+expandFallingFactorialZ (FF k) = Uni $ ZMod.fromList+  [ (U p,     c) | (p,c) <- assocs (signedStirling1stArray k) ]++--------------------------------------------------------------------------------