hspray
Simple multivariate polynomials in Haskell.
This package deals with multivariate polynomials over a commutative ring,
fractions of multivariate polynomials over a commutative field, and
multivariate polynomials with symbolic parameters in their coefficients.
The main type provided by this package is Spray a
.
An object of type Spray a
represents a multivariate polynomial whose
coefficients are represented by the objects of type a
. For example:
import Math.Algebra.Hspray
x = lone 1 :: Spray Double
y = lone 2 :: Spray Double
z = lone 3 :: Spray Double
poly = (2 *^ (x^**^3 ^*^ y ^*^ z) ^+^ x^**^2) ^*^ (4 *^ (x ^*^ y ^*^ z))
putStrLn $ prettyNumSpray poly
-- 8.0*x^4.y^2.z^2 + 4.0*x^3.y.z
This is the easiest way to construct a spray: first introduce the polynomial
variables with the lone
function, and then combine them with arithmetic
operations.
There are numerous functions to print a spray. If you don't like the letters
x
, y
, z
in the output of prettyNumSpray
, you can use prettyNumSprayXYZ
to change them to whatever you want:
putStrLn $ prettyNumSprayXYZ ["A","B","C"] poly
-- 8.0*A^4.B^2.C^2 + 4.0*A^3.B.C
Note that this function does not throw an error if you don't provide enough
letters; in such a situation, it takes the first given letter and it appends
it with the digit i
to denote the i
-th variable:
putStrLn $ prettyNumSprayXYZ ["A","B"] poly
-- 8.0*A1^4.A2^2.A3^2 + 4.0*A1^3.A2.A3
This is the same output as the one of prettyNumSprayX1X2X3 "A" poly
.
More generally, one can use the type Spray a
as long as the type a
has
the instances Eq
and Algebra.Ring
(defined in the numeric-prelude
library). For example a = Rational
:
import Math.Algebra.Hspray
import Data.Ratio
x = lone 1 :: QSpray -- QSpray = Spray Rational
y = lone 2 :: QSpray
z = lone 3 :: QSpray
poly = ((2%3) *^ (x^**^3 ^*^ y ^*^ z) ^-^ x^**^2) ^*^ ((7%4) *^ (x ^*^ y ^*^ z))
putStrLn $ prettyQSpray poly
-- (7/6)*x^4.y^2.z^2 - (7/4)*x^3.y.z
Or a = Spray Double
:
import Math.Algebra.Hspray
alpha = lone 1 :: Spray Double
x = lone 1 :: Spray (Spray Double)
y = lone 2 :: Spray (Spray Double)
poly = ((alpha *^ x) ^+^ (alpha *^ y))^**^2
showSprayXYZ' (prettyNumSprayXYZ ["alpha"]) ["x","y"] poly
-- (alpha^2)*x^2 + (2.0*alpha^2)*x.y + (alpha^2)*y^2
We will come back to these sprays of type Spray (Spray a)
. They can
be used to represent parametric polynomials.
Evaluation of a spray:
import Math.Algebra.Hspray
x = lone 1 :: Spray Double
y = lone 2 :: Spray Double
z = lone 3 :: Spray Double
spray = 2 *^ (x ^*^ y ^*^ z)
-- evaluate spray at x=2, y=1, z=2
evalSpray spray [2, 1, 2]
-- 8.0
Partial evaluation:
import Math.Algebra.Hspray
import Data.Ratio
x1 = lone 1 :: Spray Rational
x2 = lone 2 :: Spray Rational
x3 = lone 3 :: Spray Rational
spray = x1^**^2 ^+^ x2 ^+^ x3 ^-^ unitSpray
putStrLn $ prettyQSprayX1X2X3 "x" spray
-- x1^2 + x2 + x3 - 1
--
-- substitute x1 -> 2 and x3 -> 3, and don't substitute x2
spray' = substituteSpray [Just 2, Nothing, Just 3] spray
putStrLn $ prettyQSprayX1X2X3 "x" spray'
-- x2 + 6
Differentiation of a spray:
import Math.Algebra.Hspray
x = lone 1 :: Spray Double
y = lone 2 :: Spray Double
z = lone 3 :: Spray Double
spray = 2 *^ (x ^*^ y ^*^ z) ^+^ (3 *^ x^**^2)
putStrLn $ prettyNumSpray spray
-- 3.0*x^2 + 2.0*x.y.z
--
-- derivative with respect to x
putStrLn $ prettyNumSpray $ derivative 1 spray
-- 6.0*x + 2.0*y.z"
Gröbner bases
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:
import Math.Algebra.Hspray
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:
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
):
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:
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:
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\)):
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:
evaluateAt xyz implicitEquation == 0
-- True
In the unit tests,
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
To construct a spray using the ordinary symbols +
, -
, *
and ^
,
one can hide these operators from Prelude and import them from the
numeric-prelude library; constructing a spray in this context is easier:
import Prelude hiding ((+), (-), (*), (^), (*>), (<*))
import qualified Prelude as P
import Algebra.Additive
import Algebra.Module
import Algebra.Ring
import Math.Algebra.Hspray
import Data.Ratio
x = lone 1 :: QSpray
y = lone 2 :: QSpray
z = lone 3 :: QSpray
spray = ((2%3) *^ (x^**^3 ^*^ y ^*^ z) ^-^ x^**^2) ^*^ ((7%4) *^ (x ^*^ y ^*^ z))
spray' = ((2%3) *^ (x^3 * y * z) - x^2) * ((7%4) *^ (x * y * z))
spray == spray'
-- True
Note that *>
could be used instead of *^
but running lambda *> spray
possibly throws an "ambiguous type" error regarding the type of lambda
.
Maybe better (I didn't try yet), follow the "Usage" section on the
Hackage page
of numeric-prelude.
Symbolic parameters in the coefficients
Assume you have the polynomial a * (x² + y²) + 2b/3 * z
,
where a
and b
are symbolic rational numbers. You can represent this
polynomial by a Spray (Spray Rational)
spray as follows:
import Prelude hiding ((*), (+), (-), (^))
import qualified Prelude as P
import Algebra.Additive
import Algebra.Ring
import Math.Algebra.Hspray
x = lone 1 :: Spray (Spray Rational)
y = lone 2 :: Spray (Spray Rational)
z = lone 3 :: Spray (Spray Rational)
a = lone 1 :: Spray Rational
b = lone 2 :: Spray Rational
spray = a *^ (x^2 + y^2) + ((2 *^ b) /^ 3) *^ z
putStrLn $
showSprayXYZ' (prettyQSprayXYZ ["a","b"]) ["X","Y","Z"] spray
-- (a)*X^2 + (a)*Y^2 + ((2/3)*b)*Z
You can extract the powers and the coefficients as follows:
l = toList spray
map fst l
-- [[0,0,1],[2],[0,2]]
map toList $ map snd l
-- [[([0,1],2 % 3)],[([1],1 % 1)],[([1],1 % 1)]]
These Spray (Spray a)
sprays can be very useful. They represent polynomials
whose coefficients polynomially depend on some parameters.
Actually there is a type alias of Spray (Spray a)
in hspray, namely
SimpleParametricSpray a
, and there are some convenient functions to deal
with sprays of this type. There is also a type alias of
SimpleParametricSpray Rational
, namely SimpleParametricQSpray
.
For example we can print our SimpleParametricQSpray
spray spray
as follows:
putStrLn $
prettySimpleParametricQSprayABCXYZ ["a","b"] ["X","Y","Z"] spray
-- { a }*X^2 + { a }*Y^2 + { (2/3)*b }*Z
The
Gegenbauer polynomials
are a real-life example of polynomials that can be represented by
SimpleParametricQSpray
sprays. They are univariate polynomials whose
coefficients polynomially depend on a parameter \(\alpha\) (the polynomial
dependency is clearly visible from the recurrence relation given on
Wikipedia). Here is their recursive implementation in hspray:
gegenbauerPolynomial :: Int -> SimpleParametricQSpray
gegenbauerPolynomial n
| 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)
where
x = lone 1 :: SimpleParametricQSpray
a = lone 1 :: QSpray
n' = toRational n
n'' = constantSpray (n' - 1)
Let's try it:
n = 3
g = gegenbauerPolynomial n
putStrLn $
prettySimpleParametricQSprayABCXYZ ["alpha"] ["X"] g
-- { (4/3)*alpha^3 + 4*alpha^2 + (8/3)*alpha }*X^3 + { -2*alpha^2 - 2*alpha }*X
Let's check the differential equation given in the Wikipedia article:
g' = derivative 1 g
g'' = derivative 1 g'
alpha = lone 1 :: QSpray
x = lone 1 :: SimpleParametricQSpray
nAsSpray = constantSpray (toRational n)
shouldBeZero =
(unitSpray ^-^ x^**^2) ^*^ g''
^-^ (2.^alpha ^+^ unitSpray) *^ (x ^*^ g')
^+^ n.^(nAsSpray ^+^ 2.^alpha) *^ g
putStrLn $ prettySpray shouldBeZero
-- 0
Now, how to substitute a value to the parameter \(\alpha\)? For example, it is
said in the Wikipedia article that this yields the Legendre polynomials for
\(\alpha = 1/2\). The package provides the function substituteParameters
to
perform this task:
import Data.Ratio (%)
putStrLn $
prettyQSpray'' $ substituteParameters g [1%2]
-- (5/2)*X^3 - (3/2)*X
This is a Spray Rational
spray.
The Wikipedia article also provides the value at \(1\) of the Gegenbauer
polynomials in function of \(\alpha\). We can get this value with
evalParametricSpray
:
putStrLn $
prettyQSprayXYZ ["alpha"] $ evalParametricSpray g [1]
-- (4/3)*alpha^3 + 2*alpha^2 + (2/3)*alpha
This is also a Spray Rational
spray.
Ratios of sprays and general parametric sprays
Since you have just seen that the type Spray (Spray a)
is named
SimpleParametricSpray a
, you probably guessed there is also a more general
type named ParametricSpray a
. Yes, and this is an alias of
Spray (RatioOfSprays a)
, where the type RatioOfSprays a
has not been
discussed yet. The objects of this type represent fractions of multivariate
polynomials and so this type is a considerable enlargment of the Spray a
type. Thus the Spray (RatioOfSprays a)
sprays can represent multivariate
polynomials whose coefficients depend on some parameters, with a dependence
described by a fraction of polynomials in these parameters. Let's start with
a short presentation of the ratios of sprays.
The RatioOfSprays
type
The type RatioOfSprays a
, whose objects represent ratios of sprays, has
been introduced in version 0.2.7.0.
To construct a ratio of sprays, apply %//%
between its numerator and
its denominator:
import Math.Algebra.Hspray
x = qlone 1 -- shortcut for lone 1 :: Spray Rational
y = qlone 2
rOS = (x ^-^ y) %//% (x^**^2 ^-^ y^**^2)
putStrLn $ prettyRatioOfQSprays rOS
-- [ 1 ] %//% [ x + y ]
The %//%
operator always returns an irreducible fraction. If you are
sure that your numerator and your denominator are coprime, you can use
the %:%
instead, to gain some efficiency. But if they are not coprime, this
can have unfortunate consequences.
The RatioOfSprays a
type makes sense when a
has a field instance, and then
it has a field instance too. To use the field operations, import the necessary
modules from numeric-prelude, and hide these operations from the Prelude
module; then you can also use the numeric-prelude operations for sprays,
instead of using ^+^
, ^-^
, ^*^
, ^**^
:
import Prelude hiding ((+), (-), (*), (/), (^), (*>), (<*))
import qualified Prelude as P
import Algebra.Additive
import Algebra.Module
import Algebra.RightModule
import Algebra.Ring
import Algebra.Field
import Math.Algebra.Hspray
x = qlone 1
y = qlone 2
p = x^2 - 3*^(x * y) + y^3
q = x - y
rOS1 = p^2 %//% q
rOS2 = rOS1 + unitRatioOfSprays
rOS = rOS1^2 + rOS1*rOS2 - rOS1/rOS2 + rOS2 -- slow!
(rOS1 + rOS2) * (rOS1 - rOS2) == rOS1^2 - rOS2^2
-- True
rOS / rOS == unitRatioOfSprays
-- True
Actually, as of version 0.5.0.0, it is possible to apply the operations
^+^
, ^-^
, ^*^
, ^**^
and *^
to the ratios of sprays.
The RatioOfSprays a
type also has left and right module instances over a
and over Spray a
as well. That means you can multiply a ratio of sprays by
a scalar and by a spray, by using, depending on the side, either *>
or <*
:
import Data.Ratio ( (%) )
rOS' = (3%4::Rational) *> rOS^2 + p *> rOS
You can also divide a ratio of sprays by a spray with %/%
:
p *> (rOS' %/% p) == rOS'
-- True
rOS1 %/% p == p %//% q
-- True
When a
has a field instance, both a Spray a
spray and a RatioOfSprays a
ratio of sprays can be divided by a scalar with the />
operator:
k = 3 :: Rational
(p /> k) *> rOS == p *> (rOS /> k)
-- True
Use evalRatioOfSprays
to evaluate a ratio of sprays:
import Data.Ratio ( (%) )
f :: Algebra.Field.C a => a -> a -> a
f u v = u^2 + u*v - u/v + v
rOS == f rOS1 rOS2
-- True
values = [2%3, 7%4]
r1 = evalRatioOfSprays rOS1 values
r2 = evalRatioOfSprays rOS2 values
evalRatioOfSprays rOS values == f r1 r2
-- True
The ParametricSpray
type
Recall that SimpleParametricSpray a = Spray (Spray a)
and
ParametricSpray a = Spray (RatioOfSprays a)
, and we have the aliases
SimpleParametricQSpray = SimpleParametricSpray Rational
and
ParametricQSpray = ParametricSpray Rational
.
The functions substituteParameters
and evalParametricSpray
, that we
previously applied to a SimpleParametricSpray a
spray, are also applicable
to a ParametricSpray a
spray. We didn't mention the function
changeParameters
yet, which is also applicable to these two types of
parametric sprays. This function performs some polynomial transformations of
the parameters of a parametric spray. For example, consider the
Jacobi polynomials.
They are univariate polynomials with two parameters \(\alpha\) and \(\beta\).
They are implemented in hspray as ParametricQSpray
sprays. In fact
it seems that the coefficients of the Jacobi polynomials polynomially
depend on \(\alpha\) and \(\beta\), and if this is true one could implement them
as SimpleParametricQSpray
sprays. I will come back to this point later. The
recurrence relation defining the Jacobi polynomials involves a division which
makes the type ParametricQSpray
necessary anyway.
The changeParameters
function is useful to derive the Gegenbauer polynomials
from the Jacobi polynomials. Indeed, as asserted in the Wikipedia article,
the Gegenbauer polynomials coincide, up to a factor, with the Jacobi
polynomials with parameters \(\alpha - 1/2\) and \(\alpha - 1/2\). Here is how
to apply the changeParameters
function to get this special case of Jacobi
polynomials:
import Data.Ratio ( (%) )
j = jacobiPolynomial 3
alpha = qlone 1
alpha' = alpha ^-^ constantSpray (1%2)
j' = changeParameters j [alpha', alpha']
Now let's come back to the conjecture claiming that the coefficients of the
Jacobi polynomials polynomially depend on \(\alpha\) and \(\beta\), and thus
these polynomials can be represented by SimpleParametricQSpray
sprays.
Maybe this can be deduced from a formula given in the Wikipedia article, I
didn't spend some time on this problem. I made this conjecture because I
observed this fact for some small values of \(n\), and I tried the function
canCoerceToSimpleParametricSpray
for other values, which always returned
True
. One can apply the function asSimpleParametricSpray
to perform the
coercion.
The OneParameterSpray
type
There is a third type of parametric sprays in the package, namely the
OneParameterSpray
sprays. The objects of this type represent
multivariate polynomials whose coefficients are fractions
of polynomials in only one variable (the parameter). So they are less
general than the ParametricSpray
sprays.
These sprays are no longer very useful. They have been introduced in version
0.2.5.0 and this is the first type of parametric sprays that has been provided
by the package. When the more general ParametricSpray
sprays have been
introduced, I continued to develop the OneParameterSpray
sprays because they
were more efficient than the univariate ParametricSpray
sprays. But as of
version 0.4.0.0, this is no longer the case. This is what I concluded from
some benchmarks on the Jack polynomials, implemented in the
jackpolynomials package.
These are sprays of type Spray (RatioOfPolynomials a)
, where the type
RatioOfPolynomials a
deals with objects that represent fractions of
univariate polynomials. The type of these univariate polynomials is
Polynomial a
.
The functions we have seen for the simple parametric sprays and the parametric
sprays, substituteParameters
, evalParametricSpray
, and changeParameters
,
are also applicable to the one-parameter sprays.
The OneParameterSpray
sprays were used in the
jackpolynomials package to
represent the Jack polynomials with a symbolic Jack parameter but they have
been replaced with the ParametricSpray
sprays.
Other features
Other features offered by the package include: resultant and subresultants of
two polynomials, Sturm-Habicht sequence of a polynomial, number of real roots
of a univariate polynomial, and greatest common divisor of two polynomials with
coefficients in a field.