diff --git a/CHANGELOG.md b/CHANGELOG.md --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -75,4 +75,17 @@ * 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. +* New function `gcdQX`, to compute the greatest common divisor of two univariate sprays with rational coefficients. + + +## 0.2.4.0 - 2024-03-30 + +* Flipped the order of the arguments in `permuteVariables` and `swapVariables`. + +* New function `gcdSpray`, to compute the greatest common divisor of two sprays with coefficients in a field. + +* The function `gcdQX` has been removed since `gcdSpray` is more general. + +* The function `sprayDivision` has been renamed to `sprayDivisionRemainder`. + +* New function `sprayDivision`, returning the quotient and the remainder of the division of two sprays. diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -5,7 +5,7 @@ [![Stack-nightly](https://github.com/stla/hspray/actions/workflows/Stack-nightly.yml/badge.svg)](https://github.com/stla/hspray/actions/workflows/Stack-nightly.yml) -Simple multivariate polynomials in Haskell. +*Simple multivariate polynomials in Haskell.* ___ @@ -47,7 +47,7 @@ -- "((1.0) * a^(2)) * X^(0, 2) + ((2.0) * a^(2)) * X^(1, 1) + ((1.0) * a^(2)) * X^(2)" ``` -Evaluation: +#### Evaluation: ```haskell import Math.Algebra.Hspray @@ -60,7 +60,7 @@ -- 8.0 ``` -Partial evaluation: +#### Partial evaluation: ```haskell import Math.Algebra.Hspray @@ -78,7 +78,7 @@ -- "(6 % 1) + (1 % 1) x2" ``` -Differentiation: +#### Differentiation: ```haskell import Math.Algebra.Hspray @@ -174,3 +174,9 @@ map toList $ map snd l -- [[([0,1],2 % 3)],[([1],1 % 1)],[([1],1 % 1)]] ``` + + +## Other features + +Resultant and subresultants of two polynomials, and greatest common divisor of two polynomials +with coefficients in a field. diff --git a/benchmarks/Main.hs b/benchmarks/Main.hs --- a/benchmarks/Main.hs +++ b/benchmarks/Main.hs @@ -4,12 +4,27 @@ fibo :: Int -> Integer fibo n = if n < 2 then toInteger n else fibo (n - 1) + fibo (n - 2) +combn2 :: Int -> [(Int, (Int, Int))] +combn2 n = zip range (zip row1 row2) + where + range = [0 .. n-2] + row1 = concat $ scanl1 (++) (map (: []) range) + row2 = concatMap (\i -> replicate i i) [1 .. n-1] + +combn2' :: Int -> [(Int, (Int, Int))] -- winner +combn2' n = zip [0 .. n-2] (zip row1 row2) + where + range = [1 .. n-1] + row1 = concatMap (\i -> [0 .. i-1]) range + row2 = concatMap (\i -> replicate i i) range + main :: IO () main = defaultMain - [ bgroup "Fibonacci numbers" - [ bench "fifth" $ nf fibo 5 - , bench "tenth" $ nf fibo 10 - , bench "twentieth" $ nf fibo 20 + [ bgroup "combn2" + [ bench "combn2 - 20" $ nf combn2 20 + , bench "combn2' - 20" $ nf combn2' 20 + , bench "combn2 - 100" $ nf combn2 100 + , bench "combn2' - 100" $ nf combn2' 100 ] ] diff --git a/hspray.cabal b/hspray.cabal --- a/hspray.cabal +++ b/hspray.cabal @@ -1,7 +1,7 @@ name: hspray -version: 0.2.3.0 +version: 0.2.4.0 synopsis: Multivariate polynomials. -description: Manipulation of multivariate polynomials on a ring, Gröbner basis, resultant and subresultants. +description: Manipulation of multivariate polynomials on a commutative ring, Gröbner basis, resultant and subresultants, and greatest common divisor. homepage: https://github.com/stla/hspray#readme license: GPL-3 license-file: LICENSE diff --git a/src/Math/Algebra/Hspray.hs b/src/Math/Algebra/Hspray.hs --- a/src/Math/Algebra/Hspray.hs +++ b/src/Math/Algebra/Hspray.hs @@ -50,6 +50,7 @@ , swapVariables -- * Division of a spray , sprayDivision + , sprayDivisionRemainder -- * Gröbner basis , groebner , reduceGroebnerBasis @@ -62,7 +63,7 @@ , subresultants , subresultants1 -- * Greatest common divisor - , gcdQX + , gcdSpray -- * Miscellaneous , fromList , toList @@ -156,7 +157,8 @@ instance Eq Powers where (==) :: Powers -> Powers -> Bool pows1 == pows2 = exponents pows1' == exponents pows2' - where (pows1', pows2') = harmonize (pows1, pows2) + where + (pows1', pows2') = harmonize (pows1, pows2) instance Hashable Powers where hashWithSalt :: Int -> Powers -> Int @@ -167,16 +169,22 @@ type Monomial a = (Powers, a) instance (AlgAdd.C a, Eq a) => AlgAdd.C (Spray a) where - p + q = addSprays p q + (+) :: 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) => AlgRing.C (Spray a) where + (*) :: Spray a -> Spray a -> Spray a p * q = multSprays p q - one = lone 0 + one :: Spray a + one = lone 0 {- instance (AlgRing.C a, Eq a) => Num (Spray a) where p + q = addSprays p q @@ -201,7 +209,9 @@ -- | Power of a spray (^**^) :: (AlgRing.C a, Eq a) => Spray a -> Int -> Spray a -(^**^) p n = AlgRing.product (replicate n p) +(^**^) p n = if n >= 0 + then AlgRing.product (replicate n p) + else error "(^**^): negative power of a spray is not allowed." -- | Scale spray by a scalar (*^) :: (AlgRing.C a, Eq a) => a -> Spray a -> Spray a @@ -218,7 +228,8 @@ -- | drop trailing zeros simplifyPowers :: Powers -> Powers simplifyPowers pows = Powers s (S.length s) - where s = dropWhileR (== 0) (exponents pows) + where + s = dropWhileR (== 0) (exponents pows) -- | drop trailing zeros in the powers of a spray simplifySpray :: Spray a -> Spray a @@ -231,7 +242,8 @@ -- | addition of two sprays addSprays :: (AlgAdd.C a, Eq a) => Spray a -> Spray a -> Spray a addSprays p q = cleanSpray $ HM.foldlWithKey' f p q - where f s powers coef = HM.insertWith (AlgAdd.+) powers coef s + where + f s powers coef = HM.insertWith (AlgAdd.+) powers coef s -- | opposite spray negateSpray :: AlgAdd.C a => Spray a -> Spray a @@ -258,8 +270,8 @@ derivSpray :: (AlgRing.C a, Eq a) => Int -- ^ index of the variable of differentiation (starting at 1) - -> Spray a -- ^ the spray - -> Spray a + -> Spray a -- ^ the spray to be derivated + -> Spray a -- ^ the derivated spray derivSpray i p = if i >= 1 then cleanSpray $ HM.fromListWith (AlgAdd.+) monomials else error "derivSpray: invalid index." @@ -333,7 +345,7 @@ getCoefficient expnts spray = fromMaybe AlgAdd.zero (HM.lookup powers spray) where expnts' = S.dropWhileR (== 0) (S.fromList expnts) - powers = Powers expnts' (S.length expnts') + powers = Powers expnts' (S.length expnts') -- | number of variables in a spray numberOfVariables :: Spray a -> Int @@ -346,7 +358,8 @@ evalMonomial :: AlgRing.C a => [a] -> Monomial a -> a evalMonomial xyz (powers, coeff) = coeff AlgRing.* AlgRing.product (zipWith (AlgRing.^) xyz pows) - where pows = DF.toList (fromIntegral <$> exponents powers) + where + pows = DF.toList (fromIntegral <$> exponents powers) -- | Evaluates a spray -- @@ -368,14 +381,14 @@ substituteMonomial :: AlgRing.C a => [Maybe a] -> Monomial a -> Monomial a substituteMonomial subs (powers, coeff) = (powers'', coeff') where - pows = exponents powers - n = nvariables powers - indices = findIndices isJust (take n 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 + pows = exponents powers + n = nvariables powers + indices = findIndices isJust (take n 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'' = simplifyPowers $ Powers pows'' n -- | Substitutes some variables in a spray @@ -392,9 +405,10 @@ then spray' else error "substituteSpray: incorrect length of the substitutions list." where - n = numberOfVariables spray + 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) @@ -410,13 +424,13 @@ -- >>> q = composeSpray p [z, x ^+^ y ^+^ z] -- >>> putStrLn $ prettySprayXYZ q -- (1) X + (1) Y + (2) Z -composeSpray :: (AlgRing.C a, Eq a) => Spray a -> [Spray a] -> Spray a +composeSpray :: forall a. (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 :: Spray a -> Spray (Spray a) identify = HM.map constantSpray --- | Creates a spray from list of terms +-- | Creates a spray from a 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 @@ -430,30 +444,31 @@ -- >>> x3 = lone 3 :: Spray Rational -- >>> p = f x1 x2 x3 -- --- prop> permuteVariables p [3, 1, 2] == f x3 x1 x2 -permuteVariables :: Spray a -> [Int] -> Spray a -permuteVariables spray permutation = +-- prop> permuteVariables [3, 1, 2] p == f x3 x1 x2 +permuteVariables :: [Int] -> Spray a -> Spray a +permuteVariables permutation spray = if n' >= n && isPermutation permutation then spray' else error "permuteVariables: invalid permutation." where - n = numberOfVariables spray + n = numberOfVariables spray n' = maximum permutation isPermutation pmtn = minimum pmtn == 1 && length (nub pmtn) == n' - intmap = IM.fromList (zip permutation [1 .. 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 + permuteSeq x = + S.mapWithIndex (\i _ -> x `index` (invpermutation !! i - 1)) x (powers, coeffs) = unzip (HM.toList spray) - expnts = map exponents powers + expnts = map exponents powers expnts' = map (permuteSeq . growSequence' n') expnts powers' = map (\exps -> simplifyPowers (Powers exps n')) expnts' - spray' = HM.fromList (zip powers' coeffs) + spray' = HM.fromList (zip powers' coeffs) -- | Swaps two variables of a spray -- --- prop> swapVariables p (1, 3) == permuteVariables p [3, 2, 1] -swapVariables :: Spray a -> (Int, Int) -> Spray a -swapVariables spray (i, j) = +-- prop> swapVariables (1, 3) p == permuteVariables [3, 2, 1] p +swapVariables :: (Int, Int) -> Spray a -> Spray a +swapVariables (i, j) spray = if i>=1 && j>=1 then spray' else error "swapVariables: invalid indices." @@ -463,12 +478,13 @@ | k == j = i | otherwise = k transposition = map f [1 .. n] - permuteSeq x = S.mapWithIndex (\ii _ -> x `index` (transposition !! ii - 1)) x + permuteSeq x = + S.mapWithIndex (\ii _ -> x `index` (transposition !! ii - 1)) x (powers, coeffs) = unzip (HM.toList spray) - expnts = map exponents powers + expnts = map exponents powers expnts' = map (permuteSeq . growSequence' n) expnts powers' = map (\exps -> simplifyPowers (Powers exps n)) expnts' - spray' = HM.fromList (zip powers' coeffs) + spray' = HM.fromList (zip powers' coeffs) -- pretty stuff --------------------------------------------------------------- @@ -495,8 +511,9 @@ -> String prettySpray prettyCoef var p = unpack $ intercalate (pack " + ") stringTerms where - stringTerms = map stringTerm (sortBy (flip compare `on` fexpts) (HM.toList p)) - fexpts term = exponents $ fst term + stringTerms = + map stringTerm (sortBy (flip compare `on` fexpts) (HM.toList p)) + fexpts term = exponents $ fst term stringTerm term = append (snoc (snoc (cons '(' $ snoc stringCoef ')') ' ') '*') (prettyPowers var pows) @@ -510,8 +527,8 @@ where n = S.length pows f i p - | p == 0 = "" - | p == 1 = "x" ++ show i + | p == 0 = "" + | p == 1 = "x" ++ show i | otherwise = "x" ++ show i ++ "^" ++ show p x1x2x3 = concatMap (\i -> f i (pows `index` (i-1))) [1 .. n] @@ -526,14 +543,15 @@ prettySpray' :: (Show a) => Spray a -> String prettySpray' spray = unpack $ intercalate (pack " + ") terms where - terms = map stringTerm (sortBy (flip compare `on` fexpts) (HM.toList spray)) - fexpts term = exponents $ fst term + terms = map stringTerm + (sortBy (flip compare `on` fexpts) (HM.toList spray)) + fexpts term = exponents $ fst term stringTerm term = append stringCoef'' (prettyPowers' pows) where - pows = exponents (fst term) - constant = S.length pows == 0 - stringCoef = pack $ show (snd term) - stringCoef' = cons '(' $ snoc stringCoef ')' + pows = exponents (fst term) + constant = S.length pows == 0 + stringCoef = pack $ show (snd term) + stringCoef' = cons '(' $ snoc stringCoef ')' stringCoef'' = if constant then stringCoef' else snoc stringCoef' ' ' -- | prettyPowersXYZ [1, 2, 1] = XY^2Z @@ -542,15 +560,15 @@ then pack xyz else error "there is more than three variables" where - n = S.length pows + n = S.length pows gpows = growSequence pows n 3 f letter p - | p == 0 = "" - | p == 1 = letter + | p == 0 = "" + | p == 1 = letter | otherwise = letter ++ "^" ++ show p - x = f "X" (gpows `index` 0) - y = f "Y" (gpows `index` 1) - z = f "Z" (gpows `index` 2) + x = f "X" (gpows `index` 0) + y = f "Y" (gpows `index` 1) + z = f "Z" (gpows `index` 2) xyz = x ++ y ++ z -- | Pretty form of a spray having at more three variables @@ -599,14 +617,14 @@ -- | index of the maximum of a list maxIndex :: Ord a => [a] -> Int -maxIndex = fst . maximumBy (comparing snd) . zip [0..] +maxIndex = fst . maximumBy (comparing snd) . zip [0 .. ] -- | Leading term of a spray leadingTerm :: Spray a -> Monomial a leadingTerm p = (biggest, p HM.! biggest) where - powers = HM.keys p - i = maxIndex $ map exponents powers + powers = HM.keys p + i = maxIndex $ map exponents powers biggest = powers !! i -- | whether a monomial divides another monomial @@ -615,26 +633,26 @@ where expntsP = exponents powsP expntsQ = exponents powsQ - lower = DF.all (\(x, y) -> x <= y) (S.zip expntsP expntsQ) + lower = DF.all (uncurry (<=)) (S.zip expntsP expntsQ) -- | quotient of monomial Q by monomial p, assuming P divides Q quotient :: AlgField.C a => Monomial a -> Monomial a -> Monomial a quotient (powsQ, coeffQ) (powsP, coeffP) = (pows, coeff) where (powsP', powsQ') = harmonize (powsP, powsQ) - expntsP = exponents powsP' - expntsQ = exponents powsQ' - expnts = S.zipWith (-) expntsQ expntsP - n = nvariables powsP' - pows = Powers expnts n - coeff = coeffQ AlgField./ coeffP + expntsP = exponents powsP' + expntsQ = exponents powsQ' + expnts = S.zipWith (-) expntsQ expntsP + n = nvariables powsP' + pows = Powers expnts n + coeff = coeffQ AlgField./ coeffP -- | Remainder of the division of a spray by a list of divisors, -- using the lexicographic ordering of the monomials -sprayDivision :: forall a. (Eq a, AlgField.C a) => Spray a -> [Spray a] -> Spray a -sprayDivision p qs = +sprayDivisionRemainder :: forall a. (Eq a, AlgField.C a) => Spray a -> [Spray a] -> Spray a +sprayDivisionRemainder p qs = if n == 0 - then error "sprayDivision: the list of divisors is empty." + then error "sprayDivisionRemainder: the list of divisors is empty." else snd $ ogo p AlgAdd.zero where n = length qs @@ -646,27 +664,59 @@ go :: Monomial 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 + | i == n = g lts s r + | otherwise = go lts news r (i+1) newdivoccured where - (q, ltq) = qsltqs !! i + (q, ltq) = qsltqs !! i newdivoccured = divides ltq lts - news = if newdivoccured + news = if newdivoccured then s ^-^ (fromMonomial (quotient lts ltq) ^*^ q) else s ogo :: Spray a -> Spray a -> (Spray a, Spray a) ogo !s !r | s == AlgAdd.zero = (s, r) - | otherwise = ogo s' r' + | otherwise = ogo s' r' where (s', r') = go (leadingTerm s) s r 0 False +-- | Division of a spray by a spray +sprayDivision :: forall a. (Eq a, AlgField.C a) + => Spray a -- ^ dividend + -> Spray a -- ^ divisor + -> (Spray a, Spray a) -- ^ (quotient, remainder) +sprayDivision sprayA sprayB = + if sprayB == AlgAdd.zero + then error "sprayDivision: division by zero." + else ogo sprayA AlgAdd.zero AlgAdd.zero + where + go :: Monomial a -> Spray a -> Spray a -> Spray a -> Int -> Bool -> (Spray a, Spray a, Spray a) + go ltp !p !q r !i !divoccured + | divoccured = (p, q, r) + | i == 1 = (p ^-^ ltpspray, q, r ^+^ ltpspray) + | otherwise = go ltp newp newq r 1 newdivoccured + where + ltpspray = fromMonomial ltp + ltB = leadingTerm sprayB + newdivoccured = divides ltB ltp + (newp, newq) = if newdivoccured + then (p ^-^ (qtnt ^*^ sprayB), q ^+^ qtnt) + else (p, q) + where + qtnt = fromMonomial $ quotient ltp ltB + ogo :: Spray a -> Spray a -> Spray a -> (Spray a, Spray a) + ogo !p !q !r + | p == AlgAdd.zero = (q, r) + | otherwise = ogo p' q' r' + where + (p', q', r') = go (leadingTerm p) p q r 0 False + -- Groebner stuff ------------------------------------------------------------- --- | slight modification of `sprayDivision` to speed up groebner00 -sprayDivision' :: forall a. (Eq a, AlgField.C a) => Spray a -> HashMap Int (Spray a, Monomial a) -> Spray a -sprayDivision' p qsltqs = snd $ ogo p AlgAdd.zero +-- | slight modification of `sprayDivisionRemainder` to speed up groebner00 +sprayDivisionRemainder' :: forall a. (Eq a, AlgField.C a) + => Spray a -> HashMap Int (Spray a, Monomial a) -> Spray a +sprayDivisionRemainder' p qsltqs = snd $ ogo p AlgAdd.zero where n = HM.size qsltqs g :: Monomial a -> Spray a -> Spray a -> (Spray a, Spray a) @@ -676,10 +726,10 @@ go :: Monomial 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 + | i == n = g lts s r + | otherwise = go lts news r (i+1) newdivoccured where - (q, ltq) = qsltqs HM.! i + (q, ltq) = qsltqs HM.! i newdivoccured = divides ltq lts news = if newdivoccured then s ^-^ (fromMonomial (quotient lts ltq) ^*^ q) @@ -687,30 +737,34 @@ ogo :: Spray a -> Spray a -> (Spray a, Spray a) ogo !s !r | s == AlgAdd.zero = (s, r) - | otherwise = ogo 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 .. n-2] (zip row1 row2)) +combn2 n s = HM.fromList (zip range0 (zip row1 row2)) where - row1 = drop s $ concatMap (\i -> [0 .. (i-1)]) [1 .. n-1] - row2 = drop s $ concatMap (\i -> replicate i i) [1 .. n-1] + range0 = [0 .. n-2] + 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, Monomial a) -> (Spray a, Monomial a) -> Spray a sPolynomial pltp qltq = wp ^*^ p ^-^ wq ^*^ q where - p = fst pltp - q = fst qltq - (lpowsP, lcoefP) = snd pltp - (lpowsQ, lcoefQ) = snd qltq + p = fst pltp + q = fst qltq + (lpowsP, lcoefP) = snd pltp + (lpowsQ, lcoefQ) = snd qltq (lpowsP', lpowsQ') = harmonize (lpowsP, lpowsQ) - lexpntsP = exponents lpowsP' - lexpntsQ = exponents lpowsQ' + lexpntsP = exponents lpowsP' + lexpntsQ = exponents lpowsQ' gamma = S.zipWith max lexpntsP lexpntsQ betaP = S.zipWith (-) gamma lexpntsP betaQ = S.zipWith (-) gamma lexpntsQ - n = nvariables lpowsP' + n = nvariables lpowsP' wp = fromMonomial (Powers betaP n, AlgField.recip lcoefP) wq = fromMonomial (Powers betaQ n, AlgField.recip lcoefQ) @@ -718,19 +772,19 @@ groebner00 :: forall a. (Eq a, AlgField.C a) => [Spray a] -> [Spray a] groebner00 sprays = go 0 j0 combins0 spraysMap where - j0 = length sprays + j0 = length sprays combins0 = combn2 j0 0 - ltsprays = map leadingTerm sprays + ltsprays = map leadingTerm sprays spraysltsprays = zip sprays ltsprays - spraysMap = HM.fromList (zip [0 .. j0-1] spraysltsprays) + spraysMap = HM.fromList (zip [0 .. j0-1] spraysltsprays) go :: Int -> Int -> HashMap Int (Int, Int) -> HashMap Int (Spray a, Monomial a) -> [Spray a] go !i !j !combins !gpolysMap | i == length combins = map fst (HM.elems gpolysMap) - | otherwise = go i' j' combins' gpolysMap' + | otherwise = go i' j' combins' gpolysMap' where - (k, l) = combins HM.! i - sfg = sPolynomial (gpolysMap HM.! k) (gpolysMap HM.! l) - sbarfg = sprayDivision' sfg gpolysMap + (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 sbarfg == AlgAdd.zero then @@ -747,20 +801,20 @@ groebner0 sprays = if n <= 1 then sprays else [basis00 !! k | k <- [0 .. n-1] \\ discard] where - n = length basis00 + n = length basis00 basis00 = groebner00 sprays go :: Int -> [Int] -> [Int] go !i toRemove - | i == n = toRemove + | i == n = toRemove | otherwise = go (i+1) toRemove' where - ltf = leadingTerm (basis00 !! i) + ltf = leadingTerm (basis00 !! i) toDrop = toRemove ++ [i] igo :: Int -> Bool igo !j - | j == n = False + | j == n = False | j `elem` toDrop = igo (j+1) - | otherwise = ok || igo (j+1) + | otherwise = ok || igo (j+1) where ok = divides (leadingTerm (basis00 !! j)) ltf toRemove' = if igo 0 then toDrop else toRemove @@ -769,22 +823,24 @@ -- | Reduces a Groebner 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 + if length gbasis >= 2 + then map reduction [0 .. n-1] + else ngbasis where normalize :: Spray a -> Spray a normalize spray = AlgField.recip coef *^ spray where (_, coef) = leadingTerm spray ngbasis = map normalize gbasis - n = length ngbasis + n = length ngbasis reduction :: Int -> Spray a - reduction i = sprayDivision (ngbasis !! i) rest + reduction i = sprayDivisionRemainder (ngbasis !! i) rest where rest = [ngbasis !! k | k <- [0 .. n-1] \\ [i]] -- | Groebner basis (always minimal and possibly reduced) -- --- prop> groebner ps True = reduceGroebnerBasis (groebner ps False) +-- prop> groebner ps True == reduceGroebnerBasis (groebner ps False) groebner :: forall a. (Eq a, AlgField.C a) => [Spray a] -- ^ list of sprays @@ -816,12 +872,12 @@ | otherwise = map (q:) cs ++ run (n-1) i qs (drop dc cs) where f :: [a] -> (a, [a]) - f [] = error "should not happen" + f [] = error "combinationsOf: should not happen." f (b:bs) = (b, bs) - (q, qs) = f (take (n-i+1) ys) - dc = product [(n-k+1) .. (n-1)] `div` product [1 .. i-1] + (q, qs) = f (take (n-i+1) ys) + dc = product [(n-k+1) .. (n-1)] `div` product [1 .. i-1] --- | generate all permutations of a binary sequence +-- | generates all permutations of a binary sequence permutationsBinarySequence :: Int -> Int -> [Seq Int] permutationsBinarySequence nzeros nones = let n = nzeros + nones in @@ -844,8 +900,8 @@ -> Spray a esPolynomial n k | k <= 0 || n <= 0 = error "esPolynomial: both arguments must be positive integers." - | k > n = AlgAdd.zero - | otherwise = simplifySpray spray + | k > n = AlgAdd.zero + | otherwise = simplifySpray spray where perms = permutationsBinarySequence (n-k) k spray = HM.fromList $ map (\expts -> (Powers expts n, AlgRing.one)) perms @@ -856,15 +912,13 @@ where n = numberOfVariables spray indices = [1 .. n] - esPolys = map (\i -> esPolynomial n i :: Spray a) indices - yPolys = map (\i -> lone (n + i) :: Spray a) indices - gPolys = zipWith (^-^) esPolys yPolys - gbasis = groebner0 gPolys - g = sprayDivision spray gbasis + gPolys = map (\i -> esPolynomial n i ^-^ lone (n + i)) indices + gbasis = groebner0 gPolys + g = sprayDivisionRemainder spray gbasis gpowers = HM.keys g - check1 = minimum (map nvariables gpowers) > n - expnts = map exponents gpowers - check2 = DF.all (DF.all (0 ==)) (map (S.take n) expnts) + check1 = minimum (map nvariables gpowers) > n + expnts = map exponents gpowers + check2 = DF.all (DF.all (0 ==)) (map (S.take n) expnts) -- | Whether a spray can be written as a polynomial of a given list of sprays -- (the sprays in the list must belong to the same polynomial ring as the spray); @@ -883,21 +937,21 @@ n = numberOfVariables spray n' = maximum $ map numberOfVariables sprays result - | n > n' = (False, Nothing) - | n < n' = error "not enough variables in the spray" + | n > n' = (False, Nothing) + | n < n' = error "not enough variables in the spray" | otherwise = (checks, poly) where - m = length sprays - yPolys = map (\i -> lone (n + i) :: Spray a) [1 .. m] - gPolys = zipWith (^-^) sprays yPolys + m = length sprays + yPolys = map (\i -> lone (n + i) :: Spray a) [1 .. m] + gPolys = zipWith (^-^) sprays yPolys gbasis0 = groebner0 gPolys - g = sprayDivision spray gbasis0 + g = sprayDivisionRemainder spray gbasis0 gpowers = HM.keys g - check1 = minimum (map nvariables gpowers) > n - expnts = map exponents gpowers - check2 = DF.all (DF.all (0 ==)) (map (S.take n) expnts) - checks = check1 && check2 - poly = if checks + check1 = minimum (map nvariables gpowers) > n + expnts = map exponents gpowers + check2 = DF.all (DF.all (0 ==)) (map (S.take n) expnts) + checks = check1 && check2 + poly = if checks then Just $ dropXis g else Nothing dropXis = HM.mapKeys f @@ -912,28 +966,34 @@ 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]] + 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) + then fromLists [[AlgRing.one]] -- plays the role of the empty matrix: 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]] + 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]] + else suml1 [negateIf i (times (m DM.! (i,1)) (detLaplace (minorMatrix i 1 m))) + | i <- [1 .. nrows m]] where - suml1 = foldl1' (AlgAdd.+) + 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 @@ -947,18 +1007,19 @@ where n = numberOfVariables spray (powers, coeffs) = unzip (HM.toList spray) - expnts = map exponents powers + 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 + (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 permutation = [2 .. n] ++ [1] sprays = [ - permuteVariables (fromMaybe AlgAdd.zero (IM.lookup i imap')) permutation + permuteVariables permutation (fromMaybe AlgAdd.zero (IM.lookup i imap')) | i <- [0 .. maximum xpows] ] @@ -974,12 +1035,18 @@ 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 + pcoeffs = 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 + else [fromMaybe AlgAdd.zero (HM.lookup (Powers (S.singleton i) 1) p) + | i <- [maxp, maxp-1 .. 1]] ++ [p0] + where + maxp = maximum pexpnts + qcoeffs = if null qexpnts then [q0] - else q0 : [fromMaybe AlgAdd.zero (HM.lookup (Powers (S.singleton i) 1) q) | i <- [1 .. maximum qexpnts]] + else [fromMaybe AlgAdd.zero (HM.lookup (Powers (S.singleton i) 1) q) + | i <- [maxq, maxq-1 .. 1]] ++ [q0] + where + maxq = maximum qexpnts -- | Subresultants of two /univariate/ sprays subresultants1 :: (Eq a, AlgRing.C a) => Spray a -> Spray a -> [a] @@ -992,12 +1059,18 @@ 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 + pcoeffs = 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 + else [fromMaybe AlgAdd.zero (HM.lookup (Powers (S.singleton i) 1) p) + | i <- [maxp, maxp-1 .. 1]] ++ [p0] + where + maxp = maximum pexpnts + qcoeffs = if null qexpnts then [q0] - else q0 : [fromMaybe AlgAdd.zero (HM.lookup (Powers (S.singleton i) 1) q) | i <- [1 .. maximum qexpnts]] + else [fromMaybe AlgAdd.zero (HM.lookup (Powers (S.singleton i) 1) q) + | i <- [maxq, maxq-1 .. 1]] ++ [q0] + where + maxq = maximum qexpnts d = length pcoeffs e = length qcoeffs @@ -1009,15 +1082,16 @@ -> Spray a resultant var p q = if var >= 1 && var <= n - then permuteVariables det permutation' + then permuteVariables permutation' det else error "resultant: invalid variable index." where n = max (numberOfVariables p) (numberOfVariables q) - permutation = var : [1 .. var-1] ++ [var+1 .. n] + 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') + p' = permuteVariables permutation p + q' = permuteVariables permutation q + det = detLaplace $ + sylvesterMatrix (sprayCoefficients p') (sprayCoefficients q') -- | Subresultants of two sprays subresultants :: (Eq a, AlgRing.C a) @@ -1028,7 +1102,8 @@ 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] + | otherwise = map (permute' . detLaplace . sylvesterMatrix' pcoeffs qcoeffs) + [0 .. min d e - 1] where pcoeffs = sprayCoefficients p' qcoeffs = sprayCoefficients q' @@ -1036,92 +1111,158 @@ 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 + permute = permuteVariables permutation + p' = permute p + q' = permute q permutation' = [2 .. var] ++ (1 : [var+1 .. n]) - permute spray = permuteVariables spray permutation' + permute' = permuteVariables permutation' --- GCD ------------------------------------------------------------------------ +-- GCD stuff ------------------------------------------------------------------ --- 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 +-- the coefficients of a spray as a univariate spray in x_n with spray coefficients +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 - expnts = map exponents $ HM.keys spray - expnts' = filter (not . S.null) expnts + permutation = [2 .. n] ++ [1] + spray' = permuteVariables permutation 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' + 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 + deg = maximum xpows + sprays = [ + fromMaybe AlgAdd.zero (IM.lookup i imap') + | i <- [deg, deg-1 .. 0] + ] --- 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) +-- the degree of a spray as a univariate spray in x_n with spray coefficients +degree :: (Eq a, AlgAdd.C a) => Int -> Spray a -> Int +degree n spray + | numberOfVariables spray == 0 = + if spray == zeroSpray then minBound else 0 + | numberOfVariables spray /= n = 0 + | otherwise = maximum xpows + where + permutation = [2 .. n] ++ [1] + spray' = permuteVariables permutation spray + 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_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 + else 0, + constantSpray constantTerm + ) + | numberOfVariables spray /= n = (0, spray) + | otherwise = (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) + permutation = [2 .. n] ++ [1] + spray' = permuteVariables permutation 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] + deg = maximum xpows + is = elemIndices deg xpows + expnts'' = [S.deleteAt 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 +-- 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 sprayA sprayB = (ellB ^**^ delta , go sprayA zeroSpray delta) +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 sprayA - (degB, ellB) = degreeAndLeadingCoefficient sprayB - delta = degA - degB + 1 + degA = degree n sprayA + (degB, ellB) = degreeAndLeadingCoefficient n sprayB + delta = degA - degB + 1 go sprayR sprayQ e = - if degR < degB + if degR < degB || sprayR == zeroSpray then (q ^*^ sprayQ, q ^*^ sprayR) - else go (ellB ^*^ sprayR ^-^ sprayS ^*^ sprayB) (ellB ^*^ sprayQ ^+^ sprayS) (e - 1) + 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) + (degR, ellR) = degreeAndLeadingCoefficient n sprayR + q = ellB ^**^ e + sprayXn = lone n + sprayS = ellR ^*^ sprayXn ^**^ (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 +-- 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 + | sprayB == zeroSpray = sprayA + | otherwise = go sprayA' sprayB' unitSpray unitSpray 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 + 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 snd division == zeroSpray + then fst division + else error "exactDivisionBy: should not happen." + where + division = sprayDivision a b + reduceSpray :: Spray a -> Spray a + reduceSpray spray = exactDivisionBy cntnt spray + where + coeffs = sprayCoefficients' n' spray + cntnt = foldl1' gcdKX1dotsXm coeffs + 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 - | 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)) + | sprayR == zeroSpray = d ^*^ reduceSpray sprayB'' + | numberOfVariables sprayR == 0 = d + | otherwise = go sprayB'' + (exactDivisionBy (g ^*^ h^**^delta) sprayR) + ellA'' + (exactDivisionBy (h^**^delta) (h ^*^ g^**^delta)) where - (_, (_, sprayR)) = pseudoDivision sprayA'' sprayB'' - (degA'', ellA'') = degreeAndLeadingCoefficient sprayA'' - ellAq'' = getCoefficient [] ellA'' - delta = degA'' - degree sprayB'' + (_, (_, sprayR)) = pseudoDivision n' sprayA'' sprayB'' + (degA'', ellA'') = degreeAndLeadingCoefficient n' sprayA'' + degB'' = degree 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) diff --git a/tests/Main.hs b/tests/Main.hs --- a/tests/Main.hs +++ b/tests/Main.hs @@ -30,7 +30,8 @@ subresultants, resultant1, subresultants1, - gcdQX + sprayDivision, + gcdSpray ) import Test.Tasty ( defaultMain , testGroup @@ -153,7 +154,7 @@ x2 = lone 2 :: Spray Rational x3 = lone 3 :: Spray Rational p = f x1 x2 x3 - p' = permuteVariables p [3, 1, 2] + p' = permuteVariables [3, 1, 2] p assertEqual "" p' (f x3 x1 x2), testCase "swapVariables" $ do @@ -162,8 +163,8 @@ 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)), + p' = permuteVariables [3, 2, 1] p + assertEqual "" p' (swapVariables (1, 3) p), testCase "resultant w.r.t x" $ do let @@ -232,16 +233,16 @@ r1 = resultant1 f g assertEqual "" r1 (getCoefficient [] r), - testCase "gcdQX" $ do + testCase "gcdSpray - univariate example" $ 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), + sprayGCD = gcdSpray sprayA sprayB + assertEqual "" sprayGCD (9 *^ sprayD), - testCase "gcdQX with a constant spray" $ do + testCase "gcdSpray with a constant spray" $ do let x = lone 1 :: Spray Rational sprayA = 3 *^ x^**^4 ^-^ x @@ -249,6 +250,36 @@ b2 = 4 :: Rational sprayB1 = constantSpray b1 sprayB2 = constantSpray b2 - assertBool "" (gcdQX sprayA sprayB1 == constantSpray 3 && gcdQX sprayA sprayB2 == constantSpray 4) + assertBool "" (gcdSpray sprayA sprayB1 == unitSpray && gcdSpray sprayA sprayB2 == unitSpray), + + testCase "sprayDivision" $ do + let + x = lone 1 :: Spray Rational + y = lone 2 :: Spray Rational + sprayB = x^**^2 ^*^ y ^-^ x ^*^ y ^+^ constantSpray 3 + sprayQ = x^**^4 ^-^ x ^+^ y^**^2 + sprayA = sprayB ^*^ sprayQ + assertEqual "" (sprayDivision sprayA sprayB) (sprayQ, zeroSpray), + + testCase "gcdSpray - bivariate example" $ do + let + x = lone 1 :: Spray Rational + y = lone 2 :: Spray Rational + sprayD = x^**^2 ^*^ y ^-^ x ^*^ y ^+^ constantSpray 3 + sprayA = sprayD ^*^ (x^**^4 ^-^ x ^+^ y^**^2) + sprayB = sprayD ^*^ y ^*^ (2*^x ^+^ unitSpray) + g = gcdSpray sprayA sprayB + assertEqual "" g ((1%3) *^ sprayD), + + testCase "gcdSpray - trivariate example" $ do + let + x = lone 1 :: Spray Rational + y = lone 2 :: Spray Rational + z = lone 3 :: Spray Rational + sprayD = x^**^2 ^*^ y ^-^ x ^*^ y ^+^ z ^+^ constantSpray 3 + sprayA = sprayD^**^1 ^*^ (x^**^4 ^-^ x ^+^ y ^+^ x ^*^ y ^*^ z^**^2) + sprayB = sprayD^**^1 ^*^ y ^*^ (2*^x ^+^ unitSpray) ^*^ z + g = gcdSpray sprayA sprayB + assertEqual "" g sprayD ]