diff --git a/CHANGELOG.md b/CHANGELOG.md --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -32,3 +32,17 @@ * Completed the README to show how to deal with symbolic coefficients. + +## 0.2.0.0 - 2024-03-14 + +* New functions `prettySpray'` and `prettySprayXYZ`. + +* New function `substituteSpray`. + +* New function `sprayDivision`, to perform the division of a spray by a list of sprays. + +* New function `groebner`, to compute a Groebner basis of a list of sprays. + +* 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. diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -60,6 +60,24 @@ -- 8.0 ``` +Partial evaluation: + +```haskell +import Math.Algebra.Hspray +import Data.Ratio +x1 = lone 1 :: Spray Rational +x2 = lone 2 :: Spray Rational +x3 = lone 3 :: Spray Rational +poly = x1^**^2 ^+^ x2 ^+^ x3 ^-^ unitSpray +prettySpray' poly +-- "((-1) % 1) + (1 % 1) x3 + (1 % 1) x2 + (1 % 1) x1^2" +-- +-- substitute x1 -> 2 and x3 -> 3 +poly' = substituteSpray [Just 2, Nothing, Just 3] p +prettySpray' poly' +-- "(6 % 1) + (1 % 1) x2" +``` + Differentiation: ```haskell @@ -72,6 +90,34 @@ prettySpray show "X" $ derivSpray 1 poly -- "(2.0) * X^(0, 1, 1) + (6.0) * X^(1)" ``` + +## Grôbner bases + +As of version 2.0.0, it is possible to compute a Gröbner basis. + +```haskell +import Math.Algebra.Hspray +import Data.Ratio +-- define the elementary monomials +o = lone 0 :: Spray Rational +x = lone 1 :: Spray Rational +y = lone 2 :: Spray Rational +z = lone 3 :: Spray Rational +-- define three polynomials +p1 = x^**^2 ^+^ y ^+^ z ^-^ o -- X² + Y + Z - 1 +p2 = x ^+^ y^**^2 ^+^ z ^-^ o -- X + Y² + Z - 1 +p3 = x ^+^ y ^+^ z^**^2 ^-^ o -- X + Y + Z² - 1 +-- compute the reduced Gröbner basis +gbasis = groebner [p1, p2, p3] True +-- show result +prettyResult = map prettySprayXYZ gbasis +mapM_ print prettyResult +-- "((-1) % 1) + (1 % 1) Z^2 + (1 % 1) Y + (1 % 1) X" +-- "(1 % 1) Z + ((-1) % 1) Z^2 + ((-1) % 1) Y + (1 % 1) Y^2" +-- "((-1) % 2) Z^2 + (1 % 2) Z^4 + (1 % 1) YZ^2" +-- "((-1) % 1) Z^2 + (4 % 1) Z^3 + ((-4) % 1) Z^4 + (1 % 1) Z^6" +``` + ## Easier usage diff --git a/hspray.cabal b/hspray.cabal --- a/hspray.cabal +++ b/hspray.cabal @@ -1,7 +1,7 @@ name: hspray -version: 0.1.3.0 +version: 0.2.0.0 synopsis: Multivariate polynomials. -description: Manipulation of multivariate polynomials on a ring. +description: Manipulation of multivariate polynomials on a ring and Gröbner bases. homepage: https://github.com/stla/hspray#readme license: GPL-3 license-file: LICENSE @@ -18,15 +18,17 @@ hs-source-dirs: src exposed-modules: Math.Algebra.Hspray build-depends: base >= 4.7 && < 5 - , containers >= 0.6.4.1 - , hashable >= 1.3.4.0 - , unordered-containers >= 0.2.17.0 - , numeric-prelude >= 0.4.4 - , text >= 1.2.5.0 + , containers >= 0.6.4.1 && < 0.8 + , hashable >= 1.3.4.0 && < 1.5 + , unordered-containers >= 0.2.17.0 && < 0.3 + , numeric-prelude >= 0.4.4 && < 0.5 + , text >= 1.2.5.0 && < 2.2 other-extensions: FlexibleInstances , FlexibleContexts , MultiParamTypeClasses , InstanceSigs + , ScopedTypeVariables + , BangPatterns default-language: Haskell2010 ghc-options: -Wall -Wcompat @@ -42,6 +44,7 @@ type: exitcode-stdio-1.0 main-is: Main.hs hs-source-dirs: tests/ + other-modules: Approx Build-Depends: base >= 4.7 && < 5 , tasty , tasty-hunit 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 @@ -1,11 +1,23 @@ +{-| +Module : Math.Algebra.Hspray +Description : Multivariate polynomials on a ring. +Copyright : (c) Stéphane Laurent, 2023 +License : GPL-3 +Maintainer : laurent_step@outlook.fr + +Deals with multivariate polynomials on a ring. See README for examples. +-} +{-# LANGUAGE BangPatterns #-} {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE MultiParamTypeClasses #-} {-# LANGUAGE InstanceSigs #-} +{-# LANGUAGE ScopedTypeVariables #-} module Math.Algebra.Hspray ( Powers (..) , Spray + , Monomial , lone , unitSpray , constantSpray @@ -19,12 +31,23 @@ , (^*^) , (^**^) , evalSpray + , substituteSpray + , fromRationalSpray , prettySpray + , prettySpray' + , prettySprayXYZ , composeSpray , bombieriSpray , derivSpray + , leadingTerm + , sprayDivision + , groebner + , esPolynomial + , isSymmetricSpray + , isPolynomialOf ) where import qualified Algebra.Additive as AlgAdd +import qualified Algebra.Field as AlgField import qualified Algebra.Module as AlgMod import qualified Algebra.Ring as AlgRing import qualified Data.Foldable as DF @@ -32,12 +55,21 @@ import Data.HashMap.Strict ( HashMap ) import qualified Data.HashMap.Strict as HM import Data.Hashable ( Hashable(hashWithSalt) ) -import Data.List ( sortBy ) +import Data.List ( sortBy + , maximumBy + , (\\) + , findIndices + ) +import Data.Maybe ( isJust + , fromJust + ) +import Data.Ord ( comparing ) import qualified Data.Sequence as S import Data.Sequence ( (><) - , Seq + , Seq (Empty, (:<|)) , dropWhileR , (|>) + , (<|) , index , adjust ) @@ -66,9 +98,11 @@ } deriving Show +-- | append trailing zeros growSequence :: Seq Int -> Int -> Int -> Seq Int growSequence s m n = s >< t where t = S.replicate (n - m) 0 +-- | append trailing zeros to get the same length harmonize :: (Powers, Powers) -> (Powers, Powers) harmonize (pows1, pows2) = (Powers e1' n, Powers e2' n) where @@ -91,6 +125,8 @@ type Spray a = HashMap Powers a +type Monomial a = (Powers, a) + instance (AlgAdd.C a, Eq a) => AlgAdd.C (Spray a) where p + q = addSprays p q zero = HM.empty @@ -129,27 +165,34 @@ then AlgAdd.sum (replicate k pol) else AlgAdd.negate $ AlgAdd.sum (replicate (-k) pol) +-- | drop trailing zeros simplifyPowers :: Powers -> Powers simplifyPowers pows = Powers s (S.length s) where s = dropWhileR (== 0) (exponents pows) +-- | drop trailing zeros in the powers of a spray simplifySpray :: Spray a -> Spray a simplifySpray = HM.mapKeys simplifyPowers +-- | simplify powers and remove zero terms cleanSpray :: (AlgAdd.C a, Eq a) => Spray a -> Spray a cleanSpray p = HM.filter (/= AlgAdd.zero) (simplifySpray p) +-- | 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 +-- | 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 = cleanSpray $ HM.map (lambda AlgRing.*) p -derivMonomial :: AlgRing.C a => Int -> (Powers, a) -> (Powers, a) +-- | derivative of a monomial +derivMonomial :: AlgRing.C a => Int -> Monomial a -> Monomial a derivMonomial i (pows, coef) = if i' >= S.length expts then (Powers S.empty 0, AlgAdd.zero) else (pows', coef') @@ -161,20 +204,26 @@ coef' = AlgAdd.sum (replicate expt_i coef) pows' = Powers expts' (nvariables pows) -derivSpray :: (AlgRing.C a, Eq a) => Int -> Spray a -> Spray a +-- | Derivative of a spray +derivSpray + :: (AlgRing.C a, Eq a) + => Int -- ^ index of the variable of differentiation + -> Spray a -- ^ the spray + -> Spray a derivSpray i p = cleanSpray $ HM.fromListWith (AlgAdd.+) monomials where p' = HM.toList p monomials = [ derivMonomial i mp | mp <- p' ] - -multMonomial :: AlgRing.C a => (Powers, a) -> (Powers, a) -> (Powers, a) +-- | multiply two monomials +multMonomial :: AlgRing.C a => Monomial a -> Monomial a -> Monomial a multMonomial (pows1, coef1) (pows2, coef2) = (pows, coef1 AlgRing.* coef2) where (pows1', pows2') = harmonize (pows1, pows2) expts = S.zipWith (+) (exponents pows1') (exponents pows2') pows = Powers expts (nvariables pows1') +-- | multiply two sprays multSprays :: (AlgRing.C a, Eq a) => Spray a -> Spray a -> Spray a multSprays p q = cleanSpray $ HM.fromListWith (AlgAdd.+) prods where @@ -182,7 +231,7 @@ q' = HM.toList q prods = [ multMonomial mp mq | mp <- p', mq <- q' ] --- | Spray corresponding to polynomial x_n +-- | Spray corresponding to the basic monomial x_n lone :: AlgRing.C a => Int -> Spray a lone n = HM.singleton pows AlgRing.one where @@ -198,15 +247,55 @@ constantSpray :: (AlgRing.C a, Eq a) => a -> Spray a constantSpray c = c *^ lone 0 -evalMonomial :: AlgRing.C a => [a] -> (Powers, a) -> a -evalMonomial xyz (powers, coeff) = coeff - AlgRing.* AlgRing.product (zipWith (AlgRing.^) xyz pows) +-- | evaluates a monomial +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) -- | Evaluate a spray evalSpray :: AlgRing.C a => Spray a -> [a] -> a evalSpray p xyz = AlgAdd.sum $ map (evalMonomial xyz) (HM.toList p) +-- | 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') + 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 + powers'' = simplifyPowers $ Powers pows'' n + +-- | Substitute some variables in a spray +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" + where + n = numberOfVariables spray + monomials = HM.toList spray + spray' = foldl1 (^+^) (map (fromMonomial . substituteMonomial subs) monomials) + +-- | Convert a spray with rational coefficients to a spray with double coefficients +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 @@ -219,6 +308,10 @@ fromList x = cleanSpray $ HM.fromList $ map (\(expts, coef) -> (Powers (S.fromList expts) (length expts), coef)) x + +-- pretty stuff --------------------------------------------------------------- + +-- | prettyPowers "x" [0, 2, 1] = x^(0, 2, 1) prettyPowers :: String -> [Int] -> Text prettyPowers var pows = append (pack x) (cons '(' $ snoc string ')') where @@ -226,7 +319,11 @@ string = intercalate (pack ", ") (map (pack . show) pows) -- | Pretty form of a spray -prettySpray :: (a -> String) -> String -> Spray a -> String +prettySpray + :: (a -> String) -- ^ function mapping a coefficient to a string, typically 'show' + -> String -- ^ a string denoting the variable, e.g. "x" + -> Spray a -- ^ the spray + -> String prettySpray prettyCoef var p = unpack $ intercalate (pack " + ") stringTerms where stringTerms = map stringTerm (sortBy (compare `on` fexpts) (HM.toList p)) @@ -238,6 +335,65 @@ pows = DF.toList $ exponents (fst term) stringCoef = pack $ prettyCoef (snd term) +-- | prettyPowers' [0, 2, 1] = "x2^2x3" +prettyPowers' :: Seq Int -> Text +prettyPowers' pows = pack x1x2x3 + where + n = S.length pows + f i p + | p == 0 = "" + | p == 1 = "x" ++ show i + | otherwise = "x" ++ show i ++ "^" ++ show p + x1x2x3 = concatMap (\i -> f i (pows `index` (i-1))) [1 .. n] + +-- | Pretty form of a spray, with monomials showed as "x1x3^2" +prettySpray' :: (Show a) => Spray a -> String +prettySpray' spray = unpack $ intercalate (pack " + ") terms + where + terms = map stringTerm (sortBy (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 ')' + stringCoef'' = if constant then stringCoef' else snoc stringCoef' ' ' + +-- | prettyPowersXYZ [1, 2, 1] = XY^2Z +prettyPowersXYZ :: Seq Int -> Text +prettyPowersXYZ pows = if n <= 3 + then pack xyz + else error "there is more than three variables" + where + n = S.length pows + gpows = growSequence pows n 3 + f letter p + | 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) + xyz = x ++ y ++ z + +-- | Pretty form of a spray having at more three variables +prettySprayXYZ :: (Show a) => Spray a -> String +prettySprayXYZ spray = unpack $ intercalate (pack " + ") terms + where + terms = map stringTerm (sortBy (compare `on` fexpts) (HM.toList spray)) + fexpts term = exponents $ fst term + stringTerm term = append stringCoef'' (prettyPowersXYZ pows) + where + 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' ' ' + + +-- misc ----------------------------------------------------------------------- + -- | Terms of a spray sprayTerms :: Spray a -> HashMap (Seq Int) a sprayTerms = HM.mapKeys exponents @@ -246,7 +402,7 @@ toList :: Spray a -> [([Int], a)] toList p = HM.toList $ HM.mapKeys (DF.toList . exponents) p --- | Bombieri spray +-- | Bombieri spray (for internal usage in the 'scubature' library) bombieriSpray :: AlgAdd.C a => Spray a -> Spray a bombieriSpray = HM.mapWithKey f where @@ -254,3 +410,293 @@ pfactorial pows = product $ DF.toList $ factorial <$> S.filter (/= 0) pows factorial n = product [1 .. n] times k x = AlgAdd.sum (replicate k x) + + +-- division stuff ------------------------------------------------------------- + +-- | index of the maximum of a list +maxIndex :: Ord a => [a] -> Int +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 + biggest = powers !! i + +-- | whether a monomial divides another monomial +divides :: Monomial a -> Monomial a -> Bool +divides (powsP, _) (powsQ, _) = S.length expntsP <= S.length expntsQ && lower + where + expntsP = exponents powsP + expntsQ = exponents powsQ + lower = DF.all (\(x, y) -> x <= y) (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 + +-- | 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 = + if n == 0 then error "the list of divisors is empty" else snd $ ogo p AlgAdd.zero + where + n = length qs + qsltqs = zip qs (map leadingTerm qs) + g :: Monomial a -> Spray a -> Spray a -> (Spray a, Spray a) + g lts s r = (s ^-^ ltsspray, r ^+^ ltsspray) + where + ltsspray = fromMonomial lts + 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 + where + (q, ltq) = qsltqs !! i + newdivoccured = divides ltq lts + 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' + where + (s', r') = go (leadingTerm s) s 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 + where + n = HM.size qsltqs + g :: Monomial a -> Spray a -> Spray a -> (Spray a, Spray a) + g lts s r = (s ^-^ ltsspray, r ^+^ ltsspray) + where + ltsspray = fromMonomial lts + 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 + where + (q, ltq) = qsltqs HM.! i + newdivoccured = divides ltq lts + 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' + where + (s', r') = go (leadingTerm s) s r 0 False + +combn2 :: Int -> Int -> HashMap Int (Int, Int) +combn2 n s = HM.fromList (zip [0 .. n-2] (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] + +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 + (lpowsP', lpowsQ') = harmonize (lpowsP, 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' + wp = fromMonomial (Powers betaP n, AlgField.recip lcoefP) + wq = fromMonomial (Powers betaQ n, 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, Monomial a) -> [Spray a] + go !i !j !combins !gpolysMap + | 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 = sprayDivision' sfg gpolysMap + ltsbarfg = leadingTerm sbarfg + (i', j', gpolysMap', combins') = if sbarfg == AlgAdd.zero + 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 [] + +-- | reduce 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 + where + normalize :: Spray a -> Spray a + normalize spray = AlgField.recip coef *^ spray + where + (_, coef) = leadingTerm spray + ngbasis = map normalize gbasis + n = length ngbasis + reduction :: Int -> Spray a + reduction i = sprayDivision (ngbasis !! i) rest + where + rest = [ngbasis !! k | k <- [0 .. n-1] \\ [i]] + +-- | Groebner basis (always minimal and possibly reduced) +groebner + :: forall a. (Eq a, AlgField.C a) + => [Spray a] -- ^ list of sprays + -> Bool -- ^ whether to return the reduced basis + -> [Spray a] +groebner sprays reduced = + if reduced then reduceGroebnerBasis gbasis0 else map normalize gbasis0 + where + gbasis0 = groebner0 sprays + normalize :: Spray a -> Spray a + normalize spray = AlgField.recip coef *^ spray + where + (_, coef) = leadingTerm spray + + +-- elementary symmetric polynomials ------------------------------------------- + +-- | generate all permutations of a binary sequence +permutationsBinarySequence :: Int -> Int -> [Seq Int] +permutationsBinarySequence nzeros nones = unfold1 next z + where + z = (><) (S.replicate nzeros False) (S.replicate nones True) + unfold1 :: (Seq Bool -> Maybe (Seq Bool)) -> Seq Bool -> [Seq Int] + unfold1 f x = case f x of + Nothing -> [x'] + Just y -> x' : unfold1 f y + where + x' = fmap fromEnum x + next :: Seq Bool -> Maybe (Seq Bool) + next xs = case findj (S.reverse xs, S.empty) of + Nothing -> Nothing + Just ( l:<|ls , rs ) -> Just $ inc l ls (S.reverse rs, S.empty) + Just ( Empty , _ ) -> error "permutationsBinarySequence: should not happen" + findj :: (Seq Bool, Seq Bool) -> Maybe (Seq Bool, Seq Bool) + findj ( xxs@(x:<|xs), yys@(_:<|_) ) = if x + then findj ( xs, True <| yys ) + else Just ( xxs, yys ) + findj ( x:<|xs, Empty ) = findj ( xs, S.singleton x ) + findj ( Empty , _ ) = Nothing + inc :: Bool -> Seq Bool -> (Seq Bool, Seq Bool) -> Seq Bool + inc !u us ( x:<|xs , yys ) = if u + then inc True us ( xs , x <| yys ) + else (><) (S.reverse (True <| us)) ((><) (S.reverse (u <| yys)) xs) + inc _ _ ( Empty , _ ) = error "permutationsBinarySequence: should not happen" + +-- | Elementary symmetric polynomial +esPolynomial + :: (AlgRing.C a, Eq a) + => Int -- ^ number of variables + -> Int -- ^ index + -> Spray a +esPolynomial n k + | k <= 0 || n <= 0 = error "both arguments must be positive integers" + | k > n = AlgAdd.zero + | otherwise = simplifySpray spray + where + perms = permutationsBinarySequence (n-k) k + spray = HM.fromList $ map (\expts -> (Powers expts n, AlgRing.one)) perms + +-- | Whether a spray is a symmetric polynomial +isSymmetricSpray :: forall a. (AlgField.C a, Eq a) => Spray a -> Bool +isSymmetricSpray spray = check1 && check2 + 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 + 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) + +-- | 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 +isPolynomialOf :: forall a. (AlgField.C a, Eq a) => Spray a -> [Spray a] -> (Bool, Maybe (Spray a)) +isPolynomialOf spray sprays = result + where + n = numberOfVariables spray + n' = maximum $ map numberOfVariables sprays + result + | 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 + gbasis0 = groebner0 gPolys + g = sprayDivision 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 + then Just $ dropXis g + else Nothing + dropXis = HM.mapKeys f + f (Powers expnnts _) = Powers (S.drop n expnnts) n diff --git a/tests/Approx.hs b/tests/Approx.hs new file mode 100644 --- /dev/null +++ b/tests/Approx.hs @@ -0,0 +1,10 @@ +module Approx (assertApproxEqual) where +import Test.Tasty.HUnit ( Assertion, assertEqual ) + +-- round x to n digits +approx :: Int -> Double -> Double +approx n x = fromInteger (round $ x * (10^n)) / (10.0^^n) + +assertApproxEqual :: String -> Int -> Double -> Double -> Assertion +assertApproxEqual prefix n x1 x2 = + assertEqual prefix (approx n x1) (approx n x2) diff --git a/tests/Main.hs b/tests/Main.hs --- a/tests/Main.hs +++ b/tests/Main.hs @@ -1,22 +1,32 @@ module Main where +import Approx ( assertApproxEqual ) import Data.Ratio ( (%) ) import Math.Algebra.Hspray ( Spray, (^+^), + (^-^), (^*^), (^**^), (*^), lone, unitSpray, evalSpray, + substituteSpray, composeSpray, fromList, toList, bombieriSpray, - derivSpray ) + derivSpray, + groebner, + fromRationalSpray, + esPolynomial, + isSymmetricSpray, + isPolynomialOf + ) import Test.Tasty ( defaultMain , testGroup ) import Test.Tasty.HUnit ( assertEqual + , assertBool , testCase ) @@ -68,6 +78,44 @@ p1' = derivSpray 1 p1 p2' = derivSpray 1 p2 q' = derivSpray 1 q - assertEqual "" q' ((p1' ^*^ p2) ^+^ (p1 ^*^ p2')) + assertEqual "" q' ((p1' ^*^ p2) ^+^ (p1 ^*^ p2')), + testCase "groebner" $ do + let + x = lone 1 :: Spray Rational + y = lone 2 :: Spray Rational + z = lone 3 :: Spray Rational + p1 = x^**^2 ^+^ y ^+^ z ^-^ unitSpray + p2 = x ^+^ y^**^2 ^+^ z ^-^ unitSpray + p3 = x ^+^ y ^+^ z^**^2 ^-^ unitSpray + g = groebner [p1, p2, p3] True + xyz = [sqrt 2 - 1, sqrt 2 - 1, sqrt 2 - 1] + gxyz = map ((`evalSpray` xyz) . fromRationalSpray) g + sumAbsValues = sum $ map abs gxyz + assertApproxEqual "" 8 sumAbsValues 0, + + testCase "symmetric polynomials" $ do + let + e2 = esPolynomial 4 2 :: Spray Rational + e3 = esPolynomial 4 3 :: Spray Rational + p = e2^**^2 ^+^ (2*^ e3) + assertBool "" (isSymmetricSpray p), + + testCase "isPolynomialOf" $ do + let + x = lone 1 :: Spray Rational + y = lone 2 :: Spray Rational + p1 = x ^+^ y + p2 = x ^-^ y + p = p1 ^*^ p2 + assertEqual "" (isPolynomialOf p [p1, p2]) (True, Just $ x ^*^ y), + + testCase "substituteSpray" $ do + let + x1 = lone 1 :: Spray Rational + x2 = lone 2 :: Spray Rational + x3 = lone 3 :: Spray Rational + p = x1^**^2 ^+^ x2 ^+^ x3 ^-^ unitSpray + p' = substituteSpray [Just 2, Nothing, Just 3] p + assertEqual "" p' (x2 ^+^ (6*^ unitSpray)) ]