hspray 0.2.0.0 → 0.2.1.0

raw patch · 5 files changed

+329/−21 lines, 5 filesdep +matrixdep +tasty-benchdep ~tastydep ~tasty-hunitPVP ok

version bump matches the API change (PVP)

Dependencies added: matrix, tasty-bench

Dependency ranges changed: tasty, tasty-hunit

API changes (from Hackage documentation)

Files

• CHANGELOG.md +9/−0
• benchmarks/Main.hs +15/−0
• hspray.cabal +16/−3
• src/Math/Algebra/Hspray.hs +217/−16
• tests/Main.hs +72/−2

@@ -46,3 +46,12 @@ * New function `isSymmetricSpray`, to check whether a spray is a symmetric polynomial.
 
 * New function `isPolynomialOf`, to check whether a spray can be expressed as a polynomial of a given list of sprays.
+
+
+## 0.2.1.0 - 2024-03-22
+
+* New functions `permuteVariables` and `swapVariables`.
+
+* New function `resultant`, to compute the resultant of two sprays.
+
+* New function `subresultants`, to compute the subresultants of two sprays.

@@ -0,0 +1,15 @@+module Main where
+import Test.Tasty.Bench ( bench, bgroup, defaultMain, nf )
+
+fibo :: Int -> Integer
+fibo n = if n < 2 then toInteger n else fibo (n - 1) + fibo (n - 2)
+
+main :: IO ()
+main = 
+  defaultMain
+    [ bgroup "Fibonacci numbers"
+      [ bench "fifth"     $ nf fibo  5
+      , bench "tenth"     $ nf fibo 10
+      , bench "twentieth" $ nf fibo 20
+      ]
+    ]

@@ -1,5 +1,5 @@ name:                hspray
-version:             0.2.0.0
+version:             0.2.1.0
 synopsis:            Multivariate polynomials.
 description:         Manipulation of multivariate polynomials on a ring and Gröbner bases.
 homepage:            https://github.com/stla/hspray#readme
@@ -23,6 +23,7 @@                      , unordered-containers >= 0.2.17.0 && < 0.3
                      , numeric-prelude >= 0.4.4 && < 0.5
                      , text >= 1.2.5.0 && < 2.2
+                     , matrix >= 0.3.6.0 && < 0.4
   other-extensions:    FlexibleInstances
                      , FlexibleContexts
                      , MultiParamTypeClasses
@@ -46,10 +47,22 @@   hs-source-dirs:       tests/
   other-modules:        Approx
   Build-Depends:        base >= 4.7 && < 5
-                      , tasty
-                      , tasty-hunit
+                      , tasty >= 1.4 && < 1.6
+                      , tasty-hunit >= 0.10 && < 0.11
                       , hspray
   Default-Language:     Haskell2010
+
+benchmark benchmarks
+  type:                 exitcode-stdio-1.0
+  main-is:              Main.hs
+  hs-source-dirs:       benchmarks/
+  Build-Depends:        base >= 4.7 && < 5
+                      , tasty-bench >= 0.3 && < 0.4
+                      , hspray
+  ghc-options:          "-with-rtsopts=-A32m"
+  if impl(ghc >= 8.6)
+    ghc-options:        -fproc-alignment=64
+  default-language:     Haskell2010
 
 source-repository head
   type:     git

@@ -20,7 +20,9 @@   , Monomial
   , lone
   , unitSpray
+  , zeroSpray
   , constantSpray
+  , getCoefficient
   , fromList
   , toList
   , sprayTerms
@@ -33,6 +35,8 @@   , evalSpray
   , substituteSpray
   , fromRationalSpray
+  , permuteVariables
+  , swapVariables
   , prettySpray
   , prettySpray'
   , prettySprayXYZ
@@ -45,6 +49,10 @@   , esPolynomial
   , isSymmetricSpray
   , isPolynomialOf
+  , resultant
+  , resultant1
+  , subresultants
+  , subresultants1
   ) where
 import qualified Algebra.Additive              as AlgAdd
 import qualified Algebra.Field                 as AlgField
@@ -55,13 +63,23 @@ 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
+                                                , nub
+                                                , foldl1'
                                                 )
