hspray 0.2.2.0 → 0.2.3.0

raw patch · 4 files changed

+192/−20 lines, 4 filesPVP ok

version bump matches the API change (PVP)

API changes (from Hackage documentation)

+ Math.Algebra.Hspray: gcdQX :: Spray Rational -> Spray Rational -> Spray Rational

Files

CHANGELOG.md +10/−1
hspray.cabal +2/−2
src/Math/Algebra/Hspray.hs +121/−13
tests/Main.hs +59/−4

CHANGELOG.md view

@@ -66,4 +66,13 @@ 
 ## 0.2.2.0 - 2024-03-26
 
-* Fixed an error in `esPolynomial`, which resulted to a bug in `isSymmetricSpray`.+* Fixed an error in `esPolynomial`, which resulted to a bug in `isSymmetricSpray`.
+
+
+## 0.2.3.0 - 2024-03-28
+
+* New unit tests.
+
+* Fixed `resultant` and `subresultants`: the variables of the sprays they return were incorrect.
+
+* New function `gcdQX`, to compute the greatest common divisor of two univariate sprays with rational coefficients.

hspray.cabal view

@@ -1,7 +1,7 @@ name:                hspray
-version:             0.2.2.0
+version:             0.2.3.0
 synopsis:            Multivariate polynomials.
-description:         Manipulation of multivariate polynomials on a ring and Gröbner bases.
+description:         Manipulation of multivariate polynomials on a ring, Gröbner basis, resultant and subresultants.
 homepage:            https://github.com/stla/hspray#readme
 license:             GPL-3
 license-file:        LICENSE

src/Math/Algebra/Hspray.hs view

@@ -61,6 +61,8 @@   , resultant1
   , subresultants
   , subresultants1
+  -- * Greatest common divisor
+  , gcdQX
   -- * Miscellaneous
   , fromList
   , toList
@@ -83,6 +85,7 @@                                                 , maximumBy 
                                                 , (\\)
                                                 , findIndices
+                                                , elemIndices
                                                 , nub
                                                 , foldl1'
                                                 )
@@ -334,7 +337,8 @@ 
 -- | number of variables in a spray
 numberOfVariables :: Spray a -> Int
-numberOfVariables spray = maximum (map nvariables powers)
+numberOfVariables spray =
+  if null powers then 0 else maximum (map nvariables powers)
   where
     powers = HM.keys spray
 
@@ -390,7 +394,7 @@   where
     n = numberOfVariables spray
     monomials = HM.toList spray
-    spray' = foldl1 (^+^) (map (fromMonomial . substituteMonomial subs) monomials)
+    spray' = foldl1' (^+^) (map (fromMonomial . substituteMonomial subs) monomials)
 
 -- | Converts a spray with rational coefficients to a spray with double coefficients
 -- (useful for evaluation)
@@ -447,7 +451,7 @@ 
 -- | Swaps two variables of a spray
 -- 
--- prop> swapVariables (1, 3) p == permuteVariables p [3, 2, 1]
+-- prop> swapVariables p (1, 3) == permuteVariables p [3, 2, 1]
 swapVariables :: Spray a -> (Int, Int) -> Spray a
 swapVariables spray (i, j) = 
   if i>=1 && j>=1  
@@ -855,7 +859,7 @@     esPolys = map (\i -> esPolynomial n i :: Spray a) indices
     yPolys = map (\i -> lone (n + i) :: Spray a) indices
     gPolys = zipWith (^-^) esPolys yPolys
-    gbasis = groebner gPolys False
+    gbasis = groebner0 gPolys
     g = sprayDivision spray gbasis
     gpowers = HM.keys g
     check1 = minimum (map nvariables gpowers) > n
@@ -934,10 +938,14 @@     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
+-- the coefficients of a spray as a univariate spray in x with spray coefficients
 sprayCoefficients :: (Eq a, AlgRing.C a) => Spray a -> [Spray a]
-sprayCoefficients spray = reverse sprays
+sprayCoefficients spray = 
+  if n == 0 
+    then [constantSpray constantTerm]
+    else reverse sprays
   where
+    n = numberOfVariables spray 
     (powers, coeffs) = unzip (HM.toList spray)
     expnts = map exponents powers
     constantTerm = fromMaybe AlgAdd.zero (HM.lookup (Powers S.empty 0) spray)
@@ -948,12 +956,20 @@     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]]
+    permutation = [2 .. n] ++ [1]
+    sprays = [
+        permuteVariables (fromMaybe AlgAdd.zero (IM.lookup i imap')) permutation 
+        | i <- [0 .. maximum xpows]
+      ]
 
--- | Resultant of two univariate sprays
+-- | Resultant of two /univariate/ sprays
 resultant1 :: (Eq a, AlgRing.C a) => Spray a -> Spray a -> a
-resultant1 p q = detLaplace $ sylvesterMatrix pcoeffs qcoeffs
+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)
     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)
