diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -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. + diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -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. diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/polynomial-algebra.cabal b/polynomial-algebra.cabal new file mode 100644 --- /dev/null +++ b/polynomial-algebra.cabal @@ -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 + +-------------------------------------------------------------------------------- + diff --git a/src/Math/Algebra/Polynomial/Class.hs b/src/Math/Algebra/Polynomial/Class.hs new file mode 100644 --- /dev/null +++ b/src/Math/Algebra/Polynomial/Class.hs @@ -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 + +-------------------------------------------------------------------------------- diff --git a/src/Math/Algebra/Polynomial/Exterior/Indexed.hs b/src/Math/Algebra/Polynomial/Exterior/Indexed.hs new file mode 100644 --- /dev/null +++ b/src/Math/Algebra/Polynomial/Exterior/Indexed.hs @@ -0,0 +1,166 @@ + +-- | Exterior algebra where the variable set +-- looks like @{x_1, x_2, ... , x_N}@ +-- +-- See + +{-# 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!" + +-------------------------------------------------------------------------------- diff --git a/src/Math/Algebra/Polynomial/FreeModule.hs b/src/Math/Algebra/Polynomial/FreeModule.hs new file mode 100644 --- /dev/null +++ b/src/Math/Algebra/Polynomial/FreeModule.hs @@ -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 + +-------------------------------------------------------------------------------- diff --git a/src/Math/Algebra/Polynomial/Misc.hs b/src/Math/Algebra/Polynomial/Misc.hs new file mode 100644 --- /dev/null +++ b/src/Math/Algebra/Polynomial/Misc.hs @@ -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] ] + +-------------------------------------------------------------------------------- diff --git a/src/Math/Algebra/Polynomial/Monomial/Compact.hs b/src/Math/Algebra/Polynomial/Monomial/Compact.hs new file mode 100644 --- /dev/null +++ b/src/Math/Algebra/Polynomial/Monomial/Compact.hs @@ -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" + +-------------------------------------------------------------------------------- diff --git a/src/Math/Algebra/Polynomial/Monomial/Exterior/Indexed.hs b/src/Math/Algebra/Polynomial/Monomial/Exterior/Indexed.hs new file mode 100644 --- /dev/null +++ b/src/Math/Algebra/Polynomial/Monomial/Exterior/Indexed.hs @@ -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) + +-------------------------------------------------------------------------------- diff --git a/src/Math/Algebra/Polynomial/Monomial/Generic.hs b/src/Math/Algebra/Polynomial/Monomial/Generic.hs new file mode 100644 --- /dev/null +++ b/src/Math/Algebra/Polynomial/Monomial/Generic.hs @@ -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 + +-------------------------------------------------------------------------------- diff --git a/src/Math/Algebra/Polynomial/Monomial/Indexed.hs b/src/Math/Algebra/Polynomial/Monomial/Indexed.hs new file mode 100644 --- /dev/null +++ b/src/Math/Algebra/Polynomial/Monomial/Indexed.hs @@ -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!" + +-------------------------------------------------------------------------------- diff --git a/src/Math/Algebra/Polynomial/Monomial/Infinite.hs b/src/Math/Algebra/Polynomial/Monomial/Infinite.hs new file mode 100644 --- /dev/null +++ b/src/Math/Algebra/Polynomial/Monomial/Infinite.hs @@ -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 + +-------------------------------------------------------------------------------- + diff --git a/src/Math/Algebra/Polynomial/Monomial/Tensor.hs b/src/Math/Algebra/Polynomial/Monomial/Tensor.hs new file mode 100644 --- /dev/null +++ b/src/Math/Algebra/Polynomial/Monomial/Tensor.hs @@ -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 + +-------------------------------------------------------------------------------- diff --git a/src/Math/Algebra/Polynomial/Monomial/Univariate.hs b/src/Math/Algebra/Polynomial/Monomial/Univariate.hs new file mode 100644 --- /dev/null +++ b/src/Math/Algebra/Polynomial/Monomial/Univariate.hs @@ -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 + +-------------------------------------------------------------------------------- diff --git a/src/Math/Algebra/Polynomial/Multivariate/Compact.hs b/src/Math/Algebra/Polynomial/Multivariate/Compact.hs new file mode 100644 --- /dev/null +++ b/src/Math/Algebra/Polynomial/Multivariate/Compact.hs @@ -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) + +-------------------------------------------------------------------------------- diff --git a/src/Math/Algebra/Polynomial/Multivariate/Generic.hs b/src/Math/Algebra/Polynomial/Multivariate/Generic.hs new file mode 100644 --- /dev/null +++ b/src/Math/Algebra/Polynomial/Multivariate/Generic.hs @@ -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) + +-------------------------------------------------------------------------------- diff --git a/src/Math/Algebra/Polynomial/Multivariate/Indexed.hs b/src/Math/Algebra/Polynomial/Multivariate/Indexed.hs new file mode 100644 --- /dev/null +++ b/src/Math/Algebra/Polynomial/Multivariate/Indexed.hs @@ -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!" + +-------------------------------------------------------------------------------- diff --git a/src/Math/Algebra/Polynomial/Multivariate/Infinite.hs b/src/Math/Algebra/Polynomial/Multivariate/Infinite.hs new file mode 100644 --- /dev/null +++ b/src/Math/Algebra/Polynomial/Multivariate/Infinite.hs @@ -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 + +-------------------------------------------------------------------------------- + diff --git a/src/Math/Algebra/Polynomial/Pretty.hs b/src/Math/Algebra/Polynomial/Pretty.hs new file mode 100644 --- /dev/null +++ b/src/Math/Algebra/Polynomial/Pretty.hs @@ -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 ++ "]" + +-------------------------------------------------------------------------------- + +-} diff --git a/src/Math/Algebra/Polynomial/Univariate.hs b/src/Math/Algebra/Polynomial/Univariate.hs new file mode 100644 --- /dev/null +++ b/src/Math/Algebra/Polynomial/Univariate.hs @@ -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) + +-------------------------------------------------------------------------------- diff --git a/src/Math/Algebra/Polynomial/Univariate/Bernoulli.hs b/src/Math/Algebra/Polynomial/Univariate/Bernoulli.hs new file mode 100644 --- /dev/null +++ b/src/Math/Algebra/Polynomial/Univariate/Bernoulli.hs @@ -0,0 +1,135 @@ + +-- | Bernoulli and Euler polynomials +-- +-- See +-- + +{-# 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 +-- +-- 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 + +-------------------------------------------------------------------------------- + diff --git a/src/Math/Algebra/Polynomial/Univariate/Chebysev.hs b/src/Math/Algebra/Polynomial/Univariate/Chebysev.hs new file mode 100644 --- /dev/null +++ b/src/Math/Algebra/Polynomial/Univariate/Chebysev.hs @@ -0,0 +1,63 @@ + +-- | Chebysev polynomials +-- +-- See + +{-# 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) + +-------------------------------------------------------------------------------- + diff --git a/src/Math/Algebra/Polynomial/Univariate/Cyclotomic.hs b/src/Math/Algebra/Polynomial/Univariate/Cyclotomic.hs new file mode 100644 --- /dev/null +++ b/src/Math/Algebra/Polynomial/Univariate/Cyclotomic.hs @@ -0,0 +1,155 @@ + +-- | Cyclotomic polynomials +-- +-- See +-- + +{-# 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" + +-------------------------------------------------------------------------------- diff --git a/src/Math/Algebra/Polynomial/Univariate/Hermite.hs b/src/Math/Algebra/Polynomial/Univariate/Hermite.hs new file mode 100644 --- /dev/null +++ b/src/Math/Algebra/Polynomial/Univariate/Hermite.hs @@ -0,0 +1,68 @@ + +-- | Hermit polynomials +-- +-- See +-- + +{-# 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)) + +-------------------------------------------------------------------------------- + diff --git a/src/Math/Algebra/Polynomial/Univariate/Lagrange.hs b/src/Math/Algebra/Polynomial/Univariate/Lagrange.hs new file mode 100644 --- /dev/null +++ b/src/Math/Algebra/Polynomial/Univariate/Lagrange.hs @@ -0,0 +1,54 @@ + +-- | Lagrange interpolation +-- +-- See + +{-# 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 ) + ] + +-------------------------------------------------------------------------------- diff --git a/src/Math/Algebra/Polynomial/Univariate/Legendre.hs b/src/Math/Algebra/Polynomial/Univariate/Legendre.hs new file mode 100644 --- /dev/null +++ b/src/Math/Algebra/Polynomial/Univariate/Legendre.hs @@ -0,0 +1,55 @@ + +-- | Legendre polynomials +-- +-- See +-- + +{-# 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)) + +-------------------------------------------------------------------------------- + diff --git a/src/Math/Algebra/Polynomial/Univariate/Pochhammer.hs b/src/Math/Algebra/Polynomial/Univariate/Pochhammer.hs new file mode 100644 --- /dev/null +++ b/src/Math/Algebra/Polynomial/Univariate/Pochhammer.hs @@ -0,0 +1,86 @@ + +-- | Rising and falling factorials +-- +-- See + +{-# 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) ] + +--------------------------------------------------------------------------------