+import           Data.Matrix                    ( Matrix 
+                                                , fromLists
+                                                , minorMatrix
+                                                , nrows
+                                                , submatrix
+                                                )
+import qualified Data.Matrix                   as DM
 import           Data.Maybe                     ( isJust
-                                                , fromJust 
+                                                , fromJust, fromMaybe
                                                 )
 import           Data.Ord                       ( comparing )
 import qualified Data.Sequence                 as S
@@ -102,6 +120,9 @@ growSequence :: Seq Int -> Int -> Int -> Seq Int
 growSequence s m n = s >< t where t = S.replicate (n - m) 0
 
+growSequence' :: Int -> Seq Int -> Seq Int
+growSequence' n s = growSequence s (S.length s) n
+
 -- | append trailing zeros to get the same length
 harmonize :: (Powers, Powers) -> (Powers, Powers)
 harmonize (pows1, pows2) = (Powers e1' n, Powers e2' n)
@@ -139,6 +160,15 @@   p * q = multSprays p q
   one = lone 0
 
+{- instance (AlgRing.C a, Eq a) => Num (Spray a) where
+  p + q = addSprays p q
+  negate = negateSpray
+  p * q = multSprays p q
+  fromInteger n = fromInteger n .^ AlgRing.one
+  abs _ = error "Prelude.Num.abs: inappropriate abstraction"
+  signum _ = error "Prelude.Num.signum: inappropriate abstraction"
+ -} 
+
 -- | Addition of two sprays
 (^+^) :: (AlgAdd.C a, Eq a) => Spray a -> Spray a -> Spray a
 (^+^) p q = p AlgAdd.+ q
@@ -243,10 +273,27 @@ unitSpray :: AlgRing.C a => Spray a
 unitSpray = lone 0
 
+-- | The null spray
+zeroSpray :: (Eq a, AlgAdd.C a) => Spray a
+zeroSpray = AlgAdd.zero
+
 -- | Constant spray
 constantSpray :: (AlgRing.C a, Eq a) => a -> Spray a
 constantSpray c = c *^ lone 0
 
+-- | Get coefficient of a term of a spray 
+getCoefficient :: AlgAdd.C a => [Int] -> Spray a -> a
+getCoefficient expnts spray = fromMaybe AlgAdd.zero (HM.lookup powers spray)
+  where
+    expnts' = S.dropWhileR (== 0) (S.fromList expnts)
+    powers = Powers expnts' (S.length expnts')
+
+-- | number of variables in a spray
+numberOfVariables :: Spray a -> Int
+numberOfVariables spray = maximum (map nvariables powers)
+  where
+    powers = HM.keys spray
+
 -- | evaluates a monomial
 evalMonomial :: AlgRing.C a => [a] -> Monomial a -> a
 evalMonomial xyz (powers, coeff) = 
@@ -255,18 +302,14 @@ 
 -- | Evaluate a spray
 evalSpray :: AlgRing.C a => Spray a -> [a] -> a
-evalSpray p xyz = AlgAdd.sum $ map (evalMonomial xyz) (HM.toList p)
+evalSpray p xyz = if length xyz >= numberOfVariables p
+  then AlgAdd.sum $ map (evalMonomial xyz) (HM.toList p)
+  else error "evalSpray: not enough values provided."
 
 -- | spray from monomial
 fromMonomial :: Monomial a -> Spray a
 fromMonomial (pows, coeff) = HM.singleton pows coeff
 