@@ -965,10 +981,13 @@       then [q0]
       else q0 : [fromMaybe AlgAdd.zero (HM.lookup (Powers (S.singleton i) 1) q) | i <- [1 .. maximum qexpnts]]
 
--- | Subresultants of two univariate sprays
+-- | 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]
+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)
     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)
@@ -990,13 +1009,15 @@   -> Spray a
 resultant var p q = 
   if var >= 1 && var <= n 
-    then detLaplace $ sylvesterMatrix (sprayCoefficients p') (sprayCoefficients q')
+    then permuteVariables det permutation'
     else error "resultant: invalid variable index."
   where
     n = max (numberOfVariables p) (numberOfVariables q)
     permutation = var : [1 .. var-1] ++ [var+1 .. n]
+    permutation' = [2 .. var] ++ (1 : [var+1 .. n])
     p' = permuteVariables p permutation
     q' = permuteVariables q permutation
+    det = detLaplace $ sylvesterMatrix (sprayCoefficients p') (sprayCoefficients q')
 
 -- | Subresultants of two sprays
 subresultants :: (Eq a, AlgRing.C a) 
@@ -1007,7 +1028,7 @@ 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]
+  | otherwise = map (permute . detLaplace . sylvesterMatrix' pcoeffs qcoeffs) [0 .. min d e - 1]
   where
     pcoeffs = sprayCoefficients p'
     qcoeffs = sprayCoefficients q'
@@ -1017,3 +1038,90 @@     permutation = var : [1 .. var-1] ++ [var+1 .. n]
     p' = permuteVariables p permutation
     q' = permuteVariables q permutation
+    permutation' = [2 .. var] ++ (1 : [var+1 .. n])
+    permute spray = permuteVariables spray permutation'
+
+
+-- GCD ------------------------------------------------------------------------
+
+-- the degree of a spray as a univariate spray in x with spray coefficients
+degree :: Eq a => Spray a -> Int
+degree spray = 
+  if numberOfVariables spray == 0 
+    then 
+      if spray == HM.empty then minBound else 0
+    else maximum xpows
+  where
+    expnts = map exponents $ HM.keys spray
+    expnts' = filter (not . S.null) expnts
+    xpows = map (`index` 0) expnts'
+
+-- the degree and the leading coefficient of a spray as a univariate spray in x with spray coefficients
+degreeAndLeadingCoefficient :: (Eq a, AlgRing.C a) => Spray a -> (Int, Spray a)
+degreeAndLeadingCoefficient spray = 
+  if n == 0 
+    then (if constantTerm == AlgAdd.zero then minBound else 0, constantSpray constantTerm)
+    else (deg, leadingCoeff)
+  where
+    n = numberOfVariables spray
+    (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,_) -> not $ S.null s) (zip expnts coeffs)
+    xpows = map (`index` 0) expnts'
+    deg = maximum xpows
+    is = elemIndices deg xpows
+    expnts'' = [S.update 0 0 (expnts' !! i) | i <- is]
+    powers'' = map (\s -> Powers s (S.length s)) expnts''
+    coeffs'' = [coeffs' !! i | i <- is]
+    leadingCoeff = foldl1' (^+^) (zipWith (curry fromMonomial) powers'' coeffs'')
+
+-- pseudo-division of two sprays (assuming degA >= degB >= 0)
+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 = (ellB ^**^ delta , go sprayA zeroSpray delta)
+  where
+    degA = degree sprayA
+    (degB, ellB) = degreeAndLeadingCoefficient sprayB
+    delta = degA - degB + 1
+    go sprayR sprayQ e = 
+      if degR < degB
+        then (q ^*^ sprayQ, q ^*^ sprayR)
+        else go (ellB ^*^ sprayR ^-^ sprayS ^*^ sprayB) (ellB ^*^ sprayQ ^+^ sprayS) (e - 1)
+      where
+        (degR, ellR) = degreeAndLeadingCoefficient sprayR
+        q = ellB ^**^ e
+        sprayX = lone 1
+        sprayS = ellR ^*^ sprayX ^**^ (degR - degB)
+
+-- The GCD of two /univariate/ sprays with rational coefficients
+gcdQX :: Spray Rational -> Spray Rational -> Spray Rational
+gcdQX sprayA sprayB
+  | n >= 2                        = error "gcdQX: the sprays are not univariate." 
+  | degree sprayB > degree sprayA = gcdQX sprayB sprayA 
+  | sprayB == zeroSpray           = sprayA
+  | otherwise                     = go sprayA' sprayB' 1 1
+  where
+    n = max (numberOfVariables sprayA) (numberOfVariables sprayB)
+    coefficients :: Spray Rational -> [Rational]
+    coefficients spray = map (getCoefficient []) (sprayCoefficients spray)
+    cont :: Spray Rational -> Rational
+    cont spray = maximum $ coefficients spray
+    reduce :: Spray Rational -> Spray Rational
+    reduce spray = HM.map (/ cont spray) spray
+    a = cont sprayA
+    b = cont sprayB
+    d = max a b
+    sprayA' = HM.map (/ a) sprayA
+    sprayB' = HM.map (/ b) sprayB
+    go sprayA'' sprayB'' g h 
+      | sprayR == zeroSpray           = d *^ reduce sprayB''
+      | numberOfVariables sprayR == 0 = constantSpray d
+      | otherwise = go sprayB'' (HM.map (/ (g*h^delta)) sprayR) ellAq'' (g^delta / h^(delta-1))
+        where
+          (_, (_, sprayR)) = pseudoDivision sprayA'' sprayB''
+          (degA'', ellA'') = degreeAndLeadingCoefficient sprayA''
+          ellAq'' = getCoefficient [] ellA''
+          delta = degA'' - degree sprayB''

tests/Main.hs view

@@ -29,7 +29,8 @@                                                   resultant,
                                                   subresultants,
                                                   resultant1,
-                                                  subresultants1
+                                                  subresultants1,
+                                                  gcdQX
                                                 )
 import           Test.Tasty                     ( defaultMain
                                                 , testGroup
@@ -164,15 +165,33 @@         p' = permuteVariables p [3, 2, 1]
       assertEqual "" p' (swapVariables p (1, 3)),
 
-    testCase "resultant" $ do
+    testCase "resultant w.r.t x" $ 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)),
+      assertEqual "" r (y^**^4 ^-^ (8*^ y^**^6) ^+^ (16*^ y^**^8)),
 
+    testCase "resultant w.r.t y" $ 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 2 p q
+      assertEqual "" r (16*^x^**^8 ^-^ 32*^x^**^7 ^+^ 24*^x^**^6 ^-^ 8*^x^**^5 ^+^ x^**^4),
+
+    testCase "resultant product rule" $ do
+      let
+        x = lone 1 :: Spray Rational
+        y = lone 2 :: Spray Rational
+        f = x^**^4 ^-^ x^**^3 ^+^ x^**^2 ^-^ 2*^ (x ^*^ y^**^2) ^+^ y^**^4 
+        g = x ^-^ (2*^ y^**^2)
+        h = x^**^2 ^*^ y ^+^ y^**^3 ^+^ unitSpray
+      assertEqual "" (resultant 1 (f^*^g) h) (resultant 1 f h ^*^ resultant 1 g h),
+
     testCase "subresultants" $ do
       let
         x = lone 1 :: Spray Rational
@@ -189,11 +208,47 @@         q = x^**^2 ^-^ 3*^x ^+^ constantSpray 2 
       assertEqual "" (resultant1 p q) (0%1),
 
+    testCase "resultant1 product rule" $ do
+      let
+        x = lone 1 :: Spray Rational
+        f = x^**^2 ^-^ 5*^x ^+^ constantSpray 6 
+        g = x^**^2 ^-^ 3*^x ^+^ constantSpray 2 
+        h = x^**^3 ^+^ x ^-^ constantSpray 3
+      assertEqual "" (resultant1 (f^*^g) h) (resultant1 f h * resultant1 g h),
+
     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]
