{-| Module : Math.Algebra.Hspray Description : Multivariate polynomials on a ring. Copyright : (c) Stéphane Laurent, 2022-2024 License : GPL-3 Maintainer : laurent_step@outlook.fr Deals with multivariate polynomials over a commutative ring. See README for examples. -} {-# LANGUAGE BangPatterns #-} {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE MultiParamTypeClasses #-} {-# LANGUAGE InstanceSigs #-} {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE GeneralizedNewtypeDeriving #-} {-# LANGUAGE TypeFamilies #-} module Math.Algebra.Hspray ( -- * Classes FunctionLike (..) , isConstant , isUnivariate , isBivariate , isTrivariate -- * Main types , Exponents , Powers (..) , Spray , QSpray , QSpray' , Term -- * Basic sprays , lone , qlone , lone' , qlone' , monomial , qmonomial , unitSpray , zeroSpray , constantSpray -- * Showing a spray , prettySpray , prettySpray' , prettySpray'' , prettySprayXYZ , prettySprayX1X2X3 , showSpray , showSprayXYZ , showSprayXYZ' , showSprayX1X2X3 , showSprayX1X2X3' , showNumSpray , showQSpray , showQSpray' , prettyNumSprayX1X2X3 , prettyQSprayX1X2X3 , prettyQSprayX1X2X3' , prettyNumSprayXYZ , prettyQSprayXYZ , prettyQSprayXYZ' , prettyNumSpray , prettyNumSpray' , prettyQSpray , prettyQSpray'' , prettyQSpray' , prettyQSpray''' -- * Univariate polynomials and fractions of univariate polynomials -- -- | The univariate polynomials and the fractions of univariate polynomials -- are used to defined the `OneParameterSpray` sprays, which represent -- multivariate polynomials depending on a unique parameter. These sprays -- lost their interest in version 0.4.0.0 (see CHANGELOG or README). , A (..) , Rational' , Polynomial , QPolynomial , RatioOfPolynomials , RatioOfQPolynomials , prettyRatioOfPolynomials , prettyRatioOfQPolynomials , constPoly , polyFromCoeffs , soleParameter , constQPoly , qpolyFromCoeffs , qsoleParameter , evalRatioOfPolynomials -- * One-parameter sprays -- -- | The `OneParameterSpray` sprays represent multivariate polynomials with -- coefficients depending on a unique parameter. These sprays lost their -- interest in version 0.4.0.0 (see CHANGELOG or README). One can use the -- more general `ParametricSpray` sprays instead. , OneParameterSpray , OneParameterQSpray , prettyOneParameterSprayX1X2X3 , prettyOneParameterSprayXYZ , prettyOneParameterSpray , prettyOneParameterSpray' , prettyOneParameterQSprayX1X2X3 , prettyOneParameterQSprayXYZ , prettyOneParameterQSpray , prettyOneParameterQSpray' , evalOneParameterSpray , substituteTheParameter , evalOneParameterSpray' , evalOneParameterSpray'' -- * Ratios of sprays -- -- | An object of type `RatioOfSprays` represents a fraction of two -- multivariate polynomials. , RatioOfSprays (..) , RatioOfQSprays , (%:%) , (%//%) , (^/^) , (%/%) , isConstantRatioOfSprays , isPolynomialRatioOfSprays , zeroRatioOfSprays , zeroROS , unitRatioOfSprays , unitROS , constantRatioOfSprays , asRatioOfSprays , evalRatioOfSprays , substituteRatioOfSprays , fromRatioOfPolynomials , fromRatioOfQPolynomials , showRatioOfSprays , showRatioOfNumSprays , showRatioOfQSprays , showRatioOfSpraysXYZ , showRatioOfSpraysXYZ' , showRatioOfSpraysX1X2X3 , showRatioOfSpraysX1X2X3' , prettyRatioOfQSpraysXYZ , prettyRatioOfQSpraysX1X2X3 , prettyRatioOfQSprays , prettyRatioOfQSprays' , prettyRatioOfNumSpraysXYZ , prettyRatioOfNumSpraysX1X2X3 , prettyRatioOfNumSprays , prettyRatioOfNumSprays' -- * Parametric sprays -- -- | There are three types of parametric sprays: @OneParameterSpray@, -- @SimpleParametricSpray@ and @ParametricSpray@. These are sprays of -- type @Spray b@ where @b@ has the class @FunctionLike@. When we say -- \"parametric spray\" in the documentation, we mean either -- such a spray or more precisely a @ParametricSpray@ spray. , SimpleParametricSpray , SimpleParametricQSpray , ParametricSpray , ParametricQSpray , canCoerceToSimpleParametricSpray , asSimpleParametricSprayUnsafe , asSimpleParametricSpray , fromOneParameterSpray , fromOneParameterQSpray , parametricSprayToOneParameterSpray , parametricQSprayToOneParameterQSpray , gegenbauerPolynomial , jacobiPolynomial , numberOfParameters , changeParameters , substituteParameters , evalParametricSpray , evalParametricSpray' , prettyParametricQSprayABCXYZ , prettyParametricQSpray , prettyParametricNumSprayABCXYZ , prettyParametricNumSpray , prettySimpleParametricQSprayABCXYZ , prettySimpleParametricQSpray , prettySimpleParametricNumSprayABCXYZ , prettySimpleParametricNumSpray -- * Queries on a spray , getCoefficient , getConstantTerm , isZeroSpray , isConstantSpray , isHomogeneousSpray , allExponents , allCoefficients -- * Evaluation of a spray , evalSpray , substituteSpray , composeSpray , evalSpraySpray -- * Division of sprays , sprayDivision , sprayDivisionRemainder , pseudoDivision -- * Gröbner basis , groebnerBasis , reduceGroebnerBasis -- * Symmetric polynomials , esPolynomial , psPolynomial , isSymmetricSpray -- * Resultant, subresultants, and Sturm-Habicht sequence , resultant , resultant' , resultant1 , subresultants , subresultants1 , polynomialSubresultants , sturmHabichtSequence , principalSturmHabichtSequence -- * Number of real roots of a univariate spray. These functions can be very -- slow if the degree of the spray is not small. , numberOfRealRoots , numberOfRealRoots' , numberOfRealRootsInOpenInterval , numberOfRealRootsInOpenInterval' , numberOfRealRootsInClosedInterval , numberOfRealRootsInClosedInterval' -- * Greatest common divisor , gcdSpray -- * Matrices , detLaplace , detLaplace' , characteristicPolynomial -- * Miscellaneous , sumOfSprays , productOfSprays , (.^) , (/>) , (/^) , fromList , toList , fromRationalSpray , isPolynomialOf , bombieriSpray , collinearSprays ) where import qualified Algebra.Absolute as AlgAbs import qualified Algebra.Additive as AlgAdd import qualified Algebra.Differential as AlgDiff import qualified Algebra.Field as AlgField import qualified Algebra.Module as AlgMod import qualified Algebra.RightModule as AlgRightMod import qualified Algebra.Ring as AlgRing import qualified Algebra.ZeroTestable as AlgZT import qualified Data.Foldable as DF import Data.Function ( on ) import Data.HashMap.Strict ( HashMap ) import qualified Data.HashMap.Strict as HM import Data.Hashable ( Hashable ( hashWithSalt ) ) import qualified Data.IntMap.Strict as IM import Data.List ( sortBy , maximumBy , (\\) , findIndices , elemIndices , nub , foldl1' , uncons ) import Data.List.Extra ( allSame , unsnoc ) import Data.Matrix ( Matrix , fromLists , minorMatrix , nrows , ncols , submatrix ) import qualified Data.Matrix as DM import Data.Maybe ( isJust , isNothing , fromJust , fromMaybe ) import Data.Ord ( comparing ) import qualified Data.Ratio as DR import qualified GHC.Real as DR import qualified Data.Sequence as S import Data.Sequence ( (><) , Seq , dropWhileR , (|>) , index , adjust , fromFunction ) import Data.Text ( Text , append , cons , intercalate , pack , snoc , unpack ) import Data.Tuple.Extra ( both, first ) import qualified MathObj.Matrix as MathMatrix import qualified MathObj.Polynomial as MathPol import Number.Ratio ( T ( (:%) ), (%) ) import qualified Number.Ratio as NumberRatio -- import qualified Algebra.PrincipalIdealDomain as AlgPID -- import qualified Algebra.Units as AlgUnits -- import qualified Algebra.IntegralDomain as AlgID -- Classes -------------------------------------------------------------------- -- | A spray represents a multivariate polynomial so it like a function. We -- introduce a class because it will be assigned to the ratios of sprays too. class FunctionLike b where -- | Number of variables in a function-like object numberOfVariables :: b -> Int -- | Permutes the variables of a function-like object -- -- >>> f :: Spray Rational -> Spray Rational -> Spray Rational -> Spray Rational -- >>> f p1 p2 p3 = p1^**^4 ^+^ (2*^p2^**^3) ^+^ (3*^p3^**^2) ^-^ (4*^unitSpray) -- >>> x1 = lone 1 :: Spray Rational -- >>> x2 = lone 2 :: Spray Rational -- >>> x3 = lone 3 :: Spray Rational -- >>> spray = f x1 x2 x3 -- -- prop> permuteVariables [3, 1, 2] spray == f x3 x1 x2 permuteVariables :: [Int] -- ^ permutation -> b -- ^ function-like object whose variables will be permuted -> b -- ^ the function-like object with permuted variables -- | Swaps two variables of a function-like object -- -- prop> swapVariables (1, 3) x == permuteVariables [3, 2, 1] x swapVariables :: (Int, Int) -- ^ the indices of the variables to be swapped (starting at 1) -> b -- ^ function-like object whose variables will be swapped -> b -- ^ the function-like object with swapped variables -- Whether a variable is involved in a function-like object -- -- prop> involvesVariable (qlone 1 ^+^ qlone 3) 2 == False involvesVariable :: b -- ^ function-like object -> Int -- ^ index of the variable -> Bool -- | Drops a given number of leading variables in a function-like object; -- __very unsafe__, @dropVariables n x@ should /not/ be used if -- @involvesVariable x i@ is @True@ for some @i@ in @1, ... n@ -- -- prop> dropVariables 1 (qlone 2 ^+^ qlone 3) == qlone 1 ^+^ qlone 2 dropVariables :: Int -- ^ number of leading variables to drop -> b -- ^ a function-like object -> b -- | Derivative of a function-like object -- -- >>> x = lone 1 :: Spray Int -- >>> y = lone 2 :: Spray Int -- >>> spray = 2*^x ^-^ 3*^y^**^8 -- >>> spray' = derivative 1 spray -- >>> putStrLn $ prettyNumSpray spray' -- 2 derivative :: Int -- ^ index of the variable of differentiation (starting at 1) -> b -- ^ the object to be derivated -> b -- ^ the derivated object -- | The type of the coefficients (this is @a@ for both @Spray a@ and @RatioOfSprays a@) type family BaseRing b -- | The type of the variables (this is @Spray a@ for both @Spray a@ and @RatioOfSprays a@) type family VariablesType b infixl 6 ^+^ -- | Addition of two function-like objects (^+^) :: (AlgAdd.C b) => b -> b -> b (^+^) = (AlgAdd.+) infixl 6 ^-^ -- | Substraction of two function-like objects (^-^) :: (AlgAdd.C b) => b -> b -> b (^-^) = (AlgAdd.-) infixl 7 ^*^ -- | Multiplication of two function-like objects (^*^) :: (AlgRing.C b) => b -> b -> b (^*^) = (AlgRing.*) infixr 8 ^**^ -- | Power of a function-like object (^**^) :: (AlgRing.C b) => b -> Int -> b (^**^) p k = p AlgRing.^ fromIntegral k infixr 7 *^ -- | Multiply a function-like object by a scalar (*^) :: BaseRing b -> b -> b infixl 6 +> -- | Add a function-like object to a constant -- -- prop> x +> spray == constantSpray x ^+^ spray (+>) :: BaseRing b -> b -> b infixr 6 <+ -- | Add a constant to a function-like object -- -- prop> object <+ x == x +> object (<+) :: b -> BaseRing b -> b (<+) = flip (+>) -- | Evaluation (replacing the variables with some values) of a -- function-like object -- -- >>> x = lone 1 :: Spray Int -- >>> y = lone 2 :: Spray Int -- >>> spray = 2*^x^**^2 ^-^ 3*^y -- >>> evaluate spray [2, 1] -- 5 evaluate :: b -- ^ function-like object to be evaluated, e.g. a spray -> [BaseRing b] -- ^ list of values to be substituted to its variables -> BaseRing b -- | Flipped version of @evaluate@ -- -- prop> evaluateAt [2, 1] spray == evaluate spray [2, 1] evaluateAt :: [BaseRing b] -> b -> BaseRing b evaluateAt = flip evaluate -- | Partial evaluation of a function-like object (replace some variables -- with some values) -- -- >>> x1 = lone 1 :: Spray Int -- >>> x2 = lone 2 :: Spray Int -- >>> x3 = lone 3 :: Spray Int -- >>> spray = x1^**^2 ^-^ x2 ^+^ x3 ^-^ unitSpray -- >>> spray' = substitute [Just 2, Nothing, Just 3] spray -- >>> putStrLn $ prettyNumSprayX1X2X3 "x" spray' -- -x2 + 6 substitute :: [Maybe (BaseRing b)] -- ^ @Just x@ to replace the variable with @x@, @Nothing@ for no replacement -> b -- ^ function-like object to be partially evaluated -> b -- | Polynomial change of variables of a function-like object -- -- >>> x = lone 1 :: Spray Int -- >>> y = lone 2 :: Spray Int -- >>> spray = x ^*^ y -- >>> spray' = changeVariables spray [x ^+^ y, x ^-^ y] -- >>> putStrLn $ prettyNumSpray' spray' -- X^2 - Y^2 changeVariables :: b -- ^ function-like object such as a spray -> [VariablesType b] -- ^ list of new variables -> b -- | Whether a function-like object has a constant value isConstant :: FunctionLike b => b -> Bool isConstant f = numberOfVariables f == 0 -- | Whether a function-like object represents an univariate function; it is considered -- that it is univariate if it is constant isUnivariate :: FunctionLike b => b -> Bool isUnivariate f = numberOfVariables f <= 1 -- | Whether a function-like object represents a bivariate function; it is considered -- that it is bivariate if it is univariate isBivariate :: FunctionLike b => b -> Bool isBivariate f = numberOfVariables f <= 2 -- | Whether a function-like object represents a trivariate function; it is considered -- that it is trivariate if it is bivariate isTrivariate :: FunctionLike b => b -> Bool isTrivariate f = numberOfVariables f <= 3 -- Additional operations to 'numeric-prelude' infixr 7 /> -- | Divides by a scalar in a module over a field (/>) :: (AlgField.C k, AlgMod.C k a) => a -> k -> a x /> lambda = AlgField.recip lambda AlgMod.*> x infixr 7 .^ -- | Scale by an integer (I do not find this operation in __numeric-prelude__) -- -- prop> 3 .^ x == x Algebra.Additive.+ x Algebra.Additive.+ x (.^) :: (AlgAdd.C a, Eq a) => Int -> a -> a k .^ x = if k >= 0 then powerOperation (AlgAdd.+) AlgAdd.zero x k else (.^) (-k) (AlgAdd.negate x) where powerOperation op = let go acc _ 0 = acc go acc a n = go (if even n then acc else op acc a) (op a a) (div n 2) in go -- Univariate polynomials and ratios of univariate polynomials ---------------- -- | The new type @A a@ is used to attribute some instances to the -- type @Polynomial a@; it is needed to avoid orphan instances. newtype A a = A a deriving (Eq, Show, AlgAdd.C, AlgRing.C, AlgField.C) -- | The type @Rational'@ is used to introduce the univariate polynomials -- with rational coefficients (`QPolynomial`). It is similar to the well-known -- type @Rational@ (actually these two types are the same but @Rational'@ has -- more instances and we need them for the univariate polynomials). type Rational' = NumberRatio.Rational -- | The type @Polynomial a@ is used to represent univariate polynomials. type Polynomial a = MathPol.T (A a) type QPolynomial = Polynomial Rational' -- The type @RatioOfPolynomials a@ is used to represent fractions of -- univariate polynomials. type RatioOfPolynomials a = NumberRatio.T (Polynomial a) type RatioOfQPolynomials = RatioOfPolynomials Rational' instance (Eq a, AlgField.C a) => FunctionLike (Polynomial a) where -- numberOfVariables :: Polynomial a -> Int numberOfVariables p = case MathPol.degree p of Nothing -> 0 Just d -> min 1 d -- type BaseRing (Polynomial a) = a -- type VariablesType (Polynomial a) = Polynomial a -- (*^) :: a -> Polynomial a -> Polynomial a (*^) lambda pol = constPoly lambda AlgRing.* pol -- (+>) :: a -> Polynomial a -> Polynomial a (+>) x pol = constPoly x AlgAdd.+ pol -- evaluate :: Polynomial a -> [a] -> a evaluate p xs = get (MathPol.evaluate p (A (xs !! 0))) where get (A x) = x -- substitute :: [Maybe a] -> Polynomial a -> Polynomial a substitute x p = if isNothing (x !! 0) then p else constPoly (evaluate p [fromJust $ x !! 0]) -- permuteVariables :: [Int] -> Polynomial a -> Polynomial a permuteVariables = error "permuteVariables: there is only one variable." -- swapVariables :: (Int, Int) -> Polynomial a -> Polynomial a swapVariables = error "swapVariables: there is only one variable." -- involvesVariable :: Polynomial a -> Int -> Bool involvesVariable = error "involvesVariable: not available for `Polynomial`" -- dropVariables :: Int -> Polynomial a -> Polynomial a dropVariables = error "dropVariables: not available for `Polynomial`" -- derivative :: Int -> Polynomial a -> Polynomial a derivative i p = if i == 1 then AlgDiff.differentiate p else constPoly AlgAdd.zero -- changeVariables :: Polynomial a -> [Polynomial a] -> Polynomial a changeVariables p ps = MathPol.compose p (ps !! 0) instance (Eq a, AlgField.C a) => FunctionLike (RatioOfPolynomials a) where numberOfVariables :: RatioOfPolynomials a -> Int numberOfVariables (p :% q) = max (numberOfVariables p) (numberOfVariables q) -- type BaseRing (RatioOfPolynomials a) = a -- type VariablesType (RatioOfPolynomials a) = Polynomial a -- (*^) :: a -> RatioOfPolynomials a -> RatioOfPolynomials a (*^) lambda rOP = A lambda AlgMod.*> rOP -- (+>) :: a -> RatioOfPolynomials a -> RatioOfPolynomials a (+>) x rOP = NumberRatio.fromValue (constPoly x) AlgAdd.+ rOP -- evaluate :: RatioOfPolynomials a -> [a] -> a evaluate r xs = evaluate (NumberRatio.numerator r) xs AlgField./ evaluate (NumberRatio.denominator r) xs -- substitute :: [Maybe a] -> RatioOfPolynomials a -> RatioOfPolynomials a substitute x r = if isNothing (x !! 0) then r else substitute x (NumberRatio.numerator r) % substitute x (NumberRatio.denominator r) -- permuteVariables :: [Int] -> RatioOfPolynomials a -> RatioOfPolynomials a permuteVariables = error "permuteVariables: there is only one variable." -- swapVariables :: (Int, Int) -> RatioOfPolynomials a -> RatioOfPolynomials a swapVariables = error "swapVariables: there is only one variable." -- involvesVariable :: RatioOfPolynomials a -> Int -> Bool involvesVariable = error "involvesVariable: not available for `RatioOfPolynomials`" -- dropVariables :: Int -> RatioOfPolynomials a -> RatioOfPolynomials a dropVariables = error "dropVariables: not available for `RatioOfPolynomials`" -- derivative :: Int -> RatioOfPolynomials a -> RatioOfPolynomials a derivative i r = if i == 1 then (p' AlgRing.* q AlgAdd.- p AlgRing.* q') % q AlgRing.^ 2 else constPoly AlgAdd.zero :% constPoly AlgRing.one where p = NumberRatio.numerator r q = NumberRatio.denominator r p' = AlgDiff.differentiate p q' = AlgDiff.differentiate q -- changeVariables :: RatioOfPolynomials a -> [Polynomial a] -> RatioOfPolynomials a changeVariables r ps = changeVariables (NumberRatio.numerator r) ps % changeVariables (NumberRatio.denominator r) ps instance (Eq a, AlgField.C a) => AlgZT.C (A a) where isZero :: A a -> Bool isZero (A r) = r == AlgAdd.zero instance (Eq a, AlgField.C a) => AlgMod.C (A a) (RatioOfPolynomials a) where (*>) :: A a -> RatioOfPolynomials a -> RatioOfPolynomials a lambda *> rop = NumberRatio.scale (MathPol.const lambda) rop instance (Eq a, AlgField.C a) => AlgRightMod.C (A a) (RatioOfPolynomials a) where (<*) :: RatioOfPolynomials a -> A a -> RatioOfPolynomials a rop <* lambda = lambda AlgMod.*> rop instance (Eq a, AlgField.C a) => AlgMod.C a (RatioOfPolynomials a) where (*>) :: a -> RatioOfPolynomials a -> RatioOfPolynomials a lambda *> rop = A lambda AlgMod.*> rop instance (Eq a, AlgField.C a) => AlgRightMod.C a (RatioOfPolynomials a) where (<*) :: RatioOfPolynomials a -> a -> RatioOfPolynomials a rop <* lambda = lambda AlgMod.*> rop instance (Eq a, AlgField.C a) => AlgMod.C (Polynomial a) (RatioOfPolynomials a) where (*>) :: Polynomial a -> RatioOfPolynomials a -> RatioOfPolynomials a p *> r = NumberRatio.scale p r instance (Eq a, AlgField.C a) => AlgRightMod.C (Polynomial a) (RatioOfPolynomials a) where (<*) :: RatioOfPolynomials a -> Polynomial a -> RatioOfPolynomials a r <* p = p AlgMod.*> r -- | Constant univariate polynomial constPoly :: a -> Polynomial a constPoly x = MathPol.const (A x) -- | Univariate polynomial from its coefficients (ordered by increasing degrees) polyFromCoeffs :: [a] -> Polynomial a polyFromCoeffs as = MathPol.fromCoeffs (map A as) -- | The variable of a univariate polynomial; it is called \"soleParameter\" because -- this it represents the parameter of a `OneParameterSpray` spray soleParameter :: AlgRing.C a => Polynomial a soleParameter = polyFromCoeffs [AlgAdd.zero, AlgRing.one] -- | Constant rational univariate polynomial -- -- >>> import Number.Ratio ( (%) ) -- >>> constQPoly (2 % 3) -- -- prop> constQPoly (2 % 3) == qpolyFromCoeffs [2 % 3] constQPoly :: Rational' -> QPolynomial constQPoly = constPoly -- | Rational univariate polynomial from coefficients -- -- >>> import Number.Ratio ( (%) ) -- >>> qpolyFromCoeffs [2 % 3, 5, 7 % 4] qpolyFromCoeffs :: [Rational'] -> QPolynomial qpolyFromCoeffs = polyFromCoeffs -- | The variable of a univariate rational polynomial; it is called \"qsoleParameter\" -- because it represents the parameter of a `OneParameterQSpray` spray -- -- prop> qsoleParameter == qpolyFromCoeffs [0, 1] qsoleParameter :: QPolynomial qsoleParameter = qpolyFromCoeffs [0, 1] {- -- show a ratio, helper function showQ :: (Eq a, Num a, Show a) => NumberRatio.T a -> String showQ q = if d == 1 then show n else show n ++ "/" ++ show d where n = NumberRatio.numerator q d = NumberRatio.denominator q -} -- | identify a `Polynomial a` to a `Spray a`, in order to apply the show spray -- functions to the univariate polynomials polynomialToSpray :: forall a. (Eq a, AlgRing.C a) => Polynomial a -> Spray a polynomialToSpray pol = sumOfSprays terms where coeffs = MathPol.coeffs pol indices = findIndices (/= A AlgAdd.zero) coeffs get :: A a -> a get (A x) = x terms = [get (coeffs !! i) *^ lone' 1 i | i <- indices] qPolynomialToQSpray :: QPolynomial -> QSpray qPolynomialToQSpray pol = sumOfSprays terms where coeffs = MathPol.coeffs pol indices = findIndices (/= A 0) coeffs get :: A Rational' -> Rational get (A x) = NumberRatio.numerator x DR.:% NumberRatio.denominator x terms = [get (coeffs !! i) *^ qlone' 1 i | i <- indices] -- helper function; it encloses a string between two given delimiters bracify :: (String, String) -> String -> String bracify (lbrace, rbrace) x = lbrace ++ x ++ rbrace -- | helper function for prettyRatioOfPolynomials (and prettyOneParameterSpray) showRatioOfPolynomials :: forall a. (Eq a, AlgField.C a) => (Spray a -> String) -> RatioOfPolynomials a -> String showRatioOfPolynomials sprayShower polysRatio = numeratorString ++ denominatorString where numerator = NumberRatio.numerator polysRatio denominator = NumberRatio.denominator polysRatio brackets = denominator /= MathPol.const (A AlgRing.one) enclose = bracify ("[ ", " ]") numeratorString = if brackets then enclose (sprayShower (polynomialToSpray numerator)) else sprayShower (polynomialToSpray numerator) denominatorString = if not brackets then "" else " %//% " ++ enclose (sprayShower (polynomialToSpray denominator)) -- | Pretty form of a ratio of univariate polynomials with rational coefficients prettyRatioOfQPolynomials :: String -- ^ a string to denote the variable, e.g. @\"a\"@ -> RatioOfQPolynomials -> String prettyRatioOfQPolynomials var = showRatioOfPolynomials (prettyQSprayXYZ' [var]) -- | helper function for prettyRatioOfPolynomials (and prettyOneParameterSpray) showQpol :: forall a. (Eq a, AlgField.C a) => Polynomial a -> String -> (a -> String) -> Bool -> String showQpol pol variable showCoeff brackets = if brackets then "[ " ++ polyString ++ " ]" else polyString where showCoeff' :: Int -> A a -> String showCoeff' i (A coeff) = case i of 0 -> (bracify ("(", ")") . showCoeff) coeff _ -> if coeff == AlgRing.one then "" else (bracify ("(", ")") . showCoeff) coeff coeffs = MathPol.coeffs pol nonzeros = findIndices (/= A AlgAdd.zero) coeffs terms = map (pack . showTerm) nonzeros where showTerm i = case i of 0 -> showCoeff' 0 (coeffs !! 0) 1 -> showCoeff' 1 (coeffs !! 1) ++ variable _ -> showCoeff' i (coeffs !! i) ++ variable ++ "^" ++ show i polyString = unpack (intercalate (pack " + ") terms) -- | helper function for prettyRatioOfPolynomials (and prettyOneParameterSpray) showQpolysRatio :: forall a. (Eq a, AlgField.C a) => String -> (a -> String) -> RatioOfPolynomials a -> String showQpolysRatio var showCoeff polysRatio = numeratorString ++ denominatorString where denominator = NumberRatio.denominator polysRatio brackets = denominator /= constPoly AlgRing.one numeratorString = showQpol (NumberRatio.numerator polysRatio) var showCoeff brackets denominatorString = if not brackets then "" else " %//% " ++ showQpol denominator var showCoeff True -- | Pretty form of a ratio of univariate polynomials prettyRatioOfPolynomials :: (Eq a, AlgField.C a, Show a) => String -- ^ string (usually a single letter) to denote the variable, e.g. @\"a\"@ -> RatioOfPolynomials a -> String prettyRatioOfPolynomials var = showQpolysRatio var show {- -- | Pretty form of a ratio of univariate qpolynomials prettyRatioOfQPolynomials' :: String -- ^ a string to denote the variable, e.g. @\"a\"@ -> RatioOfQPolynomials -> String prettyRatioOfQPolynomials' var = showQpolysRatio var showQ -} -- | Evaluates a ratio of univariate polynomials evalRatioOfPolynomials :: AlgField.C a => a -- ^ the value at which the evaluation is desired -> RatioOfPolynomials a -> a evalRatioOfPolynomials value polysRatio = resultNumerator AlgField./ resultDenominator where A resultNumerator = MathPol.evaluate (NumberRatio.numerator polysRatio) (A value) A resultDenominator = MathPol.evaluate (NumberRatio.denominator polysRatio) (A value) -- One-parameter sprays ------------------------------------------------------- type OneParameterSpray a = Spray (RatioOfPolynomials a) type OneParameterQSpray = OneParameterSpray Rational' {- -- | simplifies a ratio of polynomials (simply by multiplying it by one) simplifyRatioOfPolynomials :: (Eq a, AlgField.C a) => RatioOfPolynomials a -> RatioOfPolynomials a simplifyRatioOfPolynomials = (AlgRing.*) AlgRing.one -- | Simplifies the coefficients (the fractions of univariate polynomials) of a -- one-parameter spray simplifyOneParameterSpray :: (Eq a, AlgField.C a) => OneParameterSpray a -> OneParameterSpray a simplifyOneParameterSpray = HM.map simplifyRatioOfPolynomials -} instance (Eq a, AlgField.C a) => AlgMod.C (Polynomial a) (OneParameterSpray a) where (*>) :: Polynomial a -> OneParameterSpray a -> OneParameterSpray a p *> spray = (p NumberRatio.:% AlgRing.one) *^ spray instance (Eq a, AlgField.C a) => AlgRightMod.C (Polynomial a) (OneParameterSpray a) where (<*) :: OneParameterSpray a -> Polynomial a -> OneParameterSpray a spray <* p = p AlgMod.*> spray instance (Eq a, AlgField.C a) => AlgMod.C a (OneParameterSpray a) where (*>) :: a -> OneParameterSpray a -> OneParameterSpray a lambda *> spray = MathPol.const (A lambda) AlgMod.*> spray instance (Eq a, AlgField.C a) => AlgRightMod.C a (OneParameterSpray a) where (<*) :: OneParameterSpray a -> a -> OneParameterSpray a spray <* lambda = lambda AlgMod.*> spray -- | Pretty form of a one-parameter spray, using a string (typically a letter) -- followed by an index to denote the variables prettyOneParameterSprayX1X2X3 :: (Eq a, Show a, AlgField.C a) => String -- ^ string to denote the parameter of the spray, e.g. @\"a\"@ -> String -- ^ typically a letter, to denote the non-indexed variables -> OneParameterSpray a -- ^ a one-parameter spray; note that this function does not simplify it -> String prettyOneParameterSprayX1X2X3 a = showSprayX1X2X3 (prettyRatioOfPolynomials a) ("{ ", " }") -- | Pretty form of a one-parameter spray, using some given strings (typically some -- letters) to denote the variables if possible, i.e. if enough letters are -- provided; otherwise this function behaves exactly like -- @prettyOneParameterSprayX1X2X3 a@ where @a@ is the first provided letter prettyOneParameterSprayXYZ :: (Eq a, Show a, AlgField.C a) => String -- ^ string to denote the parameter of the spray, e.g. @\"a\"@ -> [String] -- ^ typically some letters, to denote the main variables -> OneParameterSpray a -- ^ a one-parameter spray; note that this function does not simplify it -> String prettyOneParameterSprayXYZ a = showSprayXYZ (prettyRatioOfPolynomials a) ("{ ", " }") -- | Pretty form of a one-parameter spray; see the definition below and see -- `prettyOneParameterSprayXYZ` -- -- prop> prettyOneParameterSpray a spray == prettyOneParameterSprayXYZ a ["x","y","z"] spray prettyOneParameterSpray :: (Eq a, Show a, AlgField.C a) => String -- ^ string to denote the parameter of the spray, e.g. @\"a\"@ -> OneParameterSpray a -- ^ a one-parameter spray; note that this function does not simplify it -> String prettyOneParameterSpray a = prettyOneParameterSprayXYZ a ["x", "y", "z"] -- | Pretty form of a one-parameter spray; see the definition below and see -- `prettyOneParameterSprayXYZ` -- -- prop> prettyOneParameterSpray' a spray == prettyOneParameterSprayXYZ a ["X","Y","Z"] spray prettyOneParameterSpray' :: (Eq a, Show a, AlgField.C a) => String -- ^ string to denote the parameter of the spray, e.g. @\"a\"@ -> OneParameterSpray a -- ^ a one-parameter spray; note that this function does not simplify it -> String prettyOneParameterSpray' a = prettyOneParameterSprayXYZ a ["X", "Y", "Z"] -- | Pretty form of a one-parameter rational spray, using a string (typically a letter) -- followed by an index to denote the variables prettyOneParameterQSprayX1X2X3 :: String -- ^ usually a letter, to denote the parameter of the spray, e.g. @\"a\"@ -> String -- ^ usually a letter, to denote the non-indexed variables of the spray -> OneParameterQSpray -- ^ a one-parameter rational spray; note that this function does not simplify it -> String prettyOneParameterQSprayX1X2X3 a x = showSpray (prettyRatioOfQPolynomials a) ("{ ", " }") (showMonomialsX1X2X3 x) -- | Pretty form of a one-parameter rational spray, using some given strings (typically some -- letters) to denote the variables if possible, i.e. if enough letters are -- provided; otherwise this function behaves exactly like -- @prettyOneParameterQSprayX1X2X3 a@ where @a@ is the first provided letter prettyOneParameterQSprayXYZ :: String -- ^ usually a letter, to denote the parameter of the spray, e.g. @\"a\"@ -> [String] -- ^ usually some letters, to denote the variables of the spray -> OneParameterQSpray -- ^ a one-parameter rational spray; note that this function does not simplify it -> String prettyOneParameterQSprayXYZ a letters = showSpray (prettyRatioOfQPolynomials a) ("{ ", " }") (showMonomialsXYZ letters) -- | Pretty form of a one-parameter rational spray, using @\"x\"@, @\"y\"@ and @\"z\"@ for the variables -- if possible; i.e. if the spray does not have more than three variables, otherwise -- @\"x1\"@, @\"x2\"@, ... are used to denote the variables -- -- prop> prettyOneParameterQSpray a == prettyOneParameterQSprayXYZ a ["x","y","z"] prettyOneParameterQSpray :: String -- ^ usually a letter, to denote the parameter of the spray, e.g. @\"a\"@ -> OneParameterQSpray -- ^ the one-parameter rational spray to be printed; note that this function does not simplify it -> String prettyOneParameterQSpray a = prettyOneParameterQSprayXYZ a ["x", "y", "z"] -- | Pretty form of a one-parameter rational spray, using @\"X\"@, @\"Y\"@ and @\"Z\"@ for the variables -- if possible; i.e. if the spray does not have more than three variables, otherwise -- @\"X1\"@, @\"X2\"@, ... are used -- -- prop> prettyOneParameterQSpray' a == prettyOneParameterQSprayXYZ a ["X","Y","Z"] prettyOneParameterQSpray' :: String -- ^ usually a letter, to denote the parameter of the spray, e.g. @\"a\"@ -> OneParameterQSpray -- ^ the one-parameter rational spray to be printed; note that this function does not simplify it -> String prettyOneParameterQSpray' a = prettyOneParameterQSprayXYZ a ["X", "Y", "Z"] -- | Substitutes a value to the parameter of a one-parameter spray -- (the variable occurring in its coefficients) -- -- prop> evalOneParameterSpray spray x == substituteParameters spray [x] evalOneParameterSpray :: (Eq a, AlgField.C a) => OneParameterSpray a -> a -> Spray a evalOneParameterSpray spray x = removeZeroTerms $ HM.map (evalRatioOfPolynomials x) spray -- | Substitutes a value to the parameter of a one-parameter spray; -- same as `evalOneParameterSpray` -- -- prop> substituteTheParameter spray x == substituteParameters spray [x] substituteTheParameter :: (Eq a, AlgField.C a) => OneParameterSpray a -> a -> Spray a substituteTheParameter spray x = substituteParameters spray [x] -- | Substitutes a value to the parameter of a one-parameter spray as well -- as some values to the variables of this spray -- -- prop> evalOneParameterSpray' spray a xs == evalParametricSpray' spray [a] xs evalOneParameterSpray' :: (Eq a, AlgField.C a) => OneParameterSpray a -- ^ one-parameter spray to be evaluated -> a -- ^ a value for the parameter -> [a] -- ^ some values for the variables -> a evalOneParameterSpray' spray x xs = if length xs >= numberOfVariables spray then evalSpray (evalOneParameterSpray spray x) xs else error "evalOneParameterSpray': not enough values provided." -- | helper function for evalOneParameterSpray'' evalOneParameterTerm :: (Eq a, AlgField.C a) => [a] -> Term (RatioOfPolynomials a) -> RatioOfPolynomials a evalOneParameterTerm xs (powers, coeff) = AlgRing.product (zipWith (AlgRing.^) xs pows) *^ coeff where pows = DF.toList (fromIntegral <$> exponents powers) -- | Substitutes some values to the variables of a one-parameter spray; same -- as `evalParametricSpray` evalOneParameterSpray'' :: (Eq a, AlgField.C a) => OneParameterSpray a -> [a] -> RatioOfPolynomials a evalOneParameterSpray'' spray xs = if length xs >= numberOfVariables spray then AlgAdd.sum $ map (evalOneParameterTerm xs) (HM.toList spray) else error "evalOneParameterSpray'': not enough values provided." -- Sprays --------------------------------------------------------------------- type Exponents = Seq Int -- | The type @Powers@ is used to represent the exponents of the monomial -- occurring in a term of a spray. The integer in the field @nvariables@ is -- the number of variables involved in this monomial (it is 3, not 2, for a -- monomial such as @x^2.z^3@, because the exponents of this monomial is the -- sequence @(2, 0, 3)@). Actually this integer is always the length of the -- sequence in the field @exponents@. The reason of the presence of the field -- @nvariables@ is that I thought that it was necessary when I started to -- develop the package, but now I think it is useless. The type @Powers@ will -- possibly be abandoned in a future version of the package. However we cannot -- simply use the type `Exponents` to represent the exponents, because two -- sequences of exponents that differ only by some trailing zeros must be -- considered as identical, and they are considered as such with the type -- @Powers@ thanks to its @Eq@ instance. Instead of @Powers@, a new type -- encapsulating the `Exponents` type with such an @Eq@ instance should be -- enough to represent the exponents. data Powers = Powers { exponents :: Exponents , nvariables :: Int } deriving Show nullPowers :: Powers nullPowers = Powers S.empty 0 powerize :: Exponents -> Powers powerize expnts = Powers expnts (S.length expnts) instance Eq Powers where (==) :: Powers -> Powers -> Bool pows1 == pows2 = expts1' == expts2' where (expts1', expts2') = harmonize (pows1, pows2) instance Hashable Powers where hashWithSalt :: Int -> Powers -> Int hashWithSalt k pows = hashWithSalt k (exponents pows, nvariables pows) -- | append trailing zeros growSequence :: Exponents -> Int -> Int -> Exponents growSequence s m n = s >< S.replicate (n - m) 0 growSequence' :: Int -> Exponents -> Exponents growSequence' n s = growSequence s (S.length s) n -- | append trailing zeros to get the same length harmonize :: (Powers, Powers) -> (Exponents, Exponents) harmonize (pows1, pows2) = (e1', e2') where e1 = exponents pows1 e2 = exponents pows2 n1 = nvariables pows1 n2 = nvariables pows2 (e1', e2') = if n1 < n2 then (growSequence e1 n1 n2, e2) else (e1, growSequence e2 n2 n1) makePowers :: Exponents -> Powers makePowers expnts = powerize s where s = dropWhileR (== 0) expnts -- | An object of type @Spray a@ represents a multivariate polynomial whose -- coefficients are represented by the objects of type @a@, which must have -- a ring instance in order that we can add and multiply two polynomials. type Spray a = HashMap Powers a -- | Most often, one deals with sprays with rational coefficients, so we -- dedicate a type alias for such sprays. type QSpray = Spray Rational -- | The type `Rational'` is helpful when dealing with `OneParameterSpray` -- sprays, but this type of sprays lost its interest in version 0.4.0.0 -- (see CHANGELOG or README). type QSpray' = Spray Rational' type SafeSpray a = HashMap Exponents a -- An object of type @Term a@ represents a term of a @Spray a@ spray. Applying -- @Data.HashMap.Strict.toList@ to a @Spray a@ spray yields a list of @Term a@ -- objects. This type has probably no interest for the user, it is exported -- because it possibly has an interest for internal usage in a package using -- __hspray__. type Term a = (Powers, a) instance (AlgRing.C a, Eq a) => FunctionLike (Spray a) where type BaseRing (Spray a) = a -- type VariablesType (Spray a) = Spray a -- (*^) :: a -> Spray a -> Spray a (*^) lambda pol = lambda AlgMod.*> pol -- (+>) :: a -> Spray a -> Spray a (+>) x spray = if x == AlgAdd.zero then spray else addTerm spray (nullPowers, x) -- evaluate :: Spray a -> [a] -> a evaluate spray xyz = if length xyz >= numberOfVariables spray then evalSprayHelper xyz spray else error "evaluate: not enough values provided." -- substitute :: [Maybe a] -> Spray a -> Spray a substitute subs spray = if length subs >= n then spray' else error "substitute: incorrect length of the substitutions list." where n = numberOfVariables spray terms = HM.toList spray spray' = sumTerms (map substituteTerm terms) substituteTerm :: Term a -> Term a substituteTerm (powers, coeff) = (powers'', coeff') where pows = exponents powers nv = nvariables powers indices = findIndices isJust (take nv subs) pows' = [fromIntegral (pows `index` i) | i <- indices] xyz = [fromJust (subs !! i) | i <- indices] coeff' = coeff AlgRing.* AlgRing.product (zipWith (AlgRing.^) xyz pows') f i a = if i `elem` indices then 0 else a pows'' = S.mapWithIndex f pows powers'' = makePowers pows'' -- changeVariables :: Spray a -> [Spray a] -> Spray a changeVariables = composeSpray -- numberOfVariables :: Spray a -> Int numberOfVariables spray = if null powers then 0 else maximum (map nvariables powers) where powers = HM.keys spray -- permuteVariables :: [Int] -> Spray a -> Spray a permuteVariables permutation spray = if isPermutation permutation && n' >= n then spray' else error "permuteVariables: invalid permutation." where spray' = if isConstant spray then spray else HM.fromList (zip powers' coeffs) n = numberOfVariables spray n' = maximum permutation isPermutation pmtn = (not . null) pmtn && minimum pmtn == 1 && length (nub pmtn) == n' intmap = IM.fromList (zip permutation [1 .. n']) invpermutation = [intmap IM.! i | i <- [1 .. n']] permuteSeq x = S.mapWithIndex (\i _ -> x `index` (invpermutation !! i - 1)) x (powers, coeffs) = unzip (HM.toList spray) f pows = let expnts = (permuteSeq . growSequence' n') (exponents pows) in makePowers expnts powers' = map f powers -- swapVariables :: (Int, Int) -> Spray a -> Spray a swapVariables (i, j) spray = if i>=1 && j>=1 then spray' else error "swapVariables: invalid indices." where spray' = if isConstant spray then spray else HM.fromList (zip powers' coeffs) n = maximum [numberOfVariables spray, i, j] f k | k == i = j | k == j = i | otherwise = k transposition = map f [1 .. n] permuteSeq x = S.mapWithIndex (\ii _ -> x `index` (transposition !! ii - 1)) x (powers, coeffs) = unzip (HM.toList spray) g pows = let expnts = (permuteSeq . growSequence' n) (exponents pows) in makePowers expnts powers' = map g powers -- involvesVariable :: Spray a -> Int -> Bool involvesVariable spray i = any f (allExponents spray) where f expnts = let p = S.lookup (i - 1) expnts in isJust p && p /= Just 0 -- dropVariables :: Int -> Spray a -> Spray a dropVariables n = HM.mapKeys f where f (Powers exps nv) = Powers (S.drop n exps) (nv - n) -- derivative :: Int -> Spray a -> Spray a derivative i p = if i >= 1 then removeZeroTerms $ HM.fromListWith (AlgAdd.+) terms else error "derivative: invalid index." where terms = [ derivTerm term | term <- HM.toList p ] derivTerm :: Term a -> Term a derivTerm (pows, coef) = if i' >= S.length expts then (nullPowers, AlgAdd.zero) else (pows', coef') where i' = i - 1 expts = exponents pows expt_i = expts `index` i' expts' = adjust (subtract 1) i' expts coef' = expt_i .^ coef pows' = makePowers expts' -- | addition of two sprays addSprays :: (AlgAdd.C a, Eq a) => Spray a -> Spray a -> Spray a addSprays p q = removeZeroTerms $ HM.unionWith (AlgAdd.+) p q -- HM.foldlWithKey' f p q -- where -- f s powers coef = HM.insertWith (AlgAdd.+) powers coef s -- | addition of a term to a spray addTerm :: (AlgAdd.C a, Eq a) => Spray a -> Term a -> Spray a addTerm spray (powers, coeff) = if getCoefficient' powers spray AlgAdd.+ coeff == AlgAdd.zero then HM.delete powers spray else HM.insertWith (AlgAdd.+) powers coeff spray -- | sum list of terms sumTerms :: (Eq a, AlgAdd.C a) => [Term a] -> Spray a sumTerms = removeZeroTerms . HM.fromListWith (AlgAdd.+) -- | opposite spray negateSpray :: AlgAdd.C a => Spray a -> Spray a negateSpray = HM.map AlgAdd.negate -- | scale a spray by a scalar scaleSpray :: (AlgRing.C a, Eq a) => a -> Spray a -> Spray a scaleSpray lambda p = removeZeroTerms $ HM.map (lambda AlgRing.*) p -- | multiply two terms multTerm :: AlgRing.C a => Term a -> Term a -> Term a multTerm (pows1, coef1) (pows2, coef2) = (pows, coef1 AlgRing.* coef2) where (expts1', expts2') = harmonize (pows1, pows2) expts = S.zipWith (+) expts1' expts2' pows = makePowers expts -- | multiply a spray by a term multSprayByTerm :: (Eq a, AlgRing.C a) => Spray a -> Term a -> Spray a multSprayByTerm spray term = removeZeroTerms $ HM.fromList prods where prods = [multTerm trm term | trm <- HM.toList spray] -- | multiply two sprays multSprays :: (AlgRing.C a, Eq a) => Spray a -> Spray a -> Spray a multSprays p q = removeZeroTerms $ HM.fromListWith (AlgAdd.+) prods where p' = HM.toList p q' = HM.toList q prods = [ multTerm mp mq | mp <- p', mq <- q' ] instance (AlgAdd.C a, Eq a) => AlgAdd.C (Spray a) where (+) :: Spray a -> Spray a -> Spray a p + q = addSprays p q zero :: Spray a zero = HM.empty negate :: Spray a -> Spray a negate = negateSpray instance (AlgRing.C a, Eq a) => AlgMod.C a (Spray a) where (*>) :: a -> Spray a -> Spray a lambda *> p = scaleSpray lambda p instance (AlgRing.C a, Eq a) => AlgRightMod.C a (Spray a) where (<*) :: Spray a -> a -> Spray a p <* lambda = scaleSpray lambda p instance (AlgRing.C a, Eq a) => AlgRing.C (Spray a) where (*) :: Spray a -> Spray a -> Spray a p * q = multSprays p q one :: Spray a one = HM.singleton (Powers S.empty 0) AlgRing.one infixr 7 /^ -- | Divides a spray by a scalar; you can equivalently use `(/>)` if the type -- of the scalar is not ambiguous (/^) :: (AlgField.C a, Eq a) => Spray a -> a -> Spray a (/^) spray lambda = spray /> lambda -- | remove zero terms of a spray removeZeroTerms :: (AlgAdd.C a, Eq a) => Spray a -> Spray a removeZeroTerms = HM.filter (/= AlgAdd.zero) -- | helper function for lone and lone' lonePower :: Int -> Int -> Powers lonePower n p = if n == 0 then nullPowers else Powers (S.replicate (n - 1) 0 |> p) n -- | The @n@-th polynomial variable @x_n@ as a spray; one usually builds a -- spray by introducing these variables and combining them with the arithmetic -- operations -- -- >>> x = lone 1 :: Spray Int -- >>> y = lone 2 :: Spray Int -- >>> spray = 2*^x^**^2 ^-^ 3*^y -- >>> putStrLn $ prettyNumSpray spray -- 2*x^2 - 3*y -- -- prop> lone 0 == unitSpray lone :: AlgRing.C a => Int -> Spray a lone n = if n >= 0 then HM.singleton (lonePower n 1) AlgRing.one else error "lone: invalid index." -- | The @n@-th polynomial variable for rational sprays; this is just a -- specialization of `lone` qlone :: Int -> QSpray qlone = lone -- | The spray @x_n^p@; more efficient than exponentiating @lone n@ -- -- prop> lone' 2 10 = lone 2 ^**^ 10 lone' :: AlgRing.C a => Int -- ^ index -> Int -- ^ exponent -> Spray a lone' n p | n < 0 = error "lone': invalid index." | p < 0 = error "lone': invalid exponent" | otherwise = HM.singleton (lonePower n p) AlgRing.one -- | The rational spray @x_n^p@ qlone' :: Int -- ^ index -> Int -- ^ exponent -> QSpray qlone' = lone' loneTerm' :: AlgRing.C a => Int -> Int -> Term a loneTerm' n p = (lonePower n p, AlgRing.one) -- | Monomial spray, e.g. @monomial [(1,4),(3,2)]@ is @x^4.z^2@; indices -- and exponents must be positive but this is not checked -- prop> monomial [(1, 4), (3, 2)] == (lone 1 ^**^ 4) ^*^ (lone 3 ^**^ 2) monomial :: AlgRing.C a => [(Int, Int)] -- ^ list of (index, exponent); duplicates are deleted -> Spray a monomial nps = if null nps then unitSpray else HM.singleton (powerize expnts) AlgRing.one where nps' = nub nps nv = maximum (map fst nps') expnts = S.fromList $ [fromMaybe 0 (lookup i nps') | i <- [1 .. nv]] -- | Monomial rational spray, a specialization of 'monomial' -- -- prop> qmonomial [(1, 4), (3, 2)] == (qlone 1 ^**^ 4) ^*^ (qlone 3 ^**^ 2) qmonomial :: [(Int, Int)] -> QSpray qmonomial = monomial -- | The unit spray -- -- prop> spray ^*^ unitSpray == spray unitSpray :: AlgRing.C a => Spray a unitSpray = HM.singleton nullPowers AlgRing.one -- | The null spray -- -- prop> spray ^+^ zeroSpray == spray zeroSpray :: (Eq a, AlgAdd.C a) => Spray a zeroSpray = AlgAdd.zero -- | Whether a spray is the zero spray isZeroSpray :: Spray a -> Bool isZeroSpray = HM.null -- | Constant spray -- -- prop> constantSpray 3 == 3 *^ unitSpray constantSpray :: (Eq a, AlgAdd.C a) => a -> Spray a constantSpray c = if c == AlgAdd.zero then HM.empty else HM.singleton nullPowers c -- | Get coefficient of a term of a spray -- -- >>> x = lone 1 :: Spray Int -- >>> y = lone 2 :: Spray Int -- >>> z = lone 3 :: Spray Int -- >>> p = 2 *^ (2 *^ (x^**^3 ^*^ y^**^2)) ^+^ 4*^z ^+^ 5*^unitSpray -- >>> getCoefficient [3, 2] p -- coefficient of x^3.y^2 -- 4 -- >>> getCoefficient [3, 2, 0] p -- same as getCoefficient [3, 2] p -- 4 -- >>> getCoefficient [0, 4] p -- coefficient of y^4 -- 0 getCoefficient :: AlgAdd.C a => [Int] -- ^ sequence of exponents; trailing zeros are dropped -> Spray a -- ^ spray -> a -- ^ coefficient of the monomial given by the sequence of exponents getCoefficient expnts = getCoefficient' powers where powers = makePowers (S.fromList expnts) getCoefficient' :: AlgAdd.C a => Powers -> Spray a -> a getCoefficient' powers spray = fromMaybe AlgAdd.zero (HM.lookup powers spray) -- | Get the constant term of a spray -- -- prop> getConstantTerm p == getCoefficient [] p getConstantTerm :: AlgAdd.C a => Spray a -> a getConstantTerm = getCoefficient' nullPowers -- | remove the constant term of a spray removeConstantTerm :: Spray a -> Spray a removeConstantTerm = HM.delete nullPowers -- | Whether a spray is constant; same as `isConstant` isConstantSpray :: (Eq a, AlgRing.C a) => Spray a -> Bool isConstantSpray = isConstant -- | helper function to unify evalSpray and evalSpraySpray evalSprayHelper :: forall a. AlgRing.C a => [a] -> Spray a -> a evalSprayHelper xyz spray = AlgAdd.sum $ map evalTerm (HM.toList spray) where evalTerm :: Term a -> a evalTerm (powers, coeff) = coeff AlgRing.* AlgRing.product (zipWith (AlgRing.^) xyz pows) where pows = DF.toList (fromIntegral <$> exponents powers) -- | Evaluates a spray; same as `evaluate` -- -- >>> x = lone 1 :: Spray Int -- >>> y = lone 2 :: Spray Int -- >>> spray = 2*^x^**^2 ^-^ 3*^y -- >>> evalSpray spray [2, 1] -- 5 evalSpray :: (Eq a, AlgRing.C a) => Spray a -> [a] -> a evalSpray = evaluate -- | Evaluates the coefficients of a spray with spray coefficients; -- same as `substituteParameters` evalSpraySpray :: (Eq a, AlgRing.C a) => Spray (Spray a) -> [a] -> Spray a evalSpraySpray spray xyz = if length xyz >= n then removeZeroTerms $ HM.map (evalSprayHelper xyz) spray else error "evalSpraySpray: not enough values provided." where n = maximum (HM.elems $ HM.map numberOfVariables spray) -- | spray from term fromTerm :: Term a -> Spray a fromTerm (pows, coeff) = HM.singleton pows coeff -- | Substitutes some values to some variables of a spray; same as `substitute` -- -- >>> x1 = lone 1 :: Spray Int -- >>> x2 = lone 2 :: Spray Int -- >>> x3 = lone 3 :: Spray Int -- >>> p = x1^**^2 ^-^ x2 ^+^ x3 ^-^ unitSpray -- >>> p' = substituteSpray [Just 2, Nothing, Just 3] p -- >>> putStrLn $ prettyNumSprayX1X2X3 "x" p' -- -x2 + 6 substituteSpray :: (Eq a, AlgRing.C a) => [Maybe a] -> Spray a -> Spray a substituteSpray = substitute -- | Converts a spray with rational coefficients to a spray with double -- coefficients (useful for evaluation) fromRationalSpray :: Spray Rational -> Spray Double fromRationalSpray = HM.map fromRational -- | Sustitutes the variables of a spray with some sprays; same as `changeVariables` -- -- >>> x = lone 1 :: Spray Int -- >>> y = lone 2 :: Spray Int -- >>> z = lone 3 :: Spray Int -- >>> p = x ^+^ y -- >>> q = composeSpray p [x ^+^ y ^+^ z, z] -- >>> putStrLn $ prettyNumSpray' q -- X + Y + 2*Z composeSpray :: forall a. (AlgRing.C a, Eq a) => Spray a -> [Spray a] -> Spray a composeSpray p = removeZeroTerms . evalSpray (identify p) where identify :: Spray a -> Spray (Spray a) identify = HM.map constantSpray -- | Creates a spray from a list of terms fromList :: (AlgRing.C a, Eq a) => [([Int], a)] -- ^ list of (exponents, coefficient) -> Spray a fromList x = removeZeroTerms $ HM.fromListWith (AlgAdd.+) $ map (first (makePowers . S.fromList)) x -- pretty stuff --------------------------------------------------------------- -- | Prints a spray; this function is exported for possible usage in other packages showSpray :: (a -> String) -- ^ function mapping a coefficient to a string, typically 'show' -> (String, String) -- ^ pair of braces to enclose the coefficients -> ([Exponents] -> [String]) -- ^ function mapping a list of exponents to a list of strings representing the monomials corresponding to these exponents -> Spray a -- ^ the spray to be printed -> String showSpray showCoef braces showMonomials spray = if isZeroSpray spray then "0" else unpack $ intercalate (pack " + ") stringTerms where terms = orderedTerms spray (powers, coeffs) = unzip terms stringMonomials = showMonomials (map exponents powers) stringTerms = zipWith f coeffs stringMonomials f coeff smonomial | smonomial == "" = pack scoeff' | scoeff == "" = pack smonomial | otherwise = pack $ scoeff' ++ "*" ++ smonomial where scoeff = showCoef coeff scoeff' = bracify braces scoeff -- | Prints a spray, with monomials shown as "x.z^2", and with -- a user-defined showing function for the coefficients showSprayXYZ :: (a -> String) -- ^ function mapping a coefficient to a string, typically 'show' -> (String, String) -- ^ used to enclose the coefficients, usually a pair of braces -> [String] -- ^ typically some letters, to print the variables -> Spray a -- ^ the spray to be printed -> String showSprayXYZ showCoef braces letters spray = if null letters then error "showSprayXYZ: empty list of strings." else showSpray showCoef braces (showMonomialsXYZ letters) spray -- | Prints a spray, with monomials shown as @\"x.z^2\"@, and with -- a user-defined showing function for the coefficients; this is the same as -- the function `showSprayXYZ` with the pair of braces @("(", ")")@ showSprayXYZ' :: (a -> String) -- ^ function mapping a coefficient to a string, typically 'show' -> [String] -- ^ typically some letters, to print the variables -> Spray a -- ^ the spray to be printed -> String showSprayXYZ' showCoef = showSprayXYZ showCoef ("(", ")") -- | Pretty form of a spray with monomials displayed in the style of @\"x.z^2\"@; -- you should rather use `prettyNumSprayXYZ` or `prettyQSprayXYZ` if your -- coefficients are numeric -- -- >>> x = lone 1 :: Spray Int -- >>> y = lone 2 :: Spray Int -- >>> z = lone 3 :: Spray Int -- >>> p = 2*^x ^+^ 3*^y^**^2 ^-^ 4*^z^**^3 -- >>> putStrLn $ prettySprayXYZ ["X", "Y", "Z"] p -- (2)*X + (3)*Y^2 + (-4)*Z^3 -- >>> putStrLn $ prettySprayXYZ ["X", "Y"] p -- (2)*X1 + (3)*X2^2 + (-4)*X3^3 prettySprayXYZ :: (Show a) => [String] -- ^ typically some letters, to print the variables -> Spray a -- ^ the spray to be printed -> String prettySprayXYZ = showSprayXYZ' show -- | Pretty form of a spray, with monomials shown as "x1.x3^2", and with -- a user-defined showing function for the coefficients showSprayX1X2X3 :: (a -> String) -- ^ function mapping a coefficient to a string, typically 'show' -> (String, String) -- ^ used to enclose the coefficients -> String -- ^ typically a letter, to print the non-indexed variables -> Spray a -- ^ the spray to be printed -> String showSprayX1X2X3 showCoef braces letter = showSpray showCoef braces (showMonomialsX1X2X3 letter) -- | Pretty form of a spray, with monomials shown as "x1.x3^2", and with -- a user-defined showing function for the coefficients; this is the same as -- the function `showSprayX1X2X3` with the pair of braces @("(", ")")@ used to -- enclose the coefficients showSprayX1X2X3' :: (a -> String) -- ^ function mapping a coefficient to a string, e.g. 'show' -> String -- ^ typically a letter, to print the non-indexed variables -> Spray a -- ^ the spray to be printed -> String showSprayX1X2X3' showCoef = showSprayX1X2X3 showCoef ("(", ")") -- | Pretty form of a spray with monomials displayed in the style of @\"x1.x3^2\"@; -- you should rather use `prettyNumSprayX1X2X3` or `prettyQSprayX1X2X3` if your -- coefficients are numeric -- -- >>> x = lone 1 :: Spray Int -- >>> y = lone 2 :: Spray Int -- >>> z = lone 3 :: Spray Int -- >>> spray = 2*^x ^+^ 3*^y^**^2 ^-^ 4*^z^**^3 -- >>> putStrLn $ prettySprayX1X2X3 "X" spray -- (2)*X1 + (3)*X2^2 + (-4)*X3^3 prettySprayX1X2X3 :: Show a => String -- ^ typically a letter, to print the non-indexed variables -> Spray a -- ^ the spray to be printed -> String prettySprayX1X2X3 = showSprayX1X2X3' show -- | Pretty form of a spray with monomials displayed in the style of @\"x.z^2\"@; -- you should rather use `prettyNumSpray` or `prettyQSpray` if you deal with -- sprays with numeric coefficients -- -- >>> x = lone 1 :: Spray Int -- >>> y = lone 2 :: Spray Int -- >>> z = lone 3 :: Spray Int -- >>> p = 2*^x ^+^ 3*^y^**^2 ^-^ 4*^z^**^3 -- >>> putStrLn $ prettySpray p -- (2)*x + (3)*y^2 + (-4)*z^3 -- >>> putStrLn $ prettySpray (p ^+^ lone 4) -- (2)*x1 + (3)*x2^2 + (-4)*x3^3 + x4 -- -- prop> prettySpray spray == prettySprayXYZ ["x", "y", "z"] spray prettySpray :: (Show a) => Spray a -> String prettySpray = prettySprayXYZ ["x", "y", "z"] -- | Pretty form of a spray, with monomials shown as @\"x1.x3^2\"@; use -- `prettySprayX1X2X3` to change the letter (or `prettyNumSprayX1X2X3` -- or `prettyQSprayX1X2X3` if the coefficients are numeric) -- -- >>> x = lone 1 :: Spray Int -- >>> y = lone 2 :: Spray Int -- >>> z = lone 3 :: Spray Int -- >>> p = 2*^x ^+^ 3*^y^**^2 ^-^ 4*^z^**^3 -- >>> putStrLn $ prettySpray' p -- (2)*x1 + (3)*x2^2 + (-4)*x3^3 prettySpray' :: Show a => Spray a -> String prettySpray' = prettySprayX1X2X3 "x" -- | showMonomialOld "x" [0, 2, 1] = x^(0, 2, 1) showMonomialsOld :: String -> [Exponents] -> [String] showMonomialsOld var = map (showMonomialOld var) where showMonomialOld :: String -> Exponents -> String showMonomialOld a pows = unpack $ append (pack x) (cons '(' $ snoc string ')') where x = a ++ "^" string = intercalate (pack ", ") (map (pack . show) (DF.toList pows)) -- | Pretty form of a spray; you will probably prefer `prettySpray` or `prettySpray'` -- -- >>> x = lone 1 :: Spray Int -- >>> y = lone 2 :: Spray Int -- >>> z = lone 3 :: Spray Int -- >>> p = 2*^x ^+^ 3*^y^**^2 ^-^ 4*^z^**^3 -- >>> putStrLn $ prettySpray'' "x" p -- (2)*x^(1) + (3)*x^(0, 2) + (-4)*x^(0, 0, 3) prettySpray'' :: Show a => String -- ^ a string denoting the variables, e.g. \"x\" -> Spray a -- ^ the spray -> String prettySpray'' var = showSpray show ("(", ")") (showMonomialsOld var) -- | Show a spray with numeric coefficients; this function is exported for -- possible usage in other packages showNumSpray :: (Num a, Ord a) => ([Exponents] -> [String]) -- ^ function mapping a list of monomial exponents to a list of strings representing the monomials -> (a -> String) -- ^ function mapping a positive coefficient to a string -> Spray a -> String showNumSpray showMonomials showCoeff spray = if isZeroSpray spray then "0" else concat $ zipWith (++) stringSigns stringTerms where terms = orderedTerms spray coeffs = map snd terms (firstCoeff, otherCoeffs) = fromJust (uncons coeffs) firstSign = if firstCoeff > 0 then "" else "-" otherSigns = map (\x -> if x > 0 then " + " else " - ") otherCoeffs stringSigns = firstSign : otherSigns absCoeffs = map abs coeffs powers = map (exponents . fst) terms stringMonomials = showMonomials powers stringTerms = zipWith f absCoeffs stringMonomials f acoeff smonomial | smonomial == "" = showCoeff acoeff | scoeff == "" = smonomial | otherwise = scoeff ++ "*" ++ smonomial where scoeff = if acoeff == 1 then "" else showCoeff acoeff -- | showMonomialX1X2X3 "X" [0, 2, 1] = "X2^2.X3" showMonomialX1X2X3 :: String -> Exponents -> Text showMonomialX1X2X3 x pows = x1x2x3 where f i p | p == 0 = pack "" | p == 1 = pack $ x ++ show i | otherwise = pack $ x ++ show i ++ "^" ++ show p indices = S.findIndicesL (/= 0) pows x1x2x3 = intercalate (pack ".") [f (i+1) (pows `index` i) | i <- indices] -- | showMonomialsX1X2X3 "X" [[0, 2, 1], [1, 2]] = ["X2^2.X3", "X1.X2"] showMonomialsX1X2X3 :: String -> [Exponents] -> [String] showMonomialsX1X2X3 x = map (unpack . showMonomialX1X2X3 x) -- | showMonomialXYZ ["X", "Y", "Z"] 3 [1, 2, 1] = X.Y^2.Z -- showMonomialXYZ ["X", "Y", "Z"] 3 [1, 2, 1, 2] = X1.X2^2.X3.X4^2 showMonomialXYZ :: [String] -> Int -> Exponents -> Text showMonomialXYZ letters n pows = if n <= length letters then xyz else showMonomialX1X2X3 (letters !! 0) pows where f letter p | p == 0 = pack "" | p == 1 = pack letter | otherwise = pack $ letter ++ "^" ++ show p indices = S.findIndicesL (/= 0) pows xyz = intercalate (pack ".") [f (letters !! i) (pows `index` i) | i <- indices] -- | showMonomialsXYZ ["X", "Y", "Z"] [[0, 2, 1], [1, 2]] = ["Y^2.Z", "X.Y^2"] showMonomialsXYZ :: [String] -> [Exponents] -> [String] showMonomialsXYZ letters powers = map (unpack . showMonomialXYZ letters n) powers where n = maximum (map S.length powers) -- | Pretty form of a spray with numeric coefficients, printing monomials as @\"x1.x3^2\"@ -- -- >>> x = lone 1 :: Spray Int -- >>> y = lone 2 :: Spray Int -- >>> z = lone 3 :: Spray Int -- >>> p = 2*^x ^+^ 3*^y^**^2 ^-^ 4*^z^**^3 -- >>> putStrLn $ prettyNumSprayX1X2X3 "x" p -- 2*x1 + 3*x2^2 - 4*x3^3 prettyNumSprayX1X2X3 :: (Num a, Ord a, Show a) => String -- ^ usually a letter such as @\"x\"@ to denote the non-indexed variables -> Spray a -> String prettyNumSprayX1X2X3 x = showNumSpray (showMonomialsX1X2X3 x) show -- | Pretty form of a spray with numeric coefficients, printing monomials as @\"x.z^2\"@ -- if possible, i.e. if enough letters are provided, otherwise as @\"x1.x3^2\"@ -- -- >>> x = lone 1 :: Spray Int -- >>> y = lone 2 :: Spray Int -- >>> z = lone 3 :: Spray Int -- >>> w = lone 4 :: Spray Int -- >>> p = 2*^x ^+^ 3*^y^**^2 ^-^ 4*^z^**^3 -- >>> putStrLn $ prettyNumSprayXYZ ["x","y","z"] p -- 2*x + 3*y^2 - 4*z^3 -- >>> putStrLn $ prettyNumSprayXYZ ["x","y","z"] (p ^+^ w) -- 2*x1 + 3*x2^2 - 4*x3^3 + x4 -- >>> putStrLn $ prettyNumSprayXYZ ["a","b","c"] (p ^+^ w) -- 2*a1 + 3*a2^2 - 4*a3^3 + a4 prettyNumSprayXYZ :: (Num a, Ord a, Show a) => [String] -- ^ usually some letters, denoting the variables -> Spray a -> String prettyNumSprayXYZ letters = showNumSpray (showMonomialsXYZ letters) show -- | helper function for showQSpray showRatio :: Rational -> String showRatio q = if d == 1 then show n else "(" ++ show n ++ "/" ++ show d ++ ")" where n = DR.numerator q d = DR.denominator q -- | helper function for showQSpray' showRatio' :: (Eq a, Num a, Show a) => NumberRatio.T a -> String showRatio' q = if d == 1 then show n else "(" ++ show n ++ "/" ++ show d ++ ")" where n = NumberRatio.numerator q d = NumberRatio.denominator q -- | Prints a `QSpray`; for internal usage but exported for usage in other packages showQSpray :: ([Exponents] -> [String]) -- ^ function printing monomials -> QSpray -> String showQSpray showMonomials = showNumSpray showMonomials showRatio -- | Prints a `QSpray'`; for internal usage but exported for usage in other packages showQSpray' :: ([Exponents] -> [String]) -- ^ function mapping a list of monomials exponents to a list of strings -> QSpray' -> String showQSpray' showMonomials = showNumSpray showMonomials showRatio' -- | Pretty form of a spray with rational coefficients, printing monomials in -- the style of @\"x1.x3^2\"@ -- -- >>> x = lone 1 :: QSpray -- >>> y = lone 2 :: QSpray -- >>> z = lone 3 :: QSpray -- >>> p = 2*^x ^+^ 3*^y^**^2 ^-^ (4%3)*^z^**^3 -- >>> putStrLn $ prettyQSprayX1X2X3 "x" p -- 2*x1 + 3*x2^2 - (4/3)*x3^3 prettyQSprayX1X2X3 :: String -- ^ usually a letter such as @\"x\"@, to denote the non-indexed variables -> QSpray -> String prettyQSprayX1X2X3 x = showQSpray (showMonomialsX1X2X3 x) -- | Same as `prettyQSprayX1X2X3` but for a `QSpray'` spray prettyQSprayX1X2X3' :: String -- ^ usually a letter such as @\"x\"@, to denote the non-indexed variables -> QSpray' -> String prettyQSprayX1X2X3' x = showQSpray' (showMonomialsX1X2X3 x) -- | Pretty form of a spray with rational coefficients, printing monomials in -- the style of @\"x.z^2\"@ with the provided letters if possible, i.e. if enough -- letters are provided, otherwise in the style @\"x1.x3^2\"@, taking the first -- provided letter to denote the non-indexed variables -- -- >>> x = lone 1 :: QSpray -- >>> y = lone 2 :: QSpray -- >>> z = lone 3 :: QSpray -- >>> p = 2*^x ^+^ 3*^y^**^2 ^-^ (4%3)*^z^**^3 -- >>> putStrLn $ prettyQSprayXYZ ["x","y","z"] p -- 2*x + 3*y^2 - (4/3)*z^3 -- >>> putStrLn $ prettyQSprayXYZ ["x","y"] p -- 2*x1 + 3*x2^2 - (4%3)*x3^3 -- >>> putStrLn $ prettyQSprayXYZ ["a","b"] p -- 2*a1 + 3*a2^2 - (4/3)*a3^3 prettyQSprayXYZ :: [String] -- ^ usually some letters, to denote the variables -> QSpray -> String prettyQSprayXYZ letters = showQSpray (showMonomialsXYZ letters) -- | Same as `prettyQSprayXYZ` but for a `QSpray'` spray prettyQSprayXYZ' :: [String] -- ^ usually some letters, to denote the variables -> QSpray' -> String prettyQSprayXYZ' letters = showQSpray' (showMonomialsXYZ letters) -- | Pretty printing of a spray with rational coefficients -- prop> prettyQSpray == prettyQSprayXYZ ["x", "y", "z"] prettyQSpray :: QSpray -> String prettyQSpray = prettyQSprayXYZ ["x", "y", "z"] -- | Pretty printing of a spray with rational coefficients -- prop> prettyQSpray'' == prettyQSprayXYZ ["X", "Y", "Z"] prettyQSpray'' :: QSpray -> String prettyQSpray'' = prettyQSprayXYZ ["X", "Y", "Z"] -- | Pretty printing of a spray with rational coefficients -- prop> prettyQSpray' == prettyQSprayXYZ' ["x", "y", "z"] prettyQSpray' :: QSpray' -> String prettyQSpray' = prettyQSprayXYZ' ["x", "y", "z"] -- | Pretty printing of a spray with rational coefficients -- prop> prettyQSpray''' == prettyQSprayXYZ' ["X", "Y", "Z"] prettyQSpray''' :: QSpray' -> String prettyQSpray''' = prettyQSprayXYZ' ["X", "Y", "Z"] -- | Pretty printing of a spray with numeric coefficients -- prop> prettyNumSpray == prettyNumSprayXYZ ["x", "y", "z"] prettyNumSpray :: (Num a, Ord a, Show a) => Spray a -> String prettyNumSpray = prettyNumSprayXYZ ["x", "y", "z"] -- | Pretty printing of a spray with numeric coefficients -- prop> prettyNumSpray' == prettyNumSprayXYZ ["X", "Y", "Z"] prettyNumSpray' :: (Num a, Ord a, Show a) => Spray a -> String prettyNumSpray' = prettyNumSprayXYZ ["X", "Y", "Z"] -- misc ----------------------------------------------------------------------- -- | spray as safe spray safeSpray :: Spray a -> SafeSpray a safeSpray = HM.mapKeys exponents getCoefficient'' :: AlgAdd.C a => Exponents -> SafeSpray a -> a getCoefficient'' powers spray' = fromMaybe AlgAdd.zero (HM.lookup powers spray') getConstantTerm'' :: AlgAdd.C a => SafeSpray a -> a getConstantTerm'' = getCoefficient'' S.empty removeConstantTerm'' :: SafeSpray a -> SafeSpray a removeConstantTerm'' = HM.delete S.empty -- | Sum of sprays sumOfSprays :: (Eq a, AlgAdd.C a) => [Spray a] -> Spray a sumOfSprays = AlgAdd.sum -- | Product of sprays productOfSprays :: (Eq a, AlgRing.C a) => [Spray a] -> Spray a productOfSprays = AlgRing.product -- | ordered terms of a spray orderedTerms :: Spray a -> [Term a] orderedTerms spray = sortBy (flip compare `on` (exponents . fst)) (HM.toList spray) -- | Spray as a list toList :: Spray a -> [([Int], a)] toList p = HM.toList $ HM.mapKeys (DF.toList . exponents) p -- | Bombieri spray (for internal usage in the \'__scubature__\' package) bombieriSpray :: (Eq a, AlgAdd.C a) => Spray a -> Spray a bombieriSpray = HM.mapWithKey f where f pows = times (pfactorial $ exponents pows) pfactorial pows = product $ DF.toList $ factorial <$> S.filter (/= 0) pows factorial n = product [2 .. n] times k x = k .^ x -- | Whether two sprays are equal up to a scalar factor collinearSprays :: (Eq a, AlgRing.C a) => Spray a -> Spray a -> Bool collinearSprays spray1 spray2 = isZeroSpray spray1 && isZeroSpray spray2 || (not . isZeroSpray) spray1 && (not . isZeroSpray) spray2 && snd (leadingTerm spray1) *^ spray2 == snd (leadingTerm spray2) *^ spray1 -- | Checks whether the multivariate polynomial defined by a spray is homogeneous -- and also returns the degree in case this is true isHomogeneousSpray :: (Eq a, AlgRing.C a) => Spray a -> (Bool, Maybe Int) isHomogeneousSpray spray | isConstant spray = (True, Just 0) | getConstantTerm spray /= AlgAdd.zero = (False, Nothing) | otherwise = (check, deg) where degrees = map (sum . exponents) (HM.keys spray) check = allSame degrees deg = if check then Just (degrees !! 0) else Nothing -- | Get all the exponents of a spray allExponents :: Spray a -> [Exponents] allExponents spray = map exponents (HM.keys spray) -- | Get all the coefficients of a spray allCoefficients :: Spray a -> [a] allCoefficients = HM.elems -- division stuff ------------------------------------------------------------- -- | index of the maximum of a list -- maxWithIndex :: Ord a => [a] -> (Int, a) -- maxWithIndex = maximumBy (comparing snd) . zip [0 .. ] -- | leading term of a spray leadingTerm :: Spray a -> Term a leadingTerm p = (biggest, p HM.! biggest) where powers = HM.keys p biggest = maximumBy (comparing exponents) powers -- (i, biggest) = maxWithIndex powers -- biggest = powers !! i -- | whether a term divides another term divides :: Term a -> Term a -> Bool divides (powsP, _) (powsQ, _) = nvP <= nvQ && lower where nvP = nvariables powsP nvQ = nvariables powsQ expntsP = exponents powsP expntsQ = exponents powsQ lower = DF.all (uncurry (<=)) (S.zip expntsP expntsQ) -- | quotient of term Q by term P, assuming P divides Q quotient :: AlgField.C a => Term a -> Term a -> Term a quotient (powsQ, coeffQ) (powsP, coeffP) = (pows, coeff) where (expntsP, expntsQ) = harmonize (powsP, powsQ) expnts = S.zipWith (-) expntsQ expntsP pows = makePowers expnts coeff = coeffQ AlgField./ coeffP -- | Remainder of the division of a spray by a list of divisors, -- using the lexicographic ordering of the monomials sprayDivisionRemainder :: forall a. (Eq a, AlgField.C a) => Spray a -> [Spray a] -> Spray a sprayDivisionRemainder p qs = if n == 0 then error "sprayDivisionRemainder: the list of divisors is empty." else ogo p zeroSpray where n = length qs qsltqs = zip qs (map leadingTerm qs) g :: Term a -> Spray a -> Spray a -> (Spray a, Spray a) g lts s r = (HM.delete (fst lts) s, addTerm r lts) go :: Term a -> Spray a -> Spray a -> Int -> Bool -> (Spray a, Spray a) go lts !s r !i !divoccured | divoccured = (s, r) | i == n = g lts s r | otherwise = go lts news r (i+1) newdivoccured where (q, ltq) = qsltqs !! i newdivoccured = divides ltq lts news = if newdivoccured then s ^-^ multSprayByTerm q (quotient lts ltq) else s ogo :: Spray a -> Spray a -> Spray a ogo !s !r | isZeroSpray s = r | otherwise = ogo s' r' where (s', r') = go (leadingTerm s) s r 0 False -- | Division of a spray by a spray sprayDivision :: (Eq a, AlgField.C a) => Spray a -- ^ dividand -> Spray a -- ^ divisor -> (Spray a, Spray a) -- ^ (quotient, remainder) sprayDivision sprayA sprayB = if isConstant sprayB then if isZeroSpray sprayB then error "sprayDivision: division by zero." else let c = getConstantTerm sprayB in (sprayA /> c, zeroSpray) else sprayDivision0 sprayA sprayB sprayDivision0 :: forall a. (Eq a, AlgField.C a) => Spray a -- ^ dividand -> Spray a -- ^ divisor -> (Spray a, Spray a) -- ^ (quotient, remainder) sprayDivision0 sprayA sprayB = ogo sprayA zeroSpray zeroSpray where ltB = leadingTerm sprayB ogo :: Spray a -> Spray a -> Spray a -> (Spray a, Spray a) ogo !p !q !r | isZeroSpray p = (q, r) | otherwise = ogo p' q' r' where ltp = leadingTerm p (p', q', r') = if divides ltB ltp then (newp, newq, r) else (HM.delete (fst ltp) p, q, addTerm r ltp) qtnt = quotient ltp ltB newp = p ^-^ multSprayByTerm sprayB qtnt newq = addTerm q qtnt -- | division of univariate sprays with degree(dividend) >= degree(divisor) univariateSprayDivision :: forall a. (Eq a, AlgField.C a) => Spray a -- ^ dividand -> Spray a -- ^ divisor -> (Spray a, Spray a) -- ^ (quotient, remainder) univariateSprayDivision sprayA sprayB = if isConstant sprayB then let c = getConstantTerm sprayB in (sprayA /> c, zeroSpray) else ogo sprayA zeroSpray zeroSpray where (powsLTB, coeffLTB) = leadingTerm sprayB degB = exponents powsLTB expntLTB = degB `index` 0 ogo :: Spray a -> Spray a -> Spray a -> (Spray a, Spray a) ogo !p !q !r | isZeroSpray p = (q, r) | otherwise = ogo p' q' r' where (powsLTP, coeffLTP) = leadingTerm p degP = exponents powsLTP (p', q', r') = if degB <= degP then (newp, newq, r) else (zeroSpray, q, r ^+^ p) newp = p ^-^ multSprayByTerm sprayB qtnt newq = addTerm q qtnt qtnt = (pows, coeff) expntLTP = degP `index` 0 pows = if expntLTP == expntLTB then nullPowers else Powers (S.singleton (expntLTP - expntLTB)) 1 coeff = coeffLTP AlgField./ coeffLTB -- Groebner stuff ------------------------------------------------------------- -- | slight modification of `sprayDivisionRemainder` to speed up groebner00 sprayDivisionRemainder' :: forall a. (Eq a, AlgField.C a) => Spray a -> HashMap Int (Spray a, Term a) -> Spray a sprayDivisionRemainder' p qsltqs = ogo p zeroSpray where n = HM.size qsltqs g :: Term a -> Spray a -> Spray a -> (Spray a, Spray a) g lts s r = (HM.delete (fst lts) s, addTerm r lts) go :: Term a -> Spray a -> Spray a -> Int -> Bool -> (Spray a, Spray a) go lts !s r !i !divoccured | divoccured = (s, r) | i == n = g lts s r | otherwise = go lts news r (i+1) newdivoccured where (q, ltq) = qsltqs HM.! i newdivoccured = divides ltq lts news = if newdivoccured then s ^-^ multSprayByTerm q (quotient lts ltq) else s ogo :: Spray a -> Spray a -> Spray a ogo !s !r | isZeroSpray s = r | otherwise = ogo s' r' where (s', r') = go (leadingTerm s) s r 0 False -- combinations of two among n combn2 :: Int -> Int -> HashMap Int (Int, Int) combn2 n s = HM.fromList (zip [0 .. ] (zip row1 row2)) where range1 = [1 .. n-1] row1 = drop s $ concatMap (\i -> [0 .. i-1]) range1 row2 = drop s $ concatMap (\i -> replicate i i) range1 -- the "S polynomial" sPolynomial :: (Eq a, AlgField.C a) => (Spray a, Term a) -> (Spray a, Term a) -> Spray a sPolynomial (p, (lpowsP, lcoefP)) (q, (lpowsQ, lcoefQ)) = multSprayByTerm p wp ^-^ multSprayByTerm q wq where (lexpntsP, lexpntsQ) = harmonize (lpowsP, lpowsQ) gamma = S.zipWith max lexpntsP lexpntsQ betaP = S.zipWith (-) gamma lexpntsP betaQ = S.zipWith (-) gamma lexpntsQ wp = (makePowers betaP, AlgField.recip lcoefP) wq = (makePowers betaQ, AlgField.recip lcoefQ) -- | groebner basis, not minimal and not reduced groebner00 :: forall a. (Eq a, AlgField.C a) => [Spray a] -> [Spray a] groebner00 sprays = go 0 j0 combins0 spraysMap where j0 = length sprays combins0 = combn2 j0 0 ltsprays = map leadingTerm sprays spraysltsprays = zip sprays ltsprays spraysMap = HM.fromList (zip [0 .. j0-1] spraysltsprays) go :: Int -> Int -> HashMap Int (Int, Int) -> HashMap Int (Spray a, Term a) -> [Spray a] go !i !j !combins !gpolysMap | j == 100 = error "groebnerBasis: stopped because reached the limit; please fill an issue." | i == length combins = map fst (HM.elems gpolysMap) | otherwise = go i' j' combins' gpolysMap' where (k, l) = combins HM.! i sfg = sPolynomial (gpolysMap HM.! k) (gpolysMap HM.! l) sbarfg = sprayDivisionRemainder' sfg gpolysMap ltsbarfg = leadingTerm sbarfg (i', j', gpolysMap', combins') = if isZeroSpray sbarfg then (i + 1, j, gpolysMap, combins) else ( 0 , j+1 , HM.insert j (sbarfg, ltsbarfg) gpolysMap , combn2 (j + 1) (i + 1) ) -- | groebner basis, minimal but not reduced groebner0 :: forall a. (Eq a, AlgField.C a) => [Spray a] -> [Spray a] groebner0 sprays = if n <= 1 then sprays else [basis00 !! k | k <- [0 .. n-1] \\ discard] where n = length basis00 basis00 = groebner00 sprays go :: Int -> [Int] -> [Int] go !i toRemove | i == n = toRemove | otherwise = go (i+1) toRemove' where ltf = leadingTerm (basis00 !! i) toDrop = toRemove ++ [i] igo :: Int -> Bool igo !j | j == n = False | j `elem` toDrop = igo (j+1) | otherwise = ok || igo (j+1) where ok = divides (leadingTerm (basis00 !! j)) ltf toRemove' = if igo 0 then toDrop else toRemove discard = go 0 [] -- | Reduces a Gröbner basis reduceGroebnerBasis :: forall a. (Eq a, AlgField.C a) => [Spray a] -> [Spray a] reduceGroebnerBasis gbasis = if length gbasis >= 2 then map reduction [0 .. n-1] else ngbasis where normalize :: Spray a -> Spray a normalize spray = spray /> coef where (_, coef) = leadingTerm spray ngbasis = map normalize gbasis n = length ngbasis reduction :: Int -> Spray a reduction i = sprayDivisionRemainder (ngbasis !! i) rest where rest = [ngbasis !! k | k <- [0 .. i-1] ++ [i+1 .. n-1]] -- | Gröbner basis, always minimal and possibly reduced -- -- prop> groebnerBasis sprays True == reduceGroebnerBasis (groebnerBasis sprays False) groebnerBasis :: forall a. (Eq a, AlgField.C a) => [Spray a] -- ^ list of sprays -> Bool -- ^ whether to return the reduced basis -> [Spray a] groebnerBasis sprays reduced = if reduced then reduceGroebnerBasis gbasis0 else map normalize gbasis0 where gbasis0 = groebner0 sprays normalize :: Spray a -> Spray a normalize spray = spray /> coef where (_, coef) = leadingTerm spray -- elementary symmetric polynomials ------------------------------------------- -- | combinations of k elements among a list combinationsOf :: Int -> [a] -> [[a]] combinationsOf _ [] = error "combinationsOf: should not happen." combinationsOf 1 as = map pure as combinationsOf k as@(_:xs) = run (l-1) (k-1) as $ combinationsOf (k-1) xs where l = length as run :: Int -> Int -> [a] -> [[a]] -> [[a]] run n i ys cs | n == i = map (ys ++) cs | otherwise = map (q:) cs ++ run (n-1) i qs (drop dc cs) where f :: [a] -> (a, [a]) f [] = error "combinationsOf: should not happen." f (b:bs) = (b, bs) (q, qs) = f (take (n-i+1) ys) dc = product [n-i+1 .. n-1] `div` product [1 .. i-1] -- | generates all permutations of a binary sequence permutationsBinarySequence :: Int -> Int -> [Seq Int] permutationsBinarySequence nzeros nones = let n = nzeros + nones in map (binarySequence n) (combinationsOf nones [0 .. n-1]) where binarySequence :: Int -> [Int] -> Seq Int binarySequence n combo = fromFunction n f where f :: Int -> Int f i = fromEnum (i `elem` combo) -- | Elementary symmetric polynomial -- -- >>> putStrLn $ prettySpray' (esPolynomial 3 2) -- (1)*x1.x2 + (1)*x1.x3 + (1)*x2.x3 esPolynomial :: (AlgRing.C a, Eq a) => Int -- ^ number of variables -> Int -- ^ index -> Spray a esPolynomial n k | k < 0 || n < 0 = error "esPolynomial: both arguments must be positive integers." | k > n = zeroSpray | k == 0 = unitSpray | otherwise = spray where perms = permutationsBinarySequence (n-k) k spray = HM.fromList $ map (\expts -> (makePowers expts, AlgRing.one)) perms -- | Power sum polynomial psPolynomial :: forall a. (AlgRing.C a, Eq a) => Int -- ^ number of variables -> Int -- ^ power -> Spray a psPolynomial n k | k < 0 || n < 0 = error "psPolynomial: both arguments must be positive integers." | k > n = AlgAdd.zero | k == 0 = n .^ unitSpray | otherwise = spray where spray = HM.fromList $ map f [1 .. n] f :: Int -> (Powers, a) f j = (Powers expts j, AlgRing.one) where expts = S.replicate (j-1) 0 |> k -- | Whether a spray is a symmetric polynomial, an inefficient algorithm -- (use the function with the same name in the __jackpolynomials__ package -- if you need efficiency) isSymmetricSpray :: forall a. (AlgField.C a, Eq a) => Spray a -> Bool isSymmetricSpray spray = check where n = numberOfVariables spray indices = [1 .. n] gPolys = [esPolynomial n i ^-^ lone (n + i) | i <- indices] gbasis = groebner0 gPolys spray' = removeConstantTerm spray g = sprayDivisionRemainder spray' gbasis check = not $ any (involvesVariable g) [1 .. n] -- | Whether a spray can be written as a polynomial of a given list of sprays; -- this polynomial is returned if this is true -- -- >>> x = lone 1 :: Spray Rational -- >>> y = lone 2 :: Spray Rational -- >>> p1 = x ^+^ y -- >>> p2 = x ^-^ y -- >>> p = p1 ^*^ p2 -- -- prop> isPolynomialOf p [p1, p2] == (True, Just $ x ^*^ y) isPolynomialOf :: forall a. (AlgField.C a, Eq a) => Spray a -> [Spray a] -> (Bool, Maybe (Spray a)) isPolynomialOf spray sprays | isConstant spray = (True, Just spray) | null sprays = error "isPolynomialOf: the list of sprays is empty." | otherwise = result where nov = numberOfVariables spray n = maximum $ map numberOfVariables sprays result | nov > n = (False, Nothing) | otherwise = (check, poly) where m = length sprays yPolys = [lone (n + i) | i <- [1 .. m]] gPolys = zipWith (^-^) sprays yPolys gbasis0 = groebner0 gPolys constantTerm = getConstantTerm spray spray' = removeConstantTerm spray g = sprayDivisionRemainder spray' gbasis0 check = not $ any (involvesVariable g) [1 .. n] poly = if check then Just (constantTerm +> dropVariables n g) else Nothing -- resultant ------------------------------------------------------------------ -- | sylvester matrix sylvesterMatrix :: AlgAdd.C a => [a] -> [a] -> Matrix a sylvesterMatrix x y = fromLists (xrows ++ yrows) where m = length x - 1 n = length y - 1 xrows = [replicate i AlgAdd.zero ++ x ++ replicate (n-i-1) AlgAdd.zero | i <- [0 .. n-1]] yrows = [replicate i AlgAdd.zero ++ y ++ replicate (m-i-1) AlgAdd.zero | i <- [0 .. m-1]] -- | "truncated" Sylvester matrix sylvesterMatrix' :: AlgRing.C a => [a] -> [a] -> Int -> Matrix a sylvesterMatrix' x y k = if s == 0 then fromLists [[AlgRing.one]] -- plays the role of the empty matrix: -- the point to get is determinant=1 -- (because the empty matrix is not allowed -- in the matrix package) else submatrix 1 s 1 s $ fromLists (xrows ++ yrows) where m = length x - 1 n = length y - 1 s = m + n - 2*k xrows = [replicate i AlgAdd.zero ++ x ++ replicate (n-i-1) AlgAdd.zero | i <- [0 .. n-1-k]] yrows = [replicate i AlgAdd.zero ++ y ++ replicate (m-i-1) AlgAdd.zero | i <- [0 .. m-1-k]] -- | the coefficients of a spray as a univariate spray in x_1 with -- spray coefficients sprayCoefficients :: (Eq a, AlgRing.C a) => Spray a -> [Spray a] sprayCoefficients spray = if n == 0 then [constantTerm] else sprays where n = numberOfVariables spray spray'' = safeSpray spray spray' = removeConstantTerm'' spray'' (expnts', coeffs') = unzip (HM.toList spray') constantTerm = (constantSpray . getConstantTerm'') spray'' xpows = map (`index` 0) expnts' powers'' = map (powerize . S.deleteAt 0) expnts' sprays'' = zipWith (curry fromTerm) powers'' coeffs' imap = IM.fromListWith (^+^) (zip xpows sprays'') imap' = IM.insertWith (^+^) 0 constantTerm imap permutation = [2 .. n] ++ [1] deg = maximum xpows sprays = [ permuteVariables permutation (fromMaybe zeroSpray (IM.lookup i imap')) | i <- [deg, deg-1 .. 0] ] -- | Resultant of two /univariate/ sprays resultant1 :: (Eq a, AlgRing.C a) => Spray a -> Spray a -> a resultant1 p q = if n <= 1 then detLaplace $ sylvesterMatrix pcoeffs qcoeffs else error "resultant1: the two sprays must be univariate." where n = max (numberOfVariables p) (numberOfVariables q) p'' = safeSpray p q'' = safeSpray q pexpnts = map (`index` 0) $ HM.keys $ removeConstantTerm'' p'' qexpnts = map (`index` 0) $ HM.keys $ removeConstantTerm'' q'' p0 = getConstantTerm'' p'' q0 = getConstantTerm'' q'' pcoeffs = if null pexpnts then [p0] else [getCoefficient'' (S.singleton i) p'' | i <- [maxp, maxp-1 .. 1]] ++ [p0] where maxp = maximum pexpnts qcoeffs = if null qexpnts then [q0] else [getCoefficient'' (S.singleton i) q'' | i <- [maxq, maxq-1 .. 1]] ++ [q0] where maxq = maximum qexpnts -- | Subresultants of two /univariate/ sprays. This function makes several calls -- to `detLaplace` and then it can be slow. subresultants1 :: (Eq a, AlgRing.C a) => Spray a -> Spray a -> [a] subresultants1 p q = if n <= 1 then map (detLaplace . sylvesterMatrix' pcoeffs qcoeffs) [0 .. min d e - 1] else error "subresultants1: the two sprays must be univariate." where n = max (numberOfVariables p) (numberOfVariables q) p'' = safeSpray p q'' = safeSpray q pexpnts = map (`index` 0) $ HM.keys $ removeConstantTerm'' p'' qexpnts = map (`index` 0) $ HM.keys $ removeConstantTerm'' q'' p0 = getConstantTerm'' p'' q0 = getConstantTerm'' q'' pcoeffs = if null pexpnts then [p0] else [getCoefficient'' (S.singleton i) p'' | i <- [maxp, maxp-1 .. 1]] ++ [p0] where maxp = maximum pexpnts qcoeffs = if null qexpnts then [q0] else [getCoefficient'' (S.singleton i) q''| i <- [maxq, maxq-1 .. 1]] ++ [q0] where maxq = maximum qexpnts d = length pcoeffs e = length qcoeffs -- | Resultant of two sprays resultant :: (Eq a, AlgRing.C a) => Int -- ^ indicator of the variable with respect to which the resultant is desired (e.g. 1 for x) -> Spray a -> Spray a -> Spray a resultant var p q = if var >= 1 && var <= n then permuteVariables permutation' det else error "resultant: invalid variable index." where n = max (numberOfVariables p) (numberOfVariables q) permutation = [n-var+2 .. n] ++ [1 .. n-var+1] permutation' = [var .. n] ++ [1 .. var-1] p' = permuteVariables permutation p q' = permuteVariables permutation q det = detLaplace $ sylvesterMatrix (sprayCoefficients p') (sprayCoefficients q') -- | Subresultants of two sprays (the /principal/ subresultants, while the -- `polynomialSubresultants` function returns the polynomial subresultants). -- This function makes several calls to `detLaplace` and then it can be slow. subresultants :: (Eq a, AlgRing.C a) => Int -- ^ indicator of the variable with respect to which the subresultants are desired (e.g. 1 for x) -> Spray a -> Spray a -> [Spray a] subresultants var p q | var < 1 = error "subresultants: invalid variable index." | var > n = error "subresultants: too large variable index." | otherwise = map (permute' . detLaplace . sylvesterMatrix' pcoeffs qcoeffs) [0 .. min d e - 1] where pcoeffs = sprayCoefficients p' qcoeffs = sprayCoefficients q' d = length pcoeffs e = length qcoeffs n = max (numberOfVariables p) (numberOfVariables q) permutation = var : [1 .. var-1] ++ [var+1 .. n] permute = permuteVariables permutation p' = permute p q' = permute q permutation' = [2 .. var] ++ (1 : [var+1 .. n]) permute' = permuteVariables permutation' -- | Resultant of two sprays with coefficients in a field; this function is more -- efficient than the function `resultant` resultant' :: forall a. (Eq a, AlgField.C a) => Int -- ^ indicator of the variable with respect to which the resultant is desired (e.g. 1 for x) -> Spray a -> Spray a -> Spray a resultant' var sprayA sprayB | var < 1 || var > n = error "resultant': invalid variable index." | isZeroSpray sprayA || isZeroSpray sprayB = zeroSpray | otherwise = permuteVariables permutation' $ go unitSpray unitSpray s0 p0 q0 where n = max (numberOfVariables sprayA) (numberOfVariables sprayB) permutation = [n-var+1 .. n] ++ [1 .. n-var] permutation' = [var+1 .. n] ++ [1 .. var] sprayA' = permuteVariables permutation sprayA sprayB' = permuteVariables permutation sprayB degA = degree n sprayA' degB = degree n sprayB' content :: Spray a -> Spray a content spray = foldl1' gcdSpray (sprayCoefficients' n spray) exactDivisionBy :: Spray a -> Spray a -> Spray a exactDivisionBy b a = if isZeroSpray remainder then quo else error "exactDivisionBy: should not happen." where (quo, remainder) = sprayDivision a b contA = content sprayA' contB = content sprayB' sprayA'' = exactDivisionBy contA sprayA' sprayB'' = exactDivisionBy contB sprayB' t = contA^**^degB ^*^ contB^**^degA s0 = if degA < degB && odd degA && odd degB then AlgAdd.negate unitSpray :: Spray a else unitSpray (p0, q0) = if degA >= degB then (sprayA'', sprayB'') else (sprayB'', sprayA'') go :: Spray a -> Spray a -> Spray a -> Spray a -> Spray a -> Spray a go g h s p q = if degq' == 0 then s' ^*^ t ^*^ h'' else go g' h' s' p' q' where degp = degree n p degq = degree n q delta = degp - degq s' = if odd degp && odd degq then AlgAdd.negate s else s (_, (_, r)) = pseudoDivision' n p q p' = q q' = exactDivisionBy (g ^*^ h^**^delta) r (degp', ellp') = degreeAndLeadingCoefficient n p' (degq', ellq') = degreeAndLeadingCoefficient n q' g' = ellp' h' = exactDivisionBy (h^**^delta) (h ^*^ g'^**^delta) h'' = exactDivisionBy (h'^**^degp') (h' ^*^ ellq'^**^degp') -- | Polynomial subresultants of two sprays (the `subresultants` function -- computes the /principal/ subresultants). -- This function makes several calls to `detLaplace` and then it can be slow. polynomialSubresultants :: (Eq a, AlgRing.C a) => Int -- ^ index of the variable with respect to which the subresultants will be computed (e.g. 2 for @y@) -> Spray a -> Spray a -> [Spray a] polynomialSubresultants var p q | var < 1 || var > d = error "polynomialSubresultants: invalid variable index." | otherwise = [sres i | i <- [0 .. min n m]] where d = max (numberOfVariables p) (numberOfVariables q) permutation = [d-var+1 .. d] ++ [1 .. d-var] permutation' = [var+1 .. d] ++ [1 .. var] p' = permuteVariables permutation p q' = permuteVariables permutation q pcoeffs = reverse $ sprayCoefficients' d p' qcoeffs = reverse $ sprayCoefficients' d q' pcoeff k = if k < 0 then zeroSpray else pcoeffs !! k qcoeff k = if k < 0 then zeroSpray else qcoeffs !! k n = length pcoeffs - 1 m = length qcoeffs - 1 x = loneTerm' d prow i k = replicate k zeroSpray ++ [pcoeff j | j <- [n, n-1 .. 2*i - m + k + 2]] ++ [multSprayByTerm p' (x (m-i-1-k))] prows i = [prow i k | k <- [0 .. m - i - 1]] qrow i k = replicate k zeroSpray ++ [qcoeff j | j <- [m, m-1 .. 2*i - n + k + 2]] ++ [multSprayByTerm q' (x (n-i-1-k))] qrows i = [qrow i k | k <- [0 .. n - i - 1]] matrix i = fromLists (prows i ++ qrows i) sres i = if n == m && i == n then q else permuteVariables permutation' $ detLaplace (matrix i) -- | Sturm-Habicht sequence of a spray. This function calls `polynomialSubresultants` -- and then it can be slow. sturmHabichtSequence :: (Eq a, AlgRing.C a) => Int -- ^ index of the variable with respect to which the Sturm-Habicht sequence will be computed (e.g. 2 for @y@) -> Spray a -> [Spray a] sturmHabichtSequence var spray | var < 1 || var > d = error "sturmHabichtSequence: invalid variable index." | otherwise = signedSubresultants ++ [spray', spray] where d = numberOfVariables spray spray' = derivative var spray spraySubresultants = polynomialSubresultants var spray spray' n = length spraySubresultants delta = [odd $ ((n - k - 1)*(n - k)) `div` 2 | k <- [0 .. n - 2]] negateIf test = if test then AlgAdd.negate else id signedSubresultants = zipWith negateIf delta spraySubresultants -- | Principal Sturm-Habicht sequence of a spray. This function calls `sturmHabicht` -- sequence and then it can be slow. principalSturmHabichtSequence :: (Eq a, AlgRing.C a) => Int -- ^ index of the variable with respect to which the Sturm-Habicht sequence will be computed (e.g. 2 for @y@) -> Spray a -> [Spray a] principalSturmHabichtSequence var spray | var < 1 || var > d = error "principalSturmHabichtSequence: invalid variable index." | otherwise = [jcoeff j | j <- [0 .. length sHS - 1]] where d = numberOfVariables spray sHS = sturmHabichtSequence var spray permutation = [d - var + 1 .. d] ++ [1 .. d - var] permutation' = [var + 1 .. d] ++ [1 .. var] jcoeff j = if isZeroSpray sHSj || j >= length coeffs then zeroSpray else permuteVariables permutation' (coeffs !! j) where sHSj = sHS !! j p = permuteVariables permutation sHSj coeffs = reverse $ sprayCoefficients' d p -- what if numberOfVariables p /= d ? -- Number of real roots in interval ------------------------------------------- _signVariations :: (a -> Char) -> [a] -> Int _signVariations signFunc as = v1 + v2 + 2*v3 where count x xs = sum (map (fromEnum . (== x)) xs) l = length as signs = map signFunc as chunks2 = [(signs !! i, signs !! (i+1)) | i <- [0 .. l-2]] v1 = count ('+', '-') chunks2 + count ('-', '+') chunks2 chunks3 = [(signs !! i, signs !! (i+1), signs !! (i+2)) | i <- [0 .. l-3]] chunks4 = [(signs !! i, signs !! (i+1), signs !! (i+2), signs !! (i+3)) | i <- [0 .. l-4]] v2 = count ('-', '0', '+') chunks3 + count ('+', '0', '-') chunks3 + count ('+', '0', '0', '-') chunks4 + count ('-', '0', '0', '+') chunks4 v3 = count ('+', '0', '0', '+') chunks4 + count ('-', '0', '0', '-') chunks4 signVariations :: (Eq a, Num a) => [a] -> Int signVariations = _signVariations signFunc where signFunc a | signum a == 0 = '0' | signum a == 1 = '+' | otherwise = '-' signVariations' :: (Eq a, AlgAbs.C a) => [a] -> Int signVariations' = _signVariations signFunc where signFunc a | AlgAbs.signum a == AlgAdd.zero = '0' | AlgAbs.signum a == AlgRing.one = '+' | otherwise = '-' _numberOfRealRootsInOpenInterval :: (AlgRing.C a, Ord a) => ([a] -> Int) -> Spray a -> Maybe a -> Maybe a -> Int _numberOfRealRootsInOpenInterval signVariationsFunc spray alpha beta | isConstantSpray spray = if isZeroSpray spray then error "numberOfRealRootsInOpenInterval: the spray is null." else 0 | alphaIsJust && betaIsJust = if alpha' == beta' then 0 else if alpha > beta then error "numberOfRealRootsInOpenInterval: the bounds are not ordered." else if isZeroAtBeta then svDiff - 1 else svDiff | alphaIsJust = svDiff | betaIsJust = if isZeroAtBeta then svDiff - 1 else svDiff | otherwise = svDiff where alphaIsJust = isJust alpha betaIsJust = isJust beta alpha' = if alphaIsJust then fromJust alpha else undefined beta' = if betaIsJust then fromJust beta else undefined evaluateAtMinusInfinity p = if isConstantSpray p then getConstantTerm p else let (deg, lCoeff) = first ((`S.index` 0) . exponents) (leadingTerm p) in if even deg then lCoeff else AlgAdd.negate lCoeff evaluateAtPlusInfinity p = if isConstantSpray p then getConstantTerm p else (snd . leadingTerm) p (ginit, glast) = fromJust $ unsnoc $ filter (not . isZeroSpray) (sturmHabichtSequence 1 spray) eval1 = if alphaIsJust then evaluateAt [alpha'] else evaluateAtMinusInfinity eval2 = if betaIsJust then evaluateAt [beta'] else evaluateAtPlusInfinity sprayAtAlpha = eval1 glast sprayAtBeta = eval2 glast isZeroAtBeta = sprayAtBeta == AlgAdd.zero galpha = map eval1 ginit ++ [sprayAtAlpha] gbeta = map eval2 ginit ++ [sprayAtBeta] svAtAlpha = signVariationsFunc galpha svAtBeta = signVariationsFunc gbeta svDiff = svAtAlpha - svAtBeta _numberOfRealRootsInClosedInterval :: (AlgRing.C a, Ord a) => ([a] -> Int) -> Spray a -> Maybe a -> Maybe a -> Int _numberOfRealRootsInClosedInterval signVariationsFunc spray alpha beta = _numberOfRealRootsInOpenInterval signVariationsFunc spray alpha beta + toAdd where alphaIsJust = isJust alpha betaIsJust = isJust beta alpha' = if alphaIsJust then fromJust alpha else undefined beta' = if betaIsJust then fromJust beta else undefined isZeroAtAlpha = alphaIsJust && evaluateAt [alpha'] spray == AlgAdd.zero isZeroAtBeta = betaIsJust && evaluateAt [beta'] spray == AlgAdd.zero toAdd = if alpha == beta then fromEnum isZeroAtAlpha else fromEnum isZeroAtAlpha + fromEnum isZeroAtBeta -- | Number of real roots of a spray in an open interval (that makes sense -- only for a spray on a ring embeddable in the real numbers). numberOfRealRootsInOpenInterval :: (Num a, AlgRing.C a, Ord a) => Spray a -- ^ a spray -> Maybe a -- ^ lower bound of the interval; use @Just@ for a finite bound, and @Nothing@ for minus infinity -> Maybe a -- ^ upper bound of the interval; use @Just@ for a finite bound, and @Nothing@ for minus infinity -> Int numberOfRealRootsInOpenInterval spray = if isUnivariate spray then _numberOfRealRootsInOpenInterval signVariations spray else error "numberOfRealRootsInOpenInterval: the spray is not univariate." -- | Number of real roots of a spray in a closed interval (that makes sense -- only for a spray on a ring embeddable in the real numbers). The roots are -- not counted with their multiplicity. numberOfRealRootsInClosedInterval :: (Num a, AlgRing.C a, Ord a) => Spray a -- ^ a spray -> Maybe a -- ^ lower bound of the interval; use @Just@ for a finite bound, and @Nothing@ for minus infinity -> Maybe a -- ^ upper bound of the interval; use @Just@ for a finite bound, and @Nothing@ for minus infinity -> Int numberOfRealRootsInClosedInterval spray = if isUnivariate spray then _numberOfRealRootsInClosedInterval signVariations spray else error "numberOfRealRootsInClosedInterval: the spray is not univariate." -- | Number of real roots of a spray in an open interval (that makes sense -- only for a spray on a ring embeddable in the real numbers). numberOfRealRootsInOpenInterval' :: (AlgAbs.C a, Ord a) => Spray a -- ^ a spray -> Maybe a -- ^ lower bound of the interval; use @Just@ for a finite bound, and @Nothing@ for minus infinity -> Maybe a -- ^ upper bound of the interval; use @Just@ for a finite bound, and @Nothing@ for minus infinity -> Int numberOfRealRootsInOpenInterval' spray = if isUnivariate spray then _numberOfRealRootsInOpenInterval signVariations' spray else error "numberOfRealRootsInOpenInterval': the spray is not univariate." -- | Number of real roots of a spray in a closed interval (that makes sense -- only for a spray on a ring embeddable in the real numbers). The roots are -- not counted with their multiplicity. numberOfRealRootsInClosedInterval' :: (AlgAbs.C a, Ord a) => Spray a -- ^ a spray -> Maybe a -- ^ lower bound of the interval; use @Just@ for a finite bound, and @Nothing@ for minus infinity -> Maybe a -- ^ upper bound of the interval; use @Just@ for a finite bound, and @Nothing@ for minus infinity -> Int numberOfRealRootsInClosedInterval' spray = if isUnivariate spray then _numberOfRealRootsInClosedInterval signVariations' spray else error "numberOfRealRootsInClosedInterval': the spray is not univariate." -- Number of real roots ------------------------------------------------------- runLengthEncoding :: Eq a => [a] -> [(a,Int)] runLengthEncoding = foldr code [] where code c [] = [(c,1)] code c ((x,n):ts) | c == x = (x,n+1):ts | otherwise = (c,1):(x,n):ts _signPermanencesAndVariations :: (a -> Char) -> [a] -> (Int, Int) _signPermanencesAndVariations signFunc as = (permanences, variations) where signs = map signFunc as rle = runLengthEncoding signs l = length rle lengths = map snd rle permanences = sum lengths - l variations = l - 1 signPermanencesAndVariations :: (Eq a, Num a) => [a] -> (Int, Int) signPermanencesAndVariations = _signPermanencesAndVariations signFunc where signFunc a = if signum a == 1 then '+' else '-' signPermanencesAndVariations' :: (Eq a, AlgAbs.C a) => [a] -> (Int, Int) signPermanencesAndVariations' = _signPermanencesAndVariations signFunc where signFunc a = if AlgAbs.signum a == AlgRing.one then '+' else '-' distributionsOfZeros :: (Eq a, AlgAdd.C a) => [a] -> ([Int], [Int], [Int]) distributionsOfZeros as = (i_, k_, ik_) where symbol a = if a == AlgAdd.zero then '0' else '*' symbols = map symbol as ++ ['0'] lengths = map snd (runLengthEncoding symbols) l = length lengths cumulativeLengths = scanl (+) 0 lengths range = [0 .. l-1] ik_ = [cumulativeLengths !! i | i <- range, even i] i_ = [cumulativeLengths !! i - 1 | i <- range, odd i] k_ = [lengths !! i | i <- [0 .. l-2], odd i] ++ [lengths !! (l-1) - 1] _blocksAndEpsilons :: (Eq a, AlgAdd.C a) => (a -> Int) -> [a] -> ([[a]], [Int]) _blocksAndEpsilons signFunc as = (blocks, epsilons) where (i_, k_, ik_) = distributionsOfZeros as t = length i_ blocks = [[as !! m | m <- [ik_ !! n .. i_ !! n]] | n <- [0 .. t-1]] epsilon s = let ks = k_ !! s in if odd ks then 0 else (if even (ks `div` 2) then 1 else -1) * signFunc (as !! (ik_ !! (s+1))) * signFunc (as !! (i_ !! s)) epsilons = [epsilon n | n <- [0 .. t-2]] blocksAndEpsilons :: (Eq a, AlgAdd.C a, Num a) => [a] -> ([[a]], [Int]) blocksAndEpsilons = _blocksAndEpsilons signFunc where signFunc a = if signum a == 1 then 1 else -1 blocksAndEpsilons' :: (Eq a, AlgAbs.C a) => [a] -> ([[a]], [Int]) blocksAndEpsilons' = _blocksAndEpsilons signFunc where signFunc a = if AlgAbs.signum a == AlgRing.one then 1 else -1 _C :: ([a] -> (Int, Int)) -> ([[a]], [Int]) -> Int _C signPermanencesAndVariationsFunc (blocks, epsilons) = sum pvs - sum epsilons where permanencesMinusVariations as = let (p, v) = signPermanencesAndVariationsFunc as in p - v pvs = map permanencesMinusVariations blocks -- | Number of real roots of a spray (that makes sense only for a spray on a -- ring embeddable in the real numbers). The roots are not counted with their -- multiplicity. numberOfRealRoots :: (Eq a, AlgRing.C a, Num a) => Spray a -> Int numberOfRealRoots spray = if isUnivariate spray then if isConstantSpray spray then if isZeroSpray spray then error "numberOfRealRoots: the spray is null." else 0 else _C signPermanencesAndVariations (blocksAndEpsilons $ reverse as) else error "numberOfRealRoots: the spray is not univariate." where as = map getConstantTerm (principalSturmHabichtSequence 1 spray) -- | Number of real roots of a spray (that makes sense only for a spray on a -- ring embeddable in the real numbers). The roots are not counted with their -- multiplicity. numberOfRealRoots' :: (Eq a, AlgAbs.C a) => Spray a -> Int numberOfRealRoots' spray = if isUnivariate spray then if isConstantSpray spray then if isZeroSpray spray then error "numberOfRealRoots: the spray is null." else 0 else _C signPermanencesAndVariations' (blocksAndEpsilons' $ reverse as) else error "numberOfRealRoots': the spray is not univariate." where as = map getConstantTerm (principalSturmHabichtSequence 1 spray) -- GCD stuff ------------------------------------------------------------------ -- | the spray coefficients of a spray as a univariate spray in x_n sprayCoefficients' :: (Eq a, AlgRing.C a) => Int -> Spray a -> [Spray a] sprayCoefficients' n spray | numberOfVariables spray /= n = [spray] | n == 0 = [constantSpray constantTerm] | otherwise = sprays where permutation = [2 .. n] ++ [1] spray' = safeSpray $ permuteVariables permutation spray spray'' = removeConstantTerm'' spray' (expnts', coeffs') = unzip (HM.toList spray'') constantTerm = getConstantTerm'' spray' xpows = map (`index` 0) expnts' powers'' = map (powerize . S.deleteAt 0) expnts' sprays'' = zipWith (curry fromTerm) powers'' coeffs' imap = IM.fromListWith (^+^) (zip xpows sprays'') imap' = IM.insertWith (^+^) 0 (constantSpray constantTerm) imap deg = maximum xpows sprays = [ fromMaybe zeroSpray (IM.lookup i imap') | i <- [deg, deg-1 .. 0] ] -- | the degree of a spray as a univariate spray in x_n with spray coefficients degree :: (Eq a, AlgRing.C a) => Int -> Spray a -> Int degree n spray | isConstant spray = if isZeroSpray spray then minBound -- (should not happen) else 0 | numberOfVariables spray /= n = 0 | otherwise = maximum xpows where permutation = [2 .. n] ++ [1] spray' = permuteVariables permutation spray powers' = HM.keys $ removeConstantTerm spray' xpows = map ((`index` 0) . exponents) powers' -- | the degree and the leading coefficient of a spray as a univariate spray -- in x_n with spray coefficients degreeAndLeadingCoefficient :: (Eq a, AlgRing.C a) => Int -> Spray a -> (Int, Spray a) degreeAndLeadingCoefficient n spray | n == 0 = ( if constantTerm == AlgAdd.zero then minBound -- (should not happen) else 0, constantSpray constantTerm ) | numberOfVariables spray /= n = (0, spray) | otherwise = (deg, leadingCoeff) where permutation = [2 .. n] ++ [1] spray' = permuteVariables permutation spray constantTerm = getConstantTerm spray' spray'' = removeConstantTerm spray' (powers', coeffs') = unzip (HM.toList spray'') expnts' = map exponents powers' xpows = map (`index` 0) expnts' deg = maximum xpows is = elemIndices deg xpows powers'' = [powerize (S.deleteAt 0 (expnts' !! i)) | i <- is] coeffs'' = [coeffs' !! i | i <- is] leadingCoeff = sumTerms (zip powers'' coeffs'') -- | Pseudo-division of two sprays @A@ and @B@ such that @deg(A) >= deg(B)@ -- where @deg@ is the degree with respect to the outermost variable. pseudoDivision :: (Eq a, AlgRing.C a) => Spray a -- ^ A -> Spray a -- ^ B -> (Spray a, (Spray a, Spray a)) -- ^ (C, (Q, R)) such that C^*^A = B^*^Q ^+^ R pseudoDivision sprayA sprayB = pseudoDivision' (max (numberOfVariables sprayA) (numberOfVariables sprayB)) sprayA sprayB -- | pseudo-division of two sprays, assuming degA >= degB >= 0 pseudoDivision' :: (Eq a, AlgRing.C a) => Int -- ^ number of variables -> Spray a -- ^ A -> Spray a -- ^ B -> (Spray a, (Spray a, Spray a)) -- ^ (C, (Q, R)) such that C^*^A = B^*^Q ^+^ R pseudoDivision' n sprayA sprayB | degB == minBound = error "pseudoDivision: pseudo-division by 0." | degA < degB = error "pseudoDivision: degree(A) < degree(B)." | otherwise = (ellB ^**^ delta , go sprayA zeroSpray delta) where degA = degree n sprayA (degB, ellB) = degreeAndLeadingCoefficient n sprayB delta = degA - degB + 1 go sprayR sprayQ e = if isZeroSpray sprayR || degR < degB then (q ^*^ sprayQ, q ^*^ sprayR) else go (ellB ^*^ sprayR ^-^ sprayS ^*^ sprayB) (ellB ^*^ sprayQ ^+^ sprayS) (e - 1) where (degR, ellR) = degreeAndLeadingCoefficient n sprayR q = ellB ^**^ e sprayS = multSprayByTerm ellR (loneTerm' n (degR - degB)) {- spray1,spray2 :: QSpray spray1 = let x = qlone 1 y = qlone 2 z = qlone 3 in (-5)*^x^**^3 ^-^ 11*^(x^**^2^*^y) ^-^ 11*^(x^*^y^**^2) ^-^ 8*^(x^*^y^*^z) ^+^ 20*^(x^*^z^**^2) ^-^ 5*^y^**^3 ^+^ 20*^(y^*^z^**^2) spray2 = let x = qlone 1 y = qlone 2 z = qlone 3 in x^*^y ^-^ x^*^z ^-^ y^*^z ^+^ z^**^2 -} -- | recursive GCD function gcdKX1dotsXn :: forall a. (Eq a, AlgField.C a) => Int -> Spray a -> Spray a -> Spray a gcdKX1dotsXn n sprayA sprayB | n == 0 = constantSpray $ gcdKX0 sprayA sprayB | degB > degA = gcdKX1dotsXn n sprayB sprayA | isZeroSpray sprayB = sprayA | otherwise = go sprayA' sprayB' unitSpray unitSpray where gcdKX0 :: Spray a -> Spray a -> a gcdKX0 = const $ const AlgRing.one n' = max (numberOfVariables sprayA) (numberOfVariables sprayB) degA = degree n' sprayA degB = degree n' sprayB gcdKX1dotsXm = gcdKX1dotsXn (n-1) content :: Spray a -> Spray a content spray = foldl1' gcdKX1dotsXm (sprayCoefficients' n' spray) exactDivisionBy :: Spray a -> Spray a -> Spray a exactDivisionBy b a = if isZeroSpray remainder then quo else error "exactDivisionBy: should not happen." where (quo, remainder) = sprayDivision a b reduceSpray :: Spray a -> Spray a reduceSpray spray = exactDivisionBy (content spray) spray contA = content sprayA contB = content sprayB d = gcdKX1dotsXm contA contB sprayA' = exactDivisionBy contA sprayA sprayB' = exactDivisionBy contB sprayB go :: Spray a -> Spray a -> Spray a -> Spray a -> Spray a go sprayA'' sprayB'' g h | isZeroSpray sprayR = d ^*^ reduceSpray sprayB'' | numberOfVariables sprayR == 0 = d | otherwise = go sprayB'' (exactDivisionBy (g ^*^ h^**^delta) sprayR) ellB'' ---- ellA'' (exactDivisionBy (h^**^delta) (h ^*^ ellB''^**^delta)) --- g^**^delta where (_, (_, sprayR)) = pseudoDivision' n' sprayA'' sprayB'' -- (degA'', ellA'') = degreeAndLeadingCoefficient n' sprayA'' -- degB'' = degree n' sprayB'' degA'' = degree n' sprayA'' (degB'', ellB'') = degreeAndLeadingCoefficient n' sprayB'' delta = degA'' - degB'' -- | Greatest common divisor of two sprays with coefficients in a field gcdSpray :: forall a. (Eq a, AlgField.C a) => Spray a -> Spray a -> Spray a gcdSpray sprayA sprayB = gcdKX1dotsXn n sprayA sprayB where n = max (numberOfVariables sprayA) (numberOfVariables sprayB) -- Matrices ------------------------------------------------------------------- -- | Determinant of a matrix with entries in a ring by using Laplace -- expansion (this is slow); the __numeric-prelude__ package provides some -- stuff to deal with matrices over a ring but it does not provide the -- determinant detLaplace :: forall a. (Eq a, AlgRing.C a) => Matrix a -> a detLaplace b = if nrows b == ncols b then detUnsafe b else error "detLaplace: the matrix is not square." where detUnsafe m = if nrows m == 1 then m DM.! (1,1) else suml1 [negateIf i (times (m DM.! (i,1)) (detUnsafe (minorMatrix i 1 m))) | i <- [1 .. nrows m]] suml1 = foldl1' (AlgAdd.+) negateIf i = if even i then AlgAdd.negate else id times :: a -> a -> a times x y = if x == AlgAdd.zero then AlgAdd.zero else x AlgRing.* y -- | Determinant of a matrix over a ring by using Laplace expansion; this is -- the same as `detLaplace` but for a matrix from the __numeric-prelude__ -- package detLaplace' :: forall a. (Eq a, AlgRing.C a) => MathMatrix.T a -> a detLaplace' m = detLaplace (DM.fromLists $ MathMatrix.rows m) -- | Characteristic polynomial of a square matrix -- -- >>> import Data.Matrix (Matrix, fromLists) -- >>> m = fromLists [ [12, 16, 4] -- >>> , [16, 2, 8] -- >>> , [8, 18, 10] ] :: Matrix Int -- >>> spray = characteristicPolynomial m -- >>> putStrLn $ prettyNumSpray spray -- -x^3 + 24*x^2 + 268*x - 1936 characteristicPolynomial :: (Eq a, AlgRing.C a) => Matrix a -> Spray a characteristicPolynomial m = if nrows m /= ncols m then error "characteristicPolynomial: the matrix is not square." else detLaplace m' where m' = DM.mapPos f m f (i, j) mij = if i == j then constantSpray mij ^-^ x else constantSpray mij x = lone 1 -- Ratios of sprays ----------------------------------------------------------- -- | A @RatioOfSprays a@ object represents a fraction of two multivariate -- polynomials whose coefficients are of type @a@, which represents a field. -- These two polynomials are represented by two @Spray a@ objects. Generally -- we do not use this constructor to build a ratio of sprays: we use the `%//%` -- operator instead, because it always returns an irreducible ratio of sprays, -- meaning that its corresponding fraction of polynomials is irreducible, i.e. -- its numerator and its denominator are coprime. You can use this constructor -- if you are /sure/ that the numerator and the denominator are coprime. This -- can save some computation time, but unfortunate consequences can occur if -- the numerator and the denominator are not coprime. An arithmetic operation -- on ratios of sprays always returns an irreducible ratio of sprays under the -- condition that the ratios of sprays it involves are irreducible. Moreover, -- it never returns a ratio of sprays with a constant denominator other than -- the unit spray. If you use this constructor with a constant denominator, -- always set this denominator to the unit spray (by dividing the numerator -- by the constant value of the denominator). data RatioOfSprays a = RatioOfSprays { _numerator :: Spray a , _denominator :: Spray a } deriving Show type RatioOfQSprays = RatioOfSprays Rational instance (Eq a, AlgField.C a) => FunctionLike (RatioOfSprays a) where type BaseRing (RatioOfSprays a) = a -- type VariablesType (RatioOfSprays a) = Spray a -- (*^) :: a -> RatioOfSprays a -> RatioOfSprays a (*^) lambda rOS = lambda AlgMod.*> rOS -- (+>) :: a -> RatioOfSprays a -> RatioOfSprays a (+>) x rOS = rOS AlgAdd.+ asRatioOfSprays (constantSpray x) -- substitute :: [Maybe a] -> RatioOfSprays a -> RatioOfSprays a substitute subs (RatioOfSprays p q) = substitute subs p %//% substitute subs q -- evaluate :: RatioOfSprays a -> [a] -> a evaluate (RatioOfSprays p q) xyz = evaluate p xyz AlgField./ evaluate q xyz -- changeVariables :: RatioOfSprays a -> [Spray a] -> RatioOfSprays a changeVariables rOS newVariables = if length newVariables < numberOfVariables rOS then error "changeVariables: not enough new variables provided." else changeVariables (_numerator rOS) newVariables %//% changeVariables (_denominator rOS) newVariables -- numberOfVariables :: RatioOfSprays a -> Int numberOfVariables (RatioOfSprays p q) = max (numberOfVariables p) (numberOfVariables q) -- permuteVariables :: [Int] -> RatioOfSprays a -> RatioOfSprays a permuteVariables permutation (RatioOfSprays p q) = permuteVariables permutation p %//% permuteVariables permutation q -- swapVariables :: (Int, Int) -> RatioOfSprays a -> RatioOfSprays a swapVariables (i, j) (RatioOfSprays p q) = swapVariables (i, j) p %//% swapVariables (i, j) q -- involvesVariable :: RatioOfSprays a -> Int -> Bool involvesVariable (RatioOfSprays p q) i = involvesVariable p i || involvesVariable q i -- dropVariables :: Int -> RatioOfSprays a -> RatioOfSprays a dropVariables n (RatioOfSprays p q) = dropVariables n p %//% dropVariables n q -- derivative :: Int -> RatioOfSprays a -> RatioOfSprays a derivative i (RatioOfSprays p q) = (p' ^*^ q ^-^ p ^*^ q') %//% (q ^*^ q) where p' = derivative i p q' = derivative i q -- | quotients of two univariate sprays by their gcd -- we use `sprayDivision0` because this function is called -- (by `irreducibleFraction`) with non-constant sprays only quotientsByGCD :: (Eq a, AlgField.C a) => Spray a -> Spray a -> (Spray a, Spray a) quotientsByGCD sprayA sprayB = if isUnivariate sprayA && isUnivariate sprayB then go sprayA sprayB unitSpray zeroSpray zeroSpray unitSpray else (exactDivision sprayA g, exactDivision sprayB g) where exactDivision p q = fst (sprayDivision0 p q) g = gcdSpray sprayA sprayB go oldr r olds s oldt t | isZeroSpray r = (c *^ AlgAdd.negate t, c *^ s) -- monic denominator | otherwise = go r remainder s (olds ^-^ quo ^*^ s) t (oldt ^-^ quo ^*^ t) where (quo, remainder) = univariateSprayDivision oldr r c = AlgField.recip (snd $ leadingTerm s) -- | irreducible fraction of sprays irreducibleFraction :: (Eq a, AlgField.C a) => Spray a -> Spray a -> RatioOfSprays a irreducibleFraction p q = adjustFraction rOS where rOS = if isConstant p || isConstant q then RatioOfSprays p q else let (a, b) = quotientsByGCD p q in RatioOfSprays a b -- | set denominator to 1 if it is constant adjustFraction :: (Eq a, AlgField.C a) => RatioOfSprays a -> RatioOfSprays a adjustFraction (RatioOfSprays p q) = if isConstant q then RatioOfSprays (p /^ getConstantTerm q) unitSpray else RatioOfSprays p q instance (AlgRing.C a, Eq a) => Eq (RatioOfSprays a) where (==) :: RatioOfSprays a -> RatioOfSprays a -> Bool (==) (RatioOfSprays p q) (RatioOfSprays p' q') = isZeroSpray (p ^*^ q' ^-^ p' ^*^ q) instance (AlgField.C a, Eq a) => AlgAdd.C (RatioOfSprays a) where (+) :: RatioOfSprays a -> RatioOfSprays a -> RatioOfSprays a (+) (RatioOfSprays p q) (RatioOfSprays p' q') = irreducibleFraction (p ^*^ q' ^+^ p' ^*^ q) (q ^*^ q') zero :: RatioOfSprays a zero = RatioOfSprays zeroSpray unitSpray negate :: RatioOfSprays a -> RatioOfSprays a negate (RatioOfSprays p q) = RatioOfSprays (negateSpray p) q instance (AlgField.C a, Eq a) => AlgMod.C a (RatioOfSprays a) where (*>) :: a -> RatioOfSprays a -> RatioOfSprays a lambda *> (RatioOfSprays p q) = RatioOfSprays (lambda *^ p) q instance (AlgField.C a, Eq a) => AlgRightMod.C a (RatioOfSprays a) where (<*) :: RatioOfSprays a -> a -> RatioOfSprays a rOS <* lambda = lambda AlgMod.*> rOS instance (AlgField.C a, Eq a) => AlgMod.C (Spray a) (RatioOfSprays a) where (*>) :: Spray a -> RatioOfSprays a -> RatioOfSprays a spray *> (RatioOfSprays p q) = irreducibleFraction (spray ^*^ p) q instance (AlgField.C a, Eq a) => AlgRightMod.C (Spray a) (RatioOfSprays a) where (<*) :: RatioOfSprays a -> Spray a -> RatioOfSprays a rOS <* spray = spray AlgMod.*> rOS instance (AlgField.C a, Eq a) => AlgRing.C (RatioOfSprays a) where (*) :: RatioOfSprays a -> RatioOfSprays a -> RatioOfSprays a (*) (RatioOfSprays p q) (RatioOfSprays p' q') = irreducibleFraction (p ^*^ p') (q ^*^ q') (^) :: RatioOfSprays a -> Integer -> RatioOfSprays a (^) (RatioOfSprays p q) n = if n >= 0 then RatioOfSprays (p AlgRing.^ n) (q AlgRing.^ n) else RatioOfSprays (q AlgRing.^ (-n)) (p AlgRing.^ (-n)) one :: RatioOfSprays a one = RatioOfSprays unitSpray unitSpray instance (AlgField.C a, Eq a) => AlgField.C (RatioOfSprays a) where recip :: RatioOfSprays a -> RatioOfSprays a recip (RatioOfSprays p q) = RatioOfSprays q p infixl 7 %:% -- | Ratio of sprays from numerator and denominator, -- __without reducing the fraction__ (%:%) :: Spray a -> Spray a -> RatioOfSprays a (%:%) = RatioOfSprays infixl 7 %//% -- | Irreducible ratio of sprays from numerator and denominator; alias of @(^/^)@ (%//%) :: (Eq a, AlgField.C a) => Spray a -> Spray a -> RatioOfSprays a (%//%) = irreducibleFraction infixl 7 ^/^ -- | Irreducible ratio of sprays from numerator and denominator; alias of @(%//%)@ (^/^) :: (Eq a, AlgField.C a) => Spray a -> Spray a -> RatioOfSprays a (^/^) = irreducibleFraction infixl 7 %/% -- | Division of a ratio of sprays by a spray; the result is an -- irreducible fraction (%/%) :: (Eq a, AlgField.C a) => RatioOfSprays a -> Spray a -> RatioOfSprays a (%/%) rOS spray = rOS AlgRing.* RatioOfSprays unitSpray spray -- | Whether a ratio of sprays is constant; same as `isConstant` isConstantRatioOfSprays :: (Eq a, AlgField.C a) => RatioOfSprays a -> Bool isConstantRatioOfSprays = isConstant -- | Whether a ratio of sprays actually is polynomial, that is, whether its -- denominator is a constant spray (and then it should be the unit spray) -- -- >>> x = qlone 1 -- >>> y = qlone 2 -- >>> p = x^**^4 ^-^ y^**^4 -- >>> q = x ^-^ y -- >>> isPolynomialRatioOfSprays $ p %//% q -- True -- >>> isPolynomialRatioOfSprays $ p %:% q -- False isPolynomialRatioOfSprays :: (Eq a, AlgRing.C a) => RatioOfSprays a -> Bool isPolynomialRatioOfSprays = isConstant . _denominator -- | The null ratio of sprays zeroRatioOfSprays, zeroROS :: (AlgField.C a, Eq a) => RatioOfSprays a zeroRatioOfSprays = AlgAdd.zero zeroROS = AlgAdd.zero -- | The unit ratio of sprays unitRatioOfSprays, unitROS :: (AlgField.C a, Eq a) => RatioOfSprays a unitRatioOfSprays = AlgRing.one unitROS = AlgRing.one -- | Constant ratio of sprays constantRatioOfSprays :: (Eq a, AlgRing.C a) => a -> RatioOfSprays a constantRatioOfSprays x = asRatioOfSprays (constantSpray x) -- | Evaluates a ratio of sprays; same as `evaluate` evalRatioOfSprays :: (Eq a, AlgField.C a) => RatioOfSprays a -> [a] -> a evalRatioOfSprays = evaluate -- | Substitutes some values to some variables of a ratio of sprays; same as `substitute` substituteRatioOfSprays :: (Eq a, AlgField.C a) => [Maybe a] -> RatioOfSprays a -> RatioOfSprays a substituteRatioOfSprays = substitute -- | Coerces a spray to a ratio of sprays asRatioOfSprays :: AlgRing.C a => Spray a -> RatioOfSprays a asRatioOfSprays spray = RatioOfSprays spray unitSpray -- | Converts a ratio of polynomials to a ratio of sprays fromRatioOfPolynomials :: (Eq a, AlgRing.C a) => RatioOfPolynomials a -> RatioOfSprays a fromRatioOfPolynomials rop = RatioOfSprays (polynomialToSpray $ NumberRatio.numerator rop) (polynomialToSpray $ NumberRatio.denominator rop) -- | Converts a ratio of rational polynomials to a ratio of rational sprays; -- this is not a specialization of `fromRatioOfPolynomials` because -- @RatioOfQPolynomials@ is @RatioOfPolynomials a@ with -- @a = Rational'@, not with @a = Rational@ fromRatioOfQPolynomials :: RatioOfQPolynomials -> RatioOfQSprays fromRatioOfQPolynomials rop = RatioOfSprays (qPolynomialToQSpray $ NumberRatio.numerator rop) (qPolynomialToQSpray $ NumberRatio.denominator rop) -- | General function to print a `RatioOfSprays` object showRatioOfSprays :: (Eq a, AlgRing.C a) => ((Spray a, Spray a) -> (String, String)) -- ^ function which prints a pair of sprays that will be applied to the numerator and the denominator -> (String, String) -- ^ pair of braces to enclose the numerator and the denominator -> String -- ^ represents the quotient bar -> RatioOfSprays a -> String showRatioOfSprays spraysShower braces quotientBar (RatioOfSprays p q) = numeratorString ++ denominatorString where enclose = bracify braces (pString, qString) = spraysShower (p, q) numeratorString = enclose pString denominatorString = if q == unitSpray then "" else quotientBar ++ enclose qString showTwoSpraysXYZ :: (Eq a, AlgRing.C a) => (a -> String) -- ^ function mapping a coefficient to a string, typically 'show' -> (String, String) -- ^ used to enclose the coefficients, usually a pair of braces -> [String] -- ^ typically some letters, to print the variables -> (Spray a, Spray a) -- ^ the two sprays to be printed -> (String, String) showTwoSpraysXYZ showCoef braces letters (spray1, spray2) = both (showSpray showCoef braces showMonomials) (spray1, spray2) where n = max (numberOfVariables spray1) (numberOfVariables spray2) showMonomials = map (unpack . showMonomialXYZ letters n) showTwoSpraysX1X2X3 :: (a -> String) -- ^ function mapping a coefficient to a string, typically 'show' -> (String, String) -- ^ used to enclose the coefficients, usually a pair of braces -> String -- ^ typically a letter, to print the non-indexed variables -> (Spray a, Spray a) -- ^ the two sprays to be printed -> (String, String) showTwoSpraysX1X2X3 showCoef braces letter (spray1, spray2) = both (showSpray showCoef braces showMonomials) (spray1, spray2) where showMonomials = showMonomialsX1X2X3 letter showTwoNumSprays :: (Num a, Ord a) => (a -> String) -- ^ function mapping a positive coefficient to a string -> ([Exponents] -> [String]) -- ^ prints the monomials -> (Spray a, Spray a) -- ^ the two sprays to be printed -> (String, String) showTwoNumSprays showPositiveCoef showMonomials = both (showNumSpray showMonomials showPositiveCoef) showTwoQSprays :: ([Exponents] -> [String]) -- ^ prints the monomials -> (QSpray, QSpray) -- ^ the two sprays to be printed -> (String, String) showTwoQSprays = showTwoNumSprays showRatio showTwoNumSpraysXYZ :: (AlgRing.C a, Num a, Ord a) => (a -> String) -- ^ function mapping a positive coefficient to a string -> [String] -- ^ typically some letters, to print the variables -> (Spray a, Spray a) -- ^ the two sprays to be printed -> (String, String) showTwoNumSpraysXYZ showPositiveCoef letters (spray1, spray2) = showTwoNumSprays showPositiveCoef showMonomials (spray1, spray2) where n = max (numberOfVariables spray1) (numberOfVariables spray2) showMonomials = map (unpack . showMonomialXYZ letters n) showTwoQSpraysXYZ :: [String] -- ^ typically some letters, to print the variables -> (QSpray, QSpray) -- ^ the two sprays to be printed -> (String, String) showTwoQSpraysXYZ = showTwoNumSpraysXYZ showRatio showTwoNumSpraysX1X2X3 :: (Num a, Ord a) => (a -> String) -- ^ function mapping a positive coefficient to a string -> String -- ^ typically a letter, to print the non-indexed variable -> (Spray a, Spray a) -- ^ the two sprays to be printed -> (String, String) showTwoNumSpraysX1X2X3 showPositiveCoef letter (spray1, spray2) = showTwoNumSprays showPositiveCoef showMonomials (spray1, spray2) where showMonomials = showMonomialsX1X2X3 letter showTwoQSpraysX1X2X3 :: String -- ^ typically a letter, to print the non-indexed variables -> (QSpray, QSpray) -- ^ the two sprays to be printed -> (String, String) showTwoQSpraysX1X2X3 = showTwoNumSpraysX1X2X3 showRatio -- | Prints a ratio of sprays with numeric coefficients showRatioOfNumSprays :: (Num a, Ord a, AlgRing.C a) => (a -> String) -- ^ function mapping a positive coefficient to a string -> ([Exponents] -> [String]) -- ^ prints the monomials -> (String, String) -- ^ pair of braces to enclose the numerator and the denominator -> String -- ^ represents the quotient bar -> RatioOfSprays a -> String showRatioOfNumSprays showPositiveCoef showMonomials = showRatioOfSprays (showTwoNumSprays showPositiveCoef showMonomials) -- | Prints a ratio of sprays with rational coefficients showRatioOfQSprays :: ([Exponents] -> [String]) -- ^ prints the monomials -> (String, String) -- ^ pair of braces to enclose the numerator and the denominator -> String -- ^ represents the quotient bar -> RatioOfQSprays -> String showRatioOfQSprays showMonomials = showRatioOfSprays (showTwoQSprays showMonomials) -- | Prints a ratio of sprays with numeric coefficients showRatioOfNumSpraysXYZ :: (Num a, Ord a, AlgRing.C a) => (a -> String) -- ^ function mapping a positive coefficient to a string -> [String] -- ^ typically some letters, to print the variables -> (String, String) -- ^ pair of braces to enclose the numerator and the denominator -> String -- ^ represents the quotient bar -> RatioOfSprays a -> String showRatioOfNumSpraysXYZ showPositiveCoef letters = showRatioOfSprays (showTwoNumSpraysXYZ showPositiveCoef letters) -- | Prints a ratio of sprays with numeric coefficients showRatioOfNumSpraysX1X2X3 :: (Num a, Ord a, AlgRing.C a) => (a -> String) -- ^ function mapping a positive coefficient to a string -> String -- ^ typically a letter, to print the variables -> (String, String) -- ^ pair of braces to enclose the numerator and the denominator -> String -- ^ represents the quotient bar -> RatioOfSprays a -> String showRatioOfNumSpraysX1X2X3 showPositiveCoef letter = showRatioOfSprays (showTwoNumSpraysX1X2X3 showPositiveCoef letter) -- | Prints a ratio of sprays with rational coefficients showRatioOfQSpraysXYZ :: [String] -- ^ typically some letters, to print the variables -> (String, String) -- ^ pair of braces to enclose the numerator and the denominator -> String -- ^ represents the quotient bar -> RatioOfQSprays -> String showRatioOfQSpraysXYZ letters = showRatioOfSprays (showTwoQSpraysXYZ letters) -- | Prints a ratio of sprays with rational coefficients showRatioOfQSpraysX1X2X3 :: String -- ^ typically a letter, to print the variables -> (String, String) -- ^ pair of braces to enclose the numerator and the denominator -> String -- ^ represents the quotient bar -> RatioOfQSprays -> String showRatioOfQSpraysX1X2X3 letter = showRatioOfSprays (showTwoQSpraysX1X2X3 letter) -- | Prints a ratio of sprays showRatioOfSpraysXYZ :: forall a. (Eq a, AlgField.C a) => [String] -- ^ typically some letters, to represent the variables -> (a -> String) -- ^ function mapping a coefficient to a string, typically 'show' -> (String, String) -- ^ used to enclose the coefficients, usually a pair of braces -> (String, String) -- ^ pair of braces to enclose the numerator and the denominator -> String -- ^ represents the quotient bar -> RatioOfSprays a -> String showRatioOfSpraysXYZ letters showCoef coeffBraces = showRatioOfSprays (showTwoSpraysXYZ showCoef coeffBraces letters) -- | Prints a ratio of sprays showRatioOfSpraysXYZ' :: (Eq a, AlgField.C a) => [String] -- ^ typically some letters, to represent the variables -> (a -> String) -- ^ function mapping a coefficient to a string, typically 'show' -> RatioOfSprays a -> String showRatioOfSpraysXYZ' letters showCoef = showRatioOfSpraysXYZ letters showCoef ("(", ")") ("[ ", " ]") " %//% " -- | Prints a ratio of sprays showRatioOfSpraysX1X2X3 :: forall a. (Eq a, AlgField.C a) => String -- ^ typically a letter, to represent the variables -> (a -> String) -- ^ function mapping a coefficient to a string, typically 'show' -> (String, String) -- ^ used to enclose the coefficients, usually a pair of braces -> (String, String) -- ^ pair of braces to enclose the numerator and the denominator -> String -- ^ represents the quotient bar -> RatioOfSprays a -> String showRatioOfSpraysX1X2X3 letter showCoef coeffBraces = showRatioOfSprays (showTwoSpraysX1X2X3 showCoef coeffBraces letter) -- | Prints a ratio of sprays showRatioOfSpraysX1X2X3' :: (Eq a, AlgField.C a) => String -- ^ typically a letter, to represent the variables -> (a -> String) -- ^ function mapping a coefficient to a string, typically 'show' -> RatioOfSprays a -> String showRatioOfSpraysX1X2X3' letter showCoef = showRatioOfSpraysX1X2X3 letter showCoef ("(", ")") ("[ ", " ]") " %//% " -- | Prints a ratio of sprays with rational coefficients prettyRatioOfQSpraysXYZ :: [String] -- ^ typically some letters, to represent the variables -> RatioOfQSprays -> String prettyRatioOfQSpraysXYZ letters = showRatioOfQSpraysXYZ letters ("[ ", " ]") " %//% " -- | Prints a ratio of sprays with rational coefficients -- -- prop> prettyRatioOfQSprays rOS == prettyRatioOfQSpraysXYZ ["x","y","z"] rOS prettyRatioOfQSprays :: RatioOfQSprays -> String prettyRatioOfQSprays = prettyRatioOfQSpraysXYZ ["x", "y", "z"] -- | Prints a ratio of sprays with rational coefficients -- -- prop> prettyRatioOfQSprays' rOS == prettyRatioOfQSpraysXYZ ["X","Y","Z"] rOS prettyRatioOfQSprays' :: RatioOfQSprays -> String prettyRatioOfQSprays' = prettyRatioOfQSpraysXYZ ["X", "Y", "Z"] -- | Prints a ratio of sprays with rational coefficients, printing the monomials -- in the style of @\"x1^2.x2.x3^3\"@ prettyRatioOfQSpraysX1X2X3 :: String -- ^ typically a letter, to represent the non-indexed variables -> RatioOfQSprays -> String prettyRatioOfQSpraysX1X2X3 letter = showRatioOfQSpraysX1X2X3 letter ("[ ", " ]") " %//% " -- | Prints a ratio of sprays with numeric coefficients prettyRatioOfNumSpraysXYZ :: (Num a, Ord a, AlgRing.C a, Show a) => [String] -- ^ typically some letters, to represent the variables -> RatioOfSprays a -> String prettyRatioOfNumSpraysXYZ letters = showRatioOfNumSpraysXYZ show letters ("[ ", " ]") " %//% " -- | Prints a ratio of sprays with numeric coefficients -- -- prop> prettyRatioOfNumSprays rOS == prettyRatioOfNumSpraysXYZ ["x","y","z"] rOS prettyRatioOfNumSprays :: (Num a, Ord a, AlgRing.C a, Show a) => RatioOfSprays a -> String prettyRatioOfNumSprays = prettyRatioOfNumSpraysXYZ ["x", "y", "z"] -- | Prints a ratio of sprays with numeric coefficients -- -- prop> prettyRatioOfNumSprays' rOS == prettyRatioOfNumSpraysXYZ ["X","Y","Z"] rOS prettyRatioOfNumSprays' :: (Num a, Ord a, AlgRing.C a, Show a) => RatioOfSprays a -> String prettyRatioOfNumSprays' = prettyRatioOfNumSpraysXYZ ["X", "Y", "Z"] -- | Prints a ratio of sprays with numeric coefficients, printing the monomials -- in the style of @\"x1^2.x2.x3^3\"@ prettyRatioOfNumSpraysX1X2X3 :: (Num a, Ord a, AlgRing.C a, Show a) => String -- ^ typically a letter, to represent the variables -> RatioOfSprays a -> String prettyRatioOfNumSpraysX1X2X3 letter = showRatioOfNumSpraysX1X2X3 show letter ("[ ", " ]") " %//% " -- Parametric sprays ---------------------------------------------------------- type SimpleParametricSpray a = Spray (Spray a) type SimpleParametricQSpray = SimpleParametricSpray Rational type ParametricSpray a = Spray (RatioOfSprays a) type ParametricQSpray = ParametricSpray Rational instance (Eq a, AlgRing.C a) => AlgMod.C a (SimpleParametricSpray a) where (*>) :: a -> SimpleParametricSpray a -> SimpleParametricSpray a lambda *> pspray = HM.map (lambda AlgMod.*>) pspray instance (Eq a, AlgRing.C a) => AlgRightMod.C a (SimpleParametricSpray a) where (<*) :: SimpleParametricSpray a -> a -> SimpleParametricSpray a pspray <* lambda = HM.map (AlgRightMod.<* lambda) pspray instance (Eq a, AlgField.C a) => AlgMod.C a (ParametricSpray a) where (*>) :: a -> ParametricSpray a -> ParametricSpray a lambda *> pspray = HM.map (lambda AlgMod.*>) pspray instance (Eq a, AlgField.C a) => AlgRightMod.C a (ParametricSpray a) where (<*) :: ParametricSpray a -> a -> ParametricSpray a pspray <* lambda = HM.map (AlgRightMod.<* lambda) pspray instance (Eq a, AlgField.C a) => AlgMod.C (Spray a) (ParametricSpray a) where (*>) :: Spray a -> ParametricSpray a -> ParametricSpray a spray *> pspray = asRatioOfSprays spray *^ pspray instance (Eq a, AlgField.C a) => AlgRightMod.C (Spray a) (ParametricSpray a) where (<*) :: ParametricSpray a -> Spray a -> ParametricSpray a pspray <* spray = asRatioOfSprays spray *^ pspray -- | Number of parameters in a parametric spray -- -- >>> numberOfParameters (jacobiPolynomial 4) -- 2 numberOfParameters :: FunctionLike b => Spray b -> Int numberOfParameters pspray = if isZeroSpray pspray then 0 else maximum (map numberOfVariables (allCoefficients pspray)) -- | Apply polynomial transformations to the parameters of a parametric spray; -- e.g. you have a two-parameters polynomial \(P_{a, b}(X, Y, Z)\) and you want -- to get \(P_{a^2, b^2}(X, Y, Z)\), or the one-parameter polynomial -- \(P_{a, a}(X, Y, Z)\) -- -- >>> jp = jacobiPolynomial 4 -- >>> a = qlone 1 -- >>> b = qlone 2 -- >>> changeParameters jp [a^**^2, b^**^2] changeParameters :: (FunctionLike b, Eq b, AlgAdd.C b) => Spray b -- ^ @OneParameterSpray a@, @SimpleParametricSpray a@, or @ParametricSpray a@ -> [VariablesType b] -- ^ @[Polynomial a]@ or @[Spray a]@, the new variables -> Spray b changeParameters pspray newParameters = if length newParameters < numberOfParameters pspray then error "changeParameters: not enough new parameters provided." else removeZeroTerms $ HM.map (`changeVariables` newParameters) pspray -- | Substitutes some values to the parameters of a parametric spray -- -- >>> jacobi3 = jacobiPolynomial 3 -- >>> legendre3 = substituteParameters jp [0, 0] substituteParameters :: (FunctionLike b, Eq (BaseRing b), AlgAdd.C (BaseRing b)) => Spray b -- ^ @OneParameterSpray a@, @SimpleParametricSpray a@, or @ParametricSpray a@ -> [BaseRing b] -- ^ values of type @a@ to be substituted to the parameters -> Spray (BaseRing b) -- ^ output: a @Spray a@ spray substituteParameters pspray values = if length values < numberOfParameters pspray then error "substituteParameters: not enough values provided." else removeZeroTerms $ HM.map (evaluateAt values) pspray -- | helper function for evalParametricSpray evalTerm' :: (AlgMod.C (BaseRing b) b) => [BaseRing b] -> Term b -> b evalTerm' xs (powers, coeff) = AlgRing.product (zipWith (AlgRing.^) xs pows) AlgMod.*> coeff where pows = DF.toList (fromIntegral <$> exponents powers) -- | Substitutes some values to the variables of a parametric spray evalParametricSpray :: (Eq b, AlgMod.C (BaseRing b) b, AlgRing.C b) => Spray b -- ^ @OneParameterSpray a@, @SimpleParametricSpray a@, or @ParametricSpray a@ -> [BaseRing b] -- ^ values of type @a@ to be substituted to the variables -> b evalParametricSpray spray xs = if length xs >= numberOfVariables spray then AlgAdd.sum $ map (evalTerm' xs) (HM.toList spray) else error "evalParametricSpray: not enough values provided." -- | Substitutes some values to the parameters of a parametric spray as well as -- some values to its variables evalParametricSpray' :: (FunctionLike b, Eq (BaseRing b), AlgMod.C (BaseRing b) b) => Spray b -- ^ @OneParameterSpray a@, @SimpleParametricSpray a@, or @ParametricSpray a@ -> [BaseRing b] -- ^ values of type @a@ to be substituted to the parameters -> [BaseRing b] -- ^ values of type @a@ to be substituted to the variables -> BaseRing b -- ^ result: a value of type @a@ evalParametricSpray' spray as xs = evaluateAt xs (substituteParameters spray as) -- | Whether the coefficients of a parametric spray polynomially -- depend on their parameters; I do not know why, but it seems to be the case -- for the Jacobi polynomials -- -- >>> canCoerceToSimpleParametricSpray (jacobiPolynomial 8) -- True canCoerceToSimpleParametricSpray :: (Eq a, AlgRing.C a) => ParametricSpray a -> Bool canCoerceToSimpleParametricSpray spray = all isPolynomialRatioOfSprays (allCoefficients spray) -- | Coerces a parametric spray to a simple parametric spray, without -- checking this makes sense with `canCoerceToSimpleParametricSpray` asSimpleParametricSprayUnsafe :: ParametricSpray a -> SimpleParametricSpray a asSimpleParametricSprayUnsafe = HM.map _numerator -- | Coerces a parametric spray to a simple parametric spray, after -- checking this makes sense with `canCoerceToSimpleParametricSpray` asSimpleParametricSpray :: (Eq a, AlgRing.C a) => ParametricSpray a -> SimpleParametricSpray a asSimpleParametricSpray spray = if canCoerceToSimpleParametricSpray spray then asSimpleParametricSprayUnsafe spray else error $ "asSimpleParametricSpray: this parametric spray is not coercible" ++ " to a simple parametric spray." -- | Converts a `OneParameterSpray a` spray to a `ParametricSpray a` fromOneParameterSpray :: (Eq a, AlgRing.C a) => OneParameterSpray a -> ParametricSpray a fromOneParameterSpray = HM.map fromRatioOfPolynomials -- | Converts a `OneParameterQSpray` spray to a `ParametricQSpray` fromOneParameterQSpray :: OneParameterQSpray -> ParametricQSpray fromOneParameterQSpray = HM.map fromRatioOfQPolynomials -- | Converts a `SimpleParametricSpray a` spray to a `ParametricSpray a` fromSimpleParametricSpray :: AlgRing.C a => SimpleParametricSpray a -> ParametricSpray a fromSimpleParametricSpray = HM.map asRatioOfSprays -- | Converts a parametric spray to a one-parameter spray, without checking -- the conversion makes sense parametricSprayToOneParameterSpray :: forall a. (Eq a, AlgField.C a) => ParametricSpray a -> OneParameterSpray a parametricSprayToOneParameterSpray = HM.map toRatioOfPolynomials where toRatioOfPolynomials :: RatioOfSprays a -> RatioOfPolynomials a toRatioOfPolynomials (RatioOfSprays p q) = toPolynomial p % toPolynomial q where toPolynomial :: Spray a -> Polynomial a toPolynomial spray = polyFromCoeffs coeffs where spray'' = removeConstantTerm'' (safeSpray spray) coeffs = getConstantTerm spray : [getCoefficient'' (S.singleton i) spray'' | i <- [1 .. deg]] deg = maximum (0 : expnts) expnts = map (`index` 0) (HM.keys spray'') {- -- | division of two univariate sprays longDivision :: (Eq a, AlgField.C a) => Spray a -> Spray a -> (Spray a, Spray a) longDivision sprayA sprayB = both fromCoeffs (polydiv coeffsA coeffsB) where fromCoeffs as = if S.null as then zeroSpray else sumTerms terms where l = S.length as terms = (Powers S.empty 0, as `index` (l-1)) : map (\i -> (Powers (S.singleton i) 1, as `index` (l-1-i))) [1 .. l-1] shift n l = l >< S.replicate n AlgAdd.zero pad n l = if n > 0 then S.replicate n AlgAdd.zero >< l else l zipWith' op xs ys = S.zipWith op (pad (-d) xs) (pad d ys) where d = S.length xs - S.length ys coeffsA = coefficientsUnivariateSpray sprayA coeffsB = coefficientsUnivariateSpray sprayB coefficientsUnivariateSpray spray = coeffs where coeffs = S.fromList [getCoefficient' (Powers (S.singleton i) 1) spray' | i <- [deg, deg-1 .. 1]] |> getConstantTerm spray deg = maximum (0 : expnts) spray' = removeConstantTerm spray expnts = map ((`index` 0) . exponents) (HM.keys spray') polydiv as bs = aux as bs S.empty where aux f s q | ddif < 0 = (q, f) | otherwise = aux f' s q' where ddif = S.length (norm f) - S.length (norm s) k = f `index` 0 AlgField./ s `index` 0 ks = (AlgRing.* k) <$> shift ddif s q' = zipWith' (AlgAdd.+) q $ shift ddif (S.singleton k) f' = norm $ S.drop 1 $ zipWith' (AlgAdd.-) f ks norm = S.dropWhileL (== AlgAdd.zero) -} -- | Converts a rational parametric spray to a rational one-parameter spray, -- without checking the conversion makes sense parametricQSprayToOneParameterQSpray :: ParametricQSpray -> OneParameterQSpray parametricQSprayToOneParameterQSpray = HM.map toRatioOfQPolynomials where toRatioOfQPolynomials :: RatioOfQSprays -> RatioOfQPolynomials toRatioOfQPolynomials (RatioOfSprays p q) = toQPolynomial p % toQPolynomial q where toQPolynomial :: QSpray -> QPolynomial toQPolynomial spray = polyFromCoeffs coeffs' where spray'' = removeConstantTerm'' (safeSpray spray) coeffs' = f (getConstantTerm spray) : [f $ getCoefficient'' (S.singleton i) spray'' | i <- [1 .. deg]] f :: Rational -> Rational' f r = DR.numerator r :% DR.denominator r deg = maximum (0 : map (`index` 0) expnts) expnts = HM.keys spray'' -- | [Gegenbauer polynomials](https://en.wikipedia.org/wiki/Gegenbauer_polynomials); -- we mainly provide them to give an example of the @SimpleParametricSpray@ type -- -- >>> gp = gegenbauerPolynomial 3 -- >>> putStrLn $ prettySimpleParametricQSpray gp -- { (4/3)*a^3 + 4*a^2 + (8/3)*a }*X^3 + { -2*a^2 - 2*a }*X -- >>> putStrLn $ prettyQSpray'' $ substituteParameters gp [1] -- 8*X^3 - 4*X gegenbauerPolynomial :: Int -> SimpleParametricQSpray gegenbauerPolynomial n | n == 0 = unitSpray | n == 1 = (2.^a) *^ x | otherwise = (2.^((n'-1) +> a) /^ n') *^ (x ^*^ gegenbauerPolynomial (n - 1)) ^-^ (((n'-1) +> 2.^a ^-^ unitSpray) /^ n') *^ gegenbauerPolynomial (n - 2) where x = lone 1 :: SimpleParametricQSpray a = lone 1 :: QSpray n' = toRational n -- | [Jacobi polynomial](https://en.wikipedia.org/wiki/Jacobi_polynomials); -- the @n@-th Jacobi polynomial is a univariate polynomial of degree @n@ with -- two parameters, except for the case @n=0@ where it has no parameter -- -- >>> jP = jacobiPolynomial 1 -- >>> putStrLn $ prettyParametricQSprayABCXYZ ["alpha", "beta"] ["X"] jP -- { [ (1/2)*alpha + (1/2)*beta + 1 ] }*X + { [ (1/2)*alpha - (1/2)*beta ] } jacobiPolynomial :: Int -> ParametricQSpray jacobiPolynomial n | n < 0 = error "jacobiPolynomial: `n` must be positive." | n == 0 = unitSpray | n == 1 = fromSimpleParametricSpray $ (((gamma0 <+ 2) /^ 2) *^ (x <+ AlgAdd.negate AlgRing.one)) <+ (alpha0 <+ 1) | otherwise = spray ^*^ jacobiPolynomial (n-1) ^-^ rOS *^ jacobiPolynomial (n-2) where alpha0 = qlone 1 beta0 = qlone 2 gamma0 = alpha0 ^+^ beta0 x = lone 1 :: SimpleParametricQSpray n' = toRational n a0 = alpha0 <+ (n' - 1) b0 = beta0 <+ (n' - 1) c0 = gamma0 <+ (2 * n') c0' = c0 <+ (-1) c0'' = c0 <+ (-2) divisor = (n' *^ (c0 <+ (-n'))) ^*^ c0'' divisor' = 2 .^ divisor divide = (`RatioOfSprays` divisor') spray = HM.fromList [ ( nullPowers , divide $ c0' ^*^ (alpha0 ^-^ beta0) ^*^ gamma0 ), ( Powers (S.singleton 1) 1 , divide $ c0' ^*^ c0 ^*^ c0'' ) ] rOS = RatioOfSprays (a0 ^*^ b0 ^*^ c0) divisor -- | Pretty form of a numeric parametric spray, using some given strings (typically some -- letters) to denote the parameters and some given strings (typically some letters) to -- denote the variables; rather use `prettyParametricQSprayABCXYZ` for a rational -- parametric spray prettyParametricNumSprayABCXYZ :: (Num a, Ord a, Show a, AlgField.C a) => [String] -- ^ usually some letters, to denote the parameters of the spray -> [String] -- ^ usually some letters, to denote the variables of the spray -> ParametricSpray a -- ^ a parametric numeric spray -> String prettyParametricNumSprayABCXYZ abc xyz spray = showSpray rOSShower ("{ ", " }") (showMonomialsXYZ xyz) spray where rOSShower = if numberOfParameters spray <= length abc then prettyRatioOfNumSpraysXYZ abc else prettyRatioOfNumSpraysX1X2X3 (abc !! 0) -- | Pretty form of a numeric parametric spray; rather use `prettyParametricQSpray` for -- a rational parametric spray -- -- prop> prettyParametricNumSpray == prettyParametricNumSprayABCXYZ ["a"] ["X","Y","Z"] prettyParametricNumSpray :: (Num a, Ord a, Show a, AlgField.C a) => ParametricSpray a -- ^ a parametric numeric spray -> String prettyParametricNumSpray = prettyParametricNumSprayABCXYZ ["a"] ["X", "Y", "Z"] -- | Pretty form of a parametric rational spray, using some given strings (typically some -- letters) to denote the parameters and some given strings (typically some letters) to -- denote the variables -- -- >>> type PQS = ParametricQSpray -- >>> :{ -- >>> f :: (QSpray, QSpray) -> (PQS, PQS, PQS) -> PQS -- >>> f (a, b) (x, y, z) = -- >>> (a %:% (a ^+^ b)) *^ x^**^2 ^+^ (b %:% (a ^+^ b)) *^ (y ^*^ z) -- >>> :} -- >>> a = qlone 1 -- >>> b = qlone 2 -- >>> x = lone 1 :: PQS -- >>> y = lone 2 :: PQS -- >>> z = lone 3 :: PQS -- >>> pqs = f (a, b) (x, y, z) -- >>> putStrLn $ prettyParametricQSprayABCXYZ ["a","b"] ["X","Y","Z"] pqs -- { [ a ] %//% [ a + b ] }*X^2 + { [ b ] %//% [ a + b ] }*Y.Z prettyParametricQSprayABCXYZ :: [String] -- ^ usually some letters, to denote the parameters of the spray -> [String] -- ^ usually some letters, to denote the variables of the spray -> ParametricQSpray -- ^ a parametric rational spray -> String prettyParametricQSprayABCXYZ abc xyz spray = showSpray rOSShower ("{ ", " }") (showMonomialsXYZ xyz) spray where rOSShower = if numberOfParameters spray <= length abc then prettyRatioOfQSpraysXYZ abc else prettyRatioOfQSpraysX1X2X3 (abc !! 0) -- | Pretty form of a parametric rational spray -- -- prop> prettyParametricQSpray == prettyParametricQSprayABCXYZ ["a"] ["X","Y","Z"] prettyParametricQSpray :: ParametricQSpray -> String prettyParametricQSpray = prettyParametricQSprayABCXYZ ["a"] ["X", "Y", "Z"] -- | Pretty form of a numeric simple parametric spray, using some given strings (typically some -- letters) to denote the parameters and some given strings (typically some letters) to -- denote the variables; rather use `prettySimpleParametricQSprayABCXYZ` for a rational -- simple parametric spray prettySimpleParametricNumSprayABCXYZ :: (Num a, Ord a, Show a, AlgRing.C a) => [String] -- ^ usually some letters, to denote the parameters of the spray -> [String] -- ^ usually some letters, to denote the variables of the spray -> SimpleParametricSpray a -- ^ a numeric simple parametric spray -> String prettySimpleParametricNumSprayABCXYZ abc xyz spray = showSpray rOSShower ("{ ", " }") (showMonomialsXYZ xyz) spray where rOSShower = if numberOfParameters spray <= length abc then prettyNumSprayXYZ abc else prettyNumSprayX1X2X3 (abc !! 0) -- | Pretty form of a numeric simple parametric spray; rather use -- `prettySimpleParametricQSpray` for a numeric simple parametric spray -- -- prop> prettySimpleParametricNumSpray == prettySimpleParametricNumSprayABCXYZ ["a"] ["X","Y","Z"] prettySimpleParametricNumSpray :: (Num a, Ord a, Show a, AlgRing.C a) => SimpleParametricSpray a -- ^ a numeric simple parametric spray -> String prettySimpleParametricNumSpray = prettySimpleParametricNumSprayABCXYZ ["a"] ["X", "Y", "Z"] -- | Pretty form of a simple parametric rational spray, using some given strings (typically some -- letters) to denote the parameters and some given strings (typically some letters) to -- denote the variables -- -- >>> type SPQS = SimpleParametricQSpray -- >>> :{ -- >>> f :: (QSpray, QSpray) -> (SPQS, SPQS, SPQS) -> SPQS -- >>> f (a, b) (x, y, z) = -- >>> (a ^+^ b) *^ x^**^2 ^+^ (a^**^2 ^+^ b^**^2) *^ (y ^*^ z) -- >>> :} -- >>> a = qlone 1 -- >>> b = qlone 2 -- >>> x = lone 1 :: SPQS -- >>> y = lone 2 :: SPQS -- >>> z = lone 3 :: SPQS -- >>> spqs = f (a, b) (x, y, z) -- >>> putStrLn $ prettySimpleParametricQSprayABCXYZ ["a","b"] ["X","Y","Z"] spqs -- { a + b }*X^2 + { a^2 + b^2 }*Y.Z prettySimpleParametricQSprayABCXYZ :: [String] -- ^ usually some letters, to denote the parameters of the spray -> [String] -- ^ usually some letters, to denote the variables of the spray -> SimpleParametricQSpray -- ^ a parametric rational spray -> String prettySimpleParametricQSprayABCXYZ abc xyz spray = showSpray sprayShower ("{ ", " }") (showMonomialsXYZ xyz) spray where sprayShower = if numberOfParameters spray <= length abc then prettyQSprayXYZ abc else prettyQSprayX1X2X3 (abc !! 0) -- | Pretty form of a simple parametric rational spray -- -- prop> prettySimpleParametricQSpray == prettySimpleParametricQSprayABCXYZ ["a"] ["X","Y","Z"] prettySimpleParametricQSpray :: SimpleParametricQSpray -> String prettySimpleParametricQSpray = prettySimpleParametricQSprayABCXYZ ["a"] ["X", "Y", "Z"]