--- | number of variables in a spray
-numberOfVariables :: Spray a -> Int
-numberOfVariables spray = maximum (map nvariables powers)
-  where
-    powers = HM.keys spray
-
 -- | substitute some variables in a monomial
 substituteMonomial :: AlgRing.C a => [Maybe a] -> Monomial a -> Monomial a
 substituteMonomial subs (powers, coeff) = (powers'', coeff')
@@ -285,7 +328,7 @@ substituteSpray :: (Eq a, AlgRing.C a) => [Maybe a] -> Spray a -> Spray a
 substituteSpray subs spray = if length subs == n 
   then spray'
-  else error "incorrect length of the substitutions list"
+  else error "substituteSpray: incorrect length of the substitutions list."
   where
     n = numberOfVariables spray
     monomials = HM.toList spray
@@ -295,20 +338,57 @@ fromRationalSpray :: Spray Rational -> Spray Double
 fromRationalSpray = HM.map fromRational
 
--- | helper for `composeSpray`
-identify :: (AlgRing.C a, Eq a) => Spray a -> Spray (Spray a)
-identify = HM.map constantSpray
-
 -- | Compose a spray with a change of variables
 composeSpray :: (AlgRing.C a, Eq a) => Spray a -> [Spray a] -> Spray a
 composeSpray p = evalSpray (identify p)
+  where 
+    ---- identify :: (AlgRing.C a, Eq a) => Spray a -> Spray (Spray a)
+    identify = HM.map constantSpray
 
 -- | Create a spray from list of terms
 fromList :: (AlgRing.C a, Eq a) => [([Int], a)] -> Spray a
 fromList x = cleanSpray $ HM.fromList $ map
   (\(expts, coef) -> (Powers (S.fromList expts) (length expts), coef)) x
 
+-- | Permute the variables of a spray
+permuteVariables :: Spray a -> [Int] -> Spray a
+permuteVariables spray permutation = 
+  if n' >= n && isPermutation permutation  
+    then spray'
+    else error "permuteVariables: invalid permutation."
+  where
+    n = numberOfVariables spray
+    n' = maximum permutation
+    isPermutation 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)
+    expnts = map exponents powers
+    expnts' = map (permuteSeq . growSequence' n') expnts
+    powers' = map (\exps -> simplifyPowers (Powers exps n')) expnts'
+    spray' = HM.fromList (zip powers' coeffs)
 
+-- | Swap two variables of a spray
+swapVariables :: Spray a -> (Int, Int) -> Spray a
+swapVariables spray (i, j) = 
+  if i>=1 && j>=1  
+    then spray'
+    else error "swapVariables: invalid indices."
+  where
+    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)
+    expnts = map exponents powers
+    expnts' = map (permuteSeq . growSequence' n) expnts
+    powers' = map (\exps -> simplifyPowers (Powers exps n)) expnts'
+    spray' = HM.fromList (zip powers' coeffs)
+
+
 -- pretty stuff ---------------------------------------------------------------
 
 -- | prettyPowers "x" [0, 2, 1] = x^(0, 2, 1)
@@ -450,7 +530,9 @@ -- using the lexicographic ordering of the monomials
 sprayDivision :: forall a. (Eq a, AlgField.C a) => Spray a -> [Spray a] -> Spray a
 sprayDivision p qs = 
-  if n == 0 then error "the list of divisors is empty" else snd $ ogo p AlgAdd.zero
+  if n == 0 
+    then error "sprayDivision: the list of divisors is empty." 
+    else snd $ ogo p AlgAdd.zero
   where
     n = length qs
     qsltqs = zip qs (map leadingTerm qs)
@@ -650,7 +732,7 @@   -> Int -- ^ index
   -> Spray a
 esPolynomial n k
-  | k <= 0 || n <= 0 = error "both arguments must be positive integers"
+  | k <= 0 || n <= 0 = error "esPolynomial: both arguments must be positive integers."
   | k > n = AlgAdd.zero
   | otherwise = simplifySpray spray
   where
@@ -699,4 +781,123 @@             then Just $ dropXis g
             else Nothing
           dropXis = HM.mapKeys f