+      assertEqual "" (subresultants1 p q) [0%1, 2%1, 1%1],
+
+    testCase "resultant agrees with resultant1 for univariate case" $ do
+      let
+        x = lone 1 :: Spray Rational
+        f = x^**^4 ^-^ x^**^3 ^+^ x^**^2 ^-^ 2*^x 
+        g = x ^-^ (2*^ x^**^2) ^+^ constantSpray 4
+        r = resultant 1 f g
+        r1 = resultant1 f g
+      assertEqual "" r1 (getCoefficient [] r),
+
+    testCase "gcdQX" $ do
+      let
+        x = lone 1 :: Spray Rational
+        sprayD = x^**^2 ^+^ unitSpray
+        sprayA = sprayD ^*^ (x^**^4 ^-^ x) 
+        sprayB = sprayD ^*^ (2*^x ^+^ unitSpray)
+        sprayGCD = gcdQX sprayA sprayB
+      assertEqual "" sprayGCD (2 *^ sprayD),
+
+    testCase "gcdQX with a constant spray" $ do
+      let
+        x = lone 1 :: Spray Rational
+        sprayA = 3 *^ x^**^4 ^-^ x 
+        b1 = 2 :: Rational
+        b2 = 4 :: Rational
+        sprayB1 = constantSpray b1
+        sprayB2 = constantSpray b2
+      assertBool "" (gcdQX sprayA sprayB1 == constantSpray 3 && gcdQX sprayA sprayB2 == constantSpray 4) 
 
   ]