hspray 0.5.0.0 → 0.5.1.0
raw patch · 5 files changed
+260/−98 lines, 5 filesPVP ok
version bump matches the API change (PVP)
API changes (from Hackage documentation)
+ Math.Algebra.Hspray: allCoefficients :: Spray a -> [a]
+ Math.Algebra.Hspray: allExponents :: Spray a -> [Exponents]
+ Math.Algebra.Hspray: dropVariables :: FunctionLike b => Int -> b -> b
+ Math.Algebra.Hspray: involvesVariable :: FunctionLike b => b -> Int -> Bool
Files
- CHANGELOG.md +20/−0
- README.md +68/−19
- hspray.cabal +1/−1
- src/Math/Algebra/Hspray.hs +107/−71
- tests/Main.hs +64/−7
CHANGELOG.md view
@@ -260,3 +260,23 @@ * Function `isHomogeneousSpray`, to check whether a spray defines a homogeneous polynomial. + + +## 0.5.1.0 - 2024-05-03 + +* An error in an internal function resulted in an error in `groebnerBasis`. It +has been fixed. + +* A limit on the number of elements of a Gröbner basis has been set in the +algorithm performed by `groebnerBasis`. When this limit is attained, an error +is thrown. The reason of this limit is that I encountered an example of a +large Gröbner basis and the algorithm took a quite long time. + +* There was an error in `esPolynomial`. + +* The `FunctionLike` class provides two new functions: `involvesVariable`, to +check whether a variable is involved in a function-like object (a spray or a +ratio of sprays), and `dropVariables`, to drop a given number of leading +variables in a function-like object. The `dropVariables` functions is very +unsafe: if a variable is dropped while it is involved, the result can be an +invalid function-like object.
README.md view
@@ -133,28 +133,77 @@ As of version 2.0.0, it is possible to compute a Gröbner basis of the ideal generated by a list of spray polynomials. +One of the numerous applications of Gröbner bases is the *implicitization* +of a system of parametric equations. That means that they allow to get an +implicit equation equivalent to a given set of parametric equations, for +example an implicit equation of a parametrically defined surface. Let us give +an example of implicitization for an *Enneper surface*. Here are the parametric +equations of this surface: + ```haskell import Math.Algebra.Hspray -import Data.Ratio --- define the elementary monomials -o = lone 0 :: Spray Rational -- same as unitSpray -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 prettyQSpray gbasis -mapM_ print prettyResult --- "x + y + z^2 - 1" --- "y^2 - y - z^2 + z" --- "y.z^2 + (1/2)*z^4 - (1/2)*z^2" --- "z^6 - 4*z^4 + 4*z^3 - z^2" +import Data.Ratio ( (%) ) +u = qlone 1 +v = qlone 2 +x = 3*^u ^+^ 3*^(u ^*^ v^**^2) ^-^ u^**^3 +y = 3*^v ^+^ 3*^(u^**^2 ^*^ v) ^-^ v^**^3 +z = 3*^u^**^2 ^-^ 3*^v^**^2 ``` + +The first step of the implicitization is to compute a Gröbner basis of the +following list of polynomials: + +```haskell +generators = [x ^-^ qlone 3, y ^-^ qlone 4, z ^-^ qlone 5] +``` + +Once the Gröbner basis is obtained, one has to identify which of its elements +do not involve the first two variables (`u` and `v`): + +```haskell +gbasis = groebnerBasis generators True +isFreeOfUV :: QSpray -> Bool +isFreeOfUV p = not (involvesVariable p 1) && not (involvesVariable p 2) +freeOfUV = filter isFreeOfUV gbasis +``` + +Then the implicit equations (there can be several ones) are obtained by +removing the first two variables from these elements: + +```haskell +results = map (dropVariables 2) freeOfUV +``` + +Here we find only one implicit equation (i.e. `length results` is `1`). +We only partially display it because it is long: + +```haskell +implicitEquation = results !! 0 +putStrLn $ prettyQSpray implicitEquation +-- x^6 - 3*x^4.y^2 + (5/9)*x^4.z^3 + 6*x^4.z^2 - 3*x^4.z + 3*x^2.y^4 + ... +``` + +Let us check it is correct. We take a point on the Enneper surface, for +example the point corresponding to $u=1/4$ and $v=2/3$ ($u$ and $v$ range +from $0$ to $1$): + +```haskell +xyz = map (evaluateAt [1%4, 2%3]) [x, y, z] +``` + +If the implicitization is correct, then the polynomial `implicitEquation` +should be zero at any point on the Enneper surface. This is true for our point: + +```haskell +evaluateAt xyz implicitEquation == 0 +-- True +``` + +In the [unit tests](https://github.com/stla/hspray/blob/main/tests/Main.hs), +you will find a more complex example of implicitization: the implicitization +of the ellipse parametrically defined by $x = a \cos t$ and $y = b \sin t$. +It is more complex because there are the parameters $a$ and $b$ and because +one has to use the relation ${(\cos t)}^2 + {(\sin t)}^2 = 1$. ## Easier usage
hspray.cabal view
@@ -1,5 +1,5 @@ name: hspray -version: 0.5.0.0 +version: 0.5.1.0 synopsis: Multivariate polynomials and fractions of multivariate polynomials. description: Manipulation of multivariate polynomials over a commutative ring and fractions of multivariate polynomials over a commutative field, Gröbner bases, resultant and subresultants, and greatest common divisor. It is possible to deal with multivariate polynomials whose coefficients are fractions of multivariate polynomials, and they can be interpreted as parametric polynomials with symbolic parameters. homepage: https://github.com/stla/hspray#readme
src/Math/Algebra/Hspray.hs view
@@ -184,6 +184,8 @@ , getConstantTerm , isConstantSpray , isHomogeneousSpray + , allExponents + , allCoefficients -- * Evaluation of a spray , evalSpray , substituteSpray @@ -295,10 +297,10 @@ -- introduce a class because it will be assigned to the ratios of sprays too. class FunctionLike b where - -- | Number of variables + -- | Number of variables in a function-like object numberOfVariables :: b -> Int - -- | Permutes the variables + -- | Permutes the variables of a function-like object -- -- >>> f :: Spray Rational -> Spray Rational -> Spray Rational -> Spray Rational -- >>> f p1 p2 p3 = p1^**^4 ^+^ (2*^p2^**^3) ^+^ (3*^p3^**^2) ^-^ (4*^unitSpray) @@ -310,19 +312,37 @@ -- prop> permuteVariables [3, 1, 2] spray == f x3 x1 x2 permuteVariables :: [Int] -- ^ permutation - -> b -- ^ the object whose variables will be permuted - -> b -- ^ the object with permuted variables + -> b -- ^ function-like object whose variables will be permuted + -> b -- ^ the function-like object with permuted variables - -- | Swaps two variables + -- | Swaps two variables of a function-like object -- -- prop> swapVariables (1, 3) x == permuteVariables [3, 2, 1] x swapVariables :: (Int, Int) -- ^ the indices of the variables to be swapped (starting at 1) - -> b -- ^ the object whose variables will be swapped - -> b -- ^ the object with swapped variables + -> b -- ^ function-like object whose variables will be swapped + -> b -- ^ the function-like object with swapped variables - -- | Derivative + -- Whether a variable is involved in a function-like object -- + -- prop> involvesVariable (qlone 1 ^+^ qlone 3) 2 == False + involvesVariable :: + b -- ^ function-like object + -> Int -- ^ index of the variable + -> Bool + + -- | Drops a given number of leading variables in a function-like object; + -- __very unsafe__, @dropVariables n x@ should /not/ be used if + -- @involvesVariable x i@ is @True@ for some @i@ in @1, ... n@ + -- + -- prop> dropVariables 1 (qlone 2 ^+^ qlone 3) == qlone 1 ^+^ qlone 2 + dropVariables :: + Int -- ^ number of leading variables to drop + -> b -- ^ a function-like object + -> b + + -- | Derivative of a function-like object + -- -- >>> x = lone 1 :: Spray Int -- >>> y = lone 2 :: Spray Int -- >>> spray = 2*^x ^-^ 3*^y^**^8 @@ -341,43 +361,44 @@ type family VariablesType b infixl 6 ^+^ - -- | Addition + -- | Addition of two function-like objects (^+^) :: (AlgAdd.C b) => b -> b -> b (^+^) = (AlgAdd.+) infixl 6 ^-^ - -- | Substraction + -- | Substraction of two function-like objects (^-^) :: (AlgAdd.C b) => b -> b -> b (^-^) = (AlgAdd.-) infixl 7 ^*^ - -- | Multiplication + -- | Multiplication of two function-like objects (^*^) :: (AlgRing.C b) => b -> b -> b (^*^) = (AlgRing.*) infixr 8 ^**^ - -- | Power + -- | Power of a function-like object (^**^) :: (AlgRing.C b) => b -> Int -> b (^**^) p k = p AlgRing.^ fromIntegral k infixr 7 *^ - -- | Multiply by a scalar + -- | Multiply a function-like object by a scalar (*^) :: BaseRing b -> b -> b infixl 6 +> - -- | Add a constant at left + -- | Add a function-like object to a constant -- -- prop> x +> spray == constantSpray x ^+^ spray (+>) :: BaseRing b -> b -> b infixr 6 <+ - -- | Add a constant at right + -- | Add a constant to a function-like object -- -- prop> object <+ x == x +> object (<+) :: b -> BaseRing b -> b (<+) = flip (+>) - -- | Evaluation (replacing the variables with some values) + -- | Evaluation (replacing the variables with some values) of a + -- function-like object -- -- >>> x = lone 1 :: Spray Int -- >>> y = lone 2 :: Spray Int @@ -395,7 +416,8 @@ evaluateAt :: [BaseRing b] -> b -> BaseRing b evaluateAt = flip evaluate - -- | Partial evaluation + -- | Partial evaluation of a function-like object (replace some variables + -- with some values) -- -- >>> x1 = lone 1 :: Spray Int -- >>> x2 = lone 2 :: Spray Int @@ -409,7 +431,7 @@ -> b -- ^ function-like object to be partially evaluated -> b - -- | Polynomial change of variables + -- | Polynomial change of variables of a function-like object -- -- >>> x = lone 1 :: Spray Int -- >>> y = lone 2 :: Spray Int @@ -521,6 +543,12 @@ swapVariables :: (Int, Int) -> Polynomial a -> Polynomial a swapVariables = error "swapVariables: there is only one variable." -- + involvesVariable :: Polynomial a -> Int -> Bool + involvesVariable = error "involvesVariable: not available for `Polynomial`" + -- + dropVariables :: Int -> Polynomial a -> Polynomial a + dropVariables = error "dropVariables: not available for `Polynomial`" + -- derivative :: Int -> Polynomial a -> Polynomial a derivative i p = if i == 1 @@ -562,6 +590,12 @@ swapVariables :: (Int, Int) -> RatioOfPolynomials a -> RatioOfPolynomials a swapVariables = error "swapVariables: there is only one variable." -- + involvesVariable :: RatioOfPolynomials a -> Int -> Bool + involvesVariable = error "involvesVariable: not available for `RatioOfPolynomials`" + -- + dropVariables :: Int -> RatioOfPolynomials a -> RatioOfPolynomials a + dropVariables = error "dropVariables: not available for `RatioOfPolynomials`" + -- derivative :: Int -> RatioOfPolynomials a -> RatioOfPolynomials a derivative i r = if i == 1 @@ -579,12 +613,6 @@ changeVariables r ps = changeVariables (NumberRatio.numerator r) ps % changeVariables (NumberRatio.denominator r) ps -{- -- | Division of univariate polynomials; this is an application of `:%` --- followed by a simplification of the obtained fraction of the two polynomials -(^/^) :: (Eq a, AlgField.C a) - => Polynomial a -> Polynomial a -> RatioOfPolynomials a -(^/^) pol1 pol2 = simplifyRatioOfPolynomials $ pol1 :% pol2 - -} instance (Eq a, AlgField.C a) => AlgZT.C (A a) where isZero :: A a -> Bool isZero (A r) = r == AlgAdd.zero @@ -1125,6 +1153,17 @@ makePowers expnts powers' = map g powers -- + involvesVariable :: Spray a -> Int -> Bool + involvesVariable spray i = any f (allExponents spray) + where + f expnts = let p = S.lookup (i - 1) expnts in + isJust p && p /= Just 0 + -- + dropVariables :: Int -> Spray a -> Spray a + dropVariables n = HM.mapKeys f + where + f (Powers exps nv) = Powers (S.drop n exps) (nv - n) + -- derivative :: Int -> Spray a -> Spray a derivative i p = if i >= 1 then removeZeroTerms $ HM.fromListWith (AlgAdd.+) terms @@ -1212,16 +1251,7 @@ (*) :: Spray a -> Spray a -> Spray a p * q = multSprays p q one :: Spray a - 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" - -} + one = HM.singleton (Powers S.empty 0) AlgRing.one infixr 7 /^ -- | Divides a spray by a scalar; you can equivalently use `(/>)` if the type @@ -1336,16 +1366,20 @@ -- >>> y = lone 2 :: Spray Int -- >>> z = lone 3 :: Spray Int -- >>> p = 2 *^ (2 *^ (x^**^3 ^*^ y^**^2)) ^+^ 4*^z ^+^ 5*^unitSpray --- >>> getCoefficient [3, 2, 0] p +-- >>> getCoefficient [3, 2] p -- coefficient of x^3.y^2 -- 4 --- >>> getCoefficient [0, 4] p +-- >>> getCoefficient [3, 2, 0] p -- same as getCoefficient [3, 2] p +-- 4 +-- >>> getCoefficient [0, 4] p -- coefficient of y^4 -- 0 -getCoefficient :: AlgAdd.C a => [Int] -> Spray a -> a +getCoefficient :: + AlgAdd.C a + => [Int] -- ^ sequence of exponents; trailing zeros are dropped + -> Spray a -- ^ spray + -> a -- ^ coefficient of the monomial given by the sequence of exponents getCoefficient expnts = getCoefficient' powers where powers = makePowers (S.fromList expnts) --- expnts' = S.dropWhileR (== 0) (S.fromList expnts) --- powers = Powers expnts' (S.length expnts') getCoefficient' :: AlgAdd.C a => Powers -> Spray a -> a getCoefficient' powers spray = fromMaybe AlgAdd.zero (HM.lookup powers spray) @@ -1608,7 +1642,7 @@ showNumSpray :: (Num a, Ord a) => ([Exponents] -> [String]) -- ^ function mapping a list of monomial exponents to a list of strings representing the monomials - -> (a -> String) -- ^ function mapping a positive coefficient to a string + -> (a -> String) -- ^ function mapping a positive coefficient to a string -> Spray a -> String showNumSpray showMonomials showCoeff spray = @@ -1843,13 +1877,13 @@ toList :: Spray a -> [([Int], a)] toList p = HM.toList $ HM.mapKeys (DF.toList . exponents) p --- | Bombieri spray (for internal usage in the \'scubature\' library) +-- | Bombieri spray (for internal usage in the \'__scubature__\' package) bombieriSpray :: (Eq a, AlgAdd.C a) => Spray a -> Spray a bombieriSpray = HM.mapWithKey f where f pows = times (pfactorial $ exponents pows) pfactorial pows = product $ DF.toList $ factorial <$> S.filter (/= 0) pows - factorial n = product [1 .. n] + factorial n = product [2 .. n] times k x = k .^ x -- | Whether two sprays are equal up to a scalar factor @@ -1871,7 +1905,14 @@ check = allSame degrees deg = if check then Just (degrees !! 0) else Nothing +-- | Get all the exponents of a spray +allExponents :: Spray a -> [Exponents] +allExponents spray = map exponents (HM.keys spray) +-- | Get all the coefficients of a spray +allCoefficients :: Spray a -> [a] +allCoefficients = HM.elems + -- division stuff ------------------------------------------------------------- -- | index of the maximum of a list @@ -2039,9 +2080,8 @@ -- combinations of two among n combn2 :: Int -> Int -> HashMap Int (Int, Int) -combn2 n s = HM.fromList (zip range0 (zip row1 row2)) +combn2 n s = HM.fromList (zip [0 .. ] (zip row1 row2)) where - 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 @@ -2071,6 +2111,8 @@ go :: Int -> Int -> HashMap Int (Int, Int) -> HashMap Int (Spray a, Term a) -> [Spray a] go !i !j !combins !gpolysMap + | j == 100 = error + "groebnerBasis: stopped because reached the limit; please fill an issue." | i == length combins = map fst (HM.elems gpolysMap) | otherwise = go i' j' combins' gpolysMap' where @@ -2080,12 +2122,12 @@ ltsbarfg = leadingTerm sbarfg (i', j', gpolysMap', combins') = if isZeroSpray sbarfg then - (i+1, j, gpolysMap, combins) + (i + 1, j, gpolysMap, combins) else ( 0 , j+1 , HM.insert j (sbarfg, ltsbarfg) gpolysMap - , combn2 (j+1) (i+1) + , combn2 (j + 1) (i + 1) ) -- | groebner basis, minimal but not reduced @@ -2167,7 +2209,7 @@ 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] + dc = product [n-i+1 .. n-1] `div` product [1 .. i-1] -- | generates all permutations of a binary sequence permutationsBinarySequence :: Int -> Int -> [Seq Int] @@ -2223,7 +2265,7 @@ -- (use the function with the same name in the __jackpolynomials__ package -- if you need efficiency) isSymmetricSpray :: forall a. (AlgField.C a, Eq a) => Spray a -> Bool -isSymmetricSpray spray = check1 && check2 +isSymmetricSpray spray = check where n = numberOfVariables spray indices = [1 .. n] @@ -2231,10 +2273,7 @@ gbasis = groebner0 gPolys spray' = removeConstantTerm spray 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) + check = not $ any (involvesVariable g) [1 .. n] -- | Whether a spray can be written as a polynomial of a given list of sprays; -- this polynomial is returned if this is true @@ -2257,7 +2296,7 @@ n = maximum $ map numberOfVariables sprays result | nov > n = (False, Nothing) - | otherwise = (checks, poly) + | otherwise = (check, poly) where m = length sprays yPolys = [lone (n + i) | i <- [1 .. m]] @@ -2266,20 +2305,10 @@ constantTerm = getConstantTerm spray spray' = removeConstantTerm spray g = sprayDivisionRemainder spray' gbasis0 - gpowers = HM.keys g - check1 = minimum (map nvariables gpowers) > n - check2 = - DF.all (DF.all (0 ==)) (map (S.take n . exponents) gpowers) - checks = check1 && check2 - poly = if checks - then Just g'' + check = not $ any (involvesVariable g) [1 .. n] + poly = if check + then Just (constantTerm +> dropVariables n g) else Nothing - g' = dropXis g - g'' = if constantTerm == AlgAdd.zero - then g' - else addTerm g' (nullPowers, constantTerm) - dropXis = HM.mapKeys f - f (Powers expnnts nv) = Powers (S.drop n expnnts) (nv - n) -- resultant ------------------------------------------------------------------ @@ -2781,6 +2810,14 @@ swapVariables (i, j) (RatioOfSprays p q) = swapVariables (i, j) p %//% swapVariables (i, j) q -- + involvesVariable :: RatioOfSprays a -> Int -> Bool + involvesVariable (RatioOfSprays p q) i = + involvesVariable p i || involvesVariable q i + -- + dropVariables :: Int -> RatioOfSprays a -> RatioOfSprays a + dropVariables n (RatioOfSprays p q) = + dropVariables n p %//% dropVariables n q + -- derivative :: Int -> RatioOfSprays a -> RatioOfSprays a derivative i (RatioOfSprays p q) = (p' ^*^ q ^-^ p ^*^ q') %//% (q ^*^ q) where @@ -3245,7 +3282,7 @@ if isZeroSpray pspray then 0 else - maximum (map numberOfVariables (HM.elems pspray)) + maximum (map numberOfVariables (allCoefficients pspray)) -- | Apply polynomial transformations to the parameters of a parametric spray; -- e.g. you have a two-parameters polynomial \(P_{a, b}(X, Y, Z)\) and you want @@ -3322,7 +3359,7 @@ canCoerceToSimpleParametricSpray :: (Eq a, AlgRing.C a) => ParametricSpray a -> Bool canCoerceToSimpleParametricSpray spray = - all isPolynomialRatioOfSprays (HM.elems spray) + all isPolynomialRatioOfSprays (allCoefficients spray) -- | Coerces a parametric spray to a simple parametric spray, without -- checking this makes sense with `canCoerceToSimpleParametricSpray` @@ -3441,13 +3478,12 @@ | n == 0 = unitSpray | n == 1 = (2.^a) *^ x | otherwise = - (2.^(n'' ^+^ a) /^ n') *^ (x ^*^ gegenbauerPolynomial (n - 1)) ^-^ - ((n'' ^+^ 2.^a ^-^ unitSpray) /^ n') *^ gegenbauerPolynomial (n - 2) + (2.^((n'-1) +> a) /^ n') *^ (x ^*^ gegenbauerPolynomial (n - 1)) ^-^ + (((n'-1) +> 2.^a ^-^ unitSpray) /^ n') *^ gegenbauerPolynomial (n - 2) where x = lone 1 :: SimpleParametricQSpray a = lone 1 :: QSpray n' = toRational n - n'' = constantSpray (n' - 1) -- | [Jacobi polynomial](https://en.wikipedia.org/wiki/Jacobi_polynomials); -- the @n@-th Jacobi polynomial is a univariate polynomial of degree @n@ with
tests/Main.hs view
@@ -95,7 +95,8 @@ qlone', qmonomial, isHomogeneousSpray, - psPolynomial + psPolynomial, + prettyQSprayXYZ ) import MathObj.Matrix ( fromRows ) import qualified MathObj.Matrix as MathMatrix @@ -117,8 +118,25 @@ "Testing hspray" [ - testCase "power sum polynomials are homogeneous " $ do + testCase "involvesVariable" $ do let + x = qlone 1 + z = qlone 3 + p = x^**^4 ^-^ x ^+^ x ^*^ z + tests = map (involvesVariable p) [0 .. 4] + assertEqual "" tests [False, True, False, True, False] + + , testCase "dropVariables" $ do + let + x = qlone 1 + y = qlone 2 + z = qlone 3 + p = y^**^4 ^-^ y ^+^ y ^*^ z + p' = x^**^4 ^-^ x ^+^ x ^*^ y + assertEqual "" (dropVariables 1 p) p' + + , testCase "power sum polynomials are homogeneous " $ do + let lambda = [4, 3, 2, 2] sprays = map (psPolynomial 4) lambda :: [Spray Int] spray = AlgRing.product sprays @@ -530,9 +548,9 @@ 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, + assertApproxEqual "" 8 sumAbsValues 0 - testCase "symmetric polynomial" $ do + , testCase "symmetric polynomial" $ do let e2 = esPolynomial 4 2 :: Spray Rational e3 = esPolynomial 4 3 :: Spray Rational @@ -549,7 +567,7 @@ x^**^2 ^*^ y ^*^ z^**^3 ^+^ x ^*^ y^**^3 ^*^ z^**^2 ^+^ x ^*^ y^**^2 ^*^ z^**^3 assertBool "" (isSymmetricSpray p), - + testCase "isPolynomialOf" $ do let x = lone 1 :: Spray Rational @@ -567,9 +585,48 @@ z = lone 3 :: Spray Rational assertEqual "" (isPolynomialOf x [x ^+^ y^*^z, y, z]) - (True, Just $ x ^-^ y^*^z), + (True, Just $ x ^-^ y^*^z) - testCase "power sum polynomials" $ do + , testCase "Groebner implicitization ellipse" $ do + let + cost = qlone 1 + sint = qlone 2 + n_variables = 2 + a = qlone 3 + b = qlone 4 + equations = [a ^*^ cost, b ^*^ sint] + relations = [cost^**^2 ^+^ sint^**^2 ^-^ unitSpray] + m = maximum (map numberOfVariables equations) + coordinates = [qlone (m + i) | i <- [1 .. length equations]] + generators = relations ++ zipWith (^-^) equations coordinates + gb = groebnerBasis generators True + isfree :: QSpray -> Bool + isfree spray = not $ any (involvesVariable spray) [1 .. n_variables] + results = filter isfree gb + results' = map (dropVariables n_variables) results + showResults = map (prettyQSprayXYZ ["a", "b", "x", "y"]) results' + assertEqual "" showResults ["a^2.b^2 - a^2.y^2 - b^2.x^2"] + + , testCase "Groebner implicitization Enneper" $ do + let + u = qlone 1 + v = qlone 2 + x = 3*^u ^+^ 3*^(u ^*^ v^**^2) ^-^ u^**^3 + y = 3*^v ^+^ 3*^(u^**^2 ^*^ v) ^-^ v^**^3 + z = 3*^u^**^2 ^-^ 3*^v^**^2 + generators = [x ^-^ qlone 3, y ^-^ qlone 4, z ^-^ qlone 5] + gb = groebnerBasis generators True + isfree :: QSpray -> Bool + isfree spray = + not (involvesVariable spray 1) && not (involvesVariable spray 2) + results = filter isfree gb + results' = map (dropVariables 2) results + xyz = map (evaluateAt [1%4, 2%3]) [x, y, z] + equation = results' !! 0 + shouldBeZero = evaluateAt xyz equation + assertEqual "" (length results', shouldBeZero) (1, 0) + + , testCase "power sum polynomials" $ do let x = lone 1 :: Spray Rational y = lone 2 :: Spray Rational