-          f (Powers expnnts _) = Powers (S.drop n expnnts) n+          f (Powers expnnts _) = Powers (S.drop n expnnts) n
+
+
+-- 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: determinant=1 (because the empty matrix is not allowed)
+  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]]
+
+-- determinant
+detLaplace :: forall a. (Eq a, AlgRing.C a) => Matrix a -> a
+detLaplace m = if nrows m == 1 
+  then m DM.! (1,1)
+  else suml1 [negateIf i (times (m DM.! (i,1)) (detLaplace (minorMatrix i 1 m))) | i <- [1 .. nrows m]]
+  where 
+    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
+
+-- the coefficients of a spray as a spray with univariate spray coefficients
+sprayCoefficients :: (Eq a, AlgRing.C a) => Spray a -> [Spray a]
+sprayCoefficients spray = reverse sprays
+  where
+    (powers, coeffs) = unzip (HM.toList spray)
+    expnts = map exponents powers
+    constantTerm = fromMaybe AlgAdd.zero (HM.lookup (Powers S.empty 0) spray)
+    (expnts', coeffs') = unzip $ filter (\(s,_) -> S.length s > 0) (zip expnts coeffs)
+    xpows = map (`index` 0) expnts'
+    expnts'' = map (S.deleteAt 0) expnts'
+    powers'' = map (\s -> Powers s (S.length s)) expnts''
+    sprays'' = zipWith (curry fromMonomial) powers'' coeffs'
+    imap = IM.fromListWith (^+^) (zip xpows sprays'')
+    imap' = IM.insertWith (^+^) 0 (constantSpray constantTerm) imap
+    sprays = [fromMaybe AlgAdd.zero (IM.lookup i imap') | i <- [0 .. maximum xpows]]
+
+-- | Resultant of two univariate sprays
+resultant1 :: (Eq a, AlgRing.C a) => Spray a -> Spray a -> a
+resultant1 p q = detLaplace $ sylvesterMatrix pcoeffs qcoeffs
+  where
+    pexpnts = map (`index` 0) $ filter (not . S.null) (map exponents (HM.keys p))
+    qexpnts = map (`index` 0) $ filter (not . S.null) (map exponents (HM.keys q))
+    p0 = fromMaybe AlgAdd.zero (HM.lookup (Powers S.empty 0) p)
+    q0 = fromMaybe AlgAdd.zero (HM.lookup (Powers S.empty 0) q)
+    pcoeffs = reverse $ if null pexpnts 
+      then [p0]
+      else p0 : [fromMaybe AlgAdd.zero (HM.lookup (Powers (S.singleton i) 1) p) | i <- [1 .. maximum pexpnts]]
+    qcoeffs = reverse $ if null qexpnts 
+      then [q0]
+      else q0 : [fromMaybe AlgAdd.zero (HM.lookup (Powers (S.singleton i) 1) q) | i <- [1 .. maximum qexpnts]]
+
+-- | Subresultants of two univariate sprays
+subresultants1 :: (Eq a, AlgRing.C a) => Spray a -> Spray a -> [a]
+subresultants1 p q = map (detLaplace . sylvesterMatrix' pcoeffs qcoeffs) [0 .. min d e - 1]
+  where
+    pexpnts = map (`index` 0) $ filter (not . S.null) (map exponents (HM.keys p))
+    qexpnts = map (`index` 0) $ filter (not . S.null) (map exponents (HM.keys q))
+    p0 = fromMaybe AlgAdd.zero (HM.lookup (Powers S.empty 0) p)
+    q0 = fromMaybe AlgAdd.zero (HM.lookup (Powers S.empty 0) q)
+    pcoeffs = reverse $ if null pexpnts 
+      then [p0]
+      else p0 : [fromMaybe AlgAdd.zero (HM.lookup (Powers (S.singleton i) 1) p) | i <- [1 .. maximum pexpnts]]
+    qcoeffs = reverse $ if null qexpnts 
+      then [q0]
+      else q0 : [fromMaybe AlgAdd.zero (HM.lookup (Powers (S.singleton i) 1) q) | i <- [1 .. 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 detLaplace $ sylvesterMatrix (sprayCoefficients p') (sprayCoefficients q')
+    else error "resultant: invalid variable index."
+  where
+    n = max (numberOfVariables p) (numberOfVariables q)
+    permutation = var : [1 .. var-1] ++ [var+1 .. n]
+    p' = permuteVariables p permutation
+    q' = permuteVariables q permutation
+
+-- | Subresultants of two sprays
+subresultants :: (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]
+subresultants var p q 
+  | var < 1 = error "subresultants: invalid variable index."
+  | var > n = error "subresultants: too large variable index."
+  | otherwise = map (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]
+    p' = permuteVariables p permutation
+    q' = permuteVariables q permutation

@@ -9,9 +9,14 @@                                                   (*^),
                                                   lone,
                                                   unitSpray,
+                                                  zeroSpray,
+                                                  constantSpray,
+                                                  getCoefficient,
                                                   evalSpray,
                                                   substituteSpray,
                                                   composeSpray,
+                                                  permuteVariables,
+                                                  swapVariables,
                                                   fromList,
                                                   toList,
                                                   bombieriSpray,
@@ -20,7 +25,11 @@                                                   fromRationalSpray,
                                                   esPolynomial,
                                                   isSymmetricSpray, 
-                                                  isPolynomialOf 
+                                                  isPolynomialOf,
+                                                  resultant,
+                                                  subresultants,
+                                                  resultant1,
+                                                  subresultants1
                                                 )
 import           Test.Tasty                     ( defaultMain
                                                 , testGroup
@@ -59,6 +68,14 @@         pxyz = map (`evalSpray` xyz) [px, py, pz]
       assertEqual "" (evalSpray p pxyz) (evalSpray q xyz),
 
+    testCase "getCoefficient" $ do
+      let
+        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
+      assertEqual "" (getCoefficient [3, 2, 0] p, getCoefficient [0, 4] p) (4, 0),
+
     testCase "fromList . toList = identity" $ do
       let
         x = lone 1 :: Spray Int
@@ -117,5 +134,58 @@         x3 = lone 3 :: Spray Rational
         p = x1^**^2 ^+^ x2 ^+^ x3 ^-^ unitSpray
         p' = substituteSpray [Just 2, Nothing, Just 3] p
-      assertEqual "" p' (x2 ^+^ (6*^ unitSpray))
+      assertEqual "" p' (x2 ^+^ (6*^ unitSpray)),
+
+    testCase "permuteVariables" $ do
+      let
+        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
+        p = f x1 x2 x3
+        p' = permuteVariables p [3, 1, 2]
+      assertEqual "" p' (f x3 x1 x2),
+
+    testCase "swapVariables" $ do
+      let
+        x1 = lone 1 :: Spray Rational
+        x2 = lone 2 :: Spray Rational
+        x3 = lone 3 :: Spray Rational
+        p = x1^**^4 ^+^ (2 *^ x2^**^3) ^+^ (3 *^ x3^**^2) ^-^ (4 *^ unitSpray)
+        p' = permuteVariables p [3, 2, 1]
+      assertEqual "" p' (swapVariables p (1, 3)),
+
+    testCase "resultant" $ do
+      let
+        x = lone 1 :: Spray Rational
+        y = lone 2 :: Spray Rational
+        p = x^**^4 ^-^ x^**^3 ^+^ x^**^2 ^-^ 2*^ (x ^*^ y^**^2) ^+^ y^**^4 
+        q = x ^-^ (2*^ y^**^2)
+        r = resultant 1 p q
+      assertEqual "" r (x^**^4 ^-^ (8*^ x^**^6) ^+^ (16*^ x^**^8)),
+
+    testCase "subresultants" $ do
+      let
+        x = lone 1 :: Spray Rational
+        y = lone 2 :: Spray Rational
+        p = x^**^2 ^*^ y ^*^ (y^**^2 ^-^ 5*^ x ^+^ constantSpray 6) 
+        q = x^**^2 ^*^ y ^*^ (3*^ y ^+^ constantSpray 2)
+        sx = subresultants 1 p q
+      assertBool "" (sx!!0 == zeroSpray && sx!!1 == zeroSpray && sx!!2 /= zeroSpray),
+
+    testCase "resultant1" $ do
+      let
+        x = lone 1 :: Spray Rational
+        p = x^**^2 ^-^ 5*^x ^+^ constantSpray 6 
+        q = x^**^2 ^-^ 3*^x ^+^ constantSpray 2 
+      assertEqual "" (resultant1 p q) (0%1),
+
+    testCase "subresultants1" $ do
+      let
+        x = lone 1 :: Spray Rational
+        p = x^**^2 ^-^ 5*^x ^+^ constantSpray 6 
+        q = x^**^2 ^-^ 3*^x ^+^ constantSpray 2 
+      assertEqual "" (subresultants1 p q) [0%1, 2%1, 1%1]
+
   ]