hspray 0.1.3.0 → 0.5.4.0
raw patch · 7 files changed
Files
- CHANGELOG.md +282/−0
- README.md +466/−39
- benchmarks/Main.hs +60/−0
- hspray.cabal +43/−11
- src/Math/Algebra/Hspray.hs +4067/−256
- tests/Approx.hs +10/−0
- tests/Main.hs +1171/−71
CHANGELOG.md view
@@ -32,3 +32,285 @@ * 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 Gröbner 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. + + +## 0.2.1.0 - 2024-03-22 + +* New functions `permuteVariables` and `swapVariables`. + +* New function `resultant`, to compute the resultant of two sprays. + +* New function `subresultants`, to compute the subresultants of two sprays. + + +## 0.2.1.1 - 2024-03-25 + +* Improved the documentation. + +* Flipped the order of appearance of the terms in the output of the +`prettySpray` functions. + + +## 0.2.2.0 - 2024-03-26 + +* Fixed an error in `esPolynomial`, which resulted to a bug in +`isSymmetricSpray`. + + +## 0.2.3.0 - 2024-03-28 + +* New unit tests. + +* Fixed `resultant` and `subresultants`: the variables of the sprays they +return were incorrect. + +* New function `gcdQX`, to compute the greatest common divisor of two +univariate sprays with rational coefficients. + + +## 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. + + +## 0.2.5.0 - 2024-04-02 + +* New function `resultant'` which computes the resultant of two sprays with +coefficients in a field. Thus it is less general than the function `resultant` +but it is more efficient. + +* Fixed a small mistake in `isSymmetricSpray` and `isPolynomialOf`: these +functions didn't deal with the constant term of the spray. + +* New function `psPolynomial` which computes the power sum polynomials. + +* A particular type of sprays, namely `SymbolicSpray a`, has been introduced. +The coefficients of these sprays are ratios of univariate polynomials with `a` +coefficients. There is a specialization `SymbolicQSpray` for the case when `a` +is a type of rational numbers. The necessary instances have been defined and +there is the function `prettySymbolic(Q)Spray` to display such sprays. There are +also some functions to perform evaluation of such sprays. + + +## 0.2.6.0 - 2024-04-15 + +* New function `collinearSprays` which checks whether two sprays are collinear. + +* The function `isPolynomialOf` threw an error when the number of variables in +the spray to be tested was less than the number of variables in the list of +sprays. That is me who programmed this error and this was wrong: for example, +`x = p1 - p2^*^p3` with `p1 = x + y^*^z`, `p2 = y`, and `p3 = z`. + +* New functions to print sprays with numeric coefficients, such as +`prettyNumSpray` and `prettyQSpray`. + +* The functions `prettySpray`, `prettySpray'` and `prettySpray''` have been +changed. + +* New functions to print symbolic sprays. + +* Documentation and README have been improved. + + +## 0.2.7.0 - 2024-04-19 + +* Defined `qlone`, which is the same as `lone` but always returns a rational +spray (a `QSpray` spray). + +* The function `sprayDivision` ran into an infinite loop when the divisor was +constant. This has been fixed. + +* New function `characteristicPolynomial`, to compute the characteristic +polynomial of a matrix. + +* Gegenbauer polynomials. They have been implemented mainly to provide an +illustration of the type `Spray (Spray a)` in README. + +* New function `evalSpraySpray`, to evaluate the spray coefficients of a +`Spray (Spray a)` spray, thereby yielding a `Spray a` spray. + +* New type `RatioOfSprays a`, whose objects represent ratios of polynomials +represented by two `Spray a` objects. Thus the type `Spray (RatioOfSprays a)` +allows more possibilities than the type `SymbolicSpray a` because it is not +restricted to univariate fractions of polynomials, and obviously it also +allows more possibilities than the type `Spray (Spray a)`. Instances and +arithmetic operations for these ratios of sprays have been defined. The result +of an arithmetic operation always is an irreducible fraction. See README for +examples. + +* Jacobi polynomials. They have been implemented mainly to experiment the type +`Spray (RatioOfSprays a)`. By the way, this type has been named +`ParametricSpray a`, but this is possibly temporary. + +* The class `HasVariables a` has been introduced in order to have some functions +which apply to both `Spray a` objects and `RatioOfSprays a` objects. + +* The function `derivSpray` no longer exists. To get a derivative of a spray, +use the `derivative` function, which is also applicable to a ratio of sprays +(this is a method of the class `HasVariables`). + +* A spray with coefficients in a field can now be divided by a scalar by using +the `/>` operator. This operator can also be used to divide a ratio of sprays +(with coefficients in a field) by a scalar. + + +## 0.3.0.0 - 2024-04-21 + +* The type `SymbolicSpray a` has been renamed to `OneParameterSpray a`, and +all functions names which contained the string `Symbolic` have been changed by +replacing `Symbolic` with `OnePerameter`. + +* The class `HasVariables`, which is instantiated for `Spray` and +`RatioOfSprays`, has a new method `changeVariables` allowing to perform +polynomial transformations of the variables of a spray and of a ratio of +sprays. For sprays, this is the same as the `composeSpray` function. + +* The class `HasVariables` is now also instantiated for `RatioOfPolynomials`. + +* The type alias `ParametricSpray a = Spray (RatioOfSprays a)` has been kept +and the type alias `SimpleParametricSpray a = Spray (Spray a)` has been +introduced. We say that a `Spray b` spray is parametric when `b` has the +`HasVariables` instance. So this applies to a `ParametricSpray a` spray, to +a `SimpleParametricSpray a` spray, and also to a `OneParameterSpray a` spray +(recall that `OneParameterSpray a = Spray (RatioOfPolynomials a)`). + +* Functions to print `ParametricSpray` sprays and `SimpleParametricSpray` +sprays. + +* Function `numberOfParameters`, returning the number of parameters of a +parametric spray, that is to say the number of variables occurring in the +coefficients of this spray. + +* Function `changeParameters`, to perform polynomial transformations of the +parameters of a parametric spray. + +* Function `substituteParameters`, to replace the parameters of a parametric +spray with some values. For a `SimpleParametricSpray` spray, this function is +the same as `evalSpraySpray`(which will probably disappear in the future). + +* Function `evalParametricSpray`, to replace the variables of a parametric +spray with some values. + + +## 0.4.0.0 - 2024-04-27 + +* The efficiency of the arithmetic on the `RatioOfSprays` fractions of +polynomials has been greatly improved for the univariate case. According to +some benchmarks on the Jack polynomials, the `ParametricSpray` sprays with +only one parameter are now more efficient than the `OneParameterSpray` sprays. + +* For this reason, the Jack polynomials with a symbolic Jack parameter, +implemented in the 'jackpolynomials' package, are represented by +`ParametricSpray` sprays in a new version of the package, while there were +previously represented by `OneParameterSpray` sprays. + +* Slight improvements of the code in general. + +* Function `lone'` to construct monomials like `x_n^p` more efficiently than +`lone n ^**^ p`. + +* Function `monomial` to construct monomials like `x_1^4.x_3^7`. + + +## 0.5.0.0 - 2024-05-01 + +* The class `HasVariables` has been renamed to `FunctionLike`, and the +arithmetic operations for sprays (`^+^`, `^-^`, `^*^`, `^**^`, `*^`) have +been moved to this class. Therefore it is now possible to apply these +operations to the ratios of sprays. Moreover, there are two new operators in +this class, `(+>)` and `(<+)`, which allow to add a constant to an object of +this class. For example, `x +> spray` gives the same result as +`constantSpray x ^+^ spray` but is more efficient. + +* There was an error in `gcdSpray`. + +* There was a small mistake in `collinearSprays`. + +* The functions `changeParameters` and `changeVariables` did not remove the +possibly null terms of their result. + +* The function `groebner` has been renamed to `groebnerBasis`. + +* 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. + + +## 0.5.2.0 - 2024-05-05 + +* Function `polynomialSubresultants`, to compute the polynomial subresultants +of two sprays (while the `subresultants` function computes the principal +subresultants). + +* Function `sturmHabichtSequence`, to compute the Sturm-Habicht sequence of a +spray. + +* Function `principalSturmHabichtSequence`, to compute the principal +Sturm-Habicht sequence of a spray. + +* Functions `numberOfRealRoots`, `numberOfRealRootsInOpenInterval` and +`numberOfRealRootsInClosedInterval`, to compute the total number of real +roots of a suitable univariate spray or its number of real roots in a given +interval. This can be very slow if the degree of the spray is not small. + + +## 0.5.3.0 - 2024-06-03 + +* Now the denominator of a `RatioOfSprays` is always monic, i.e. it is a +polynomial whose leading coefficient is 1. + + +## 0.5.4.0 - 2024-06-19 + +* Now the greatest common divisor of two sprays (function `gcdSpray`) is +always monic, i.e. it is a polynomial whose leading coefficient is 1. + +* Function `evalRatioOfSprays'`, to substitute the variables of a ratio of +sprays with some ratios of sprays. +
README.md view
@@ -5,10 +5,16 @@ [](https://github.com/stla/hspray/actions/workflows/Stack-nightly.yml) <!-- badges: end --> -Simple multivariate polynomials in Haskell. +***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: ```haskell import Math.Algebra.Hspray @@ -16,10 +22,34 @@ y = lone 2 :: Spray Double z = lone 3 :: Spray Double poly = (2 *^ (x^**^3 ^*^ y ^*^ z) ^+^ x^**^2) ^*^ (4 *^ (x ^*^ y ^*^ z)) -prettySpray show "X" poly --- "(4.0) * X^(3, 1, 1) + (8.0) * X^(4, 2, 2)" +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: + +```haskell +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: + +```haskell +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`: @@ -27,85 +57,198 @@ ```haskell import Math.Algebra.Hspray import Data.Ratio -x = lone 1 :: Spray Rational -y = lone 2 :: Spray Rational -z = lone 3 :: Spray Rational -poly = ((2%3) *^ (x^**^3 ^*^ y ^*^ z) ^+^ x^**^2) ^*^ ((7%4) *^ (x ^*^ y ^*^ z)) -prettySpray show "X" poly --- "(7 % 4) * X^(3, 1, 1) + (7 % 6) * X^(4, 2, 2)" +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`: ```haskell import Math.Algebra.Hspray -p = lone 1 :: Spray Double +alpha = lone 1 :: Spray Double x = lone 1 :: Spray (Spray Double) y = lone 2 :: Spray (Spray Double) -poly = ((p *^ x) ^+^ (p *^ y))^**^2 -prettySpray (prettySpray show "a") "X" poly --- "((1.0) * a^(2)) * X^(0, 2) + ((2.0) * a^(2)) * X^(1, 1) + ((1.0) * a^(2)) * X^(2)" +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 ``` -Evaluation: +We will come back to these sprays of type `Spray (Spray a)`. They can +be used to represent parametric polynomials. + +#### Evaluation of a spray: + ```haskell import Math.Algebra.Hspray x = lone 1 :: Spray Double y = lone 2 :: Spray Double z = lone 3 :: Spray Double -poly = 2 *^ (x ^*^ y ^*^ z) --- evaluate poly at x=2, y=1, z=2 -evalSpray poly [2, 1, 2] +spray = 2 *^ (x ^*^ y ^*^ z) +-- evaluate spray at x=2, y=1, z=2 +evalSpray spray [2, 1, 2] -- 8.0 ``` -Differentiation: +#### 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 +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: + +```haskell +import Math.Algebra.Hspray x = lone 1 :: Spray Double y = lone 2 :: Spray Double z = lone 3 :: Spray Double -poly = 2 *^ (x ^*^ y ^*^ z) ^+^ (3 *^ x^**^2) --- derivate with respect to x -prettySpray show "X" $ derivSpray 1 poly --- "(2.0) * X^(0, 1, 1) + (6.0) * X^(1)" +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: + +```haskell +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: + +```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 -To construct a polynomial using the ordinary symbols `+`, `*` and `-`, -one can hide these operators from **Prelude** and import them from -the **numeric-prelude** library: +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: ```haskell -import Prelude hiding ((*), (+), (-)) +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 ``` -Or, maybe better (I didn't try yet), follow the "Usage" section on the +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](https://hackage.haskell.org/package/numeric-prelude-0.4.4#usage) of **numeric-prelude**. -## Symbolic coefficients +## Symbolic parameters in the coefficients + Assume you have the polynomial `a * (x² + y²) + 2b/3 * z`, -where `a` and `b` are symbolic coefficients. -You can define this polynomial as a `Spray` as follows: +where `a` and `b` are symbolic rational numbers. You can represent this +polynomial by a `Spray (Spray Rational)` spray as follows: ```haskell -import Prelude hiding ((*), (+), (-)) +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 :: Spray (Spray Rational) y = lone 2 :: Spray (Spray Rational) @@ -113,18 +256,302 @@ a = lone 1 :: Spray Rational b = lone 2 :: Spray Rational -poly = a *^ (x*x + y*y) + ((2%3) *^ b) *^ z -prettySpray (prettySpray show "a") "X" poly --- "((2 % 3) * a^(0, 1)) * X^(0, 0, 1) + ((1 % 1) * a^(1)) * X^(0, 2) + ((1 % 1) * a^(1)) * X^(2)" +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 ``` -The `prettySpray` function shows the expansion of the polynomial. You can extract the powers and the coefficients as follows: ```haskell -l = toList poly +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: + +```haskell +putStrLn $ + prettySimpleParametricQSprayABCXYZ ["a","b"] ["X","Y","Z"] spray +-- { a }*X^2 + { a }*Y^2 + { (2/3)*b }*Z +``` + +The +[Gegenbauer polynomials](https://en.wikipedia.org/wiki/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**: + +```haskell +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: + +```haskell +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: + +```haskell +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: + +```haskell +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`: + +```haskell +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: + +```haskell +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 `^+^`, `^-^`, `^*^`, `^**^`: + +```haskell +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 `<*`: + +```haskell +import Data.Ratio ( (%) ) +rOS' = (3%4::Rational) *> rOS^2 + p *> rOS +``` + +You can also divide a ratio of sprays by a spray with `%/%`: + +```haskell +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: + +```haskell +k = 3 :: Rational +(p /> k) *> rOS == p *> (rOS /> k) +-- True +``` + +Use `evalRatioOfSprays` to evaluate a ratio of sprays: + +```haskell +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](https://en.wikipedia.org/wiki/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: + +```haskell +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](https://github.com/stla/jackpolynomials). + +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](https://github.com/stla/jackpolynomials) 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.
+ benchmarks/Main.hs view
@@ -0,0 +1,60 @@+module Main where +import Test.Tasty.Bench ( bench, bgroup, defaultMain, nf, whnf ) +import Math.Algebra.Hspray +import qualified Algebra.Additive as AlgAdd +import qualified Algebra.Ring as AlgRing +import qualified Algebra.Module as AlgMod +import Number.Ratio ( (%), T ( (:%) ) ) + +f :: Integer -> RatioOfQPolynomials +f n = ((a AlgRing.^ 8 AlgAdd.- AlgRing.one) % (a AlgAdd.- AlgRing.one)) AlgRing.^ n AlgAdd.+ + ((a AlgAdd.+ AlgRing.one) :% a) AlgAdd.+ + ((a AlgRing.^2 AlgAdd.+ 3 .^ a) :% (a AlgRing.^ 3 AlgAdd.- a AlgAdd.+ AlgRing.one)) AlgRing.* + ((a AlgAdd.+ AlgRing.one) :% (a AlgRing.^ 3 AlgAdd.- a AlgAdd.+ AlgRing.one)) + where + a = qsoleParameter + +g :: Integer -> RatioOfQSprays +g n = ((x^**^8 ^-^ unitSpray) %//% (x ^-^ unitSpray)) AlgRing.^ n AlgAdd.+ + RatioOfSprays (x ^+^ unitSpray) x AlgAdd.+ + RatioOfSprays (x^**^2 ^+^ 3 .^ x) (x^**^3 ^-^ x ^+^ unitSpray) AlgRing.* + RatioOfSprays (x ^+^ unitSpray) (x^**^3 ^-^ x ^+^ unitSpray) + where + x = qlone 1 + +h :: Integer -> RatioOfQSprays +h n = ((x^**^8 ^-^ unitSpray) %//% (x ^-^ unitSpray)) AlgRing.^ n AlgAdd.+ + RatioOfSprays (x ^+^ unitSpray) x AlgAdd.+ + RatioOfSprays (x^**^2 ^+^ 3 .^ x) (x^**^3 ^-^ x ^+^ unitSpray) AlgRing.* + RatioOfSprays (x ^+^ unitSpray) (x^**^3 ^-^ x ^+^ unitSpray) + where + x = qlone 2 + + +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 "ratios" + [ bench "f 2" $ whnf f 2 + , bench "g 2" $ whnf g 2 + , bench "f 5" $ whnf f 5 + , bench "g 5" $ whnf g 5 + ] + ]
hspray.cabal view
@@ -1,13 +1,13 @@ name: hspray -version: 0.1.3.0 -synopsis: Multivariate polynomials. -description: Manipulation of multivariate polynomials on a ring. +version: 0.5.4.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, subresultants, Sturm-Habicht sequence, 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 license: GPL-3 license-file: LICENSE author: Stéphane Laurent maintainer: laurent_step@outlook.fr -copyright: 2022 Stéphane Laurent +copyright: 2022-2024 Stéphane Laurent category: Math, Algebra build-type: Simple extra-source-files: README.md @@ -18,15 +18,21 @@ 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 + , matrix >= 0.3.6.0 && < 0.4 + , extra >= 1.7 && < 1.8 other-extensions: FlexibleInstances , FlexibleContexts , MultiParamTypeClasses , InstanceSigs + , ScopedTypeVariables + , GeneralizedNewtypeDeriving + , BangPatterns + , TypeFamilies default-language: Haskell2010 ghc-options: -Wall -Wcompat @@ -42,11 +48,37 @@ 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 + , tasty >= 1.4 && < 1.6 + , tasty-hunit >= 0.10 && < 0.11 + , numeric-prelude >= 0.4.4 && < 0.5 + , matrix >= 0.3.6.0 && < 0.4 + , unordered-containers >= 0.2.17.0 && < 0.3 , hspray Default-Language: Haskell2010 + ghc-options: -Wall + -Wcompat + -Widentities + -Wincomplete-record-updates + -Wincomplete-uni-patterns + -Wmissing-export-lists + -Wmissing-home-modules + -Wpartial-fields + -Wredundant-constraints + +benchmark benchmarks + type: exitcode-stdio-1.0 + main-is: Main.hs + hs-source-dirs: benchmarks/ + Build-Depends: base >= 4.7 && < 5 + , tasty-bench >= 0.3 && < 0.4 + , numeric-prelude >= 0.4.4 && < 0.5 + , hspray + ghc-options: "-with-rtsopts=-A32m" + if impl(ghc >= 8.6) + ghc-options: -fproc-alignment=64 + default-language: Haskell2010 source-repository head type: git
src/Math/Algebra/Hspray.hs view
@@ -1,256 +1,4067 @@-{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE MultiParamTypeClasses #-} -{-# LANGUAGE InstanceSigs #-} - -module Math.Algebra.Hspray - ( Powers (..) - , Spray - , lone - , unitSpray - , constantSpray - , fromList - , toList - , sprayTerms - , (*^) - , (.^) - , (^+^) - , (^-^) - , (^*^) - , (^**^) - , evalSpray - , prettySpray - , composeSpray - , bombieriSpray - , derivSpray - ) where -import qualified Algebra.Additive as AlgAdd -import qualified Algebra.Module as AlgMod -import qualified Algebra.Ring as AlgRing -import qualified Data.Foldable as DF -import Data.Function ( on ) -import Data.HashMap.Strict ( HashMap ) -import qualified Data.HashMap.Strict as HM -import Data.Hashable ( Hashable(hashWithSalt) ) -import Data.List ( sortBy ) -import qualified Data.Sequence as S -import Data.Sequence ( (><) - , Seq - , dropWhileR - , (|>) - , index - , adjust - ) -import Data.Text ( Text - , append - , cons - , intercalate - , pack - , snoc - , unpack - ) - - -infixr 7 *^, .^ - -infixl 6 ^+^, ^-^ - -infixl 7 ^*^ - -infixr 8 ^**^ - - -data Powers = Powers - { exponents :: Seq Int - , nvariables :: Int - } - deriving Show - -growSequence :: Seq Int -> Int -> Int -> Seq Int -growSequence s m n = s >< t where t = S.replicate (n - m) 0 - -harmonize :: (Powers, Powers) -> (Powers, Powers) -harmonize (pows1, pows2) = (Powers e1' n, Powers e2' n) - where - e1 = exponents pows1 - e2 = exponents pows2 - n1 = nvariables pows1 - n2 = nvariables pows2 - (e1', e2', n) = if n1 < n2 - then (growSequence e1 n1 n2, e2, n2) - else (e1, growSequence e2 n2 n1, n1) - -instance Eq Powers where - (==) :: Powers -> Powers -> Bool - pows1 == pows2 = exponents pows1' == exponents pows2' - where (pows1', pows2') = harmonize (pows1, pows2) - -instance Hashable Powers where - hashWithSalt :: Int -> Powers -> Int - hashWithSalt k pows = hashWithSalt k (exponents pows, nvariables pows) - -type Spray a = HashMap Powers a - -instance (AlgAdd.C a, Eq a) => AlgAdd.C (Spray a) where - p + q = addSprays p q - zero = HM.empty - negate = negateSpray - -instance (AlgRing.C a, Eq a) => AlgMod.C a (Spray a) where - lambda *> p = scaleSpray lambda p - -instance (AlgRing.C a, Eq a) => AlgRing.C (Spray a) where - p * q = multSprays p q - one = lone 0 - --- | Addition of two sprays -(^+^) :: (AlgAdd.C a, Eq a) => Spray a -> Spray a -> Spray a -(^+^) p q = p AlgAdd.+ q - --- | Substraction of two sprays -(^-^) :: (AlgAdd.C a, Eq a) => Spray a -> Spray a -> Spray a -(^-^) p q = p AlgAdd.- q - --- | Multiply two sprays -(^*^) :: (AlgRing.C a, Eq a) => Spray a -> Spray a -> Spray a -(^*^) p q = p AlgRing.* q - --- | Power of a spray -(^**^) :: (AlgRing.C a, Eq a) => Spray a -> Int -> Spray a -(^**^) p n = AlgRing.product (replicate n p) - --- | Scale spray by a scalar -(*^) :: (AlgRing.C a, Eq a) => a -> Spray a -> Spray a -(*^) lambda pol = lambda AlgMod.*> pol - --- | Scale spray by an integer -(.^) :: (AlgAdd.C a, Eq a) => Int -> Spray a -> Spray a -(.^) k pol = if k >= 0 - then AlgAdd.sum (replicate k pol) - else AlgAdd.negate $ AlgAdd.sum (replicate (-k) pol) - -simplifyPowers :: Powers -> Powers -simplifyPowers pows = Powers s (S.length s) - where s = dropWhileR (== 0) (exponents pows) - -simplifySpray :: Spray a -> Spray a -simplifySpray = HM.mapKeys simplifyPowers - -cleanSpray :: (AlgAdd.C a, Eq a) => Spray a -> Spray a -cleanSpray p = HM.filter (/= AlgAdd.zero) (simplifySpray p) - -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 - -negateSpray :: AlgAdd.C a => Spray a -> Spray a -negateSpray = HM.map AlgAdd.negate - -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) -derivMonomial i (pows, coef) = if i' >= S.length expts - then (Powers S.empty 0, AlgAdd.zero) - else (pows', coef') - where - i' = i - 1 - expts = exponents pows - expt_i = expts `index` i' - expts' = adjust (subtract 1) i' expts - coef' = AlgAdd.sum (replicate expt_i coef) - pows' = Powers expts' (nvariables pows) - -derivSpray :: (AlgRing.C a, Eq a) => Int -> Spray a -> 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) -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') - -multSprays :: (AlgRing.C a, Eq a) => Spray a -> Spray a -> Spray a -multSprays p q = cleanSpray $ HM.fromListWith (AlgAdd.+) prods - where - p' = HM.toList p - q' = HM.toList q - prods = [ multMonomial mp mq | mp <- p', mq <- q' ] - --- | Spray corresponding to polynomial x_n -lone :: AlgRing.C a => Int -> Spray a -lone n = HM.singleton pows AlgRing.one - where - pows = if n == 0 - then Powers S.empty 0 - else Powers (S.replicate (n - 1) AlgAdd.zero |> AlgRing.one) n - --- | Unit spray -unitSpray :: AlgRing.C a => Spray a -unitSpray = lone 0 - --- | Constant spray -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) - 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) - -identify :: (AlgRing.C a, Eq a) => Spray a -> Spray (Spray a) -identify = HM.map constantSpray - --- | Compose a spray with a change of variables -composeSpray :: (AlgRing.C a, Eq a) => Spray a -> [Spray a] -> Spray a -composeSpray p = evalSpray (identify p) - --- | Create a spray from 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 - -prettyPowers :: String -> [Int] -> Text -prettyPowers var pows = append (pack x) (cons '(' $ snoc string ')') - where - x = " " ++ var ++ "^" - string = intercalate (pack ", ") (map (pack . show) pows) - --- | Pretty form of a spray -prettySpray :: (a -> String) -> String -> Spray a -> String -prettySpray prettyCoef var p = unpack $ intercalate (pack " + ") stringTerms - where - stringTerms = map stringTerm (sortBy (compare `on` fexpts) (HM.toList p)) - fexpts term = exponents $ fst term - stringTerm term = append - (snoc (snoc (cons '(' $ snoc stringCoef ')') ' ') '*') - (prettyPowers var pows) - where - pows = DF.toList $ exponents (fst term) - stringCoef = pack $ prettyCoef (snd term) - --- | Terms of a spray -sprayTerms :: Spray a -> HashMap (Seq Int) a -sprayTerms = HM.mapKeys exponents - --- | Spray as list -toList :: Spray a -> [([Int], a)] -toList p = HM.toList $ HM.mapKeys (DF.toList . exponents) p - --- | Bombieri spray -bombieriSpray :: 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] - times k x = AlgAdd.sum (replicate k x) +{-| +Module : Math.Algebra.Hspray +Description : Multivariate polynomials on a ring. +Copyright : (c) Stéphane Laurent, 2022-2024 +License : GPL-3 +Maintainer : laurent_step@outlook.fr + +Deals with multivariate polynomials over a commutative ring. +See README for examples. +-} +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE InstanceSigs #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE GeneralizedNewtypeDeriving #-} +{-# LANGUAGE TypeFamilies #-} + +module Math.Algebra.Hspray + ( + -- * Classes + FunctionLike (..) + , isConstant + , isUnivariate + , isBivariate + , isTrivariate + -- * Main types + , Exponents + , Powers (..) + , Spray + , QSpray + , QSpray' + , Term + -- * Basic sprays + , lone + , qlone + , lone' + , qlone' + , monomial + , qmonomial + , unitSpray + , zeroSpray + , constantSpray + -- * Showing a spray + , prettySpray + , prettySpray' + , prettySpray'' + , prettySprayXYZ + , prettySprayX1X2X3 + , showSpray + , showSprayXYZ + , showSprayXYZ' + , showSprayX1X2X3 + , showSprayX1X2X3' + , showNumSpray + , showQSpray + , showQSpray' + , prettyNumSprayX1X2X3 + , prettyQSprayX1X2X3 + , prettyQSprayX1X2X3' + , prettyNumSprayXYZ + , prettyQSprayXYZ + , prettyQSprayXYZ' + , prettyNumSpray + , prettyNumSpray' + , prettyQSpray + , prettyQSpray'' + , prettyQSpray' + , prettyQSpray''' + -- * Univariate polynomials and fractions of univariate polynomials + -- + -- | The univariate polynomials and the fractions of univariate polynomials + -- are used to defined the `OneParameterSpray` sprays, which represent + -- multivariate polynomials depending on a unique parameter. These sprays + -- lost their interest in version 0.4.0.0 (see CHANGELOG or README). + , A (..) + , Rational' + , Polynomial + , QPolynomial + , RatioOfPolynomials + , RatioOfQPolynomials + , prettyRatioOfPolynomials + , prettyRatioOfQPolynomials + , constPoly + , polyFromCoeffs + , soleParameter + , constQPoly + , qpolyFromCoeffs + , qsoleParameter + , evalRatioOfPolynomials + -- * One-parameter sprays + -- + -- | The `OneParameterSpray` sprays represent multivariate polynomials with + -- coefficients depending on a unique parameter. These sprays lost their + -- interest in version 0.4.0.0 (see CHANGELOG or README). One can use the + -- more general `ParametricSpray` sprays instead. + , OneParameterSpray + , OneParameterQSpray + , prettyOneParameterSprayX1X2X3 + , prettyOneParameterSprayXYZ + , prettyOneParameterSpray + , prettyOneParameterSpray' + , prettyOneParameterQSprayX1X2X3 + , prettyOneParameterQSprayXYZ + , prettyOneParameterQSpray + , prettyOneParameterQSpray' + , evalOneParameterSpray + , substituteTheParameter + , evalOneParameterSpray' + , evalOneParameterSpray'' + -- * Ratios of sprays + -- + -- | An object of type `RatioOfSprays` represents a fraction of two + -- multivariate polynomials. + , RatioOfSprays (..) + , RatioOfQSprays + , (%:%) + , (%//%) + , (^/^) + , (%/%) + , isConstantRatioOfSprays + , isPolynomialRatioOfSprays + , zeroRatioOfSprays + , zeroROS + , unitRatioOfSprays + , unitROS + , constantRatioOfSprays + , asRatioOfSprays + , evalRatioOfSprays + , evalRatioOfSprays' + , substituteRatioOfSprays + , fromRatioOfPolynomials + , fromRatioOfQPolynomials + , showRatioOfSprays + , showRatioOfNumSprays + , showRatioOfQSprays + , showRatioOfSpraysXYZ + , showRatioOfSpraysXYZ' + , showRatioOfSpraysX1X2X3 + , showRatioOfSpraysX1X2X3' + , prettyRatioOfQSpraysXYZ + , prettyRatioOfQSpraysX1X2X3 + , prettyRatioOfQSprays + , prettyRatioOfQSprays' + , prettyRatioOfNumSpraysXYZ + , prettyRatioOfNumSpraysX1X2X3 + , prettyRatioOfNumSprays + , prettyRatioOfNumSprays' + -- * Parametric sprays + -- + -- | There are three types of parametric sprays: @OneParameterSpray@, + -- @SimpleParametricSpray@ and @ParametricSpray@. These are sprays of + -- type @Spray b@ where @b@ has the class @FunctionLike@. When we say + -- \"parametric spray\" in the documentation, we mean either + -- such a spray or more precisely a @ParametricSpray@ spray. + , SimpleParametricSpray + , SimpleParametricQSpray + , ParametricSpray + , ParametricQSpray + , canCoerceToSimpleParametricSpray + , asSimpleParametricSprayUnsafe + , asSimpleParametricSpray + , fromOneParameterSpray + , fromOneParameterQSpray + , parametricSprayToOneParameterSpray + , parametricQSprayToOneParameterQSpray + , gegenbauerPolynomial + , jacobiPolynomial + , numberOfParameters + , changeParameters + , substituteParameters + , evalParametricSpray + , evalParametricSpray' + , prettyParametricQSprayABCXYZ + , prettyParametricQSpray + , prettyParametricNumSprayABCXYZ + , prettyParametricNumSpray + , prettySimpleParametricQSprayABCXYZ + , prettySimpleParametricQSpray + , prettySimpleParametricNumSprayABCXYZ + , prettySimpleParametricNumSpray + -- * Queries on a spray + , getCoefficient + , getConstantTerm + , isZeroSpray + , isConstantSpray + , isHomogeneousSpray + , allExponents + , allCoefficients + -- * Evaluation of a spray + , evalSpray + , substituteSpray + , composeSpray + , evalSpraySpray + -- * Division of sprays + , sprayDivision + , sprayDivisionRemainder + , pseudoDivision + -- * Gröbner basis + , groebnerBasis + , reduceGroebnerBasis + -- * Symmetric polynomials + , esPolynomial + , psPolynomial + , isSymmetricSpray + -- * Resultant, subresultants, and Sturm-Habicht sequence + , resultant + , resultant' + , resultant1 + , subresultants + , subresultants1 + , polynomialSubresultants + , sturmHabichtSequence + , principalSturmHabichtSequence + -- * Number of real roots of a univariate spray. These functions can be very + -- slow if the degree of the spray is not small. + , numberOfRealRoots + , numberOfRealRoots' + , numberOfRealRootsInOpenInterval + , numberOfRealRootsInOpenInterval' + , numberOfRealRootsInClosedInterval + , numberOfRealRootsInClosedInterval' + -- * Greatest common divisor + , gcdSpray + -- * Matrices + , detLaplace + , detLaplace' + , characteristicPolynomial + -- * Miscellaneous + , sumOfSprays + , productOfSprays + , (.^) + , (/>) + , (/^) + , fromList + , toList + , fromRationalSpray + , isPolynomialOf + , bombieriSpray + , collinearSprays + ) where +import qualified Algebra.Absolute as AlgAbs +import qualified Algebra.Additive as AlgAdd +import qualified Algebra.Differential as AlgDiff +import qualified Algebra.Field as AlgField +import qualified Algebra.Module as AlgMod +import qualified Algebra.RightModule as AlgRightMod +import qualified Algebra.Ring as AlgRing +import qualified Algebra.ZeroTestable as AlgZT +import qualified Data.Foldable as DF +import Data.Function ( on ) +import Data.HashMap.Strict ( HashMap ) +import qualified Data.HashMap.Strict as HM +import Data.Hashable ( Hashable ( hashWithSalt ) ) +import qualified Data.IntMap.Strict as IM +import Data.List ( sortBy + , maximumBy + , (\\) + , findIndices + , elemIndices + , nub + , foldl1' + , foldl' + , uncons + ) +import Data.List.Extra ( allSame + , unsnoc + ) +import Data.Matrix ( Matrix + , fromLists + , minorMatrix + , nrows + , ncols + , submatrix + ) +import qualified Data.Matrix as DM +import Data.Maybe ( isJust + , isNothing + , fromJust + , fromMaybe + ) +import Data.Ord ( comparing ) +import qualified Data.Ratio as DR +import qualified GHC.Real as DR +import qualified Data.Sequence as S +import Data.Sequence ( (><) + , Seq + , dropWhileR + , (|>) + , index + , adjust + , fromFunction + ) +import Data.Text ( Text + , append + , cons + , intercalate + , pack + , snoc + , unpack + ) +import Data.Tuple.Extra ( both, first ) +import qualified MathObj.Matrix as MathMatrix +import qualified MathObj.Polynomial as MathPol +import Number.Ratio ( T ( (:%) ), (%) ) +import qualified Number.Ratio as NumberRatio +-- import qualified Algebra.PrincipalIdealDomain as AlgPID +-- import qualified Algebra.Units as AlgUnits +-- import qualified Algebra.IntegralDomain as AlgID + + +-- Classes -------------------------------------------------------------------- + +-- | A spray represents a multivariate polynomial so it like a function. We +-- introduce a class because it will be assigned to the ratios of sprays too. +class FunctionLike b where + + -- | Number of variables in a function-like object + numberOfVariables :: b -> Int + + -- | 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) + -- >>> x1 = lone 1 :: Spray Rational + -- >>> x2 = lone 2 :: Spray Rational + -- >>> x3 = lone 3 :: Spray Rational + -- >>> spray = f x1 x2 x3 + -- + -- prop> permuteVariables [3, 1, 2] spray == f x3 x1 x2 + permuteVariables :: + [Int] -- ^ permutation + -> b -- ^ function-like object whose variables will be permuted + -> b -- ^ the function-like object with permuted 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 -- ^ function-like object whose variables will be swapped + -> b -- ^ the function-like object with swapped variables + + -- 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 + -- >>> spray' = derivative 1 spray + -- >>> putStrLn $ prettyNumSpray spray' + -- 2 + derivative :: + Int -- ^ index of the variable of differentiation (starting at 1) + -> b -- ^ the object to be derivated + -> b -- ^ the derivated object + + -- | The type of the coefficients (this is @a@ for both @Spray a@ and @RatioOfSprays a@) + type family BaseRing b + + -- | The type of the variables (this is @Spray a@ for both @Spray a@ and @RatioOfSprays a@) + type family VariablesType b + + infixl 6 ^+^ + -- | Addition of two function-like objects + (^+^) :: (AlgAdd.C b) => b -> b -> b + (^+^) = (AlgAdd.+) + + infixl 6 ^-^ + -- | Substraction of two function-like objects + (^-^) :: (AlgAdd.C b) => b -> b -> b + (^-^) = (AlgAdd.-) + + infixl 7 ^*^ + -- | Multiplication of two function-like objects + (^*^) :: (AlgRing.C b) => b -> b -> b + (^*^) = (AlgRing.*) + + infixr 8 ^**^ + -- | Power of a function-like object + (^**^) :: (AlgRing.C b) => b -> Int -> b + (^**^) p k = p AlgRing.^ fromIntegral k + + infixr 7 *^ + -- | Multiply a function-like object by a scalar + (*^) :: BaseRing b -> b -> b + + infixl 6 +> + -- | Add a function-like object to a constant + -- + -- prop> x +> spray == constantSpray x ^+^ spray + (+>) :: BaseRing b -> b -> b + + infixr 6 <+ + -- | 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) of a + -- function-like object + -- + -- >>> x = lone 1 :: Spray Int + -- >>> y = lone 2 :: Spray Int + -- >>> spray = 2*^x^**^2 ^-^ 3*^y + -- >>> evaluate spray [2, 1] + -- 5 + evaluate :: + b -- ^ function-like object to be evaluated, e.g. a spray + -> [BaseRing b] -- ^ list of values to be substituted to its variables + -> BaseRing b + + -- | Flipped version of @evaluate@ + -- + -- prop> evaluateAt [2, 1] spray == evaluate spray [2, 1] + evaluateAt :: [BaseRing b] -> b -> BaseRing b + evaluateAt = flip evaluate + + -- | Partial evaluation of a function-like object (replace some variables + -- with some values) + -- + -- >>> x1 = lone 1 :: Spray Int + -- >>> x2 = lone 2 :: Spray Int + -- >>> x3 = lone 3 :: Spray Int + -- >>> spray = x1^**^2 ^-^ x2 ^+^ x3 ^-^ unitSpray + -- >>> spray' = substitute [Just 2, Nothing, Just 3] spray + -- >>> putStrLn $ prettyNumSprayX1X2X3 "x" spray' + -- -x2 + 6 + substitute :: + [Maybe (BaseRing b)] -- ^ @Just x@ to replace the variable with @x@, @Nothing@ for no replacement + -> b -- ^ function-like object to be partially evaluated + -> b + + -- | Polynomial change of variables of a function-like object + -- + -- >>> x = lone 1 :: Spray Int + -- >>> y = lone 2 :: Spray Int + -- >>> spray = x ^*^ y + -- >>> spray' = changeVariables spray [x ^+^ y, x ^-^ y] + -- >>> putStrLn $ prettyNumSpray' spray' + -- X^2 - Y^2 + changeVariables :: + b -- ^ function-like object such as a spray + -> [VariablesType b] -- ^ list of new variables + -> b + +-- | Whether a function-like object has a constant value +isConstant :: FunctionLike b => b -> Bool +isConstant f = numberOfVariables f == 0 + +-- | Whether a function-like object represents an univariate function; it is considered +-- that it is univariate if it is constant +isUnivariate :: FunctionLike b => b -> Bool +isUnivariate f = numberOfVariables f <= 1 + +-- | Whether a function-like object represents a bivariate function; it is considered +-- that it is bivariate if it is univariate +isBivariate :: FunctionLike b => b -> Bool +isBivariate f = numberOfVariables f <= 2 + +-- | Whether a function-like object represents a trivariate function; it is considered +-- that it is trivariate if it is bivariate +isTrivariate :: FunctionLike b => b -> Bool +isTrivariate f = numberOfVariables f <= 3 + + +-- Additional operations to 'numeric-prelude' + +infixr 7 /> +-- | Divides by a scalar in a module over a field +(/>) :: (AlgField.C k, AlgMod.C k a) => a -> k -> a +x /> lambda = AlgField.recip lambda AlgMod.*> x + +infixr 7 .^ +-- | Scale by an integer (I do not find this operation in __numeric-prelude__) +-- +-- prop> 3 .^ x == x Algebra.Additive.+ x Algebra.Additive.+ x +(.^) :: (AlgAdd.C a, Eq a) => Int -> a -> a +k .^ x = if k >= 0 + then powerOperation (AlgAdd.+) AlgAdd.zero x k + else (.^) (-k) (AlgAdd.negate x) + where + powerOperation op = + let go acc _ 0 = acc + go acc a n = go (if even n then acc else op acc a) (op a a) (div n 2) + in go + + +-- Univariate polynomials and ratios of univariate polynomials ---------------- + +-- | The new type @A a@ is used to attribute some instances to the +-- type @Polynomial a@; it is needed to avoid orphan instances. +newtype A a = A a + deriving + (Eq, Show, AlgAdd.C, AlgRing.C, AlgField.C) + +-- | The type @Rational'@ is used to introduce the univariate polynomials +-- with rational coefficients (`QPolynomial`). It is similar to the well-known +-- type @Rational@ (actually these two types are the same but @Rational'@ has +-- more instances and we need them for the univariate polynomials). +type Rational' = NumberRatio.Rational + +-- | The type @Polynomial a@ is used to represent univariate polynomials. +type Polynomial a = MathPol.T (A a) +type QPolynomial = Polynomial Rational' + +-- The type @RatioOfPolynomials a@ is used to represent fractions of +-- univariate polynomials. +type RatioOfPolynomials a = NumberRatio.T (Polynomial a) +type RatioOfQPolynomials = RatioOfPolynomials Rational' + +instance (Eq a, AlgField.C a) => FunctionLike (Polynomial a) where + -- + numberOfVariables :: Polynomial a -> Int + numberOfVariables p = case MathPol.degree p of + Nothing -> 0 + Just d -> min 1 d + -- + type BaseRing (Polynomial a) = a + -- + type VariablesType (Polynomial a) = Polynomial a + -- + (*^) :: a -> Polynomial a -> Polynomial a + (*^) lambda pol = constPoly lambda AlgRing.* pol + -- + (+>) :: a -> Polynomial a -> Polynomial a + (+>) x pol = constPoly x AlgAdd.+ pol + -- + evaluate :: Polynomial a -> [a] -> a + evaluate p xs = get (MathPol.evaluate p (A (xs !! 0))) + where + get (A x) = x + -- + substitute :: [Maybe a] -> Polynomial a -> Polynomial a + substitute x p = + if isNothing (x !! 0) + then p + else constPoly (evaluate p [fromJust $ x !! 0]) + -- + permuteVariables :: [Int] -> Polynomial a -> Polynomial a + permuteVariables = error "permuteVariables: there is only one variable." + -- + 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 + then AlgDiff.differentiate p + else constPoly AlgAdd.zero + -- + changeVariables :: Polynomial a -> [Polynomial a] -> Polynomial a + changeVariables p ps = MathPol.compose p (ps !! 0) + +instance (Eq a, AlgField.C a) => FunctionLike (RatioOfPolynomials a) where + numberOfVariables :: RatioOfPolynomials a -> Int + numberOfVariables (p :% q) = + max (numberOfVariables p) (numberOfVariables q) + -- + type BaseRing (RatioOfPolynomials a) = a + -- + type VariablesType (RatioOfPolynomials a) = Polynomial a + -- + (*^) :: a -> RatioOfPolynomials a -> RatioOfPolynomials a + (*^) lambda rOP = A lambda AlgMod.*> rOP + -- + (+>) :: a -> RatioOfPolynomials a -> RatioOfPolynomials a + (+>) x rOP = NumberRatio.fromValue (constPoly x) AlgAdd.+ rOP + -- + evaluate :: RatioOfPolynomials a -> [a] -> a + evaluate r xs = evaluate (NumberRatio.numerator r) xs AlgField./ + evaluate (NumberRatio.denominator r) xs + -- + substitute :: [Maybe a] -> RatioOfPolynomials a -> RatioOfPolynomials a + substitute x r = + if isNothing (x !! 0) + then r + else substitute x (NumberRatio.numerator r) % + substitute x (NumberRatio.denominator r) + -- + permuteVariables :: [Int] -> RatioOfPolynomials a -> RatioOfPolynomials a + permuteVariables = error "permuteVariables: there is only one variable." + -- + 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 + then + (p' AlgRing.* q AlgAdd.- p AlgRing.* q') % q AlgRing.^ 2 + else constPoly AlgAdd.zero :% constPoly AlgRing.one + where + p = NumberRatio.numerator r + q = NumberRatio.denominator r + p' = AlgDiff.differentiate p + q' = AlgDiff.differentiate q + -- + changeVariables :: + RatioOfPolynomials a -> [Polynomial a] -> RatioOfPolynomials a + changeVariables r ps = changeVariables (NumberRatio.numerator r) ps % + changeVariables (NumberRatio.denominator r) ps + +instance (Eq a, AlgField.C a) => AlgZT.C (A a) where + isZero :: A a -> Bool + isZero (A r) = r == AlgAdd.zero + +instance (Eq a, AlgField.C a) => AlgMod.C (A a) (RatioOfPolynomials a) where + (*>) :: A a -> RatioOfPolynomials a -> RatioOfPolynomials a + lambda *> rop = NumberRatio.scale (MathPol.const lambda) rop + +instance (Eq a, AlgField.C a) => AlgRightMod.C (A a) (RatioOfPolynomials a) + where + (<*) :: RatioOfPolynomials a -> A a -> RatioOfPolynomials a + rop <* lambda = lambda AlgMod.*> rop + +instance (Eq a, AlgField.C a) => AlgMod.C a (RatioOfPolynomials a) where + (*>) :: a -> RatioOfPolynomials a -> RatioOfPolynomials a + lambda *> rop = A lambda AlgMod.*> rop + +instance (Eq a, AlgField.C a) => AlgRightMod.C a (RatioOfPolynomials a) where + (<*) :: RatioOfPolynomials a -> a -> RatioOfPolynomials a + rop <* lambda = lambda AlgMod.*> rop + +instance (Eq a, AlgField.C a) + => AlgMod.C (Polynomial a) (RatioOfPolynomials a) + where + (*>) :: Polynomial a -> RatioOfPolynomials a -> RatioOfPolynomials a + p *> r = NumberRatio.scale p r + +instance (Eq a, AlgField.C a) + => AlgRightMod.C (Polynomial a) (RatioOfPolynomials a) + where + (<*) :: RatioOfPolynomials a -> Polynomial a -> RatioOfPolynomials a + r <* p = p AlgMod.*> r + +-- | Constant univariate polynomial +constPoly :: a -> Polynomial a +constPoly x = MathPol.const (A x) + +-- | Univariate polynomial from its coefficients (ordered by increasing degrees) +polyFromCoeffs :: [a] -> Polynomial a +polyFromCoeffs as = MathPol.fromCoeffs (map A as) + +-- | The variable of a univariate polynomial; it is called \"soleParameter\" because +-- this it represents the parameter of a `OneParameterSpray` spray +soleParameter :: AlgRing.C a => Polynomial a +soleParameter = polyFromCoeffs [AlgAdd.zero, AlgRing.one] + +-- | Constant rational univariate polynomial +-- +-- >>> import Number.Ratio ( (%) ) +-- >>> constQPoly (2 % 3) +-- +-- prop> constQPoly (2 % 3) == qpolyFromCoeffs [2 % 3] +constQPoly :: Rational' -> QPolynomial +constQPoly = constPoly + +-- | Rational univariate polynomial from coefficients +-- +-- >>> import Number.Ratio ( (%) ) +-- >>> qpolyFromCoeffs [2 % 3, 5, 7 % 4] +qpolyFromCoeffs :: [Rational'] -> QPolynomial +qpolyFromCoeffs = polyFromCoeffs + +-- | The variable of a univariate rational polynomial; it is called \"qsoleParameter\" +-- because it represents the parameter of a `OneParameterQSpray` spray +-- +-- prop> qsoleParameter == qpolyFromCoeffs [0, 1] +qsoleParameter :: QPolynomial +qsoleParameter = qpolyFromCoeffs [0, 1] + +{- +-- show a ratio, helper function +showQ :: (Eq a, Num a, Show a) => NumberRatio.T a -> String +showQ q = if d == 1 + then show n + else show n ++ "/" ++ show d + where + n = NumberRatio.numerator q + d = NumberRatio.denominator q + -} + +-- | identify a `Polynomial a` to a `Spray a`, in order to apply the show spray +-- functions to the univariate polynomials +polynomialToSpray :: forall a. (Eq a, AlgRing.C a) => Polynomial a -> Spray a +polynomialToSpray pol = sumOfSprays terms + where + coeffs = MathPol.coeffs pol + indices = findIndices (/= A AlgAdd.zero) coeffs + get :: A a -> a + get (A x) = x + terms = [get (coeffs !! i) *^ lone' 1 i | i <- indices] + +qPolynomialToQSpray :: QPolynomial -> QSpray +qPolynomialToQSpray pol = sumOfSprays terms + where + coeffs = MathPol.coeffs pol + indices = findIndices (/= A 0) coeffs + get :: A Rational' -> Rational + get (A x) = NumberRatio.numerator x DR.:% NumberRatio.denominator x + terms = [get (coeffs !! i) *^ qlone' 1 i | i <- indices] + +-- helper function; it encloses a string between two given delimiters +bracify :: (String, String) -> String -> String +bracify (lbrace, rbrace) x = lbrace ++ x ++ rbrace + +-- | helper function for prettyRatioOfPolynomials (and prettyOneParameterSpray) +showRatioOfPolynomials :: forall a. (Eq a, AlgField.C a) + => (Spray a -> String) -> RatioOfPolynomials a -> String +showRatioOfPolynomials sprayShower polysRatio = + numeratorString ++ denominatorString + where + numerator = NumberRatio.numerator polysRatio + denominator = NumberRatio.denominator polysRatio + brackets = denominator /= MathPol.const (A AlgRing.one) + enclose = bracify ("[ ", " ]") + numeratorString = if brackets + then enclose (sprayShower (polynomialToSpray numerator)) + else sprayShower (polynomialToSpray numerator) + denominatorString = if not brackets + then "" + else " %//% " ++ enclose (sprayShower (polynomialToSpray denominator)) + +-- | Pretty form of a ratio of univariate polynomials with rational coefficients +prettyRatioOfQPolynomials :: + String -- ^ a string to denote the variable, e.g. @\"a\"@ + -> RatioOfQPolynomials + -> String +prettyRatioOfQPolynomials var = showRatioOfPolynomials (prettyQSprayXYZ' [var]) + +-- | helper function for prettyRatioOfPolynomials (and prettyOneParameterSpray) +showQpol :: forall a. (Eq a, AlgField.C a) + => Polynomial a -> String -> (a -> String) -> Bool -> String +showQpol pol variable showCoeff brackets = if brackets + then "[ " ++ polyString ++ " ]" + else polyString + where + showCoeff' :: Int -> A a -> String + showCoeff' i (A coeff) = case i of + 0 -> (bracify ("(", ")") . showCoeff) coeff + _ -> if coeff == AlgRing.one + then "" + else (bracify ("(", ")") . showCoeff) coeff + coeffs = MathPol.coeffs pol + nonzeros = findIndices (/= A AlgAdd.zero) coeffs + terms = map (pack . showTerm) nonzeros + where + showTerm i = case i of + 0 -> showCoeff' 0 (coeffs !! 0) + 1 -> showCoeff' 1 (coeffs !! 1) ++ variable + _ -> showCoeff' i (coeffs !! i) ++ variable ++ "^" ++ show i + polyString = unpack (intercalate (pack " + ") terms) + +-- | helper function for prettyRatioOfPolynomials (and prettyOneParameterSpray) +showQpolysRatio :: forall a. (Eq a, AlgField.C a) + => String -> (a -> String) -> RatioOfPolynomials a -> String +showQpolysRatio var showCoeff polysRatio = numeratorString ++ denominatorString + where + denominator = NumberRatio.denominator polysRatio + brackets = denominator /= constPoly AlgRing.one + numeratorString = + showQpol (NumberRatio.numerator polysRatio) var showCoeff brackets + denominatorString = if not brackets + then "" + else " %//% " ++ showQpol denominator var showCoeff True + +-- | Pretty form of a ratio of univariate polynomials +prettyRatioOfPolynomials :: (Eq a, AlgField.C a, Show a) + => String -- ^ string (usually a single letter) to denote the variable, e.g. @\"a\"@ + -> RatioOfPolynomials a + -> String +prettyRatioOfPolynomials var = showQpolysRatio var show + +{- -- | Pretty form of a ratio of univariate qpolynomials +prettyRatioOfQPolynomials' + :: String -- ^ a string to denote the variable, e.g. @\"a\"@ + -> RatioOfQPolynomials + -> String +prettyRatioOfQPolynomials' var = showQpolysRatio var showQ + -} + +-- | Evaluates a ratio of univariate polynomials +evalRatioOfPolynomials :: AlgField.C a + => a -- ^ the value at which the evaluation is desired + -> RatioOfPolynomials a + -> a +evalRatioOfPolynomials value polysRatio = + resultNumerator AlgField./ resultDenominator + where + A resultNumerator = + MathPol.evaluate (NumberRatio.numerator polysRatio) (A value) + A resultDenominator = + MathPol.evaluate (NumberRatio.denominator polysRatio) (A value) + + +-- One-parameter sprays ------------------------------------------------------- + +type OneParameterSpray a = Spray (RatioOfPolynomials a) +type OneParameterQSpray = OneParameterSpray Rational' + +{- -- | simplifies a ratio of polynomials (simply by multiplying it by one) +simplifyRatioOfPolynomials :: + (Eq a, AlgField.C a) => RatioOfPolynomials a -> RatioOfPolynomials a +simplifyRatioOfPolynomials = (AlgRing.*) AlgRing.one + +-- | Simplifies the coefficients (the fractions of univariate polynomials) of a +-- one-parameter spray +simplifyOneParameterSpray :: + (Eq a, AlgField.C a) => OneParameterSpray a -> OneParameterSpray a +simplifyOneParameterSpray = HM.map simplifyRatioOfPolynomials + -} + +instance (Eq a, AlgField.C a) => AlgMod.C (Polynomial a) (OneParameterSpray a) where + (*>) :: Polynomial a -> OneParameterSpray a -> OneParameterSpray a + p *> spray = (p NumberRatio.:% AlgRing.one) *^ spray + +instance (Eq a, AlgField.C a) => AlgRightMod.C (Polynomial a) (OneParameterSpray a) where + (<*) :: OneParameterSpray a -> Polynomial a -> OneParameterSpray a + spray <* p = p AlgMod.*> spray + +instance (Eq a, AlgField.C a) => AlgMod.C a (OneParameterSpray a) where + (*>) :: a -> OneParameterSpray a -> OneParameterSpray a + lambda *> spray = MathPol.const (A lambda) AlgMod.*> spray + +instance (Eq a, AlgField.C a) => AlgRightMod.C a (OneParameterSpray a) where + (<*) :: OneParameterSpray a -> a -> OneParameterSpray a + spray <* lambda = lambda AlgMod.*> spray + +-- | Pretty form of a one-parameter spray, using a string (typically a letter) +-- followed by an index to denote the variables +prettyOneParameterSprayX1X2X3 :: + (Eq a, Show a, AlgField.C a) + => String -- ^ string to denote the parameter of the spray, e.g. @\"a\"@ + -> String -- ^ typically a letter, to denote the non-indexed variables + -> OneParameterSpray a -- ^ a one-parameter spray; note that this function does not simplify it + -> String +prettyOneParameterSprayX1X2X3 a = showSprayX1X2X3 (prettyRatioOfPolynomials a) ("{ ", " }") + +-- | Pretty form of a one-parameter spray, using some given strings (typically some +-- letters) to denote the variables if possible, i.e. if enough letters are +-- provided; otherwise this function behaves exactly like +-- @prettyOneParameterSprayX1X2X3 a@ where @a@ is the first provided letter +prettyOneParameterSprayXYZ :: + (Eq a, Show a, AlgField.C a) + => String -- ^ string to denote the parameter of the spray, e.g. @\"a\"@ + -> [String] -- ^ typically some letters, to denote the main variables + -> OneParameterSpray a -- ^ a one-parameter spray; note that this function does not simplify it + -> String +prettyOneParameterSprayXYZ a = showSprayXYZ (prettyRatioOfPolynomials a) ("{ ", " }") + +-- | Pretty form of a one-parameter spray; see the definition below and see +-- `prettyOneParameterSprayXYZ` +-- +-- prop> prettyOneParameterSpray a spray == prettyOneParameterSprayXYZ a ["x","y","z"] spray +prettyOneParameterSpray :: + (Eq a, Show a, AlgField.C a) + => String -- ^ string to denote the parameter of the spray, e.g. @\"a\"@ + -> OneParameterSpray a -- ^ a one-parameter spray; note that this function does not simplify it + -> String +prettyOneParameterSpray a = prettyOneParameterSprayXYZ a ["x", "y", "z"] + +-- | Pretty form of a one-parameter spray; see the definition below and see +-- `prettyOneParameterSprayXYZ` +-- +-- prop> prettyOneParameterSpray' a spray == prettyOneParameterSprayXYZ a ["X","Y","Z"] spray +prettyOneParameterSpray' :: + (Eq a, Show a, AlgField.C a) + => String -- ^ string to denote the parameter of the spray, e.g. @\"a\"@ + -> OneParameterSpray a -- ^ a one-parameter spray; note that this function does not simplify it + -> String +prettyOneParameterSpray' a = prettyOneParameterSprayXYZ a ["X", "Y", "Z"] + +-- | Pretty form of a one-parameter rational spray, using a string (typically a letter) +-- followed by an index to denote the variables +prettyOneParameterQSprayX1X2X3 :: + String -- ^ usually a letter, to denote the parameter of the spray, e.g. @\"a\"@ + -> String -- ^ usually a letter, to denote the non-indexed variables of the spray + -> OneParameterQSpray -- ^ a one-parameter rational spray; note that this function does not simplify it + -> String +prettyOneParameterQSprayX1X2X3 a x = + showSpray (prettyRatioOfQPolynomials a) ("{ ", " }") (showMonomialsX1X2X3 x) + +-- | Pretty form of a one-parameter rational spray, using some given strings (typically some +-- letters) to denote the variables if possible, i.e. if enough letters are +-- provided; otherwise this function behaves exactly like +-- @prettyOneParameterQSprayX1X2X3 a@ where @a@ is the first provided letter +prettyOneParameterQSprayXYZ :: + String -- ^ usually a letter, to denote the parameter of the spray, e.g. @\"a\"@ + -> [String] -- ^ usually some letters, to denote the variables of the spray + -> OneParameterQSpray -- ^ a one-parameter rational spray; note that this function does not simplify it + -> String +prettyOneParameterQSprayXYZ a letters = + showSpray (prettyRatioOfQPolynomials a) ("{ ", " }") (showMonomialsXYZ letters) + +-- | Pretty form of a one-parameter rational spray, using @\"x\"@, @\"y\"@ and @\"z\"@ for the variables +-- if possible; i.e. if the spray does not have more than three variables, otherwise +-- @\"x1\"@, @\"x2\"@, ... are used to denote the variables +-- +-- prop> prettyOneParameterQSpray a == prettyOneParameterQSprayXYZ a ["x","y","z"] +prettyOneParameterQSpray :: + String -- ^ usually a letter, to denote the parameter of the spray, e.g. @\"a\"@ + -> OneParameterQSpray -- ^ the one-parameter rational spray to be printed; note that this function does not simplify it + -> String +prettyOneParameterQSpray a = prettyOneParameterQSprayXYZ a ["x", "y", "z"] + +-- | Pretty form of a one-parameter rational spray, using @\"X\"@, @\"Y\"@ and @\"Z\"@ for the variables +-- if possible; i.e. if the spray does not have more than three variables, otherwise +-- @\"X1\"@, @\"X2\"@, ... are used +-- +-- prop> prettyOneParameterQSpray' a == prettyOneParameterQSprayXYZ a ["X","Y","Z"] +prettyOneParameterQSpray' :: + String -- ^ usually a letter, to denote the parameter of the spray, e.g. @\"a\"@ + -> OneParameterQSpray -- ^ the one-parameter rational spray to be printed; note that this function does not simplify it + -> String +prettyOneParameterQSpray' a = prettyOneParameterQSprayXYZ a ["X", "Y", "Z"] + +-- | Substitutes a value to the parameter of a one-parameter spray +-- (the variable occurring in its coefficients) +-- +-- prop> evalOneParameterSpray spray x == substituteParameters spray [x] +evalOneParameterSpray :: + (Eq a, AlgField.C a) => OneParameterSpray a -> a -> Spray a +evalOneParameterSpray spray x = + removeZeroTerms $ HM.map (evalRatioOfPolynomials x) spray + +-- | Substitutes a value to the parameter of a one-parameter spray; +-- same as `evalOneParameterSpray` +-- +-- prop> substituteTheParameter spray x == substituteParameters spray [x] +substituteTheParameter :: + (Eq a, AlgField.C a) => OneParameterSpray a -> a -> Spray a +substituteTheParameter spray x = substituteParameters spray [x] + +-- | Substitutes a value to the parameter of a one-parameter spray as well +-- as some values to the variables of this spray +-- +-- prop> evalOneParameterSpray' spray a xs == evalParametricSpray' spray [a] xs +evalOneParameterSpray' :: (Eq a, AlgField.C a) + => OneParameterSpray a -- ^ one-parameter spray to be evaluated + -> a -- ^ a value for the parameter + -> [a] -- ^ some values for the variables + -> a +evalOneParameterSpray' spray x xs = if length xs >= numberOfVariables spray + then evalSpray (evalOneParameterSpray spray x) xs + else error "evalOneParameterSpray': not enough values provided." + +-- | helper function for evalOneParameterSpray'' +evalOneParameterTerm :: (Eq a, AlgField.C a) + => [a] -> Term (RatioOfPolynomials a) -> RatioOfPolynomials a +evalOneParameterTerm xs (powers, coeff) = + AlgRing.product (zipWith (AlgRing.^) xs pows) *^ coeff + where + pows = DF.toList (fromIntegral <$> exponents powers) + +-- | Substitutes some values to the variables of a one-parameter spray; same +-- as `evalParametricSpray` +evalOneParameterSpray'' :: + (Eq a, AlgField.C a) => OneParameterSpray a -> [a] -> RatioOfPolynomials a +evalOneParameterSpray'' spray xs = if length xs >= numberOfVariables spray + then AlgAdd.sum $ map (evalOneParameterTerm xs) (HM.toList spray) + else error "evalOneParameterSpray'': not enough values provided." + + +-- Sprays --------------------------------------------------------------------- + +type Exponents = Seq Int + +-- | The type @Powers@ is used to represent the exponents of the monomial +-- occurring in a term of a spray. The integer in the field @nvariables@ is +-- the number of variables involved in this monomial (it is 3, not 2, for a +-- monomial such as @x^2.z^3@, because the exponents of this monomial is the +-- sequence @(2, 0, 3)@). Actually this integer is always the length of the +-- sequence in the field @exponents@. The reason of the presence of the field +-- @nvariables@ is that I thought that it was necessary when I started to +-- develop the package, but now I think it is useless. The type @Powers@ will +-- possibly be abandoned in a future version of the package. However we cannot +-- simply use the type `Exponents` to represent the exponents, because two +-- sequences of exponents that differ only by some trailing zeros must be +-- considered as identical, and they are considered as such with the type +-- @Powers@ thanks to its @Eq@ instance. Instead of @Powers@, a new type +-- encapsulating the `Exponents` type with such an @Eq@ instance should be +-- enough to represent the exponents. +data Powers = Powers + { exponents :: Exponents + , nvariables :: Int + } + deriving Show + +nullPowers :: Powers +nullPowers = Powers S.empty 0 + +powerize :: Exponents -> Powers +powerize expnts = Powers expnts (S.length expnts) + +instance Eq Powers where + (==) :: Powers -> Powers -> Bool + pows1 == pows2 = expts1' == expts2' + where + (expts1', expts2') = harmonize (pows1, pows2) + +instance Hashable Powers where + hashWithSalt :: Int -> Powers -> Int + hashWithSalt k pows = hashWithSalt k (exponents pows, nvariables pows) + +-- | append trailing zeros +growSequence :: Exponents -> Int -> Int -> Exponents +growSequence s m n = s >< S.replicate (n - m) 0 + +growSequence' :: Int -> Exponents -> Exponents +growSequence' n s = growSequence s (S.length s) n + +-- | append trailing zeros to get the same length +harmonize :: (Powers, Powers) -> (Exponents, Exponents) +harmonize (pows1, pows2) = (e1', e2') + where + e1 = exponents pows1 + e2 = exponents pows2 + n1 = nvariables pows1 + n2 = nvariables pows2 + (e1', e2') = if n1 < n2 + then (growSequence e1 n1 n2, e2) + else (e1, growSequence e2 n2 n1) + +makePowers :: Exponents -> Powers +makePowers expnts = powerize s + where + s = dropWhileR (== 0) expnts + +-- | An object of type @Spray a@ represents a multivariate polynomial whose +-- coefficients are represented by the objects of type @a@, which must have +-- a ring instance in order that we can add and multiply two polynomials. +type Spray a = HashMap Powers a + +-- | Most often, one deals with sprays with rational coefficients, so we +-- dedicate a type alias for such sprays. +type QSpray = Spray Rational + +-- | The type `Rational'` is helpful when dealing with `OneParameterSpray` +-- sprays, but this type of sprays lost its interest in version 0.4.0.0 +-- (see CHANGELOG or README). +type QSpray' = Spray Rational' + +type SafeSpray a = HashMap Exponents a + +-- An object of type @Term a@ represents a term of a @Spray a@ spray. Applying +-- @Data.HashMap.Strict.toList@ to a @Spray a@ spray yields a list of @Term a@ +-- objects. This type has probably no interest for the user, it is exported +-- because it possibly has an interest for internal usage in a package using +-- __hspray__. +type Term a = (Powers, a) + +instance (AlgRing.C a, Eq a) => FunctionLike (Spray a) where + type BaseRing (Spray a) = a + -- + type VariablesType (Spray a) = Spray a + -- + (*^) :: a -> Spray a -> Spray a + (*^) lambda pol = lambda AlgMod.*> pol + -- + (+>) :: a -> Spray a -> Spray a + (+>) x spray = if x == AlgAdd.zero + then spray + else addTerm spray (nullPowers, x) + -- + evaluate :: Spray a -> [a] -> a + evaluate spray xyz = if length xyz >= numberOfVariables spray + then evalSprayHelper xyz spray + else error "evaluate: not enough values provided." + -- + substitute :: [Maybe a] -> Spray a -> Spray a + substitute subs spray = if length subs >= n + then spray' + else error "substitute: incorrect length of the substitutions list." + where + n = numberOfVariables spray + terms = HM.toList spray + spray' = sumTerms (map substituteTerm terms) + substituteTerm :: Term a -> Term a + substituteTerm (powers, coeff) = (powers'', coeff') + where + pows = exponents powers + nv = nvariables powers + indices = findIndices isJust (take nv 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'' = makePowers pows'' + -- + changeVariables :: Spray a -> [Spray a] -> Spray a + changeVariables = composeSpray + -- + numberOfVariables :: Spray a -> Int + numberOfVariables spray = + if null powers then 0 else maximum (map nvariables powers) + where + powers = HM.keys spray + -- + permuteVariables :: [Int] -> Spray a -> Spray a + permuteVariables permutation spray = + if isPermutation permutation && n' >= n + then spray' + else error "permuteVariables: invalid permutation." + where + spray' = if isConstant spray + then spray + else HM.fromList (zip powers' coeffs) + n = numberOfVariables spray + n' = maximum permutation + isPermutation pmtn = + (not . null) pmtn && minimum pmtn == 1 && length (nub pmtn) == 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 + (powers, coeffs) = unzip (HM.toList spray) + f pows = let expnts = (permuteSeq . growSequence' n') (exponents pows) in + makePowers expnts + powers' = map f powers + -- + swapVariables :: (Int, Int) -> Spray a -> Spray a + swapVariables (i, j) spray = + if i>=1 && j>=1 + then spray' + else error "swapVariables: invalid indices." + where + spray' = if isConstant spray + then spray + else HM.fromList (zip powers' coeffs) + n = maximum [numberOfVariables spray, i, j] + f k | k == i = j + | k == j = i + | otherwise = k + transposition = map f [1 .. n] + permuteSeq x = + S.mapWithIndex (\ii _ -> x `index` (transposition !! ii - 1)) x + (powers, coeffs) = unzip (HM.toList spray) + g pows = let expnts = (permuteSeq . growSequence' n) (exponents pows) in + 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 + else error "derivative: invalid index." + where + terms = [ derivTerm term | term <- HM.toList p ] + derivTerm :: Term a -> Term a + derivTerm (pows, coef) = if i' >= S.length expts + then (nullPowers, AlgAdd.zero) + else (pows', coef') + where + i' = i - 1 + expts = exponents pows + expt_i = expts `index` i' + expts' = adjust (subtract 1) i' expts + coef' = expt_i .^ coef + pows' = makePowers expts' + +-- | addition of two sprays +addSprays :: (AlgAdd.C a, Eq a) => Spray a -> Spray a -> Spray a +addSprays p q = removeZeroTerms $ HM.unionWith (AlgAdd.+) p q -- HM.foldlWithKey' f p q +-- where +-- f s powers coef = HM.insertWith (AlgAdd.+) powers coef s + +-- | addition of a term to a spray +addTerm :: (AlgAdd.C a, Eq a) => Spray a -> Term a -> Spray a +addTerm spray (powers, coeff) = + if getCoefficient' powers spray AlgAdd.+ coeff == AlgAdd.zero + then + HM.delete powers spray + else + HM.insertWith (AlgAdd.+) powers coeff spray + +-- | sum list of terms +sumTerms :: (Eq a, AlgAdd.C a) => [Term a] -> Spray a +sumTerms = removeZeroTerms . HM.fromListWith (AlgAdd.+) + +-- | 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 = removeZeroTerms $ HM.map (lambda AlgRing.*) p + +-- | multiply two terms +multTerm :: AlgRing.C a => Term a -> Term a -> Term a +multTerm (pows1, coef1) (pows2, coef2) = (pows, coef1 AlgRing.* coef2) + where + (expts1', expts2') = harmonize (pows1, pows2) + expts = S.zipWith (+) expts1' expts2' + pows = makePowers expts + +-- | multiply a spray by a term +multSprayByTerm :: (Eq a, AlgRing.C a) => Spray a -> Term a -> Spray a +multSprayByTerm spray term = removeZeroTerms $ HM.fromList prods + where + prods = [multTerm trm term | trm <- HM.toList spray] + +-- | multiply two sprays +multSprays :: (AlgRing.C a, Eq a) => Spray a -> Spray a -> Spray a +multSprays p q = removeZeroTerms $ HM.fromListWith (AlgAdd.+) prods + where + p' = HM.toList p + q' = HM.toList q + prods = [ multTerm mp mq | mp <- p', mq <- q' ] + +instance (AlgAdd.C a, Eq a) => AlgAdd.C (Spray a) where + (+) :: 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) => AlgRightMod.C a (Spray a) where + (<*) :: Spray a -> a -> Spray a + p <* lambda = 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 :: Spray a + one = HM.singleton (Powers S.empty 0) AlgRing.one + +infixr 7 /^ +-- | Divides a spray by a scalar; you can equivalently use `(/>)` if the type +-- of the scalar is not ambiguous +(/^) :: (AlgField.C a, Eq a) => Spray a -> a -> Spray a +(/^) spray lambda = spray /> lambda + +-- | remove zero terms of a spray +removeZeroTerms :: (AlgAdd.C a, Eq a) => Spray a -> Spray a +removeZeroTerms = HM.filter (/= AlgAdd.zero) + +-- | helper function for lone and lone' +lonePower :: Int -> Int -> Powers +lonePower n p = if n == 0 || p == 0 + then nullPowers + else Powers (S.replicate (n - 1) 0 |> p) n + +-- | The @n@-th polynomial variable @x_n@ as a spray; one usually builds a +-- spray by introducing these variables and combining them with the arithmetic +-- operations +-- +-- >>> x = lone 1 :: Spray Int +-- >>> y = lone 2 :: Spray Int +-- >>> spray = 2*^x^**^2 ^-^ 3*^y +-- >>> putStrLn $ prettyNumSpray spray +-- 2*x^2 - 3*y +-- +-- prop> lone 0 == unitSpray +lone :: AlgRing.C a => Int -> Spray a +lone n = if n >= 0 + then HM.singleton (lonePower n 1) AlgRing.one + else error "lone: invalid index." + +-- | The @n@-th polynomial variable for rational sprays; this is just a +-- specialization of `lone` +qlone :: Int -> QSpray +qlone = lone + +-- | The spray @x_n^p@; more efficient than exponentiating @lone n@ +-- +-- prop> lone' 2 10 = lone 2 ^**^ 10 +lone' :: + AlgRing.C a + => Int -- ^ index + -> Int -- ^ exponent + -> Spray a +lone' n p + | n < 0 = error "lone': invalid index." + | p < 0 = error "lone': invalid exponent" + | otherwise = HM.singleton (lonePower n p) AlgRing.one + +-- | The rational spray @x_n^p@ +qlone' :: + Int -- ^ index + -> Int -- ^ exponent + -> QSpray +qlone' = lone' + +loneTerm' :: AlgRing.C a => Int -> Int -> Term a +loneTerm' n p = (lonePower n p, AlgRing.one) + +-- | Monomial spray, e.g. @monomial [(1,4),(3,2)]@ is @x^4.z^2@; indices +-- and exponents must be positive but this is not checked +-- prop> monomial [(1, 4), (3, 2)] == (lone 1 ^**^ 4) ^*^ (lone 3 ^**^ 2) +monomial :: + AlgRing.C a + => [(Int, Int)] -- ^ list of (index, exponent); duplicates are deleted + -> Spray a +monomial nps = if null nps + then unitSpray + else HM.singleton (powerize expnts) AlgRing.one + where + nps' = nub nps + nv = maximum (map fst nps') + expnts = S.fromList $ [fromMaybe 0 (lookup i nps') | i <- [1 .. nv]] + +-- | Monomial rational spray, a specialization of 'monomial' +-- +-- prop> qmonomial [(1, 4), (3, 2)] == (qlone 1 ^**^ 4) ^*^ (qlone 3 ^**^ 2) +qmonomial :: + [(Int, Int)] + -> QSpray +qmonomial = monomial + +-- | The unit spray +-- +-- prop> spray ^*^ unitSpray == spray +unitSpray :: AlgRing.C a => Spray a +unitSpray = HM.singleton nullPowers AlgRing.one + +-- | The null spray +-- +-- prop> spray ^+^ zeroSpray == spray +zeroSpray :: (Eq a, AlgAdd.C a) => Spray a +zeroSpray = AlgAdd.zero + +-- | Whether a spray is the zero spray +isZeroSpray :: Spray a -> Bool +isZeroSpray = HM.null + +-- | Constant spray +-- +-- prop> constantSpray 3 == 3 *^ unitSpray +constantSpray :: (Eq a, AlgAdd.C a) => a -> Spray a +constantSpray c = if c == AlgAdd.zero + then HM.empty + else HM.singleton nullPowers c + +-- | Get coefficient of a term of a spray +-- +-- >>> x = lone 1 :: Spray Int +-- >>> y = lone 2 :: Spray Int +-- >>> z = lone 3 :: Spray Int +-- >>> p = 2 *^ (2 *^ (x^**^3 ^*^ y^**^2)) ^+^ 4*^z ^+^ 5*^unitSpray +-- >>> getCoefficient [3, 2] p -- coefficient of x^3.y^2 +-- 4 +-- >>> 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] -- ^ 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) + +getCoefficient' :: AlgAdd.C a => Powers -> Spray a -> a +getCoefficient' powers spray = fromMaybe AlgAdd.zero (HM.lookup powers spray) + +-- | Get the constant term of a spray +-- +-- prop> getConstantTerm p == getCoefficient [] p +getConstantTerm :: AlgAdd.C a => Spray a -> a +getConstantTerm = getCoefficient' nullPowers + +-- | remove the constant term of a spray +removeConstantTerm :: Spray a -> Spray a +removeConstantTerm = HM.delete nullPowers + +-- | Whether a spray is constant; same as `isConstant` +isConstantSpray :: (Eq a, AlgRing.C a) => Spray a -> Bool +isConstantSpray = isConstant + +-- | helper function to unify evalSpray and evalSpraySpray +evalSprayHelper :: forall a. AlgRing.C a => [a] -> Spray a -> a +evalSprayHelper xyz spray = + AlgAdd.sum $ map evalTerm (HM.toList spray) + where + evalTerm :: Term a -> a + evalTerm (powers, coeff) = + coeff AlgRing.* AlgRing.product (zipWith (AlgRing.^) xyz pows) + where + pows = DF.toList (fromIntegral <$> exponents powers) + +-- | Evaluates a spray; same as `evaluate` +-- +-- >>> x = lone 1 :: Spray Int +-- >>> y = lone 2 :: Spray Int +-- >>> spray = 2*^x^**^2 ^-^ 3*^y +-- >>> evalSpray spray [2, 1] +-- 5 +evalSpray :: (Eq a, AlgRing.C a) => Spray a -> [a] -> a +evalSpray = evaluate + +-- | Evaluates the coefficients of a spray with spray coefficients; +-- same as `substituteParameters` +evalSpraySpray :: (Eq a, AlgRing.C a) => Spray (Spray a) -> [a] -> Spray a +evalSpraySpray spray xyz = if length xyz >= n + then removeZeroTerms $ HM.map (evalSprayHelper xyz) spray + else error "evalSpraySpray: not enough values provided." + where + n = maximum (HM.elems $ HM.map numberOfVariables spray) + +-- | spray from term +fromTerm :: Term a -> Spray a +fromTerm (pows, coeff) = HM.singleton pows coeff + +-- | Substitutes some values to some variables of a spray; same as `substitute` +-- +-- >>> x1 = lone 1 :: Spray Int +-- >>> x2 = lone 2 :: Spray Int +-- >>> x3 = lone 3 :: Spray Int +-- >>> p = x1^**^2 ^-^ x2 ^+^ x3 ^-^ unitSpray +-- >>> p' = substituteSpray [Just 2, Nothing, Just 3] p +-- >>> putStrLn $ prettyNumSprayX1X2X3 "x" p' +-- -x2 + 6 +substituteSpray :: (Eq a, AlgRing.C a) => [Maybe a] -> Spray a -> Spray a +substituteSpray = substitute + +-- | Converts a spray with rational coefficients to a spray with double +-- coefficients (useful for evaluation) +fromRationalSpray :: Spray Rational -> Spray Double +fromRationalSpray = HM.map fromRational + +-- | Sustitutes the variables of a spray with some sprays; same as `changeVariables` +-- +-- >>> x = lone 1 :: Spray Int +-- >>> y = lone 2 :: Spray Int +-- >>> z = lone 3 :: Spray Int +-- >>> p = x ^+^ y +-- >>> q = composeSpray p [x ^+^ y ^+^ z, z] +-- >>> putStrLn $ prettyNumSpray' q +-- X + Y + 2*Z +composeSpray :: + forall a. (AlgRing.C a, Eq a) => Spray a -> [Spray a] -> Spray a +composeSpray p = removeZeroTerms . evalSpray (identify p) + where + identify :: Spray a -> Spray (Spray a) + identify = HM.map constantSpray + +-- | Creates a spray from a list of terms +fromList :: + (AlgRing.C a, Eq a) + => [([Int], a)] -- ^ list of (exponents, coefficient) + -> Spray a +fromList x = removeZeroTerms $ HM.fromListWith (AlgAdd.+) $ map + (first (makePowers . S.fromList)) x + + +-- pretty stuff --------------------------------------------------------------- + +-- | Prints a spray; this function is exported for possible usage in other packages +showSpray :: + (a -> String) -- ^ function mapping a coefficient to a string, typically 'show' + -> (String, String) -- ^ pair of braces to enclose the coefficients + -> ([Exponents] -> [String]) -- ^ function mapping a list of exponents to a list of strings representing the monomials corresponding to these exponents + -> Spray a -- ^ the spray to be printed + -> String +showSpray showCoef braces showMonomials spray = + if isZeroSpray spray + then "0" + else unpack $ intercalate (pack " + ") stringTerms + where + terms = orderedTerms spray + (powers, coeffs) = unzip terms + stringMonomials = showMonomials (map exponents powers) + stringTerms = zipWith f coeffs stringMonomials + f coeff smonomial + | smonomial == "" = pack scoeff' + | scoeff == "" = pack smonomial + | otherwise = pack $ scoeff' ++ "*" ++ smonomial + where + scoeff = showCoef coeff + scoeff' = bracify braces scoeff + +-- | Prints a spray, with monomials shown as "x.z^2", and with +-- a user-defined showing function for the coefficients +showSprayXYZ :: + (a -> String) -- ^ function mapping a coefficient to a string, typically 'show' + -> (String, String) -- ^ used to enclose the coefficients, usually a pair of braces + -> [String] -- ^ typically some letters, to print the variables + -> Spray a -- ^ the spray to be printed + -> String +showSprayXYZ showCoef braces letters spray = + if null letters + then error "showSprayXYZ: empty list of strings." + else showSpray showCoef braces (showMonomialsXYZ letters) spray + +-- | Prints a spray, with monomials shown as @\"x.z^2\"@, and with +-- a user-defined showing function for the coefficients; this is the same as +-- the function `showSprayXYZ` with the pair of braces @("(", ")")@ +showSprayXYZ' :: + (a -> String) -- ^ function mapping a coefficient to a string, typically 'show' + -> [String] -- ^ typically some letters, to print the variables + -> Spray a -- ^ the spray to be printed + -> String +showSprayXYZ' showCoef = showSprayXYZ showCoef ("(", ")") + +-- | Pretty form of a spray with monomials displayed in the style of @\"x.z^2\"@; +-- you should rather use `prettyNumSprayXYZ` or `prettyQSprayXYZ` if your +-- coefficients are numeric +-- +-- >>> x = lone 1 :: Spray Int +-- >>> y = lone 2 :: Spray Int +-- >>> z = lone 3 :: Spray Int +-- >>> p = 2*^x ^+^ 3*^y^**^2 ^-^ 4*^z^**^3 +-- >>> putStrLn $ prettySprayXYZ ["X", "Y", "Z"] p +-- (2)*X + (3)*Y^2 + (-4)*Z^3 +-- >>> putStrLn $ prettySprayXYZ ["X", "Y"] p +-- (2)*X1 + (3)*X2^2 + (-4)*X3^3 +prettySprayXYZ :: + (Show a) + => [String] -- ^ typically some letters, to print the variables + -> Spray a -- ^ the spray to be printed + -> String +prettySprayXYZ = showSprayXYZ' show + +-- | Pretty form of a spray, with monomials shown as "x1.x3^2", and with +-- a user-defined showing function for the coefficients +showSprayX1X2X3 :: + (a -> String) -- ^ function mapping a coefficient to a string, typically 'show' + -> (String, String) -- ^ used to enclose the coefficients + -> String -- ^ typically a letter, to print the non-indexed variables + -> Spray a -- ^ the spray to be printed + -> String +showSprayX1X2X3 showCoef braces letter = + showSpray showCoef braces (showMonomialsX1X2X3 letter) + +-- | Pretty form of a spray, with monomials shown as "x1.x3^2", and with +-- a user-defined showing function for the coefficients; this is the same as +-- the function `showSprayX1X2X3` with the pair of braces @("(", ")")@ used to +-- enclose the coefficients +showSprayX1X2X3' :: + (a -> String) -- ^ function mapping a coefficient to a string, e.g. 'show' + -> String -- ^ typically a letter, to print the non-indexed variables + -> Spray a -- ^ the spray to be printed + -> String +showSprayX1X2X3' showCoef = showSprayX1X2X3 showCoef ("(", ")") + +-- | Pretty form of a spray with monomials displayed in the style of @\"x1.x3^2\"@; +-- you should rather use `prettyNumSprayX1X2X3` or `prettyQSprayX1X2X3` if your +-- coefficients are numeric +-- +-- >>> x = lone 1 :: Spray Int +-- >>> y = lone 2 :: Spray Int +-- >>> z = lone 3 :: Spray Int +-- >>> spray = 2*^x ^+^ 3*^y^**^2 ^-^ 4*^z^**^3 +-- >>> putStrLn $ prettySprayX1X2X3 "X" spray +-- (2)*X1 + (3)*X2^2 + (-4)*X3^3 +prettySprayX1X2X3 :: + Show a + => String -- ^ typically a letter, to print the non-indexed variables + -> Spray a -- ^ the spray to be printed + -> String +prettySprayX1X2X3 = showSprayX1X2X3' show + +-- | Pretty form of a spray with monomials displayed in the style of @\"x.z^2\"@; +-- you should rather use `prettyNumSpray` or `prettyQSpray` if you deal with +-- sprays with numeric coefficients +-- +-- >>> x = lone 1 :: Spray Int +-- >>> y = lone 2 :: Spray Int +-- >>> z = lone 3 :: Spray Int +-- >>> p = 2*^x ^+^ 3*^y^**^2 ^-^ 4*^z^**^3 +-- >>> putStrLn $ prettySpray p +-- (2)*x + (3)*y^2 + (-4)*z^3 +-- >>> putStrLn $ prettySpray (p ^+^ lone 4) +-- (2)*x1 + (3)*x2^2 + (-4)*x3^3 + x4 +-- +-- prop> prettySpray spray == prettySprayXYZ ["x", "y", "z"] spray +prettySpray :: (Show a) => Spray a -> String +prettySpray = prettySprayXYZ ["x", "y", "z"] + +-- | Pretty form of a spray, with monomials shown as @\"x1.x3^2\"@; use +-- `prettySprayX1X2X3` to change the letter (or `prettyNumSprayX1X2X3` +-- or `prettyQSprayX1X2X3` if the coefficients are numeric) +-- +-- >>> x = lone 1 :: Spray Int +-- >>> y = lone 2 :: Spray Int +-- >>> z = lone 3 :: Spray Int +-- >>> p = 2*^x ^+^ 3*^y^**^2 ^-^ 4*^z^**^3 +-- >>> putStrLn $ prettySpray' p +-- (2)*x1 + (3)*x2^2 + (-4)*x3^3 +prettySpray' :: Show a => Spray a -> String +prettySpray' = prettySprayX1X2X3 "x" + +-- | showMonomialOld "x" [0, 2, 1] = x^(0, 2, 1) +showMonomialsOld :: String -> [Exponents] -> [String] +showMonomialsOld var = map (showMonomialOld var) + where + showMonomialOld :: String -> Exponents -> String + showMonomialOld a pows = + unpack $ append (pack x) (cons '(' $ snoc string ')') + where + x = a ++ "^" + string = intercalate (pack ", ") (map (pack . show) (DF.toList pows)) + +-- | Pretty form of a spray; you will probably prefer `prettySpray` or `prettySpray'` +-- +-- >>> x = lone 1 :: Spray Int +-- >>> y = lone 2 :: Spray Int +-- >>> z = lone 3 :: Spray Int +-- >>> p = 2*^x ^+^ 3*^y^**^2 ^-^ 4*^z^**^3 +-- >>> putStrLn $ prettySpray'' "x" p +-- (2)*x^(1) + (3)*x^(0, 2) + (-4)*x^(0, 0, 3) +prettySpray'' :: + Show a + => String -- ^ a string denoting the variables, e.g. \"x\" + -> Spray a -- ^ the spray + -> String +prettySpray'' var = showSpray show ("(", ")") (showMonomialsOld var) + +-- | Show a spray with numeric coefficients; this function is exported for +-- possible usage in other packages +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 + -> Spray a + -> String +showNumSpray showMonomials showCoeff spray = + if isZeroSpray spray + then "0" + else concat $ zipWith (++) stringSigns stringTerms + where + terms = orderedTerms spray + coeffs = map snd terms + (firstCoeff, otherCoeffs) = fromJust (uncons coeffs) + firstSign = if firstCoeff > 0 then "" else "-" + otherSigns = map (\x -> if x > 0 then " + " else " - ") otherCoeffs + stringSigns = firstSign : otherSigns + absCoeffs = map abs coeffs + powers = map (exponents . fst) terms + stringMonomials = showMonomials powers + stringTerms = zipWith f absCoeffs stringMonomials + f acoeff smonomial + | smonomial == "" = showCoeff acoeff + | scoeff == "" = smonomial + | otherwise = scoeff ++ "*" ++ smonomial + where + scoeff = if acoeff == 1 then "" else showCoeff acoeff + +-- | showMonomialX1X2X3 "X" [0, 2, 1] = "X2^2.X3" +showMonomialX1X2X3 :: String -> Exponents -> Text +showMonomialX1X2X3 x pows = x1x2x3 + where + f i p + | p == 0 = pack "" + | p == 1 = pack $ x ++ show i + | otherwise = pack $ x ++ show i ++ "^" ++ show p + indices = S.findIndicesL (/= 0) pows + x1x2x3 = + intercalate (pack ".") [f (i+1) (pows `index` i) | i <- indices] + +-- | showMonomialsX1X2X3 "X" [[0, 2, 1], [1, 2]] = ["X2^2.X3", "X1.X2"] +showMonomialsX1X2X3 :: String -> [Exponents] -> [String] +showMonomialsX1X2X3 x = map (unpack . showMonomialX1X2X3 x) + +-- | showMonomialXYZ ["X", "Y", "Z"] 3 [1, 2, 1] = X.Y^2.Z +-- showMonomialXYZ ["X", "Y", "Z"] 3 [1, 2, 1, 2] = X1.X2^2.X3.X4^2 +showMonomialXYZ :: [String] -> Int -> Exponents -> Text +showMonomialXYZ letters n pows = if n <= length letters + then xyz + else showMonomialX1X2X3 (letters !! 0) pows + where + f letter p + | p == 0 = pack "" + | p == 1 = pack letter + | otherwise = pack $ letter ++ "^" ++ show p + indices = S.findIndicesL (/= 0) pows + xyz = intercalate (pack ".") + [f (letters !! i) (pows `index` i) | i <- indices] + +-- | showMonomialsXYZ ["X", "Y", "Z"] [[0, 2, 1], [1, 2]] = ["Y^2.Z", "X.Y^2"] +showMonomialsXYZ :: [String] -> [Exponents] -> [String] +showMonomialsXYZ letters powers = + map (unpack . showMonomialXYZ letters n) powers + where + n = maximum (map S.length powers) + +-- | Pretty form of a spray with numeric coefficients, printing monomials as @\"x1.x3^2\"@ +-- +-- >>> x = lone 1 :: Spray Int +-- >>> y = lone 2 :: Spray Int +-- >>> z = lone 3 :: Spray Int +-- >>> p = 2*^x ^+^ 3*^y^**^2 ^-^ 4*^z^**^3 +-- >>> putStrLn $ prettyNumSprayX1X2X3 "x" p +-- 2*x1 + 3*x2^2 - 4*x3^3 +prettyNumSprayX1X2X3 :: (Num a, Ord a, Show a) + => String -- ^ usually a letter such as @\"x\"@ to denote the non-indexed variables + -> Spray a + -> String +prettyNumSprayX1X2X3 x = showNumSpray (showMonomialsX1X2X3 x) show + +-- | Pretty form of a spray with numeric coefficients, printing monomials as @\"x.z^2\"@ +-- if possible, i.e. if enough letters are provided, otherwise as @\"x1.x3^2\"@ +-- +-- >>> x = lone 1 :: Spray Int +-- >>> y = lone 2 :: Spray Int +-- >>> z = lone 3 :: Spray Int +-- >>> w = lone 4 :: Spray Int +-- >>> p = 2*^x ^+^ 3*^y^**^2 ^-^ 4*^z^**^3 +-- >>> putStrLn $ prettyNumSprayXYZ ["x","y","z"] p +-- 2*x + 3*y^2 - 4*z^3 +-- >>> putStrLn $ prettyNumSprayXYZ ["x","y","z"] (p ^+^ w) +-- 2*x1 + 3*x2^2 - 4*x3^3 + x4 +-- >>> putStrLn $ prettyNumSprayXYZ ["a","b","c"] (p ^+^ w) +-- 2*a1 + 3*a2^2 - 4*a3^3 + a4 +prettyNumSprayXYZ :: (Num a, Ord a, Show a) + => [String] -- ^ usually some letters, denoting the variables + -> Spray a + -> String +prettyNumSprayXYZ letters = showNumSpray (showMonomialsXYZ letters) show + +-- | helper function for showQSpray +showRatio :: Rational -> String +showRatio q = if d == 1 + then show n + else "(" ++ show n ++ "/" ++ show d ++ ")" + where + n = DR.numerator q + d = DR.denominator q + +-- | helper function for showQSpray' +showRatio' :: (Eq a, Num a, Show a) => NumberRatio.T a -> String +showRatio' q = if d == 1 + then show n + else "(" ++ show n ++ "/" ++ show d ++ ")" + where + n = NumberRatio.numerator q + d = NumberRatio.denominator q + +-- | Prints a `QSpray`; for internal usage but exported for usage in other packages +showQSpray :: + ([Exponents] -> [String]) -- ^ function printing monomials + -> QSpray + -> String +showQSpray showMonomials = showNumSpray showMonomials showRatio + +-- | Prints a `QSpray'`; for internal usage but exported for usage in other packages +showQSpray' :: + ([Exponents] -> [String]) -- ^ function mapping a list of monomials exponents to a list of strings + -> QSpray' + -> String +showQSpray' showMonomials = showNumSpray showMonomials showRatio' + +-- | Pretty form of a spray with rational coefficients, printing monomials in +-- the style of @\"x1.x3^2\"@ +-- +-- >>> x = lone 1 :: QSpray +-- >>> y = lone 2 :: QSpray +-- >>> z = lone 3 :: QSpray +-- >>> p = 2*^x ^+^ 3*^y^**^2 ^-^ (4%3)*^z^**^3 +-- >>> putStrLn $ prettyQSprayX1X2X3 "x" p +-- 2*x1 + 3*x2^2 - (4/3)*x3^3 +prettyQSprayX1X2X3 :: + String -- ^ usually a letter such as @\"x\"@, to denote the non-indexed variables + -> QSpray + -> String +prettyQSprayX1X2X3 x = showQSpray (showMonomialsX1X2X3 x) + +-- | Same as `prettyQSprayX1X2X3` but for a `QSpray'` spray +prettyQSprayX1X2X3' :: + String -- ^ usually a letter such as @\"x\"@, to denote the non-indexed variables + -> QSpray' + -> String +prettyQSprayX1X2X3' x = showQSpray' (showMonomialsX1X2X3 x) + +-- | Pretty form of a spray with rational coefficients, printing monomials in +-- the style of @\"x.z^2\"@ with the provided letters if possible, i.e. if enough +-- letters are provided, otherwise in the style @\"x1.x3^2\"@, taking the first +-- provided letter to denote the non-indexed variables +-- +-- >>> x = lone 1 :: QSpray +-- >>> y = lone 2 :: QSpray +-- >>> z = lone 3 :: QSpray +-- >>> p = 2*^x ^+^ 3*^y^**^2 ^-^ (4%3)*^z^**^3 +-- >>> putStrLn $ prettyQSprayXYZ ["x","y","z"] p +-- 2*x + 3*y^2 - (4/3)*z^3 +-- >>> putStrLn $ prettyQSprayXYZ ["x","y"] p +-- 2*x1 + 3*x2^2 - (4%3)*x3^3 +-- >>> putStrLn $ prettyQSprayXYZ ["a","b"] p +-- 2*a1 + 3*a2^2 - (4/3)*a3^3 +prettyQSprayXYZ :: + [String] -- ^ usually some letters, to denote the variables + -> QSpray + -> String +prettyQSprayXYZ letters = showQSpray (showMonomialsXYZ letters) + +-- | Same as `prettyQSprayXYZ` but for a `QSpray'` spray +prettyQSprayXYZ' :: + [String] -- ^ usually some letters, to denote the variables + -> QSpray' + -> String +prettyQSprayXYZ' letters = showQSpray' (showMonomialsXYZ letters) + +-- | Pretty printing of a spray with rational coefficients +-- prop> prettyQSpray == prettyQSprayXYZ ["x", "y", "z"] +prettyQSpray :: QSpray -> String +prettyQSpray = prettyQSprayXYZ ["x", "y", "z"] + +-- | Pretty printing of a spray with rational coefficients +-- prop> prettyQSpray'' == prettyQSprayXYZ ["X", "Y", "Z"] +prettyQSpray'' :: QSpray -> String +prettyQSpray'' = prettyQSprayXYZ ["X", "Y", "Z"] + +-- | Pretty printing of a spray with rational coefficients +-- prop> prettyQSpray' == prettyQSprayXYZ' ["x", "y", "z"] +prettyQSpray' :: QSpray' -> String +prettyQSpray' = prettyQSprayXYZ' ["x", "y", "z"] + +-- | Pretty printing of a spray with rational coefficients +-- prop> prettyQSpray''' == prettyQSprayXYZ' ["X", "Y", "Z"] +prettyQSpray''' :: QSpray' -> String +prettyQSpray''' = prettyQSprayXYZ' ["X", "Y", "Z"] + +-- | Pretty printing of a spray with numeric coefficients +-- prop> prettyNumSpray == prettyNumSprayXYZ ["x", "y", "z"] +prettyNumSpray :: (Num a, Ord a, Show a) => Spray a -> String +prettyNumSpray = prettyNumSprayXYZ ["x", "y", "z"] + +-- | Pretty printing of a spray with numeric coefficients +-- prop> prettyNumSpray' == prettyNumSprayXYZ ["X", "Y", "Z"] +prettyNumSpray' :: (Num a, Ord a, Show a) => Spray a -> String +prettyNumSpray' = prettyNumSprayXYZ ["X", "Y", "Z"] + + +-- misc ----------------------------------------------------------------------- + +-- | spray as safe spray +safeSpray :: Spray a -> SafeSpray a +safeSpray = HM.mapKeys exponents + +getCoefficient'' :: AlgAdd.C a => Exponents -> SafeSpray a -> a +getCoefficient'' powers spray' = + fromMaybe AlgAdd.zero (HM.lookup powers spray') + +getConstantTerm'' :: AlgAdd.C a => SafeSpray a -> a +getConstantTerm'' = getCoefficient'' S.empty + +removeConstantTerm'' :: SafeSpray a -> SafeSpray a +removeConstantTerm'' = HM.delete S.empty + +-- | Sum of sprays +sumOfSprays :: (Eq a, AlgAdd.C a) => [Spray a] -> Spray a +sumOfSprays = AlgAdd.sum + +-- | Product of sprays +productOfSprays :: (Eq a, AlgRing.C a) => [Spray a] -> Spray a +productOfSprays = AlgRing.product + +-- | ordered terms of a spray +orderedTerms :: Spray a -> [Term a] +orderedTerms spray = + sortBy (flip compare `on` (exponents . fst)) (HM.toList spray) + +-- | Spray as a list +toList :: Spray a -> [([Int], a)] +toList p = HM.toList $ HM.mapKeys (DF.toList . exponents) p + +-- | 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 [2 .. n] + times k x = k .^ x + +-- | Whether two sprays are equal up to a scalar factor +collinearSprays :: (Eq a, AlgRing.C a) => Spray a -> Spray a -> Bool +collinearSprays spray1 spray2 = + isZeroSpray spray1 && isZeroSpray spray2 || + (not . isZeroSpray) spray1 && (not . isZeroSpray) spray2 && + snd (leadingTerm spray1) *^ spray2 == snd (leadingTerm spray2) *^ spray1 + +-- | Checks whether the multivariate polynomial defined by a spray is homogeneous +-- and also returns the degree in case this is true +isHomogeneousSpray :: (Eq a, AlgRing.C a) => Spray a -> (Bool, Maybe Int) +isHomogeneousSpray spray + | isConstant spray = (True, Just 0) + | getConstantTerm spray /= AlgAdd.zero = (False, Nothing) + | otherwise = (check, deg) + where + degrees = map (sum . exponents) (HM.keys spray) + 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 +-- maxWithIndex :: Ord a => [a] -> (Int, a) +-- maxWithIndex = maximumBy (comparing snd) . zip [0 .. ] + +-- | leading term of a spray +leadingTerm :: Spray a -> Term a +leadingTerm p = (biggest, p HM.! biggest) + where + powers = HM.keys p + biggest = maximumBy (comparing exponents) powers +-- (i, biggest) = maxWithIndex powers + -- biggest = powers !! i + +-- | whether a term divides another term +divides :: Term a -> Term a -> Bool +divides (powsP, _) (powsQ, _) = nvP <= nvQ && lower + where + nvP = nvariables powsP + nvQ = nvariables powsQ + expntsP = exponents powsP + expntsQ = exponents powsQ + lower = DF.all (uncurry (<=)) (S.zip expntsP expntsQ) + +-- | quotient of term Q by term P, assuming P divides Q +quotient :: AlgField.C a => Term a -> Term a -> Term a +quotient (powsQ, coeffQ) (powsP, coeffP) = (pows, coeff) + where + (expntsP, expntsQ) = harmonize (powsP, powsQ) + expnts = S.zipWith (-) expntsQ expntsP + pows = makePowers expnts + coeff = coeffQ AlgField./ coeffP + +-- | Remainder of the division of a spray by a list of divisors, +-- using the lexicographic ordering of the monomials +sprayDivisionRemainder :: forall a. (Eq a, AlgField.C a) + => Spray a -> [Spray a] -> Spray a +sprayDivisionRemainder p qs = + if n == 0 + then error "sprayDivisionRemainder: the list of divisors is empty." + else ogo p zeroSpray + where + n = length qs + qsltqs = zip qs (map leadingTerm qs) + g :: Term a -> Spray a -> Spray a -> (Spray a, Spray a) + g lts s r = (HM.delete (fst lts) s, addTerm r lts) + go :: Term 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 ^-^ multSprayByTerm q (quotient lts ltq) + else s + ogo :: Spray a -> Spray a -> Spray a + ogo !s !r + | isZeroSpray s = r + | otherwise = ogo s' r' + where + (s', r') = go (leadingTerm s) s r 0 False + +-- | Division of a spray by a spray +sprayDivision :: + (Eq a, AlgField.C a) + => Spray a -- ^ dividand + -> Spray a -- ^ divisor + -> (Spray a, Spray a) -- ^ (quotient, remainder) +sprayDivision sprayA sprayB = + if isConstant sprayB + then if isZeroSpray sprayB + then + error "sprayDivision: division by zero." + else + let c = getConstantTerm sprayB in (sprayA /> c, zeroSpray) + else sprayDivision0 sprayA sprayB + +sprayDivision0 :: forall a. (Eq a, AlgField.C a) + => Spray a -- ^ dividand + -> Spray a -- ^ divisor + -> (Spray a, Spray a) -- ^ (quotient, remainder) +sprayDivision0 sprayA sprayB = +-- both ((*^) coeff) (ogo sprayA' zeroSpray zeroSpray) + ogo sprayA zeroSpray zeroSpray + where + -- q'*(sprayB) + r' = sprayA * lc sprayB / lc sprayA + -- (q' * lc sprayA / lc sprayB) * sprayB + r' * lc sprayA / lc sprayB = sprayA +-- lcA = snd (leadingTerm sprayA) +-- coeff = lcA AlgField./ snd ltB +-- sprayA' = sprayA /> coeff + + ltB = leadingTerm sprayB + ogo :: Spray a -> Spray a -> Spray a -> (Spray a, Spray a) + ogo !p !q !r + | isZeroSpray p = (q, r) + | otherwise = ogo p' q' r' + where + ltp = leadingTerm p + (p', q', r') = if divides ltB ltp + then (newp, newq, r) + else (HM.delete (fst ltp) p, q, addTerm r ltp) + qtnt = quotient ltp ltB + newp = p ^-^ multSprayByTerm sprayB qtnt + newq = addTerm q qtnt + +-- | division of univariate sprays with degree(dividend) >= degree(divisor) +univariateSprayDivision :: forall a. (Eq a, AlgField.C a) + => Spray a -- ^ dividand + -> Spray a -- ^ divisor + -> (Spray a, Spray a) -- ^ (quotient, remainder) +univariateSprayDivision sprayA sprayB = + if isConstant sprayB + then + let c = getConstantTerm sprayB in (sprayA /> c, zeroSpray) + else + ogo sprayA zeroSpray zeroSpray + where + (powsLTB, coeffLTB) = leadingTerm sprayB + degB = exponents powsLTB + expntLTB = degB `index` 0 + ogo :: Spray a -> Spray a -> Spray a -> (Spray a, Spray a) + ogo !p !q !r + | isZeroSpray p = (q, r) + | otherwise = ogo p' q' r' + where + (powsLTP, coeffLTP) = leadingTerm p + degP = exponents powsLTP + (p', q', r') = if degB <= degP + then (newp, newq, r) + else (zeroSpray, q, r ^+^ p) + newp = p ^-^ multSprayByTerm sprayB qtnt + newq = addTerm q qtnt + qtnt = (pows, coeff) + expntLTP = degP `index` 0 + pows = if expntLTP == expntLTB + then nullPowers + else Powers (S.singleton (expntLTP - expntLTB)) 1 + coeff = coeffLTP AlgField./ coeffLTB + + +-- Groebner stuff ------------------------------------------------------------- + +-- | slight modification of `sprayDivisionRemainder` to speed up groebner00 +sprayDivisionRemainder' :: + forall a. (Eq a, AlgField.C a) + => Spray a -> HashMap Int (Spray a, Term a) -> Spray a +sprayDivisionRemainder' p qsltqs = ogo p zeroSpray + where + n = HM.size qsltqs + g :: Term a -> Spray a -> Spray a -> (Spray a, Spray a) + g lts s r = (HM.delete (fst lts) s, addTerm r lts) + go :: Term 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 ^-^ multSprayByTerm q (quotient lts ltq) + else s + ogo :: Spray a -> Spray a -> Spray a + ogo !s !r + | isZeroSpray 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 .. ] (zip row1 row2)) + where + 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, Term a) -> (Spray a, Term a) -> Spray a +sPolynomial (p, (lpowsP, lcoefP)) (q, (lpowsQ, lcoefQ)) = + multSprayByTerm p wp ^-^ multSprayByTerm q wq + where + (lexpntsP, lexpntsQ) = harmonize (lpowsP, lpowsQ) + gamma = S.zipWith max lexpntsP lexpntsQ + betaP = S.zipWith (-) gamma lexpntsP + betaQ = S.zipWith (-) gamma lexpntsQ + wp = (makePowers betaP, AlgField.recip lcoefP) + wq = (makePowers betaQ, 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, 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 + (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 isZeroSpray sbarfg + 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 [] + +-- | Reduces a Gröbner 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 = spray /> coef + where + (_, coef) = leadingTerm spray + ngbasis = map normalize gbasis + n = length ngbasis + reduction :: Int -> Spray a + reduction i = sprayDivisionRemainder (ngbasis !! i) rest + where + rest = [ngbasis !! k | k <- [0 .. i-1] ++ [i+1 .. n-1]] + +-- | Gröbner basis, always minimal and possibly reduced +-- +-- prop> groebnerBasis sprays True == reduceGroebnerBasis (groebnerBasis sprays False) +groebnerBasis :: + forall a. (Eq a, AlgField.C a) + => [Spray a] -- ^ list of sprays + -> Bool -- ^ whether to return the reduced basis + -> [Spray a] +groebnerBasis sprays reduced = + if reduced then reduceGroebnerBasis gbasis0 else map normalize gbasis0 + where + gbasis0 = groebner0 sprays + normalize :: Spray a -> Spray a + normalize spray = spray /> coef + where + (_, coef) = leadingTerm spray + + +-- elementary symmetric polynomials ------------------------------------------- + +-- | combinations of k elements among a list +combinationsOf :: Int -> [a] -> [[a]] +combinationsOf _ [] = error "combinationsOf: should not happen." +combinationsOf 1 as = map pure as +combinationsOf k as@(_:xs) = + run (l-1) (k-1) as $ combinationsOf (k-1) xs + where + l = length as + run :: Int -> Int -> [a] -> [[a]] -> [[a]] + run n i ys cs + | n == i = map (ys ++) cs + | otherwise = map (q:) cs ++ run (n-1) i qs (drop dc cs) + where + f :: [a] -> (a, [a]) + f [] = error "combinationsOf: should not happen." + f (b:bs) = (b, bs) + (q, qs) = f (take (n-i+1) ys) + dc = product [n-i+1 .. n-1] `div` product [1 .. i-1] + +-- | generates all permutations of a binary sequence +permutationsBinarySequence :: Int -> Int -> [Seq Int] +permutationsBinarySequence nzeros nones = + let n = nzeros + nones in + map (binarySequence n) (combinationsOf nones [0 .. n-1]) + where + binarySequence :: Int -> [Int] -> Seq Int + binarySequence n combo = fromFunction n f + where + f :: Int -> Int + f i = fromEnum (i `elem` combo) + +-- | Elementary symmetric polynomial +-- +-- >>> putStrLn $ prettySpray' (esPolynomial 3 2) +-- (1)*x1.x2 + (1)*x1.x3 + (1)*x2.x3 +esPolynomial :: + (AlgRing.C a, Eq a) + => Int -- ^ number of variables + -> Int -- ^ index + -> Spray a +esPolynomial n k + | k < 0 || n < 0 + = error "esPolynomial: both arguments must be positive integers." + | k > n = zeroSpray + | k == 0 = unitSpray + | otherwise = spray + where + perms = permutationsBinarySequence (n-k) k + spray = HM.fromList $ map (\expts -> (makePowers expts, AlgRing.one)) perms + +-- | Power sum polynomial +psPolynomial :: + forall a. (AlgRing.C a, Eq a) + => Int -- ^ number of variables + -> Int -- ^ power + -> Spray a +psPolynomial n k + | k < 0 || n < 0 + = error "psPolynomial: both arguments must be positive integers." + | k > n = AlgAdd.zero + | k == 0 = unitSpray + | otherwise = spray + where + spray = HM.fromList $ map f [1 .. n] + f :: Int -> (Powers, a) + f j = (Powers expts j, AlgRing.one) + where + expts = S.replicate (j-1) 0 |> k + +-- | Whether a spray is a symmetric polynomial, an inefficient algorithm +-- (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 = check + where + n = numberOfVariables spray + indices = [1 .. n] + gPolys = [esPolynomial n i ^-^ lone (n + i) | i <- indices] + gbasis = groebner0 gPolys + spray' = removeConstantTerm spray + g = sprayDivisionRemainder spray' gbasis + 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 +-- +-- >>> x = lone 1 :: Spray Rational +-- >>> y = lone 2 :: Spray Rational +-- >>> p1 = x ^+^ y +-- >>> p2 = x ^-^ y +-- >>> p = p1 ^*^ p2 +-- +-- prop> isPolynomialOf p [p1, p2] == (True, Just $ x ^*^ y) +isPolynomialOf :: forall a. (AlgField.C a, Eq a) + => Spray a -> [Spray a] -> (Bool, Maybe (Spray a)) +isPolynomialOf spray sprays + | isConstant spray = (True, Just spray) + | null sprays = error "isPolynomialOf: the list of sprays is empty." + | otherwise = result + where + nov = numberOfVariables spray + n = maximum $ map numberOfVariables sprays + result + | nov > n = (False, Nothing) + | otherwise = (check, poly) + where + m = length sprays + yPolys = [lone (n + i) | i <- [1 .. m]] + gPolys = zipWith (^-^) sprays yPolys + gbasis0 = groebner0 gPolys + constantTerm = getConstantTerm spray + spray' = removeConstantTerm spray + g = sprayDivisionRemainder spray' gbasis0 + check = not $ any (involvesVariable g) [1 .. n] + poly = if check + then Just (constantTerm +> dropVariables n g) + else Nothing + + +-- resultant ------------------------------------------------------------------ + +-- | sylvester matrix +sylvesterMatrix :: AlgAdd.C a => [a] -> [a] -> Matrix a +sylvesterMatrix x y = fromLists (xrows ++ yrows) + 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]] + +-- | "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: + -- the point to get is 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]] + +-- | the coefficients of a spray as a univariate spray in x_1 with +-- spray coefficients +sprayCoefficients :: (Eq a, AlgRing.C a) => Spray a -> [Spray a] +sprayCoefficients spray = + if n == 0 + then [constantTerm] + else sprays + where + n = numberOfVariables spray + spray'' = safeSpray spray + spray' = removeConstantTerm'' spray'' + (expnts', coeffs') = unzip (HM.toList spray') + constantTerm = (constantSpray . getConstantTerm'') spray'' + xpows = map (`index` 0) expnts' + powers'' = + map (powerize . S.deleteAt 0) expnts' + sprays'' = zipWith (curry fromTerm) powers'' coeffs' + imap = IM.fromListWith (^+^) (zip xpows sprays'') + imap' = IM.insertWith (^+^) 0 constantTerm imap + permutation = [2 .. n] ++ [1] + deg = maximum xpows + sprays = [ + permuteVariables permutation (fromMaybe zeroSpray (IM.lookup i imap')) + | i <- [deg, deg-1 .. 0] + ] + +-- | Resultant of two /univariate/ sprays +resultant1 :: (Eq a, AlgRing.C a) => Spray a -> Spray a -> a +resultant1 p q = + if n <= 1 + then detLaplace $ sylvesterMatrix pcoeffs qcoeffs + else error "resultant1: the two sprays must be univariate." + where + n = max (numberOfVariables p) (numberOfVariables q) + p'' = safeSpray p + q'' = safeSpray q + pexpnts = + map (`index` 0) $ HM.keys $ removeConstantTerm'' p'' + qexpnts = + map (`index` 0) $ HM.keys $ removeConstantTerm'' q'' + p0 = getConstantTerm'' p'' + q0 = getConstantTerm'' q'' + pcoeffs = if null pexpnts + then [p0] + else [getCoefficient'' (S.singleton i) p'' | i <- [maxp, maxp-1 .. 1]] + ++ [p0] + where + maxp = maximum pexpnts + qcoeffs = if null qexpnts + then [q0] + else [getCoefficient'' (S.singleton i) q'' | i <- [maxq, maxq-1 .. 1]] + ++ [q0] + where + maxq = maximum qexpnts + +-- | Subresultants of two /univariate/ sprays. This function makes several calls +-- to `detLaplace` and then it can be slow. +subresultants1 :: (Eq a, AlgRing.C a) => Spray a -> Spray a -> [a] +subresultants1 p q = if n <= 1 + then map (detLaplace . sylvesterMatrix' pcoeffs qcoeffs) [0 .. min d e - 1] + else error "subresultants1: the two sprays must be univariate." + where + n = max (numberOfVariables p) (numberOfVariables q) + p'' = safeSpray p + q'' = safeSpray q + pexpnts = + map (`index` 0) $ HM.keys $ removeConstantTerm'' p'' + qexpnts = + map (`index` 0) $ HM.keys $ removeConstantTerm'' q'' + p0 = getConstantTerm'' p'' + q0 = getConstantTerm'' q'' + pcoeffs = if null pexpnts + then [p0] + else [getCoefficient'' (S.singleton i) p'' | i <- [maxp, maxp-1 .. 1]] + ++ [p0] + where + maxp = maximum pexpnts + qcoeffs = if null qexpnts + then [q0] + else [getCoefficient'' (S.singleton i) q''| i <- [maxq, maxq-1 .. 1]] + ++ [q0] + where + maxq = maximum qexpnts + d = length pcoeffs + e = length qcoeffs + +-- | Resultant of two sprays +resultant :: (Eq a, AlgRing.C a) + => Int -- ^ indicator of the variable with respect to which the resultant is desired (e.g. 1 for x) + -> Spray a + -> Spray a + -> Spray a +resultant var p q = + if var >= 1 && var <= n + then permuteVariables permutation' det + else error "resultant: invalid variable index." + where + n = max (numberOfVariables p) (numberOfVariables q) + permutation = [n-var+2 .. n] ++ [1 .. n-var+1] + permutation' = [var .. n] ++ [1 .. var-1] + p' = permuteVariables permutation p + q' = permuteVariables permutation q + det = detLaplace $ + sylvesterMatrix (sprayCoefficients p') (sprayCoefficients q') + +-- | Subresultants of two sprays (the /principal/ subresultants, while the +-- `polynomialSubresultants` function returns the polynomial subresultants). +-- This function makes several calls to `detLaplace` and then it can be slow. +subresultants :: (Eq a, AlgRing.C a) + => Int -- ^ indicator of the variable with respect to which the subresultants are desired (e.g. 1 for x) + -> Spray a + -> Spray a + -> [Spray a] +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] + where + pcoeffs = sprayCoefficients p' + qcoeffs = sprayCoefficients q' + d = length pcoeffs + e = length qcoeffs + n = max (numberOfVariables p) (numberOfVariables q) + permutation = var : [1 .. var-1] ++ [var+1 .. n] + permute = permuteVariables permutation + p' = permute p + q' = permute q + permutation' = [2 .. var] ++ (1 : [var+1 .. n]) + permute' = permuteVariables permutation' + +-- | Resultant of two sprays with coefficients in a field; this function is more +-- efficient than the function `resultant` +resultant' :: forall a. (Eq a, AlgField.C a) + => Int -- ^ indicator of the variable with respect to which the resultant is desired (e.g. 1 for x) + -> Spray a + -> Spray a + -> Spray a +resultant' var sprayA sprayB + | var < 1 || var > n + = error "resultant': invalid variable index." + | isZeroSpray sprayA || isZeroSpray sprayB + = zeroSpray + | otherwise + = permuteVariables permutation' $ go unitSpray unitSpray s0 p0 q0 + where + n = max (numberOfVariables sprayA) (numberOfVariables sprayB) + permutation = [n-var+1 .. n] ++ [1 .. n-var] + permutation' = [var+1 .. n] ++ [1 .. var] + sprayA' = permuteVariables permutation sprayA + sprayB' = permuteVariables permutation sprayB + degA = degree n sprayA' + degB = degree n sprayB' + content :: Spray a -> Spray a + content spray = foldl1' gcdSpray (sprayCoefficients' n spray) + exactDivisionBy :: Spray a -> Spray a -> Spray a + exactDivisionBy b a = + if isZeroSpray remainder + then quo + else error "exactDivisionBy: should not happen." + where + (quo, remainder) = sprayDivision a b + contA = content sprayA' + contB = content sprayB' + sprayA'' = exactDivisionBy contA sprayA' + sprayB'' = exactDivisionBy contB sprayB' + t = contA^**^degB ^*^ contB^**^degA + s0 = if degA < degB && odd degA && odd degB + then AlgAdd.negate unitSpray :: Spray a + else unitSpray + (p0, q0) = if degA >= degB + then (sprayA'', sprayB'') + else (sprayB'', sprayA'') + go :: Spray a -> Spray a -> Spray a -> Spray a -> Spray a -> Spray a + go g h s p q = + if degq' == 0 + then s' ^*^ t ^*^ h'' + else go g' h' s' p' q' + where + degp = degree n p + degq = degree n q + delta = degp - degq + s' = if odd degp && odd degq + then AlgAdd.negate s + else s + (_, (_, r)) = pseudoDivision' n p q + p' = q + q' = exactDivisionBy (g ^*^ h^**^delta) r + (degp', ellp') = degreeAndLeadingCoefficient n p' + (degq', ellq') = degreeAndLeadingCoefficient n q' + g' = ellp' + h' = exactDivisionBy (h^**^delta) (h ^*^ g'^**^delta) + h'' = exactDivisionBy (h'^**^degp') (h' ^*^ ellq'^**^degp') + +-- | Polynomial subresultants of two sprays (the `subresultants` function +-- computes the /principal/ subresultants). +-- This function makes several calls to `detLaplace` and then it can be slow. +polynomialSubresultants :: + (Eq a, AlgRing.C a) + => Int -- ^ index of the variable with respect to which the subresultants will be computed (e.g. 2 for @y@) + -> Spray a + -> Spray a + -> [Spray a] +polynomialSubresultants var p q + | var < 1 || var > d + = error "polynomialSubresultants: invalid variable index." + | otherwise = [sres i | i <- [0 .. min n m]] + where + d = max (numberOfVariables p) (numberOfVariables q) + permutation = [d-var+1 .. d] ++ [1 .. d-var] + permutation' = [var+1 .. d] ++ [1 .. var] + p' = permuteVariables permutation p + q' = permuteVariables permutation q + pcoeffs = reverse $ sprayCoefficients' d p' + qcoeffs = reverse $ sprayCoefficients' d q' + pcoeff k = if k < 0 then zeroSpray else pcoeffs !! k + qcoeff k = if k < 0 then zeroSpray else qcoeffs !! k + n = length pcoeffs - 1 + m = length qcoeffs - 1 + x = loneTerm' d + prow i k = replicate k zeroSpray ++ + [pcoeff j | j <- [n, n-1 .. 2*i - m + k + 2]] ++ + [multSprayByTerm p' (x (m-i-1-k))] + prows i = [prow i k | k <- [0 .. m - i - 1]] + qrow i k = replicate k zeroSpray ++ + [qcoeff j | j <- [m, m-1 .. 2*i - n + k + 2]] ++ + [multSprayByTerm q' (x (n-i-1-k))] + qrows i = [qrow i k | k <- [0 .. n - i - 1]] + matrix i = fromLists (prows i ++ qrows i) + sres i = if n == m && i == n + then q + else permuteVariables permutation' $ detLaplace (matrix i) + +-- | Sturm-Habicht sequence of a spray. This function calls `polynomialSubresultants` +-- and then it can be slow. +sturmHabichtSequence :: + (Eq a, AlgRing.C a) + => Int -- ^ index of the variable with respect to which the Sturm-Habicht sequence will be computed (e.g. 2 for @y@) + -> Spray a + -> [Spray a] +sturmHabichtSequence var spray + | var < 1 || var > d + = error "sturmHabichtSequence: invalid variable index." + | otherwise = signedSubresultants ++ [spray', spray] + where + d = numberOfVariables spray + spray' = derivative var spray + spraySubresultants = polynomialSubresultants var spray spray' + n = length spraySubresultants + delta = [odd $ ((n - k - 1)*(n - k)) `div` 2 | k <- [0 .. n - 2]] + negateIf test = if test then AlgAdd.negate else id + signedSubresultants = zipWith negateIf delta spraySubresultants + +-- | Principal Sturm-Habicht sequence of a spray. This function calls `sturmHabicht` +-- sequence and then it can be slow. +principalSturmHabichtSequence :: + (Eq a, AlgRing.C a) + => Int -- ^ index of the variable with respect to which the Sturm-Habicht sequence will be computed (e.g. 2 for @y@) + -> Spray a + -> [Spray a] +principalSturmHabichtSequence var spray + | var < 1 || var > d + = error "principalSturmHabichtSequence: invalid variable index." + | otherwise = [jcoeff j | j <- [0 .. length sHS - 1]] + where + d = numberOfVariables spray + sHS = sturmHabichtSequence var spray + permutation = [d - var + 1 .. d] ++ [1 .. d - var] + permutation' = [var + 1 .. d] ++ [1 .. var] + jcoeff j = if isZeroSpray sHSj || j >= length coeffs + then zeroSpray + else permuteVariables permutation' (coeffs !! j) + where + sHSj = sHS !! j + p = permuteVariables permutation sHSj + coeffs = reverse $ sprayCoefficients' d p -- what if numberOfVariables p /= d ? + + +-- Number of real roots in interval ------------------------------------------- + +_signVariations :: (a -> Char) -> [a] -> Int +_signVariations signFunc as = v1 + v2 + 2*v3 + where + count x xs = sum (map (fromEnum . (== x)) xs) + l = length as + signs = map signFunc as + chunks2 = [(signs !! i, signs !! (i+1)) | i <- [0 .. l-2]] + v1 = count ('+', '-') chunks2 + + count ('-', '+') chunks2 + chunks3 = [(signs !! i, signs !! (i+1), signs !! (i+2)) | i <- [0 .. l-3]] + chunks4 = [(signs !! i, signs !! (i+1), signs !! (i+2), signs !! (i+3)) + | i <- [0 .. l-4]] + v2 = count ('-', '0', '+') chunks3 + + count ('+', '0', '-') chunks3 + + count ('+', '0', '0', '-') chunks4 + + count ('-', '0', '0', '+') chunks4 + v3 = count ('+', '0', '0', '+') chunks4 + + count ('-', '0', '0', '-') chunks4 + +signVariations :: (Eq a, Num a) => [a] -> Int +signVariations = _signVariations signFunc + where + signFunc a + | signum a == 0 = '0' + | signum a == 1 = '+' + | otherwise = '-' + +signVariations' :: (Eq a, AlgAbs.C a) => [a] -> Int +signVariations' = _signVariations signFunc + where + signFunc a + | AlgAbs.signum a == AlgAdd.zero = '0' + | AlgAbs.signum a == AlgRing.one = '+' + | otherwise = '-' + +_numberOfRealRootsInOpenInterval :: + (AlgRing.C a, Ord a) => ([a] -> Int) -> Spray a -> Maybe a -> Maybe a -> Int +_numberOfRealRootsInOpenInterval signVariationsFunc spray alpha beta + | isConstantSpray spray = + if isZeroSpray spray + then error "numberOfRealRootsInOpenInterval: the spray is null." + else 0 + | alphaIsJust && betaIsJust = + if alpha' == beta' + then 0 + else + if alpha > beta + then + error "numberOfRealRootsInOpenInterval: the bounds are not ordered." + else + if isZeroAtBeta then svDiff - 1 else svDiff + | alphaIsJust = svDiff + | betaIsJust = if isZeroAtBeta then svDiff - 1 else svDiff + | otherwise = svDiff + where + alphaIsJust = isJust alpha + betaIsJust = isJust beta + alpha' = if alphaIsJust then fromJust alpha else undefined + beta' = if betaIsJust then fromJust beta else undefined + evaluateAtMinusInfinity p = if isConstantSpray p + then getConstantTerm p + else let + (deg, lCoeff) = first ((`S.index` 0) . exponents) (leadingTerm p) in + if even deg then lCoeff else AlgAdd.negate lCoeff + evaluateAtPlusInfinity p = if isConstantSpray p + then getConstantTerm p + else (snd . leadingTerm) p + (ginit, glast) = + fromJust $ unsnoc $ filter (not . isZeroSpray) (sturmHabichtSequence 1 spray) + eval1 = + if alphaIsJust then evaluateAt [alpha'] else evaluateAtMinusInfinity + eval2 = + if betaIsJust then evaluateAt [beta'] else evaluateAtPlusInfinity + sprayAtAlpha = eval1 glast + sprayAtBeta = eval2 glast + isZeroAtBeta = sprayAtBeta == AlgAdd.zero + galpha = map eval1 ginit ++ [sprayAtAlpha] + gbeta = map eval2 ginit ++ [sprayAtBeta] + svAtAlpha = signVariationsFunc galpha + svAtBeta = signVariationsFunc gbeta + svDiff = svAtAlpha - svAtBeta + +_numberOfRealRootsInClosedInterval :: + (AlgRing.C a, Ord a) => ([a] -> Int) -> Spray a -> Maybe a -> Maybe a -> Int +_numberOfRealRootsInClosedInterval signVariationsFunc spray alpha beta = + _numberOfRealRootsInOpenInterval signVariationsFunc spray alpha beta + toAdd + where + alphaIsJust = isJust alpha + betaIsJust = isJust beta + alpha' = if alphaIsJust then fromJust alpha else undefined + beta' = if betaIsJust then fromJust beta else undefined + isZeroAtAlpha = alphaIsJust && evaluateAt [alpha'] spray == AlgAdd.zero + isZeroAtBeta = betaIsJust && evaluateAt [beta'] spray == AlgAdd.zero + toAdd = + if alpha == beta + then fromEnum isZeroAtAlpha + else fromEnum isZeroAtAlpha + fromEnum isZeroAtBeta + +-- | Number of real roots of a spray in an open interval (that makes sense +-- only for a spray on a ring embeddable in the real numbers). +numberOfRealRootsInOpenInterval :: + (Num a, AlgRing.C a, Ord a) + => Spray a -- ^ a spray + -> Maybe a -- ^ lower bound of the interval; use @Just@ for a finite bound, and @Nothing@ for minus infinity + -> Maybe a -- ^ upper bound of the interval; use @Just@ for a finite bound, and @Nothing@ for minus infinity + -> Int +numberOfRealRootsInOpenInterval spray = + if isUnivariate spray + then _numberOfRealRootsInOpenInterval signVariations spray + else error "numberOfRealRootsInOpenInterval: the spray is not univariate." + +-- | Number of real roots of a spray in a closed interval (that makes sense +-- only for a spray on a ring embeddable in the real numbers). The roots are +-- not counted with their multiplicity. +numberOfRealRootsInClosedInterval :: + (Num a, AlgRing.C a, Ord a) + => Spray a -- ^ a spray + -> Maybe a -- ^ lower bound of the interval; use @Just@ for a finite bound, and @Nothing@ for minus infinity + -> Maybe a -- ^ upper bound of the interval; use @Just@ for a finite bound, and @Nothing@ for minus infinity + -> Int +numberOfRealRootsInClosedInterval spray = + if isUnivariate spray + then _numberOfRealRootsInClosedInterval signVariations spray + else error "numberOfRealRootsInClosedInterval: the spray is not univariate." + +-- | Number of real roots of a spray in an open interval (that makes sense +-- only for a spray on a ring embeddable in the real numbers). +numberOfRealRootsInOpenInterval' :: + (AlgAbs.C a, Ord a) + => Spray a -- ^ a spray + -> Maybe a -- ^ lower bound of the interval; use @Just@ for a finite bound, and @Nothing@ for minus infinity + -> Maybe a -- ^ upper bound of the interval; use @Just@ for a finite bound, and @Nothing@ for minus infinity + -> Int +numberOfRealRootsInOpenInterval' spray = + if isUnivariate spray + then _numberOfRealRootsInOpenInterval signVariations' spray + else error "numberOfRealRootsInOpenInterval': the spray is not univariate." + +-- | Number of real roots of a spray in a closed interval (that makes sense +-- only for a spray on a ring embeddable in the real numbers). The roots are +-- not counted with their multiplicity. +numberOfRealRootsInClosedInterval' :: + (AlgAbs.C a, Ord a) + => Spray a -- ^ a spray + -> Maybe a -- ^ lower bound of the interval; use @Just@ for a finite bound, and @Nothing@ for minus infinity + -> Maybe a -- ^ upper bound of the interval; use @Just@ for a finite bound, and @Nothing@ for minus infinity + -> Int +numberOfRealRootsInClosedInterval' spray = + if isUnivariate spray + then _numberOfRealRootsInClosedInterval signVariations' spray + else error "numberOfRealRootsInClosedInterval': the spray is not univariate." + + +-- Number of real roots ------------------------------------------------------- + +runLengthEncoding :: Eq a => [a] -> [(a,Int)] +runLengthEncoding = foldr code [] + where + code c [] = [(c,1)] + code c ((x,n):ts) + | c == x = (x,n+1):ts + | otherwise = (c,1):(x,n):ts + +_signPermanencesAndVariations :: (a -> Char) -> [a] -> (Int, Int) +_signPermanencesAndVariations signFunc as = (permanences, variations) + where + signs = map signFunc as + rle = runLengthEncoding signs + l = length rle + lengths = map snd rle + permanences = sum lengths - l + variations = l - 1 + +signPermanencesAndVariations :: (Eq a, Num a) => [a] -> (Int, Int) +signPermanencesAndVariations = _signPermanencesAndVariations signFunc + where + signFunc a = if signum a == 1 then '+' else '-' + +signPermanencesAndVariations' :: (Eq a, AlgAbs.C a) => [a] -> (Int, Int) +signPermanencesAndVariations' = _signPermanencesAndVariations signFunc + where + signFunc a = if AlgAbs.signum a == AlgRing.one then '+' else '-' + +distributionsOfZeros :: (Eq a, AlgAdd.C a) => [a] -> ([Int], [Int], [Int]) +distributionsOfZeros as = (i_, k_, ik_) + where + symbol a = if a == AlgAdd.zero then '0' else '*' + symbols = map symbol as ++ ['0'] + lengths = map snd (runLengthEncoding symbols) + l = length lengths + cumulativeLengths = scanl (+) 0 lengths + range = [0 .. l-1] + ik_ = [cumulativeLengths !! i | i <- range, even i] + i_ = [cumulativeLengths !! i - 1 | i <- range, odd i] + k_ = [lengths !! i | i <- [0 .. l-2], odd i] ++ [lengths !! (l-1) - 1] + +_blocksAndEpsilons :: (Eq a, AlgAdd.C a) => (a -> Int) -> [a] -> ([[a]], [Int]) +_blocksAndEpsilons signFunc as = (blocks, epsilons) + where + (i_, k_, ik_) = distributionsOfZeros as + t = length i_ + blocks = [[as !! m | m <- [ik_ !! n .. i_ !! n]] | n <- [0 .. t-1]] + epsilon s = let ks = k_ !! s in if odd ks + then + 0 + else + (if even (ks `div` 2) then 1 else -1) * + signFunc (as !! (ik_ !! (s+1))) * signFunc (as !! (i_ !! s)) + epsilons = [epsilon n | n <- [0 .. t-2]] + +blocksAndEpsilons :: (Eq a, AlgAdd.C a, Num a) => [a] -> ([[a]], [Int]) +blocksAndEpsilons = _blocksAndEpsilons signFunc + where + signFunc a = if signum a == 1 then 1 else -1 + +blocksAndEpsilons' :: (Eq a, AlgAbs.C a) => [a] -> ([[a]], [Int]) +blocksAndEpsilons' = _blocksAndEpsilons signFunc + where + signFunc a = if AlgAbs.signum a == AlgRing.one then 1 else -1 + +_C :: ([a] -> (Int, Int)) -> ([[a]], [Int]) -> Int +_C signPermanencesAndVariationsFunc (blocks, epsilons) = sum pvs - sum epsilons + where + permanencesMinusVariations as = + let (p, v) = signPermanencesAndVariationsFunc as in p - v + pvs = map permanencesMinusVariations blocks + +-- | Number of real roots of a spray (that makes sense only for a spray on a +-- ring embeddable in the real numbers). The roots are not counted with their +-- multiplicity. +numberOfRealRoots :: (Eq a, AlgRing.C a, Num a) => Spray a -> Int +numberOfRealRoots spray = + if isUnivariate spray + then + if isConstantSpray spray + then if isZeroSpray spray + then error "numberOfRealRoots: the spray is null." + else 0 + else _C signPermanencesAndVariations (blocksAndEpsilons $ reverse as) + else + error "numberOfRealRoots: the spray is not univariate." + where + as = map getConstantTerm (principalSturmHabichtSequence 1 spray) + +-- | Number of real roots of a spray (that makes sense only for a spray on a +-- ring embeddable in the real numbers). The roots are not counted with their +-- multiplicity. +numberOfRealRoots' :: (Eq a, AlgAbs.C a) => Spray a -> Int +numberOfRealRoots' spray = + if isUnivariate spray + then + if isConstantSpray spray + then if isZeroSpray spray + then error "numberOfRealRoots: the spray is null." + else 0 + else _C signPermanencesAndVariations' (blocksAndEpsilons' $ reverse as) + else + error "numberOfRealRoots': the spray is not univariate." + where + as = map getConstantTerm (principalSturmHabichtSequence 1 spray) + + +-- GCD stuff ------------------------------------------------------------------ + +-- | the spray coefficients of a spray as a univariate spray in x_n +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 + permutation = [2 .. n] ++ [1] + spray' = safeSpray $ permuteVariables permutation spray + spray'' = removeConstantTerm'' spray' + (expnts', coeffs') = unzip (HM.toList spray'') + constantTerm = getConstantTerm'' spray' + xpows = map (`index` 0) expnts' + powers'' = map (powerize . S.deleteAt 0) expnts' + sprays'' = zipWith (curry fromTerm) powers'' coeffs' + imap = IM.fromListWith (^+^) (zip xpows sprays'') + imap' = IM.insertWith (^+^) 0 (constantSpray constantTerm) imap + deg = maximum xpows + sprays = [ + fromMaybe zeroSpray (IM.lookup i imap') + | i <- [deg, deg-1 .. 0] + ] + +-- | the degree of a spray as a univariate spray in x_n with spray coefficients +degree :: (Eq a, AlgRing.C a) => Int -> Spray a -> Int +degree n spray + | isConstant spray = + if isZeroSpray spray + then minBound -- (should not happen) + else 0 + | numberOfVariables spray /= n = 0 + | otherwise = maximum xpows + where + permutation = [2 .. n] ++ [1] + spray' = permuteVariables permutation spray + powers' = HM.keys $ removeConstantTerm spray' + xpows = map ((`index` 0) . exponents) powers' + +-- | 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 -- (should not happen) + else 0, + constantSpray constantTerm + ) + | numberOfVariables spray /= n = (0, spray) + | otherwise = (deg, leadingCoeff) + where + permutation = [2 .. n] ++ [1] + spray' = permuteVariables permutation spray + constantTerm = getConstantTerm spray' + spray'' = removeConstantTerm spray' + (powers', coeffs') = unzip (HM.toList spray'') + expnts' = map exponents powers' + xpows = map (`index` 0) expnts' + deg = maximum xpows + is = elemIndices deg xpows + powers'' = [powerize (S.deleteAt 0 (expnts' !! i)) | i <- is] + coeffs'' = [coeffs' !! i | i <- is] + leadingCoeff = sumTerms (zip powers'' coeffs'') + +-- | Pseudo-division of two sprays @A@ and @B@ such that @deg(A) >= deg(B)@ +-- where @deg@ is the degree with respect to the outermost variable. +pseudoDivision :: (Eq a, AlgRing.C a) + => Spray a -- ^ A + -> Spray a -- ^ B + -> (Spray a, (Spray a, Spray a)) -- ^ (C, (Q, R)) such that C^*^A = B^*^Q ^+^ R +pseudoDivision sprayA sprayB = + pseudoDivision' (max (numberOfVariables sprayA) (numberOfVariables sprayB)) + sprayA sprayB + +-- | 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' 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 n sprayA + (degB, ellB) = degreeAndLeadingCoefficient n sprayB + delta = degA - degB + 1 + go sprayR sprayQ e = + if isZeroSpray sprayR || degR < degB + then (q ^*^ sprayQ, q ^*^ sprayR) + else go (ellB ^*^ sprayR ^-^ sprayS ^*^ sprayB) + (ellB ^*^ sprayQ ^+^ sprayS) + (e - 1) + where + (degR, ellR) = degreeAndLeadingCoefficient n sprayR + q = ellB ^**^ e + sprayS = multSprayByTerm ellR (loneTerm' n (degR - degB)) + +-- | 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 + | isZeroSpray sprayB = normalize sprayA + | otherwise = go sprayA' sprayB' unitSpray unitSpray + where + gcdKX0 :: Spray a -> Spray a -> a + gcdKX0 = const $ const AlgRing.one + n' = 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 isZeroSpray remainder + then quo + else error "exactDivisionBy: should not happen." + where + (quo, remainder) = sprayDivision a b + reduceSpray :: Spray a -> Spray a + reduceSpray spray = exactDivisionBy (content spray) spray + normalize :: Spray a -> Spray a + normalize spray = + let leadingCoeff = snd $ leadingTerm spray in spray /> leadingCoeff + 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 + | isZeroSpray sprayR = d ^*^ reduceSpray (normalize sprayB'') + | numberOfVariables sprayR == 0 = d + | otherwise = go sprayB'' + (exactDivisionBy (g ^*^ h^**^delta) sprayR) + ellB'' ---- ellA'' + (exactDivisionBy (h^**^delta) (h ^*^ ellB''^**^delta)) --- g^**^delta + where + (_, (_, sprayR)) = pseudoDivision' n' sprayA'' sprayB'' + -- (degA'', ellA'') = degreeAndLeadingCoefficient n' sprayA'' + -- degB'' = degree n' sprayB'' + degA'' = degree n' sprayA'' + (degB'', ellB'') = degreeAndLeadingCoefficient 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) + + +-- Matrices ------------------------------------------------------------------- + +-- | Determinant of a matrix with entries in a ring by using Laplace +-- expansion (this is slow); the __numeric-prelude__ package provides some +-- stuff to deal with matrices over a ring but it does not provide the +-- determinant +detLaplace :: forall a. (Eq a, AlgRing.C a) => Matrix a -> a +detLaplace b = + if nrows b == ncols b + then detUnsafe b + else error "detLaplace: the matrix is not square." + where + detUnsafe m = if nrows m == 1 + then + m DM.! (1,1) + else + suml1 + [negateIf i (times (m DM.! (i,1)) (detUnsafe (minorMatrix i 1 m))) + | i <- [1 .. nrows m]] + 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 + +-- | Determinant of a matrix over a ring by using Laplace expansion; this is +-- the same as `detLaplace` but for a matrix from the __numeric-prelude__ +-- package +detLaplace' :: forall a. (Eq a, AlgRing.C a) => MathMatrix.T a -> a +detLaplace' m = detLaplace (DM.fromLists $ MathMatrix.rows m) + +-- | Characteristic polynomial of a square matrix +-- +-- >>> import Data.Matrix (Matrix, fromLists) +-- >>> m = fromLists [ [12, 16, 4] +-- >>> , [16, 2, 8] +-- >>> , [8, 18, 10] ] :: Matrix Int +-- >>> spray = characteristicPolynomial m +-- >>> putStrLn $ prettyNumSpray spray +-- -x^3 + 24*x^2 + 268*x - 1936 +characteristicPolynomial :: (Eq a, AlgRing.C a) => Matrix a -> Spray a +characteristicPolynomial m = + if nrows m /= ncols m + then error "characteristicPolynomial: the matrix is not square." + else detLaplace m' + where + m' = DM.mapPos f m + f (i, j) mij = if i == j + then constantSpray mij ^-^ x + else constantSpray mij + x = lone 1 + + +-- Ratios of sprays ----------------------------------------------------------- + +-- | A @RatioOfSprays a@ object represents a fraction of two multivariate +-- polynomials whose coefficients are of type @a@, which represents a field. +-- These two polynomials are represented by two @Spray a@ objects. Generally +-- we do not use this constructor to build a ratio of sprays: we use the `%//%` +-- operator instead, because it always returns an irreducible ratio of sprays, +-- meaning that its corresponding fraction of polynomials is irreducible, i.e. +-- its numerator and its denominator are coprime. You can use this constructor +-- if you are /sure/ that the numerator and the denominator are coprime. This +-- can save some computation time, but unfortunate consequences can occur if +-- the numerator and the denominator are not coprime. An arithmetic operation +-- on ratios of sprays always returns an irreducible ratio of sprays under the +-- condition that the ratios of sprays it involves are irreducible. Moreover, +-- it never returns a ratio of sprays with a constant denominator other than +-- the unit spray. If you use this constructor with a constant denominator, +-- always set this denominator to the unit spray (by dividing the numerator +-- by the constant value of the denominator). +data RatioOfSprays a = RatioOfSprays + { _numerator :: Spray a + , _denominator :: Spray a + } + deriving Show + +type RatioOfQSprays = RatioOfSprays Rational + +instance (Eq a, AlgField.C a) => FunctionLike (RatioOfSprays a) where + type BaseRing (RatioOfSprays a) = a + -- + type VariablesType (RatioOfSprays a) = Spray a + -- + (*^) :: a -> RatioOfSprays a -> RatioOfSprays a + (*^) lambda rOS = lambda AlgMod.*> rOS + -- + (+>) :: a -> RatioOfSprays a -> RatioOfSprays a + (+>) x rOS = rOS AlgAdd.+ asRatioOfSprays (constantSpray x) + -- + substitute :: [Maybe a] -> RatioOfSprays a -> RatioOfSprays a + substitute subs (RatioOfSprays p q) = + substitute subs p %//% substitute subs q + -- + evaluate :: RatioOfSprays a -> [a] -> a + evaluate (RatioOfSprays p q) xyz = evaluate p xyz AlgField./ evaluate q xyz + -- + changeVariables :: RatioOfSprays a -> [Spray a] -> RatioOfSprays a + changeVariables rOS newVariables = + if length newVariables < numberOfVariables rOS + then + error "changeVariables: not enough new variables provided." + else + changeVariables (_numerator rOS) newVariables + %//% changeVariables (_denominator rOS) newVariables + -- + numberOfVariables :: RatioOfSprays a -> Int + numberOfVariables (RatioOfSprays p q) = + max (numberOfVariables p) (numberOfVariables q) + -- + permuteVariables :: [Int] -> RatioOfSprays a -> RatioOfSprays a + permuteVariables permutation (RatioOfSprays p q) = + permuteVariables permutation p %//% permuteVariables permutation q + -- + swapVariables :: (Int, Int) -> RatioOfSprays a -> RatioOfSprays a + 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 + p' = derivative i p + q' = derivative i q + +-- | quotients of two univariate sprays by their gcd +-- we use `sprayDivision0` because this function is called +-- (by `irreducibleFraction`) with non-constant sprays only +quotientsByGCD :: + (Eq a, AlgField.C a) => Spray a -> Spray a -> (Spray a, Spray a) +quotientsByGCD sprayA sprayB = + if isUnivariate sprayA && isUnivariate sprayB + then + go sprayA sprayB unitSpray zeroSpray zeroSpray unitSpray + else + (exactDivision sprayA g, exactDivision sprayB g) + where + exactDivision p q = fst (sprayDivision0 p q) + g = gcdSpray sprayA sprayB + go oldr r olds s oldt t + | isZeroSpray r = (AlgAdd.negate t, s) + | otherwise = + go r remainder s (olds ^-^ quo ^*^ s) t (oldt ^-^ quo ^*^ t) + where + (quo, remainder) = univariateSprayDivision oldr r + +-- | irreducible fraction of sprays +irreducibleFraction :: + (Eq a, AlgField.C a) => Spray a -> Spray a -> RatioOfSprays a +irreducibleFraction p q = adjustFraction rOS + where + rOS = if isConstant p || isConstant q + then RatioOfSprays p q + else let (a, b) = quotientsByGCD p q in RatioOfSprays a b + +-- | make monic denominator +adjustFraction :: (Eq a, AlgField.C a) => RatioOfSprays a -> RatioOfSprays a +adjustFraction (RatioOfSprays p q) = + let leadingCoeff = snd $ leadingTerm q in + RatioOfSprays (p /^ leadingCoeff) (q /^ leadingCoeff) + +instance (AlgRing.C a, Eq a) => Eq (RatioOfSprays a) where + (==) :: RatioOfSprays a -> RatioOfSprays a -> Bool + (==) (RatioOfSprays p q) (RatioOfSprays p' q') = + isZeroSpray (p ^*^ q' ^-^ p' ^*^ q) + +instance (AlgField.C a, Eq a) => AlgAdd.C (RatioOfSprays a) where + (+) :: RatioOfSprays a -> RatioOfSprays a -> RatioOfSprays a + (+) (RatioOfSprays p q) (RatioOfSprays p' q') = + irreducibleFraction (p ^*^ q' ^+^ p' ^*^ q) (q ^*^ q') + zero :: RatioOfSprays a + zero = RatioOfSprays zeroSpray unitSpray + negate :: RatioOfSprays a -> RatioOfSprays a + negate (RatioOfSprays p q) = RatioOfSprays (negateSpray p) q + +instance (AlgField.C a, Eq a) => AlgMod.C a (RatioOfSprays a) where + (*>) :: a -> RatioOfSprays a -> RatioOfSprays a + lambda *> (RatioOfSprays p q) = RatioOfSprays (lambda *^ p) q + +instance (AlgField.C a, Eq a) => AlgRightMod.C a (RatioOfSprays a) where + (<*) :: RatioOfSprays a -> a -> RatioOfSprays a + rOS <* lambda = lambda AlgMod.*> rOS + +instance (AlgField.C a, Eq a) => AlgMod.C (Spray a) (RatioOfSprays a) where + (*>) :: Spray a -> RatioOfSprays a -> RatioOfSprays a + spray *> (RatioOfSprays p q) = irreducibleFraction (spray ^*^ p) q + +instance (AlgField.C a, Eq a) => AlgRightMod.C (Spray a) (RatioOfSprays a) where + (<*) :: RatioOfSprays a -> Spray a -> RatioOfSprays a + rOS <* spray = spray AlgMod.*> rOS + +instance (AlgField.C a, Eq a) => AlgRing.C (RatioOfSprays a) where + (*) :: RatioOfSprays a -> RatioOfSprays a -> RatioOfSprays a + (*) (RatioOfSprays p q) (RatioOfSprays p' q') = + irreducibleFraction (p ^*^ p') (q ^*^ q') + (^) :: RatioOfSprays a -> Integer -> RatioOfSprays a + (^) (RatioOfSprays p q) n = if n >= 0 + then RatioOfSprays (p AlgRing.^ n) (q AlgRing.^ n) + else RatioOfSprays (q AlgRing.^ (-n)) (p AlgRing.^ (-n)) + one :: RatioOfSprays a + one = RatioOfSprays unitSpray unitSpray + +instance (AlgField.C a, Eq a) => AlgField.C (RatioOfSprays a) where + recip :: RatioOfSprays a -> RatioOfSprays a + recip (RatioOfSprays p q) = RatioOfSprays q p + +infixl 7 %:% +-- | Ratio of sprays from numerator and denominator, +-- __without reducing the fraction__ +(%:%) :: Spray a -> Spray a -> RatioOfSprays a +(%:%) = RatioOfSprays + +infixl 7 %//% +-- | Irreducible ratio of sprays from numerator and denominator; alias of @(^/^)@ +(%//%) :: (Eq a, AlgField.C a) => Spray a -> Spray a -> RatioOfSprays a +(%//%) = irreducibleFraction + +infixl 7 ^/^ +-- | Irreducible ratio of sprays from numerator and denominator; alias of @(%//%)@ +(^/^) :: (Eq a, AlgField.C a) => Spray a -> Spray a -> RatioOfSprays a +(^/^) = irreducibleFraction + +infixl 7 %/% +-- | Division of a ratio of sprays by a spray; the result is an +-- irreducible fraction +(%/%) :: (Eq a, AlgField.C a) => RatioOfSprays a -> Spray a -> RatioOfSprays a +(%/%) rOS spray = rOS AlgRing.* RatioOfSprays unitSpray spray + +-- | Whether a ratio of sprays is constant; same as `isConstant` +isConstantRatioOfSprays :: (Eq a, AlgField.C a) => RatioOfSprays a -> Bool +isConstantRatioOfSprays = isConstant + +-- | Whether a ratio of sprays actually is polynomial, that is, whether its +-- denominator is a constant spray (and then it should be the unit spray) +-- +-- >>> x = qlone 1 +-- >>> y = qlone 2 +-- >>> p = x^**^4 ^-^ y^**^4 +-- >>> q = x ^-^ y +-- >>> isPolynomialRatioOfSprays $ p %//% q +-- True +-- >>> isPolynomialRatioOfSprays $ p %:% q +-- False +isPolynomialRatioOfSprays :: (Eq a, AlgRing.C a) => RatioOfSprays a -> Bool +isPolynomialRatioOfSprays = isConstant . _denominator + +-- | The null ratio of sprays +zeroRatioOfSprays, zeroROS :: (AlgField.C a, Eq a) => RatioOfSprays a +zeroRatioOfSprays = AlgAdd.zero +zeroROS = AlgAdd.zero + +-- | The unit ratio of sprays +unitRatioOfSprays, unitROS :: (AlgField.C a, Eq a) => RatioOfSprays a +unitRatioOfSprays = AlgRing.one +unitROS = AlgRing.one + +-- | Constant ratio of sprays +constantRatioOfSprays :: (Eq a, AlgRing.C a) => a -> RatioOfSprays a +constantRatioOfSprays x = asRatioOfSprays (constantSpray x) + +-- | Evaluates a ratio of sprays; same as `evaluate` +evalRatioOfSprays :: (Eq a, AlgField.C a) => RatioOfSprays a -> [a] -> a +evalRatioOfSprays = evaluate + +-- | Substitutes the variables of a ratio of sprays with some ratios of sprays +evalRatioOfSprays' :: + (Eq a, AlgField.C a) => RatioOfSprays a -> [RatioOfSprays a] -> RatioOfSprays a +evalRatioOfSprays' (RatioOfSprays spray1 spray2) listOfROS = + if length listOfROS >= + max (numberOfVariables spray1) (numberOfVariables spray2) + then + eval_spray spray1 AlgField./ eval_spray spray2 + else + error "evalRatioOfSprays': not enough ratios of sprays provided." + where + seqOfROS = S.fromList listOfROS + evalTermOfSpray powers coeff = + coeff AlgMod.*> + foldl' + (AlgRing.*) + unitRatioOfSprays + (S.zipWith + (\ros e -> + ros AlgRing.^ (fromIntegral e)) seqOfROS (exponents powers) + ) + eval_spray spray = + HM.foldl' + (AlgAdd.+) + zeroRatioOfSprays + (HM.mapWithKey evalTermOfSpray spray) + +-- | Substitutes some values to some variables of a ratio of sprays; same as `substitute` +substituteRatioOfSprays :: + (Eq a, AlgField.C a) => [Maybe a] -> RatioOfSprays a -> RatioOfSprays a +substituteRatioOfSprays = substitute + +-- | Coerces a spray to a ratio of sprays +asRatioOfSprays :: AlgRing.C a => Spray a -> RatioOfSprays a +asRatioOfSprays spray = RatioOfSprays spray unitSpray + +-- | Converts a ratio of polynomials to a ratio of sprays +fromRatioOfPolynomials :: + (Eq a, AlgRing.C a) => RatioOfPolynomials a -> RatioOfSprays a +fromRatioOfPolynomials rop = + RatioOfSprays + (polynomialToSpray $ NumberRatio.numerator rop) + (polynomialToSpray $ NumberRatio.denominator rop) + +-- | Converts a ratio of rational polynomials to a ratio of rational sprays; +-- this is not a specialization of `fromRatioOfPolynomials` because +-- @RatioOfQPolynomials@ is @RatioOfPolynomials a@ with +-- @a = Rational'@, not with @a = Rational@ +fromRatioOfQPolynomials :: RatioOfQPolynomials -> RatioOfQSprays +fromRatioOfQPolynomials rop = + RatioOfSprays + (qPolynomialToQSpray $ NumberRatio.numerator rop) + (qPolynomialToQSpray $ NumberRatio.denominator rop) + +-- | General function to print a `RatioOfSprays` object +showRatioOfSprays :: (Eq a, AlgRing.C a) + => ((Spray a, Spray a) -> (String, String)) -- ^ function which prints a pair of sprays that will be applied to the numerator and the denominator + -> (String, String) -- ^ pair of braces to enclose the numerator and the denominator + -> String -- ^ represents the quotient bar + -> RatioOfSprays a + -> String +showRatioOfSprays spraysShower braces quotientBar (RatioOfSprays p q) = + numeratorString ++ denominatorString + where + enclose = bracify braces + (pString, qString) = spraysShower (p, q) + numeratorString = enclose pString + denominatorString = if q == unitSpray + then "" + else quotientBar ++ enclose qString + +showTwoSpraysXYZ :: (Eq a, AlgRing.C a) + => (a -> String) -- ^ function mapping a coefficient to a string, typically 'show' + -> (String, String) -- ^ used to enclose the coefficients, usually a pair of braces + -> [String] -- ^ typically some letters, to print the variables + -> (Spray a, Spray a) -- ^ the two sprays to be printed + -> (String, String) +showTwoSpraysXYZ showCoef braces letters (spray1, spray2) = + both (showSpray showCoef braces showMonomials) (spray1, spray2) + where + n = max (numberOfVariables spray1) (numberOfVariables spray2) + showMonomials = map (unpack . showMonomialXYZ letters n) + +showTwoSpraysX1X2X3 :: + (a -> String) -- ^ function mapping a coefficient to a string, typically 'show' + -> (String, String) -- ^ used to enclose the coefficients, usually a pair of braces + -> String -- ^ typically a letter, to print the non-indexed variables + -> (Spray a, Spray a) -- ^ the two sprays to be printed + -> (String, String) +showTwoSpraysX1X2X3 showCoef braces letter (spray1, spray2) = + both (showSpray showCoef braces showMonomials) (spray1, spray2) + where + showMonomials = showMonomialsX1X2X3 letter + +showTwoNumSprays :: (Num a, Ord a) + => (a -> String) -- ^ function mapping a positive coefficient to a string + -> ([Exponents] -> [String]) -- ^ prints the monomials + -> (Spray a, Spray a) -- ^ the two sprays to be printed + -> (String, String) +showTwoNumSprays showPositiveCoef showMonomials = + both (showNumSpray showMonomials showPositiveCoef) + +showTwoQSprays :: + ([Exponents] -> [String]) -- ^ prints the monomials + -> (QSpray, QSpray) -- ^ the two sprays to be printed + -> (String, String) +showTwoQSprays = showTwoNumSprays showRatio + +showTwoNumSpraysXYZ :: (AlgRing.C a, Num a, Ord a) + => (a -> String) -- ^ function mapping a positive coefficient to a string + -> [String] -- ^ typically some letters, to print the variables + -> (Spray a, Spray a) -- ^ the two sprays to be printed + -> (String, String) +showTwoNumSpraysXYZ showPositiveCoef letters (spray1, spray2) = + showTwoNumSprays showPositiveCoef showMonomials (spray1, spray2) + where + n = max (numberOfVariables spray1) (numberOfVariables spray2) + showMonomials = map (unpack . showMonomialXYZ letters n) + +showTwoQSpraysXYZ :: + [String] -- ^ typically some letters, to print the variables + -> (QSpray, QSpray) -- ^ the two sprays to be printed + -> (String, String) +showTwoQSpraysXYZ = showTwoNumSpraysXYZ showRatio + +showTwoNumSpraysX1X2X3 :: (Num a, Ord a) + => (a -> String) -- ^ function mapping a positive coefficient to a string + -> String -- ^ typically a letter, to print the non-indexed variable + -> (Spray a, Spray a) -- ^ the two sprays to be printed + -> (String, String) +showTwoNumSpraysX1X2X3 showPositiveCoef letter (spray1, spray2) = + showTwoNumSprays showPositiveCoef showMonomials (spray1, spray2) + where + showMonomials = showMonomialsX1X2X3 letter + +showTwoQSpraysX1X2X3 :: + String -- ^ typically a letter, to print the non-indexed variables + -> (QSpray, QSpray) -- ^ the two sprays to be printed + -> (String, String) +showTwoQSpraysX1X2X3 = showTwoNumSpraysX1X2X3 showRatio + +-- | Prints a ratio of sprays with numeric coefficients +showRatioOfNumSprays :: (Num a, Ord a, AlgRing.C a) + => (a -> String) -- ^ function mapping a positive coefficient to a string + -> ([Exponents] -> [String]) -- ^ prints the monomials + -> (String, String) -- ^ pair of braces to enclose the numerator and the denominator + -> String -- ^ represents the quotient bar + -> RatioOfSprays a + -> String +showRatioOfNumSprays showPositiveCoef showMonomials = + showRatioOfSprays (showTwoNumSprays showPositiveCoef showMonomials) + +-- | Prints a ratio of sprays with rational coefficients +showRatioOfQSprays :: + ([Exponents] -> [String]) -- ^ prints the monomials + -> (String, String) -- ^ pair of braces to enclose the numerator and the denominator + -> String -- ^ represents the quotient bar + -> RatioOfQSprays + -> String +showRatioOfQSprays showMonomials = + showRatioOfSprays (showTwoQSprays showMonomials) + +-- | Prints a ratio of sprays with numeric coefficients +showRatioOfNumSpraysXYZ :: (Num a, Ord a, AlgRing.C a) + => (a -> String) -- ^ function mapping a positive coefficient to a string + -> [String] -- ^ typically some letters, to print the variables + -> (String, String) -- ^ pair of braces to enclose the numerator and the denominator + -> String -- ^ represents the quotient bar + -> RatioOfSprays a + -> String +showRatioOfNumSpraysXYZ showPositiveCoef letters = + showRatioOfSprays (showTwoNumSpraysXYZ showPositiveCoef letters) + +-- | Prints a ratio of sprays with numeric coefficients +showRatioOfNumSpraysX1X2X3 :: (Num a, Ord a, AlgRing.C a) + => (a -> String) -- ^ function mapping a positive coefficient to a string + -> String -- ^ typically a letter, to print the variables + -> (String, String) -- ^ pair of braces to enclose the numerator and the denominator + -> String -- ^ represents the quotient bar + -> RatioOfSprays a + -> String +showRatioOfNumSpraysX1X2X3 showPositiveCoef letter = + showRatioOfSprays (showTwoNumSpraysX1X2X3 showPositiveCoef letter) + +-- | Prints a ratio of sprays with rational coefficients +showRatioOfQSpraysXYZ :: + [String] -- ^ typically some letters, to print the variables + -> (String, String) -- ^ pair of braces to enclose the numerator and the denominator + -> String -- ^ represents the quotient bar + -> RatioOfQSprays + -> String +showRatioOfQSpraysXYZ letters = showRatioOfSprays (showTwoQSpraysXYZ letters) + +-- | Prints a ratio of sprays with rational coefficients +showRatioOfQSpraysX1X2X3 :: + String -- ^ typically a letter, to print the variables + -> (String, String) -- ^ pair of braces to enclose the numerator and the denominator + -> String -- ^ represents the quotient bar + -> RatioOfQSprays + -> String +showRatioOfQSpraysX1X2X3 letter = showRatioOfSprays (showTwoQSpraysX1X2X3 letter) + +-- | Prints a ratio of sprays +showRatioOfSpraysXYZ :: forall a. (Eq a, AlgField.C a) + => [String] -- ^ typically some letters, to represent the variables + -> (a -> String) -- ^ function mapping a coefficient to a string, typically 'show' + -> (String, String) -- ^ used to enclose the coefficients, usually a pair of braces + -> (String, String) -- ^ pair of braces to enclose the numerator and the denominator + -> String -- ^ represents the quotient bar + -> RatioOfSprays a + -> String +showRatioOfSpraysXYZ letters showCoef coeffBraces = + showRatioOfSprays (showTwoSpraysXYZ showCoef coeffBraces letters) + +-- | Prints a ratio of sprays +showRatioOfSpraysXYZ' :: (Eq a, AlgField.C a) + => [String] -- ^ typically some letters, to represent the variables + -> (a -> String) -- ^ function mapping a coefficient to a string, typically 'show' + -> RatioOfSprays a + -> String +showRatioOfSpraysXYZ' letters showCoef = + showRatioOfSpraysXYZ letters showCoef ("(", ")") ("[ ", " ]") " %//% " + +-- | Prints a ratio of sprays +showRatioOfSpraysX1X2X3 :: forall a. (Eq a, AlgField.C a) + => String -- ^ typically a letter, to represent the variables + -> (a -> String) -- ^ function mapping a coefficient to a string, typically 'show' + -> (String, String) -- ^ used to enclose the coefficients, usually a pair of braces + -> (String, String) -- ^ pair of braces to enclose the numerator and the denominator + -> String -- ^ represents the quotient bar + -> RatioOfSprays a + -> String +showRatioOfSpraysX1X2X3 letter showCoef coeffBraces = + showRatioOfSprays (showTwoSpraysX1X2X3 showCoef coeffBraces letter) + +-- | Prints a ratio of sprays +showRatioOfSpraysX1X2X3' :: (Eq a, AlgField.C a) + => String -- ^ typically a letter, to represent the variables + -> (a -> String) -- ^ function mapping a coefficient to a string, typically 'show' + -> RatioOfSprays a + -> String +showRatioOfSpraysX1X2X3' letter showCoef = + showRatioOfSpraysX1X2X3 letter showCoef ("(", ")") ("[ ", " ]") " %//% " + +-- | Prints a ratio of sprays with rational coefficients +prettyRatioOfQSpraysXYZ :: + [String] -- ^ typically some letters, to represent the variables + -> RatioOfQSprays + -> String +prettyRatioOfQSpraysXYZ letters = + showRatioOfQSpraysXYZ letters ("[ ", " ]") " %//% " + +-- | Prints a ratio of sprays with rational coefficients +-- +-- prop> prettyRatioOfQSprays rOS == prettyRatioOfQSpraysXYZ ["x","y","z"] rOS +prettyRatioOfQSprays :: RatioOfQSprays -> String +prettyRatioOfQSprays = prettyRatioOfQSpraysXYZ ["x", "y", "z"] + +-- | Prints a ratio of sprays with rational coefficients +-- +-- prop> prettyRatioOfQSprays' rOS == prettyRatioOfQSpraysXYZ ["X","Y","Z"] rOS +prettyRatioOfQSprays' :: RatioOfQSprays -> String +prettyRatioOfQSprays' = prettyRatioOfQSpraysXYZ ["X", "Y", "Z"] + +-- | Prints a ratio of sprays with rational coefficients, printing the monomials +-- in the style of @\"x1^2.x2.x3^3\"@ +prettyRatioOfQSpraysX1X2X3 :: + String -- ^ typically a letter, to represent the non-indexed variables + -> RatioOfQSprays + -> String +prettyRatioOfQSpraysX1X2X3 letter = + showRatioOfQSpraysX1X2X3 letter ("[ ", " ]") " %//% " + +-- | Prints a ratio of sprays with numeric coefficients +prettyRatioOfNumSpraysXYZ :: (Num a, Ord a, AlgRing.C a, Show a) + => [String] -- ^ typically some letters, to represent the variables + -> RatioOfSprays a + -> String +prettyRatioOfNumSpraysXYZ letters = + showRatioOfNumSpraysXYZ show letters ("[ ", " ]") " %//% " + +-- | Prints a ratio of sprays with numeric coefficients +-- +-- prop> prettyRatioOfNumSprays rOS == prettyRatioOfNumSpraysXYZ ["x","y","z"] rOS +prettyRatioOfNumSprays :: + (Num a, Ord a, AlgRing.C a, Show a) => RatioOfSprays a -> String +prettyRatioOfNumSprays = prettyRatioOfNumSpraysXYZ ["x", "y", "z"] + +-- | Prints a ratio of sprays with numeric coefficients +-- +-- prop> prettyRatioOfNumSprays' rOS == prettyRatioOfNumSpraysXYZ ["X","Y","Z"] rOS +prettyRatioOfNumSprays' :: + (Num a, Ord a, AlgRing.C a, Show a) => RatioOfSprays a -> String +prettyRatioOfNumSprays' = prettyRatioOfNumSpraysXYZ ["X", "Y", "Z"] + +-- | Prints a ratio of sprays with numeric coefficients, printing the monomials +-- in the style of @\"x1^2.x2.x3^3\"@ +prettyRatioOfNumSpraysX1X2X3 :: (Num a, Ord a, AlgRing.C a, Show a) + => String -- ^ typically a letter, to represent the variables + -> RatioOfSprays a + -> String +prettyRatioOfNumSpraysX1X2X3 letter = + showRatioOfNumSpraysX1X2X3 show letter ("[ ", " ]") " %//% " + + +-- Parametric sprays ---------------------------------------------------------- + +type SimpleParametricSpray a = Spray (Spray a) +type SimpleParametricQSpray = SimpleParametricSpray Rational +type ParametricSpray a = Spray (RatioOfSprays a) +type ParametricQSpray = ParametricSpray Rational + +instance (Eq a, AlgRing.C a) => AlgMod.C a (SimpleParametricSpray a) where + (*>) :: a -> SimpleParametricSpray a -> SimpleParametricSpray a + lambda *> pspray = HM.map (lambda AlgMod.*>) pspray + +instance (Eq a, AlgRing.C a) => AlgRightMod.C a (SimpleParametricSpray a) where + (<*) :: SimpleParametricSpray a -> a -> SimpleParametricSpray a + pspray <* lambda = HM.map (AlgRightMod.<* lambda) pspray + +instance (Eq a, AlgField.C a) => AlgMod.C a (ParametricSpray a) where + (*>) :: a -> ParametricSpray a -> ParametricSpray a + lambda *> pspray = HM.map (lambda AlgMod.*>) pspray + +instance (Eq a, AlgField.C a) => AlgRightMod.C a (ParametricSpray a) where + (<*) :: ParametricSpray a -> a -> ParametricSpray a + pspray <* lambda = HM.map (AlgRightMod.<* lambda) pspray + +instance (Eq a, AlgField.C a) => AlgMod.C (Spray a) (ParametricSpray a) where + (*>) :: Spray a -> ParametricSpray a -> ParametricSpray a + spray *> pspray = asRatioOfSprays spray *^ pspray + +instance (Eq a, AlgField.C a) => AlgRightMod.C (Spray a) (ParametricSpray a) where + (<*) :: ParametricSpray a -> Spray a -> ParametricSpray a + pspray <* spray = asRatioOfSprays spray *^ pspray + +-- | Number of parameters in a parametric spray +-- +-- >>> numberOfParameters (jacobiPolynomial 4) +-- 2 +numberOfParameters :: FunctionLike b => Spray b -> Int +numberOfParameters pspray = + if isZeroSpray pspray + then 0 + else + 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 +-- to get \(P_{a^2, b^2}(X, Y, Z)\), or the one-parameter polynomial +-- \(P_{a, a}(X, Y, Z)\) +-- +-- >>> jp = jacobiPolynomial 4 +-- >>> a = qlone 1 +-- >>> b = qlone 2 +-- >>> changeParameters jp [a^**^2, b^**^2] +changeParameters :: + (FunctionLike b, Eq b, AlgAdd.C b) + => Spray b -- ^ @OneParameterSpray a@, @SimpleParametricSpray a@, or @ParametricSpray a@ + -> [VariablesType b] -- ^ @[Polynomial a]@ or @[Spray a]@, the new variables + -> Spray b +changeParameters pspray newParameters = + if length newParameters < numberOfParameters pspray + then + error "changeParameters: not enough new parameters provided." + else + removeZeroTerms $ HM.map (`changeVariables` newParameters) pspray + +-- | Substitutes some values to the parameters of a parametric spray +-- +-- >>> jacobi3 = jacobiPolynomial 3 +-- >>> legendre3 = substituteParameters jp [0, 0] +substituteParameters :: + (FunctionLike b, Eq (BaseRing b), AlgAdd.C (BaseRing b)) + => Spray b -- ^ @OneParameterSpray a@, @SimpleParametricSpray a@, or @ParametricSpray a@ + -> [BaseRing b] -- ^ values of type @a@ to be substituted to the parameters + -> Spray (BaseRing b) -- ^ output: a @Spray a@ spray +substituteParameters pspray values = + if length values < numberOfParameters pspray + then + error "substituteParameters: not enough values provided." + else + removeZeroTerms $ HM.map (evaluateAt values) pspray + +-- | helper function for evalParametricSpray +evalTerm' :: + (AlgMod.C (BaseRing b) b) => [BaseRing b] -> Term b -> b +evalTerm' xs (powers, coeff) = + AlgRing.product (zipWith (AlgRing.^) xs pows) AlgMod.*> coeff + where + pows = DF.toList (fromIntegral <$> exponents powers) + +-- | Substitutes some values to the variables of a parametric spray +evalParametricSpray :: + (Eq b, AlgMod.C (BaseRing b) b, AlgRing.C b) + => Spray b -- ^ @OneParameterSpray a@, @SimpleParametricSpray a@, or @ParametricSpray a@ + -> [BaseRing b] -- ^ values of type @a@ to be substituted to the variables + -> b +evalParametricSpray spray xs = if length xs >= numberOfVariables spray + then AlgAdd.sum $ map (evalTerm' xs) (HM.toList spray) + else error "evalParametricSpray: not enough values provided." + +-- | Substitutes some values to the parameters of a parametric spray as well as +-- some values to its variables +evalParametricSpray' :: + (FunctionLike b, Eq (BaseRing b), AlgMod.C (BaseRing b) b) + => Spray b -- ^ @OneParameterSpray a@, @SimpleParametricSpray a@, or @ParametricSpray a@ + -> [BaseRing b] -- ^ values of type @a@ to be substituted to the parameters + -> [BaseRing b] -- ^ values of type @a@ to be substituted to the variables + -> BaseRing b -- ^ result: a value of type @a@ +evalParametricSpray' spray as xs = + evaluateAt xs (substituteParameters spray as) + +-- | Whether the coefficients of a parametric spray polynomially +-- depend on their parameters; I do not know why, but it seems to be the case +-- for the Jacobi polynomials +-- +-- >>> canCoerceToSimpleParametricSpray (jacobiPolynomial 8) +-- True +canCoerceToSimpleParametricSpray :: + (Eq a, AlgRing.C a) => ParametricSpray a -> Bool +canCoerceToSimpleParametricSpray spray = + all isPolynomialRatioOfSprays (allCoefficients spray) + +-- | Coerces a parametric spray to a simple parametric spray, without +-- checking this makes sense with `canCoerceToSimpleParametricSpray` +asSimpleParametricSprayUnsafe :: ParametricSpray a -> SimpleParametricSpray a +asSimpleParametricSprayUnsafe = HM.map _numerator + +-- | Coerces a parametric spray to a simple parametric spray, after +-- checking this makes sense with `canCoerceToSimpleParametricSpray` +asSimpleParametricSpray :: + (Eq a, AlgRing.C a) => ParametricSpray a -> SimpleParametricSpray a +asSimpleParametricSpray spray = + if canCoerceToSimpleParametricSpray spray + then asSimpleParametricSprayUnsafe spray + else error $ + "asSimpleParametricSpray: this parametric spray is not coercible" ++ + " to a simple parametric spray." + +-- | Converts a `OneParameterSpray a` spray to a `ParametricSpray a` +fromOneParameterSpray :: + (Eq a, AlgRing.C a) => OneParameterSpray a -> ParametricSpray a +fromOneParameterSpray = HM.map fromRatioOfPolynomials + +-- | Converts a `OneParameterQSpray` spray to a `ParametricQSpray` +fromOneParameterQSpray :: OneParameterQSpray -> ParametricQSpray +fromOneParameterQSpray = HM.map fromRatioOfQPolynomials + +-- | Converts a `SimpleParametricSpray a` spray to a `ParametricSpray a` +fromSimpleParametricSpray :: + AlgRing.C a => SimpleParametricSpray a -> ParametricSpray a +fromSimpleParametricSpray = HM.map asRatioOfSprays + +-- | Converts a parametric spray to a one-parameter spray, without checking +-- the conversion makes sense +parametricSprayToOneParameterSpray :: + forall a. (Eq a, AlgField.C a) => ParametricSpray a -> OneParameterSpray a +parametricSprayToOneParameterSpray = HM.map toRatioOfPolynomials + where + toRatioOfPolynomials :: RatioOfSprays a -> RatioOfPolynomials a + toRatioOfPolynomials (RatioOfSprays p q) = + toPolynomial p % toPolynomial q + where + toPolynomial :: Spray a -> Polynomial a + toPolynomial spray = polyFromCoeffs coeffs + where + spray'' = removeConstantTerm'' (safeSpray spray) + coeffs = getConstantTerm spray : + [getCoefficient'' (S.singleton i) spray'' | i <- [1 .. deg]] + deg = maximum (0 : expnts) + expnts = map (`index` 0) (HM.keys spray'') + +{- -- | division of two univariate sprays +longDivision :: (Eq a, AlgField.C a) => Spray a -> Spray a -> (Spray a, Spray a) +longDivision sprayA sprayB = both fromCoeffs (polydiv coeffsA coeffsB) + where + fromCoeffs as = if S.null as + then zeroSpray + else sumTerms terms + where + l = S.length as + terms = (Powers S.empty 0, as `index` (l-1)) : + map (\i -> (Powers (S.singleton i) 1, as `index` (l-1-i))) [1 .. l-1] + shift n l = l >< S.replicate n AlgAdd.zero + pad n l = if n > 0 then S.replicate n AlgAdd.zero >< l else l + zipWith' op xs ys = S.zipWith op (pad (-d) xs) (pad d ys) + where d = S.length xs - S.length ys + coeffsA = coefficientsUnivariateSpray sprayA + coeffsB = coefficientsUnivariateSpray sprayB + coefficientsUnivariateSpray spray = coeffs + where + coeffs = S.fromList [getCoefficient' (Powers (S.singleton i) 1) spray' | + i <- [deg, deg-1 .. 1]] |> getConstantTerm spray + deg = maximum (0 : expnts) + spray' = removeConstantTerm spray + expnts = map ((`index` 0) . exponents) (HM.keys spray') + polydiv as bs = aux as bs S.empty + where aux f s q | ddif < 0 = (q, f) + | otherwise = aux f' s q' + where ddif = S.length (norm f) - S.length (norm s) + k = f `index` 0 AlgField./ s `index` 0 + ks = (AlgRing.* k) <$> shift ddif s + q' = zipWith' (AlgAdd.+) q $ shift ddif (S.singleton k) + f' = norm $ S.drop 1 $ zipWith' (AlgAdd.-) f ks + norm = S.dropWhileL (== AlgAdd.zero) + -} + +-- | Converts a rational parametric spray to a rational one-parameter spray, +-- without checking the conversion makes sense +parametricQSprayToOneParameterQSpray :: ParametricQSpray -> OneParameterQSpray +parametricQSprayToOneParameterQSpray = HM.map toRatioOfQPolynomials + where + toRatioOfQPolynomials :: RatioOfQSprays -> RatioOfQPolynomials + toRatioOfQPolynomials (RatioOfSprays p q) = + toQPolynomial p % toQPolynomial q + where + toQPolynomial :: QSpray -> QPolynomial + toQPolynomial spray = polyFromCoeffs coeffs' + where + spray'' = removeConstantTerm'' (safeSpray spray) + coeffs' = f (getConstantTerm spray) : + [f $ getCoefficient'' (S.singleton i) spray'' | i <- [1 .. deg]] + f :: Rational -> Rational' + f r = DR.numerator r :% DR.denominator r + deg = maximum (0 : map (`index` 0) expnts) + expnts = HM.keys spray'' + +-- | [Gegenbauer polynomials](https://en.wikipedia.org/wiki/Gegenbauer_polynomials); +-- we mainly provide them to give an example of the @SimpleParametricSpray@ type +-- +-- >>> gp = gegenbauerPolynomial 3 +-- >>> putStrLn $ prettySimpleParametricQSpray gp +-- { (4/3)*a^3 + 4*a^2 + (8/3)*a }*X^3 + { -2*a^2 - 2*a }*X +-- >>> putStrLn $ prettyQSpray'' $ substituteParameters gp [1] +-- 8*X^3 - 4*X +gegenbauerPolynomial :: Int -> SimpleParametricQSpray +gegenbauerPolynomial n + | n == 0 = unitSpray + | n == 1 = (2.^a) *^ x + | otherwise = + (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 + +-- | [Jacobi polynomial](https://en.wikipedia.org/wiki/Jacobi_polynomials); +-- the @n@-th Jacobi polynomial is a univariate polynomial of degree @n@ with +-- two parameters, except for the case @n=0@ where it has no parameter +-- +-- >>> jP = jacobiPolynomial 1 +-- >>> putStrLn $ prettyParametricQSprayABCXYZ ["alpha", "beta"] ["X"] jP +-- { [ (1/2)*alpha + (1/2)*beta + 1 ] }*X + { [ (1/2)*alpha - (1/2)*beta ] } +jacobiPolynomial :: Int -> ParametricQSpray +jacobiPolynomial n + | n < 0 = error "jacobiPolynomial: `n` must be positive." + | n == 0 = unitSpray + | n == 1 = + fromSimpleParametricSpray $ + (((gamma0 <+ 2) /^ 2) *^ + (x <+ AlgAdd.negate AlgRing.one)) <+ (alpha0 <+ 1) + | otherwise = + spray ^*^ jacobiPolynomial (n-1) ^-^ rOS *^ jacobiPolynomial (n-2) + where + alpha0 = qlone 1 + beta0 = qlone 2 + gamma0 = alpha0 ^+^ beta0 + x = lone 1 :: SimpleParametricQSpray + n' = toRational n + a0 = alpha0 <+ (n' - 1) + b0 = beta0 <+ (n' - 1) + c0 = gamma0 <+ (2 * n') + c0' = c0 <+ (-1) + c0'' = c0 <+ (-2) + divisor = (n' *^ (c0 <+ (-n'))) ^*^ c0'' + divisor' = 2 .^ divisor + divide = (`RatioOfSprays` divisor') + spray = HM.fromList [ + ( + nullPowers + , divide $ c0' ^*^ (alpha0 ^-^ beta0) ^*^ gamma0 + ), + ( + Powers (S.singleton 1) 1 + , divide $ c0' ^*^ c0 ^*^ c0'' + ) + ] + rOS = RatioOfSprays (a0 ^*^ b0 ^*^ c0) divisor + +-- | Pretty form of a numeric parametric spray, using some given strings (typically some +-- letters) to denote the parameters and some given strings (typically some letters) to +-- denote the variables; rather use `prettyParametricQSprayABCXYZ` for a rational +-- parametric spray +prettyParametricNumSprayABCXYZ :: + (Num a, Ord a, Show a, AlgField.C a) + => [String] -- ^ usually some letters, to denote the parameters of the spray + -> [String] -- ^ usually some letters, to denote the variables of the spray + -> ParametricSpray a -- ^ a parametric numeric spray + -> String +prettyParametricNumSprayABCXYZ abc xyz spray = + showSpray rOSShower ("{ ", " }") (showMonomialsXYZ xyz) spray + where + rOSShower = if numberOfParameters spray <= length abc + then prettyRatioOfNumSpraysXYZ abc + else prettyRatioOfNumSpraysX1X2X3 (abc !! 0) + +-- | Pretty form of a numeric parametric spray; rather use `prettyParametricQSpray` for +-- a rational parametric spray +-- +-- prop> prettyParametricNumSpray == prettyParametricNumSprayABCXYZ ["a"] ["X","Y","Z"] +prettyParametricNumSpray :: + (Num a, Ord a, Show a, AlgField.C a) + => ParametricSpray a -- ^ a parametric numeric spray + -> String +prettyParametricNumSpray = prettyParametricNumSprayABCXYZ ["a"] ["X", "Y", "Z"] + +-- | Pretty form of a parametric rational spray, using some given strings (typically some +-- letters) to denote the parameters and some given strings (typically some letters) to +-- denote the variables +-- +-- >>> type PQS = ParametricQSpray +-- >>> :{ +-- >>> f :: (QSpray, QSpray) -> (PQS, PQS, PQS) -> PQS +-- >>> f (a, b) (x, y, z) = +-- >>> (a %:% (a ^+^ b)) *^ x^**^2 ^+^ (b %:% (a ^+^ b)) *^ (y ^*^ z) +-- >>> :} +-- >>> a = qlone 1 +-- >>> b = qlone 2 +-- >>> x = lone 1 :: PQS +-- >>> y = lone 2 :: PQS +-- >>> z = lone 3 :: PQS +-- >>> pqs = f (a, b) (x, y, z) +-- >>> putStrLn $ prettyParametricQSprayABCXYZ ["a","b"] ["X","Y","Z"] pqs +-- { [ a ] %//% [ a + b ] }*X^2 + { [ b ] %//% [ a + b ] }*Y.Z +prettyParametricQSprayABCXYZ :: + [String] -- ^ usually some letters, to denote the parameters of the spray + -> [String] -- ^ usually some letters, to denote the variables of the spray + -> ParametricQSpray -- ^ a parametric rational spray + -> String +prettyParametricQSprayABCXYZ abc xyz spray = + showSpray rOSShower ("{ ", " }") (showMonomialsXYZ xyz) spray + where + rOSShower = if numberOfParameters spray <= length abc + then prettyRatioOfQSpraysXYZ abc + else prettyRatioOfQSpraysX1X2X3 (abc !! 0) + +-- | Pretty form of a parametric rational spray +-- +-- prop> prettyParametricQSpray == prettyParametricQSprayABCXYZ ["a"] ["X","Y","Z"] +prettyParametricQSpray :: ParametricQSpray -> String +prettyParametricQSpray = prettyParametricQSprayABCXYZ ["a"] ["X", "Y", "Z"] + +-- | Pretty form of a numeric simple parametric spray, using some given strings (typically some +-- letters) to denote the parameters and some given strings (typically some letters) to +-- denote the variables; rather use `prettySimpleParametricQSprayABCXYZ` for a rational +-- simple parametric spray +prettySimpleParametricNumSprayABCXYZ :: + (Num a, Ord a, Show a, AlgRing.C a) + => [String] -- ^ usually some letters, to denote the parameters of the spray + -> [String] -- ^ usually some letters, to denote the variables of the spray + -> SimpleParametricSpray a -- ^ a numeric simple parametric spray + -> String +prettySimpleParametricNumSprayABCXYZ abc xyz spray = + showSpray rOSShower ("{ ", " }") (showMonomialsXYZ xyz) spray + where + rOSShower = if numberOfParameters spray <= length abc + then prettyNumSprayXYZ abc + else prettyNumSprayX1X2X3 (abc !! 0) + +-- | Pretty form of a numeric simple parametric spray; rather use +-- `prettySimpleParametricQSpray` for a numeric simple parametric spray +-- +-- prop> prettySimpleParametricNumSpray == prettySimpleParametricNumSprayABCXYZ ["a"] ["X","Y","Z"] +prettySimpleParametricNumSpray :: + (Num a, Ord a, Show a, AlgRing.C a) + => SimpleParametricSpray a -- ^ a numeric simple parametric spray + -> String +prettySimpleParametricNumSpray = + prettySimpleParametricNumSprayABCXYZ ["a"] ["X", "Y", "Z"] + +-- | Pretty form of a simple parametric rational spray, using some given strings (typically some +-- letters) to denote the parameters and some given strings (typically some letters) to +-- denote the variables +-- +-- >>> type SPQS = SimpleParametricQSpray +-- >>> :{ +-- >>> f :: (QSpray, QSpray) -> (SPQS, SPQS, SPQS) -> SPQS +-- >>> f (a, b) (x, y, z) = +-- >>> (a ^+^ b) *^ x^**^2 ^+^ (a^**^2 ^+^ b^**^2) *^ (y ^*^ z) +-- >>> :} +-- >>> a = qlone 1 +-- >>> b = qlone 2 +-- >>> x = lone 1 :: SPQS +-- >>> y = lone 2 :: SPQS +-- >>> z = lone 3 :: SPQS +-- >>> spqs = f (a, b) (x, y, z) +-- >>> putStrLn $ prettySimpleParametricQSprayABCXYZ ["a","b"] ["X","Y","Z"] spqs +-- { a + b }*X^2 + { a^2 + b^2 }*Y.Z +prettySimpleParametricQSprayABCXYZ :: + [String] -- ^ usually some letters, to denote the parameters of the spray + -> [String] -- ^ usually some letters, to denote the variables of the spray + -> SimpleParametricQSpray -- ^ a parametric rational spray + -> String +prettySimpleParametricQSprayABCXYZ abc xyz spray = + showSpray sprayShower ("{ ", " }") (showMonomialsXYZ xyz) spray + where + sprayShower = if numberOfParameters spray <= length abc + then prettyQSprayXYZ abc + else prettyQSprayX1X2X3 (abc !! 0) + +-- | Pretty form of a simple parametric rational spray +-- +-- prop> prettySimpleParametricQSpray == prettySimpleParametricQSprayABCXYZ ["a"] ["X","Y","Z"] +prettySimpleParametricQSpray :: SimpleParametricQSpray -> String +prettySimpleParametricQSpray = + prettySimpleParametricQSprayABCXYZ ["a"] ["X", "Y", "Z"] +
+ tests/Approx.hs view
@@ -0,0 +1,10 @@+module Approx (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)
tests/Main.hs view
@@ -1,73 +1,1173 @@-module Main where -import Data.Ratio ( (%) ) -import Math.Algebra.Hspray ( Spray, - (^+^), - (^*^), - (^**^), - (*^), - lone, - unitSpray, - evalSpray, - composeSpray, - fromList, - toList, - bombieriSpray, - derivSpray ) -import Test.Tasty ( defaultMain - , testGroup - ) -import Test.Tasty.HUnit ( assertEqual - , testCase - ) - -main :: IO () -main = defaultMain $ testGroup - "Testing hspray" - - [ testCase "bombieriSpray" $ do - let - x = lone 1 :: Spray Rational - y = lone 2 :: Spray Rational - z = lone 3 :: Spray Rational - poly = - (2 % 1) *^ ((2 % 1) *^ (x ^**^ 3 ^*^ y ^**^ 2)) ^+^ (4 % 1) *^ z ^+^ (5 % 1) *^ unitSpray - bpoly = - (24 % 1) *^ ((2 % 1) *^ (x ^**^ 3 ^*^ y ^**^ 2)) ^+^ (4 % 1) *^ z ^+^ (5 % 1) *^ unitSpray - assertEqual "" bpoly (bombieriSpray poly), - - testCase "composeSpray" $ do - let - x = lone 1 :: Spray Int - y = lone 2 :: Spray Int - z = lone 3 :: Spray Int - p = 2 *^ (2 *^ (x ^**^ 3 ^*^ y ^**^ 2)) ^+^ 4 *^ z ^+^ 5 *^ unitSpray - px = x ^+^ y ^+^ z - py = x ^*^ y ^*^ z - pz = y ^**^ 2 - q = composeSpray p [px, py, pz] - xyz = [2, 3, 4] - pxyz = map (`evalSpray` xyz) [px, py, pz] - assertEqual "" (evalSpray p pxyz) (evalSpray q xyz), - - testCase "fromList . toList = identity" $ do - let - x = lone 1 :: Spray Int - y = lone 2 :: Spray Int - z = lone 3 :: Spray Int - p = 2 *^ (2 *^ (x ^**^ 3 ^*^ y ^**^ 2)) ^+^ 4 *^ z ^+^ 5 *^ unitSpray - assertEqual "" p (fromList . toList $ p), - - testCase "derivSpray" $ do - let - x = lone 1 :: Spray Int - y = lone 2 :: Spray Int - z = lone 3 :: Spray Int - p1 = x ^+^ y ^*^ z ^**^ 3 - p2 = (x ^*^ y ^*^ z) ^+^ (2 *^ (x ^**^ 3 ^*^ y ^**^ 2)) - q = p1 ^*^ p2 - p1' = derivSpray 1 p1 - p2' = derivSpray 1 p2 - q' = derivSpray 1 q - assertEqual "" q' ((p1' ^*^ p2) ^+^ (p1 ^*^ p2')) +{-# LANGUAGE TupleSections #-} +module Main (main) where +import qualified Algebra.Additive as AlgAdd +import qualified Algebra.Module as AlgMod +import qualified Algebra.Ring as AlgRing +import qualified Algebra.Field as AlgField +import Approx ( approx, assertApproxEqual ) +import qualified Data.HashMap.Strict as HM +import Data.List ( sortOn ) +import Data.Matrix ( Matrix, fromLists ) +import Data.Maybe ( fromJust ) +import Data.Ratio ( (%) ) +import Math.Algebra.Hspray ( Spray, + QSpray, + FunctionLike (..), + (.^), + (/>), + (/^), + lone, + qlone, + unitSpray, + zeroSpray, + constantSpray, + getCoefficient, + getConstantTerm, + evalSpray, + substituteSpray, + composeSpray, + permuteVariables, + swapVariables, + fromList, + toList, + bombieriSpray, + collinearSprays, + derivative, + groebnerBasis, + fromRationalSpray, + esPolynomial, + psPolynomial, + isSymmetricSpray, + isPolynomialOf, + resultant, + resultant', + subresultants, + resultant1, + subresultants1, + polynomialSubresultants, + sturmHabichtSequence, + principalSturmHabichtSequence, + numberOfRealRootsInOpenInterval, + numberOfRealRootsInClosedInterval, + numberOfRealRoots, + sprayDivision, + gcdSpray, + QSpray', + Rational', + Polynomial, + OneParameterQSpray, + evalRatioOfPolynomials, + evalOneParameterSpray', + constQPoly, + evalOneParameterSpray'', + prettyQSpray, + prettyQSprayX1X2X3, + prettySpray, + prettySpray'', + qsoleParameter, + soleParameter, + constQPoly, + constPoly, + prettyOneParameterQSpray', + RatioOfSprays (..), + RatioOfQSprays, + (%//%), + (%/%), + (%:%), + unitRatioOfSprays, + constantRatioOfSprays, + isPolynomialRatioOfSprays, + evalRatioOfSprays, + evalRatioOfSprays', + substituteRatioOfSprays, + prettyRatioOfQSprays, + characteristicPolynomial, + detLaplace', + gegenbauerPolynomial, + jacobiPolynomial, + asRatioOfSprays, + SimpleParametricQSpray, + ParametricQSpray, + ParametricSpray, + zeroRatioOfSprays, + fromRatioOfQPolynomials, + numberOfParameters, + changeParameters, + substituteParameters, + evalParametricSpray, + asSimpleParametricSpray, + parametricSprayToOneParameterSpray, + prettyParametricQSprayABCXYZ, + asSimpleParametricSpray, + constPoly, + qlone', + qmonomial, + isHomogeneousSpray, + psPolynomial, + prettyQSprayXYZ, + productOfSprays + ) +import MathObj.Matrix ( fromRows ) +import qualified MathObj.Matrix as MathMatrix +import Number.Ratio ( T ( (:%) ) ) +import qualified Number.Ratio as NR +import Test.Tasty ( defaultMain + , testGroup + ) +import Test.Tasty.HUnit ( assertEqual + , assertBool + , testCase + ) + +type PQS = ParametricQSpray + +main :: IO () +main = defaultMain $ testGroup + + "Testing hspray" + + [ + 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 + homogeneous = isHomogeneousSpray spray + assertEqual "" homogeneous (True, Just $ sum lambda) + + , testCase "qmonomial and qlone'" $ do + let + assocs = [(1, 5), (3, 7), (1, 5)] + assocs0 = [(0, 4), (0, 0)] + assertEqual "" + (qmonomial assocs, qmonomial assocs0) + (qlone' 1 5 ^*^ qlone' 3 7, unitSpray) + + , testCase "asSimpleParametricSpray" $ do + let + jp = jacobiPolynomial 8 + jp' = asSimpleParametricSpray jp + jp'' = HM.map asRatioOfSprays jp' + assertEqual "" jp jp'' + + , testCase "prettyParametricQSprayABCXYZ" $ do + let + f :: (QSpray, QSpray) -> (PQS, PQS, PQS) -> PQS + f (a, b) (x, y, z) = + (a %:% (a ^+^ unitSpray)) *^ x^**^2 ^+^ + (b %:% (a ^+^ b)) *^ (y ^*^ z) + pqs = f (lone 1, lone 2) (lone 1, lone 2, lone 3) + s1 = prettyParametricQSprayABCXYZ ["a","b"] ["X","Y","Z"] pqs + s2 = prettyParametricQSprayABCXYZ ["a"] ["X","Y","Z"] pqs + s3 = prettyParametricQSprayABCXYZ ["a","b"] ["X","Y"] pqs + assertEqual "" + [ + s1, s2, s3 + ] + [ + "{ [ a ] %//% [ a + 1 ] }*X^2 + { [ b ] %//% [ a + b ] }*Y.Z", + "{ [ a1 ] %//% [ a1 + 1 ] }*X^2 + { [ a2 ] %//% [ a1 + a2 ] }*Y.Z", + "{ [ a ] %//% [ a + 1 ] }*X1^2 + { [ b ] %//% [ a + b ] }*X2.X3" + ] + + , testCase "substituteParameters and evalParametricSpray" $ do + let + jacobi = jacobiPolynomial 5 + jacobi' = substituteParameters jacobi [3, 2%7] + r1 = evaluate jacobi' [13] + r2 = evaluate (evalParametricSpray jacobi [13]) [3, 2%7] + assertEqual "" r1 r2 + + , testCase "changeParameters in Jacobi polynomial -> Gegenbauer" $ do + let + risingFactorial :: QSpray -> Int -> QSpray + risingFactorial theta n = + AlgRing.product + (map (\k -> theta ^+^ constantSpray (toRational k)) [0 .. n-1]) + m = 5 + alpha = qlone 1 + jacobi = + changeParameters (jacobiPolynomial m) + [alpha ^-^ constantSpray (1%2), alpha ^-^ constantSpray (1%2)] + factor = + risingFactorial (2 *^ alpha) m %//% + risingFactorial (alpha ^+^ constantSpray (1%2)) m + obtained = asSimpleParametricSpray (factor *^ jacobi) + gegenbauer = gegenbauerPolynomial m + assertEqual "" gegenbauer obtained + + , testCase "substituteParameters in Jacobi polynomial -> Legendre" $ do + let + x = qlone 1 + jacobi = jacobiPolynomial 5 + legendre = (63*^x^**^5 ^-^ 70*^x^**^3 ^+^ 15*^x) /^ 8 + assertEqual "" legendre (substituteParameters jacobi [0, 0]) + + , testCase "changeParameters in Jacobi polynomial" $ do + let + n = 5 + jp_n = jacobiPolynomial n + jp_nminus1 = jacobiPolynomial (n-1) + jp_n' = derivative 1 jp_n + a = qlone 1 + b = qlone 2 + spray = (a ^+^ b ^+^ constantSpray (toRational $ n + 1)) /^ 2 + rOS = asRatioOfSprays spray + rhs = rOS *^ + changeParameters jp_nminus1 [a ^+^ unitSpray, b ^+^ unitSpray] + assertEqual "" jp_n' rhs + + , testCase "numberOfParameters in Jacobi polynomial" $ do + let + jp = jacobiPolynomial 5 + assertEqual "" (numberOfParameters jp) 2 + + , testCase "changeVariables in RatioOfSprays" $ do + let + f :: QSpray -> QSpray -> RatioOfQSprays + f p1 p2 = (p1^**^2 ^+^ 2 *^ p2) %//% (p1^**^3 ^-^ unitSpray) + x = qlone 1 + y = qlone 2 + rOS = f x y + u = x ^*^ y + v = x^**^2 ^-^ y ^+^ unitSpray + rOS' = f u v + assertEqual "" rOS' (changeVariables rOS [u, v]) + + , testCase "changeVariables in constant RatioOfSprays" $ do + let + rOS = constantRatioOfSprays 3 :: RatioOfQSprays + x = qlone 1 + y = qlone 2 + u = x ^*^ y + v = x^**^2 ^-^ y ^+^ unitSpray + assertEqual "" rOS (changeVariables rOS [u, v]) + + , testCase "changeVariables in Polynomial" $ do + let + f :: Polynomial Rational -> Polynomial Rational + f p = (constPoly (2::Rational) AlgRing.* p) AlgRing.^ 3 + AlgAdd.- p AlgAdd.+ constPoly 4 + a = soleParameter + pol = f a + u = a AlgRing.^ 2 AlgAdd.- constPoly 10 + assertEqual "" (f u) (changeVariables pol [u]) + + , testCase "fromRatioOfQPolynomials" $ do + let + a = qsoleParameter + x = qlone 1 + rOP = ((a AlgRing.^ 8 AlgAdd.- AlgRing.one) NR.% + (a AlgAdd.- AlgRing.one)) AlgRing.^ 3 + AlgAdd.+ (a AlgAdd.+ AlgRing.one) :% a + rOQ = ((x^**^8 ^-^ unitSpray) %//% (x ^-^ unitSpray)) AlgRing.^ 3 + AlgAdd.+ RatioOfSprays (x ^+^ unitSpray) x + assertEqual "" rOQ (fromRatioOfQPolynomials rOP) + + , testCase "(.^)" $ do + let + x = lone 1 :: QSpray + y = lone 2 :: QSpray + rOS = (x^**^4 ^-^ y^**^4) %//% (x ^+^ y ^-^ unitSpray) + assertEqual "" (10 .^ rOS AlgAdd.+ (-10) .^ rOS) zeroRatioOfSprays + + , testCase "module `ParametricSpray a` over `a`" $ do + let + x = lone 1 :: ParametricQSpray + y = lone 2 :: ParametricQSpray + p = x^**^2 ^+^ x^*^y ^-^ unitSpray + lambda = 3 :: Rational + p' = asRatioOfSprays (lambda *^ unitSpray) *^ p + assertEqual "" (lambda AlgMod.*> p) p' + + , testCase "Jacobi polynomial" $ do + let + jp2 = jacobiPolynomial 2 + alpha0 = qlone 1 + beta0 = qlone 2 + cst :: Rational -> QSpray + cst = constantSpray + x = lone 1 :: ParametricQSpray + p = x ^-^ unitSpray + t1 = asRatioOfSprays (((alpha0 ^+^ cst 1)^*^(alpha0 ^+^ cst 2)) /^ 2) + t2 = asRatioOfSprays + (((alpha0 ^+^ cst 2)^*^(alpha0 ^+^ beta0 ^+^ cst 3)) /^ 2) + t3 = asRatioOfSprays + (((alpha0 ^+^ beta0 ^+^ cst 3)^*^(alpha0 ^+^ beta0 ^+^ cst 4)) /^ 8) + expected = t1 *^ unitSpray ^+^ t2 *^ p ^+^ t3 *^ p^**^2 + assertEqual "" jp2 expected + + , testCase "characteristic polynomial" $ do + let + m = fromLists [ [12, 16, 4] + , [16, 2, 8] + , [8, 18, 10] ] :: Matrix Int + spray = characteristicPolynomial m + x = lone 1 :: Spray Int + expected = + AlgAdd.negate x^**^3 ^+^ 24*^x^**^2 ^+^ 268*^x ^-^ constantSpray 1936 + assertEqual "" spray expected + + , testCase "determinant of product" $ do + let + m = fromRows 3 3 [ [12, 16, 4] + , [16, 2, 8] + , [8, 18, 10] ] :: MathMatrix.T Int + assertEqual "" (detLaplace' (m AlgRing.* m)) (AlgRing.sqr (detLaplace' m)) + + , testCase "ratio of sprays is irreducible" $ do + let + x = lone 1 :: QSpray + y = lone 2 :: QSpray + w = lone 4 :: QSpray + rOS = (x^**^4 ^-^ y^**^4) %//% (x ^-^ y) + rOS' = rOS %/% (x^**^4 ^-^ y^**^4) + rOS'' = w AlgMod.*> rOS' + assertEqual "" + ( + prettyRatioOfQSprays rOS + , prettyRatioOfQSprays rOS' + , prettyRatioOfQSprays rOS'' + ) + ( + "[ x^3 + x^2.y + x.y^2 + y^3 ]" + , "[ 1 ] %//% [ x - y ]" + , "[ x4 ] %//% [ x1 - x2 ]" + ) + + , testCase "isPolynomialRatioOfSprays" $ do + let + x = qlone 1 + y = qlone 2 + rOS = (x^**^4 ^-^ y^**^4) %//% (x ^-^ y) + assertBool "" (isPolynomialRatioOfSprays rOS) + + , testCase "power of ratio of sprays" $ do + let + x = qlone 1 + y = qlone 2 + z = qlone 3 + p = x^**^4 ^-^ x ^*^ y^**^4 ^+^ x ^*^ z + q = x ^-^ x ^*^ y + rOS = p %//% q + rOS' = p^**^4 %//% q^**^4 + assertEqual "" rOS' (rOS AlgRing.^ 4) + + , testCase "arithmetic on ratio of sprays" $ do + let + x = qlone 1 + y = qlone 2 + p = x^**^2 ^-^ 3*^(x ^*^ y) ^+^ y^**^3 + q = x ^-^ y + rOS1 = p^**^2 %//% q + rOS2 = rOS1 AlgAdd.+ unitRatioOfSprays + rOS = rOS1 AlgRing.^ 2 AlgAdd.+ rOS1 AlgRing.* rOS2 AlgAdd.- rOS1 + test1 = + (rOS1 AlgAdd.+ rOS2) AlgRing.* (rOS1 AlgAdd.- rOS2) == + rOS1 AlgRing.^ 2 AlgAdd.- rOS2 AlgRing.^ 2 + rOS' = + (3%4 :: Rational) AlgMod.*> rOS AlgRing.^ 2 AlgAdd.+ p AlgMod.*> rOS + test2 = p AlgMod.*> (rOS' %/% p) == rOS' + test3 = rOS1 %/% p == p %//% q + test4 = rOS' AlgField./ rOS' == unitRatioOfSprays + k = 3 :: Rational + test5 = (p /> k) AlgMod.*> rOS1 == p AlgMod.*> (rOS1 /> k) + assertEqual "" + [test1, test2, test3, test4, test5] [True, True, True, True, True] + + , testCase "arithmetic on univariate ratio of sprays" $ do + let + x = qlone 1 + y = qlone 2 + f :: QSpray -> RatioOfQSprays + f p = + ((p^**^4 ^-^ 3*^p^**^3 ^+^ p^**^2)^**^2 %//% (p^**^2 ^+^ unitSpray)) + AlgAdd.+ + ((p^**^2 ^-^ unitSpray) %//% (p^**^2 ^+^ 3*^p)^**^2) AlgRing.* + ((p^**^3 ^-^ p^**^2 ^+^ unitSpray) %//% (2*^p^**^3 ^-^ unitSpray)) + rOSx = f x + rOSy = f y + assertEqual "" rOSx (swapVariables (1, 2) rOSy) + + , testCase "evaluate ratio of sprays" $ do + let + x = qlone 1 + y = qlone 2 + p = x ^+^ y + q = x ^-^ y + rOS1 = p %//% q + rOS2 = q %//% p + f :: AlgField.C a => a -> a -> a + f u v = u AlgRing.^ 2 AlgAdd.+ u AlgRing.* v AlgAdd.- u AlgField./ v + rOS = f rOS1 rOS2 + values = [2%3, 7%4] + r1 = evalRatioOfSprays rOS1 values + r2 = evalRatioOfSprays rOS2 values + r = evalRatioOfSprays rOS values + assertEqual "" r (f r1 r2) + + , testCase "evalRatioOfSprays'" $ do + let + x = qlone 1 + y = qlone 2 + z = qlone 3 + p = x ^+^ y + q = x ^-^ y + rOS = p %//% q + rOS1 = x %//% y + rOS2 = z %//% x + obtained = evalRatioOfSprays' rOS [rOS1, rOS2] + expected = + (rOS1 AlgAdd.+ rOS2) AlgField./ (rOS1 AlgAdd.- rOS2) + assertEqual "" obtained expected + + , testCase "Gegenbauer: differential equation and Chebyshev case" $ do + let + n = 5 + g = gegenbauerPolynomial n + 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 + chebyshev = fromRationalSpray $ substituteParameters g [1] + theta = 2.5 + assertEqual "" + (shouldBeZero, approx 8 $ sin theta * evalSpray chebyshev [cos theta]) + (zeroSpray, approx 8 $ sin (fromIntegral (n+1) * theta)) + + , testCase "scale spray by integer" $ do + let + x = lone 1 :: Spray Int + y = lone 2 :: Spray Int + z = lone 3 :: Spray Int + p = 2 *^ (2 *^ (x ^**^ 3 ^*^ y ^**^ 2)) ^+^ 4 *^ z ^+^ 5 *^ unitSpray + spray = p ^+^ p ^+^ p ^+^ p ^+^ p ^+^ p ^+^ p + assertEqual "" + (7 .^ p, (-7) .^ p, 0 .^ p) (spray, AlgAdd.negate spray, zeroSpray) + + , testCase "collinearSprays" $ do + let + x = lone 1 :: Spray Rational + y = lone 2 :: Spray Rational + z = lone 3 :: Spray Rational + term1 = (2 % 1) *^ ((2 % 1) *^ (x ^**^ 3 ^*^ y ^**^ 2)) + term2 = (4 % 1) *^ z + term3 = (5 % 1) *^ unitSpray + spray1 = term1 ^+^ term2 ^+^ term3 + spray2 = 121 *^ (term3 ^+^ term1 ^+^ term2) + assertBool "" (collinearSprays spray1 spray2) + + , testCase "bombieriSpray" $ do + let + x = lone 1 :: Spray Rational + y = lone 2 :: Spray Rational + z = lone 3 :: Spray Rational + poly = + (2 % 1) *^ ((2 % 1) *^ (x ^**^ 3 ^*^ y ^**^ 2)) ^+^ + (4 % 1) *^ z ^+^ (5 % 1) *^ unitSpray + bpoly = + (24 % 1) *^ ((2 % 1) *^ (x ^**^ 3 ^*^ y ^**^ 2)) ^+^ + (4 % 1) *^ z ^+^ (5 % 1) *^ unitSpray + assertEqual "" bpoly (bombieriSpray poly), + + testCase "composeSpray" $ do + let + x = lone 1 :: Spray Int + y = lone 2 :: Spray Int + z = lone 3 :: Spray Int + p = 2 *^ (2 *^ (x ^**^ 3 ^*^ y ^**^ 2)) ^+^ 4 *^ z ^+^ 5 *^ unitSpray + px = x ^+^ y ^+^ z + py = x ^*^ y ^*^ z + pz = y ^**^ 2 + q = composeSpray p [px, py, pz] + xyz = [2, 3, 4] + pxyz = map (`evalSpray` xyz) [px, py, pz] + assertEqual "" (evalSpray p pxyz) (evalSpray q xyz), + + testCase "getCoefficient" $ do + let + x = lone 1 :: Spray Int + y = lone 2 :: Spray Int + z = lone 3 :: Spray Int + p = 2 *^ (2 *^ (x^**^3 ^*^ y^**^2)) ^+^ 4 *^ z ^+^ 5 *^ unitSpray + assertEqual "" + (getCoefficient [3, 2, 0] p, getCoefficient [0, 4] p) (4, 0), + + testCase "getConstantTerm" $ do + let + x = lone 1 :: Spray Int + y = lone 2 :: Spray Int + z = lone 3 :: Spray Int + p = 2 *^ (2 *^ (x^**^3 ^*^ y^**^2)) ^+^ 4 *^ z ^+^ 5 *^ unitSpray + assertEqual "" (getConstantTerm p) 5 + + , testCase "fromList . toList = identity" $ do + let + x = lone 1 :: Spray Int + y = lone 2 :: Spray Int + z = lone 3 :: Spray Int + p = 2 *^ (2 *^ (x ^**^ 3 ^*^ y ^**^ 2)) ^+^ 4 *^ z ^+^ 5 *^ unitSpray + assertEqual "" p (fromList . toList $ p) + + , testCase "toList . fromList cleans the list" $ do + let + l = [([], 2), ([3, 1], 4), ([3, 1], 6), ([2, 0, 0], 5), ([9, 3], 0)] + spray = fromList l :: Spray Int + l' = sortOn fst (toList spray) + assertEqual "" l' [([], 2), ([2], 5), ([3, 1], 10)] + + , testCase "derivative of spray" $ do + let + x = lone 1 :: Spray Int + y = lone 2 :: Spray Int + z = lone 3 :: Spray Int + p1 = x ^+^ y ^*^ z ^**^ 3 + p2 = (x ^*^ y ^*^ z) ^+^ (2 *^ (x ^**^ 3 ^*^ y ^**^ 2)) + q = p1 ^*^ p2 + p1' = derivative 1 p1 + p2' = derivative 1 p2 + q' = derivative 1 q + assertEqual "" q' ((p1' ^*^ p2) ^+^ (p1 ^*^ p2')) + + , testCase "derivative of a ratio of sprays" $ do + let + x = lone 1 :: QSpray + y = lone 2 :: QSpray + rOS = x %//% y + rOS' = derivative 2 rOS + expected = AlgAdd.negate x %//% y^**^2 + assertEqual "" rOS' expected + + , 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 = groebnerBasis [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 polynomial" $ do + let + e2 = esPolynomial 4 2 :: Spray Rational + e3 = esPolynomial 4 3 :: Spray Rational + p = e2^**^2 ^+^ (2*^ e3) + assertBool "" (isSymmetricSpray p), + + testCase "Schur polynomial is symmetric" $ do + let + x = lone 1 :: Spray Rational + y = lone 2 :: Spray Rational + z = lone 3 :: Spray Rational + p = x^**^3 ^*^ y^**^2 ^*^ z ^+^ x^**^3 ^*^ y ^*^ z^**^2 ^+^ + x^**^2 ^*^ y^**^3 ^*^ z ^+^ 2*^(x^**^2 ^*^ y^**^2 ^*^ z^**^2) ^+^ + 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 + y = lone 2 :: Spray Rational + p1 = x ^+^ y + p2 = x ^-^ y + p = p1 ^*^ p2 ^+^ unitSpray + assertEqual "" + (isPolynomialOf p [p1, p2]) (True, Just $ x ^*^ y ^+^ unitSpray), + + testCase "isPolynomialOf - 2" $ do + let + x = lone 1 :: Spray Rational + y = lone 2 :: Spray Rational + z = lone 3 :: Spray Rational + assertEqual "" + (isPolynomialOf x [x ^+^ y^*^z, y, z]) + (True, Just $ x ^-^ y^*^z) + + , 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 + symSpray = x^**^2 ^+^ y^**^2 ^+^ x ^+^ y + p1 = psPolynomial 2 1 :: Spray Rational + p2 = psPolynomial 2 2 :: Spray Rational + p = fromJust $ snd $ isPolynomialOf symSpray [p1, p2] + symSpray' = composeSpray p [p1, p2] + assertEqual "" symSpray symSpray' + + , 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)) + + , testCase "substituteRatioOfSprays" $ do + let + x1 = lone 1 :: Spray Rational + x2 = lone 2 :: Spray Rational + x3 = lone 3 :: Spray Rational + p = x1^**^2 ^+^ x2 ^+^ x3 ^-^ unitSpray + q = x1 ^-^ x2 + rOS = p %//% q + subs = [Just 2, Nothing, Just 3] + p' = substituteSpray subs p + q' = substituteSpray subs q + rOS' = substituteRatioOfSprays subs rOS + assertEqual "" rOS' (p' %//% q') + + , testCase "permuteVariables of a spray" $ do + let + f :: Spray Rational -> Spray Rational -> Spray Rational + -> Spray Rational + f p1 p2 p3 = + p1^**^4 ^+^ (2 *^ p2^**^3) ^+^ (3 *^ p3^**^2) ^-^ (4 *^ unitSpray) + x1 = lone 1 :: Spray Rational + x2 = lone 2 :: Spray Rational + x3 = lone 3 :: Spray Rational + p = f x1 x2 x3 + p' = permuteVariables [3, 1, 2] p + assertEqual "" p' (f x3 x1 x2) + + , testCase "permuteVariables can add a new variable" $ do + let + f :: Spray Rational -> Spray Rational -> Spray Rational + -> Spray Rational + f p1 p2 p3 = + p1^**^4 ^+^ (2 *^ p2^**^3) ^+^ (3 *^ p3^**^2) ^-^ (4 *^ unitSpray) + x1 = lone 1 :: Spray Rational + x2 = lone 2 :: Spray Rational + x3 = lone 3 :: Spray Rational + x4 = lone 4 :: Spray Rational + p = f x1 x2 x3 + p' = permuteVariables [2, 1, 4, 3] p + assertEqual "" p' (f x2 x1 x4) + + , testCase "swapVariables of a spray" $ do + let + x1 = lone 1 :: Spray Rational + x2 = lone 2 :: Spray Rational + x3 = lone 3 :: Spray Rational + p = x1^**^4 ^+^ (2 *^ x2^**^3) ^+^ (3 *^ x3^**^2) ^-^ (4 *^ unitSpray) + p' = permuteVariables [3, 2, 1] p + assertEqual "" p' (swapVariables (1, 3) p) + + , testCase "permuteVariables of a ratio of sprays" $ do + let + f :: QSpray -> QSpray -> QSpray -> RatioOfQSprays + f p1 p2 p3 = (p1^**^4 ^+^ (2 *^ p2^**^3)) %//% + ((3 *^ p3^**^2) ^-^ (4 *^ unitSpray)) + x1 = lone 1 :: QSpray + x2 = lone 2 :: QSpray + x3 = lone 3 :: QSpray + rOS = f x1 x2 x3 + rOS' = permuteVariables [3, 1, 2] rOS + assertEqual "" rOS' (f x3 x1 x2) + + , testCase "swapVariables of a ratio of sprays" $ do + let + x1 = lone 1 :: Spray Rational + x2 = lone 2 :: Spray Rational + x3 = lone 3 :: Spray Rational + rOS = (x1^**^4 ^+^ (2 *^ x2^**^3)) %//% + ((3 *^ x3^**^2) ^-^ (4 *^ unitSpray)) + rOS' = permuteVariables [3, 2, 1] rOS + assertEqual "" rOS' (swapVariables (1, 3) rOS) + + , testCase "resultant w.r.t x" $ do + let + x = lone 1 :: Spray Rational + y = lone 2 :: Spray Rational + p = x^**^4 ^-^ x^**^3 ^+^ x^**^2 ^-^ 2*^ (x ^*^ y^**^2) ^+^ y^**^4 + q = x ^-^ (2*^ y^**^2) + r = resultant 1 p q + assertEqual "" r (y^**^4 ^-^ (8*^ y^**^6) ^+^ (16*^ y^**^8)), + + testCase "resultant w.r.t y" $ do + let + x = lone 1 :: Spray Rational + y = lone 2 :: Spray Rational + p = x^**^4 ^-^ x^**^3 ^+^ x^**^2 ^-^ 2*^ (x ^*^ y^**^2) ^+^ y^**^4 + q = x ^-^ (2*^ y^**^2) + r = resultant 2 p q + assertEqual "" + r (16*^x^**^8 ^-^ 32*^x^**^7 ^+^ 24*^x^**^6 ^-^ 8*^x^**^5 ^+^ x^**^4), + + testCase "resultant product rule" $ do + let + x = lone 1 :: Spray Rational + y = lone 2 :: Spray Rational + f = x^**^4 ^-^ x^**^3 ^+^ x^**^2 ^-^ 2*^ (x ^*^ y^**^2) ^+^ y^**^4 + g = x ^-^ (2*^ y^**^2) + h = x^**^2 ^*^ y ^+^ y^**^3 ^+^ unitSpray + assertEqual "" + (resultant 1 (f^*^g) h) (resultant 1 f h ^*^ resultant 1 g h), + + testCase "subresultants" $ do + let + x = lone 1 :: Spray Rational + y = lone 2 :: Spray Rational + p = x^**^2 ^*^ y ^*^ (y^**^2 ^-^ 5*^ x ^+^ constantSpray 6) + q = x^**^2 ^*^ y ^*^ (3*^ y ^+^ constantSpray 2) + sx = subresultants 1 p q + assertBool "" + (sx!!0 == zeroSpray && sx!!1 == zeroSpray && sx!!2 /= zeroSpray), + + testCase "resultant1" $ do + let + x = lone 1 :: Spray Rational + p = x^**^2 ^-^ 5*^x ^+^ constantSpray 6 + q = x^**^2 ^-^ 3*^x ^+^ constantSpray 2 + assertEqual "" (resultant1 p q) (0%1), + + testCase "resultant1 product rule" $ do + let + x = lone 1 :: Spray Rational + f = x^**^2 ^-^ 5*^x ^+^ constantSpray 6 + g = x^**^2 ^-^ 3*^x ^+^ constantSpray 2 + h = x^**^3 ^+^ x ^-^ constantSpray 3 + assertEqual "" (resultant1 (f^*^g) h) (resultant1 f h * resultant1 g h), + + testCase "subresultants1" $ do + let + x = lone 1 :: Spray Rational + p = x^**^2 ^-^ 5*^x ^+^ constantSpray 6 + q = x^**^2 ^-^ 3*^x ^+^ constantSpray 2 + assertEqual "" (subresultants1 p q) [0%1, 2%1, 1%1], + + testCase "resultant agrees with resultant1 for univariate case" $ do + let + x = lone 1 :: Spray Rational + f = x^**^4 ^-^ x^**^3 ^+^ x^**^2 ^-^ 2*^x + g = x ^-^ (2*^ x^**^2) ^+^ constantSpray 4 + r = resultant 1 f g + r1 = resultant1 f g + assertEqual "" r1 (getCoefficient [] r), + + testCase "resultant and resultant' are in agreement" $ do + let + x = lone 1 :: Spray Rational + y = lone 2 :: Spray Rational + z = lone 3 :: Spray Rational + p = x^**^4 ^-^ x^**^3 ^+^ x^**^2 ^-^ 2*^ (x ^*^ y^**^2) ^+^ z^**^4 + q = x ^-^ (2*^ y^**^2) ^*^ z^**^2 ^*^ unitSpray + rx = resultant 1 p q + rx' = resultant' 1 p q + ry = resultant 2 p q + ry' = resultant' 2 p q + rz = resultant 3 p q + rz' = resultant' 3 p q + test1 = rx == rx' + test2 = ry == ry' + test3 = rz == rz' + assertBool "" (test1 && test2 && test3) + + , testCase "polynomialSubresultants" $ do + let + x = qlone 1 + y = qlone 2 + p = x^**^3^*^y^**^3 ^+^ x^**^2 ^+^ 3*^x ^+^ unitSpray + q = x^**^4 ^+^ 2*^x^**^3 ^-^ x^*^y^**^4 + srs1 = map prettyQSpray (polynomialSubresultants 1 p q) + srs2 = map prettyQSpray (polynomialSubresultants 2 p q) + assertEqual "" + (srs1, srs2) + ( + [ + "y^21 + 15*y^14 - 11*y^11 + y^8 + 37*y^7 + 4*y^4 - 8*y^3 - 1" + , "x.y^17 + 14*x.y^10 - 11*x.y^7 + x.y^4 + 39*x.y^3 + 5*x + 2*y^10 - y^7 + 14*y^3 + 2" + , "-5*x^2.y^3 + x^2 - x.y^10 - 7*x.y^3 + 3*x - 2*y^3 + 1" + , "x^3.y^3 + x^2 + 3*x + 1" + ], + [ + "x^24 + 6*x^23 + 12*x^22 + 8*x^21 - x^11 - 12*x^10 - 58*x^9 - 144*x^8 - 195*x^7 - 144*x^6 - 58*x^5 - 12*x^4 - x^3" + , "x^13 + 5*x^12 + 7*x^11 + 2*x^10 + x^9.y + 6*x^8.y + 11*x^7.y + 6*x^6.y + x^5.y" + , "x^10 + 2*x^9 + x^6.y + 3*x^5.y + x^4.y" + , "x^3.y^3 + x^2 + 3*x + 1" + ] + ) + + , testCase "principalSturmHabichtSequence" $ do + let + x = qlone 1 + y = qlone 2 + p = x^**^4 ^*^ y ^+^ y^**^2 + psh1 = principalSturmHabichtSequence 1 p + psh2 = principalSturmHabichtSequence 2 p + assertEqual "" + (psh1, psh2) + ( + [ + 256 *^ y^**^10 + , zeroSpray + , zeroSpray + , 4 *^ y + , y + ], + [ + x^**^8 + , constantSpray 2 + , constantSpray 1 + ] + ) + + , testCase "sturmHabichtSequence & principalSturmHabichtSequence" $ do + let + x = qlone 1 + y = qlone 2 + p = x^**^3^*^y^**^4 ^+^ x^**^2^*^y^**^2 ^+^ 3*^y ^-^ constantSpray 6 + sh1 = map prettyQSpray (sturmHabichtSequence 1 p) + sh2 = map prettyQSpray (sturmHabichtSequence 2 p) + psh1 = map prettyQSpray (principalSturmHabichtSequence 1 p) + psh2 = map prettyQSpray (principalSturmHabichtSequence 2 p) + assertEqual "" + ((sh1, psh1), (sh2, psh2)) + ( + ( + [ + "-243*y^14 + 972*y^13 - 972*y^12 - 12*y^11 + 24*y^10" + , "2*x.y^8 - 27*y^9 + 54*y^8" + , "3*x^2.y^4 + 2*x.y^2" + , "x^3.y^4 + x^2.y^2 + 3*y - 6" + ], + [ + "-243*y^14 + 972*y^13 - 972*y^12 - 12*y^11 + 24*y^10" + , "2*y^8" + , "3*y^4" + , "y^4" + ] + ), + ( + [ + "-96*x^14 - 4608*x^13 - 55332*x^12 - 7776*x^11 - 2187*x^9" + , "-8*x^12.y - 192*x^11.y - 12*x^10 - 324*x^9.y + 864*x^9" + , "-8*x^8.y^2 - 36*x^6.y + 96*x^6" + , "4*x^3.y^3 + 2*x^2.y + 3" + , "x^3.y^4 + x^2.y^2 + 3*y - 6" + ], + [ + "-96*x^14 - 4608*x^13 - 55332*x^12 - 7776*x^11 - 2187*x^9" + , "-8*x^12 - 192*x^11 - 324*x^9" + , "-8*x^8" + , "4*x^3" + , "x^3" + ] + ) + ) + + , testCase "number of real roots in interval" $ do + let + x = qlone 1 + factors = [x ^-^ constantSpray (toRational i) | i <- [1::Int .. 5]] + spray = productOfSprays factors + intervals = [ + (Just 0, Just 9) + , (Just 1, Just 6) + , (Just 2, Just 3) + , (Just 0, Just 4) + , (Just $ 2 + (1%4), Just $ 3 - (1%4)) + ] + nroots = + map (uncurry (numberOfRealRootsInClosedInterval spray)) intervals + nroots' = + map (uncurry (numberOfRealRootsInOpenInterval spray)) intervals + rightUnboundedIntervals = map ((, Nothing) . Just) [0, 1, 2, 2 + (1%4), 5] + nroots'' = map (uncurry (numberOfRealRootsInClosedInterval spray)) + rightUnboundedIntervals + nroots''' = map (uncurry (numberOfRealRootsInOpenInterval spray)) + rightUnboundedIntervals + leftUnboundedIntervals = map ((Nothing, ) . Just) [0, 1, 2, 2 + (1%4), 5] + nroots'''' = map (uncurry (numberOfRealRootsInClosedInterval spray)) + leftUnboundedIntervals + nroots''''' = map (uncurry (numberOfRealRootsInOpenInterval spray)) + leftUnboundedIntervals + ntotalroots = numberOfRealRootsInClosedInterval spray Nothing Nothing + assertEqual "" + ( nroots, nroots' + , nroots'', nroots''' + , nroots'''', nroots''''' + , ntotalroots + ) + ( [5, 5, 2, 4, 0], [5, 4, 0, 3, 0] + , [5, 5, 4, 3, 1], [5, 4, 3, 3, 0] + , [0, 1, 2, 2, 5], [0, 0, 1, 2, 4] + , 5 + ) + + , testCase "number of real roots" $ do + let + x = qlone 1 + factors = [x ^-^ constantSpray (toRational i) | i <- [1::Int .. 5]] + spray = productOfSprays factors + nroots = map numberOfRealRoots [spray, x^**^3, x^**^2 <+ 1] + assertEqual "" nroots [5, 1, 0] + + , 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 = gcdSpray sprayA sprayB + assertEqual "" sprayGCD sprayD, + + testCase "gcdSpray with a constant spray" $ do + let + x = lone 1 :: Spray Rational + sprayA = 3 *^ x^**^4 ^-^ x + b1 = 2 :: Rational + b2 = 4 :: Rational + sprayB1 = constantSpray b1 + sprayB2 = constantSpray b2 + assertBool "" + (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 "sprayDivision by constant spray" $ do + let + x = lone 1 :: Spray Rational + y = lone 2 :: Spray Rational + sprayA = 3*^(x^**^2 ^*^ y) ^-^ 3*^(x ^*^ y) ^+^ constantSpray 3 + sprayB = constantSpray 3 + expected = x^**^2 ^*^ y ^-^ x ^*^ y ^+^ unitSpray + assertEqual "" (sprayDivision sprayA sprayB) (expected, 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 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 + + , testCase "gcdSpray is commutative" $ do + let + x = lone 1 :: Spray Rational + y = lone 2 :: Spray Rational + z = lone 3 :: Spray Rational + sprayA = (-5)*^x^**^3 ^-^ 11*^(x^**^2 ^*^ y) ^-^ 11*^(x^*^ y^**^2) ^-^ + 8*^(x ^*^ y ^*^ z) ^+^ 20*^(x ^*^ z^**^2) ^-^ 5 *^ y^**^3 ^+^ + 20*^(y ^*^ z^**^2) + sprayB = x^*^y ^-^ x^*^z ^-^ y^*^z ^+^ z^**^2 + assertEqual "" (gcdSpray sprayA sprayB) (gcdSpray sprayB sprayA) + + , testCase "evaluation of one-parameter spray" $ do + let + a = qsoleParameter + p = a AlgRing.^ 2 AlgAdd.- constQPoly 4 + q1 = a AlgAdd.- constQPoly 3 + q2 = a AlgAdd.- constQPoly 2 + rop1 = p :% q1 + rop2 = p :% q2 + f :: (Eq a, AlgRing.C a) + => Spray a -> Spray a -> Spray a -> (Spray a, Spray a) + f x y z = (x^**^2 ^+^ y^**^2, z) + g :: (Eq a, AlgRing.C a) + => Spray a -> Spray a -> Spray a -> (a, a, a) -> (a, a) + g px py pz (x, y, z) = + ( + evalSpray f1 [x, y, AlgAdd.zero] + , evalSpray f2 [AlgAdd.zero, AlgAdd.zero, z] + ) + where (f1, f2) = f px py pz + (r1, r2) = g (lone 1 :: QSpray') (lone 2) (lone 3) (2, 3, 4) + r = evalRatioOfPolynomials 5 rop1 AlgRing.* r1 AlgAdd.+ + evalRatioOfPolynomials 5 rop2 AlgRing.* r2 + (f1', f2') = f (lone 1 :: OneParameterQSpray) (lone 2) (lone 3) + opSpray = rop1 *^ f1' ^+^ rop2 *^ f2' + r' = evalOneParameterSpray' opSpray 5 [2, 3, 4] + rop1' = evalOneParameterSpray'' f1' [2, 3] + rop2' = evalOneParameterSpray'' f2' [0, 0, 4] + r'' = evalRatioOfPolynomials 5 + (rop1 AlgRing.* rop1' AlgAdd.+ rop2 AlgRing.* rop2') + spray = substituteParameters opSpray [5] + r''' = evalSpray spray [2, 3, 4] + assertEqual "" ((r, r'), r'') ((r', r''), r'''), + + testCase "pretty spray" $ do + let + x = lone 1 :: Spray Rational + y = lone 2 :: Spray Rational + z = lone 3 :: Spray Rational + p1 = ((2%3) *^ x^**^3) ^*^ y ^-^ x^**^2 ^+^ + y ^*^ z ^-^ (2%3) *^ unitSpray + p2 = (3%2) *^ p1 + p3 = AlgAdd.negate $ + swapVariables (1, 3) $ + p2 ^+^ unitSpray ^-^ (x^**^3 ^*^ y ^-^ ((3%2) *^ x^**^2)) + strings = + [ + prettyQSpray (zeroSpray ^*^ p1) + , prettyQSpray p1 + , prettyQSpray (AlgAdd.negate p2) + , prettyQSpray (p2 ^+^ lone 4) + , prettyQSpray p3 + , " ---------- " + , prettyQSprayX1X2X3 "a" (zeroSpray ^*^ p1) + , prettyQSprayX1X2X3 "a" p1 + , prettyQSprayX1X2X3 "a" (p2 ^+^ lone 4) + , prettyQSprayX1X2X3 "a" p3 + , " ---------- " + , prettySpray (zeroSpray ^*^ p1) + , prettySpray p1 + , prettySpray (p2 ^+^ lone 4) + , prettySpray p3 + , " ---------- " + , prettySpray'' "w" (zeroSpray ^*^ p1) + , prettySpray'' "w" p1 + , prettySpray'' "w" (p2 ^+^ lone 4) + , prettySpray'' "w" p3 + ] + strings' = + [ + "0" + , "(2/3)*x^3.y - x^2 + y.z - (2/3)" + , "-x^3.y + (3/2)*x^2 - (3/2)*y.z + 1" + , "x1^3.x2 - (3/2)*x1^2 + (3/2)*x2.x3 + x4 - 1" + , "-(3/2)*x.y" + , " ---------- " + , "0" + , "(2/3)*a1^3.a2 - a1^2 + a2.a3 - (2/3)" + , "a1^3.a2 - (3/2)*a1^2 + (3/2)*a2.a3 + a4 - 1" + , "-(3/2)*a1.a2" + , " ---------- " + , "0" + , "(2 % 3)*x^3.y + ((-1) % 1)*x^2 + (1 % 1)*y.z + ((-2) % 3)" + , "(1 % 1)*x1^3.x2 + ((-3) % 2)*x1^2 + (3 % 2)*x2.x3 + (1 % 1)*x4 + ((-1) % 1)" + , "((-3) % 2)*x.y" + , " ---------- " + , "0" + , "(2 % 3)*w^(3, 1) + ((-1) % 1)*w^(2) + (1 % 1)*w^(0, 1, 1) + ((-2) % 3)*w^()" + , "(1 % 1)*w^(3, 1) + ((-3) % 2)*w^(2) + (3 % 2)*w^(0, 1, 1) + (1 % 1)*w^(0, 0, 0, 1) + ((-1) % 1)*w^()" + , "((-3) % 2)*w^(1, 1)" + ] + assertEqual "" strings strings', + + testCase "prettyOneParameterQSpray'" $ do + let + x = lone 1 :: OneParameterQSpray + y = lone 2 :: OneParameterQSpray + z = lone 3 :: OneParameterQSpray + a = qsoleParameter + sSpray + = ((4 NR.% 5) *^ (a :% (a AlgRing.^ 2 AlgAdd.+ AlgRing.one))) + *^ (x^**^2 ^-^ y^**^2) + ^+^ (constQPoly (2 NR.% 3) AlgRing.* a) AlgMod.*> (y ^*^ z) + string = prettyOneParameterQSpray' "a" sSpray + string' = + "{ [ (4/5)*a ] %//% [ a^2 + 1 ] }*X^2 + { [ -(4/5)*a ] %//% [ a^2 + 1 ] }*Y^2 + { (2/3)*a }*Y.Z" + assertEqual "" string string' + + , testCase "parametricSprayToOneParameterSpray" $ do + let + x = lone 1 :: OneParameterQSpray + y = lone 2 :: OneParameterQSpray + z = lone 3 :: OneParameterQSpray + a = qsoleParameter + sSpray + = ((4 NR.% 5) *^ (a :% (a AlgRing.^ 2 AlgAdd.+ AlgRing.one))) + *^ (x^**^2 ^-^ y^**^2) + ^+^ (constQPoly (2 NR.% 3) AlgRing.* a) AlgMod.*> (y ^*^ z) + x' = lone 1 :: ParametricSpray Rational' + y' = lone 2 :: ParametricSpray Rational' + z' = lone 3 :: ParametricSpray Rational' + a' = lone 1 :: Spray Rational' + spray + = ((4 NR.% 5 :: Rational') AlgMod.*> + (a' %:% (a'^**^2 ^+^ unitSpray))) *^ (x'^**^2 ^-^ y'^**^2) + ^+^ ((2 NR.% 3) *^ a') AlgMod.*> (y' ^*^ z') + assertEqual "" (parametricSprayToOneParameterSpray spray) sSpray + ]