R package to deal with multivariate polynomials withrational coefficients.
This package is strongly inspired by Robin Hankin’sspray package. The C++ implementations are verysimilar.
library(qspray)Theqspray package provides theqsprayobjects, which represent multivariate polynomials whose coefficients arerational numbers.
qspray polynomial and arithmeticThe easiest way to build a multivariate polynomial withqspray is to start by introducing the generatingvariables with the help of theqlone function and then tocombine them with arithmetic operations:
x<-qlone(1); y<-qlone(2); z<-qlone(3)( pol<-4*x^2+"1/2"*y-5*x*y*z/3 )## 4*x^2 - 5/3*x.y.z + 1/2*yI often like to use a function like this:
f<-function(x, y, z) {4*x^2+ y/2-5*x*y*z/3}f(x, y, z)## 4*x^2 - 5/3*x.y.z + 1/2*yOr maybe you prefer to define the polynomial by giving it as astring:
qsprayMaker(string ="4 x^(2) + 1/2 x^(0, 1) - 5/3 x^(1, 1, 1)")## 4*x^2 - 5/3*x.y.z + 1/2*yAs you want, but this method is not highly robust. And it is not veryeasy to figure out what is the monomial represented by a string such as"x^(i,j,k)" (this isx^i*y^j*z^k).
Some arithmetic on this polynomial:
-pol## -4*x^2 + 5/3*x.y.z - 1/2*y2* pol## 8*x^2 - 10/3*x.y.z + ypol/2## 2*x^2 - 5/6*x.y.z + 1/4*y"5/3"* pol## 20/3*x^2 - 25/9*x.y.z + 5/6*ypol+5## 4*x^2 - 5/3*x.y.z + 1/2*y + 5pol- gmp::as.bigq("2/5")## 4*x^2 - 5/3*x.y.z + 1/2*y - 2/5pol^2## 16*x^4 - 40/3*x^3.y.z + 25/9*x^2.y^2.z^2 + 4*x^2.y - 5/3*x.y^2.z + 1/4*y^2Two polynomials can be added and multiplied:
pol1<- polpol2<- polpol1+ pol2## 8*x^2 - 10/3*x.y.z + ypol1- pol2## 0pol1* pol2## 16*x^4 - 40/3*x^3.y.z + 25/9*x^2.y^2.z^2 + 4*x^2.y - 5/3*x.y^2.z + 1/4*y^2qsprayUseevalQspray to evaluate a polynomial for some valuesof the variables:
evalQspray(pol,c("1","2","3/2"))## Big Rational ('bigq') :## [1] 0Alternatively, you can convert the polynomial to a function:
g<-as.function(pol)g("1","2","3/2")## [1] "0"You can pass the strings you want as the arguments of thisfunction:
g("x","y","z")## [1] "(24*x^2-10*x*z*y+3*y)/6"g("x+1","2*x","y^2")## [1] "((-10)*x^2*y^2+12*x^2+(-10)*x*y^2+27*x+12)/3"The output ofg("x+1", "2*x", "y^2") is the expressionof a bivariate polynomial. You can get it as aqspraypolynomial with the help of the functionchangeVariables(see sectionTransforming aqspray).
If you want a function returning numerical approximations, use theoptionN=TRUE:
h<-as.function(pol,N =TRUE)h("1","2","3/2")## [1] 2e-29h("x","y","z")## expression(4 * x^2 - 1.6666666666 * x * y * z + 0.5 * y)h("x+1","2*x","y^2")## expression(-3.3333333333 * x^2 * y^2 + 4 * x^2 - 3.3333333333 *## x * y^2 + 9 * x + 4)You can also perform “partial evaluation” of aqspray,that is to say replacing only certain variables. This is done by usingthe functionsubstituteQspray and indicating the variablesto be kept withNA:
substituteQspray(pol,c("1",NA,"3/2"))## -2*y + 4f(gmp::as.bigq(1), y, gmp::as.bigq("3/2"))## -2*y + 4g("1","y","3/2")## [1] "2*(2-y)"h("1","y","3/2")## expression(-2 * y + 4)qsprayYou can control the way of printing aqspray with thehelp of the functionshowQsprayOption<-. By default, themonomials of aqspray are printed in the style ofx^2.y.z^3 if there are no more than three variables,otherwise they are printed in the style ofx1^2.x2.x3^3:
set.seed(3141)( qspray<-rQspray() )# a random qspray## -2*x^4.y^3.z^4 - 4*y^2.z^2qspray+qlone(4)^99## -2*x1^4.x2^3.x3^4 - 4*x2^2.x3^2 + x4^99If you want to always use the second way, you can do:
showQsprayOption(qspray,"x")<-"x"qspray## -2*x1^4.x2^3.x3^4 - 4*x2^2.x3^2If you want to restore the wayqspray objects wereprinted in previous versions, you can do
showQsprayOption(qspray,"showMonomial")<-showMonomialOld()qspray## -2*x^(4, 3, 4) - 4*x^(0, 2, 2)There are three possible show options that can be passed toshowQsprayOption:
"showQspray". AshowQspray function, that is to say a function appropriatefor the"showQspray" option, must be a function whichtransforms aqspray to a string. The package provides somehelper functions to built such functions, likeshowQsprayXYZ andshowQsprayX1X2X3. WithshowQsprayXYZ, you can choose the letters you want todenote the variables:f<-showQsprayXYZ(c("A","B","C"))f(qspray)## [1] "-2*A^4.B^3.C^4 - 4*B^2.C^2"WithshowQsprayX1X2X3, you choose only one letter forthe variables and they will be appended with a digit:
f<-showQsprayX1X2X3("X")f(qspray)## [1] "-2*X1^4.X2^3.X3^4 - 4*X2^2.X3^2"Once you have constructed such a function, you pass it as a showoption by doingshowQsprayOption(qspray, "showQspray") <- f.
"showMonomial", tocontrol the way the monomials are printed. Actually in the two aboveexamples ofshowQsprayXYZ andshowQsprayX1X2X3we only changed the way the monomials are printed. Indeed, these twocommands are equivalent:showQsprayOption(qspray,"showQspray")<-showQsprayXYZ(c("A","B","C"))showQsprayOption(qspray,"showMonomial")<-showMonomialXYZ(c("A","B","C"))and these two commands are equivalent as well:
showQsprayOption(qspray,"showQspray")<-showQsprayX1X2X3("X")showQsprayOption(qspray,"showMonomial")<-showMonomialX1X2X3("X")But theshowQspray functions allow finer control,e.g. they allow to control the multiplication symbol which separates acoefficient and a monomial within a term.
"x". Setting thisoption to a letterx:showQsprayOption(qspray,"x")<- xis equivalent to:
showQsprayOption(qspray,"showMonomial")<-showMonomialX1X2X3(x)ButshowMonomialX1X2X3 also allows to control the waythe individual powers are collapsed, e.g. "x^2.y.z^3" (thedefault) or"x^2*y*z^3", or"x^2yz^3". If thedot is nice for you, use the"x" option, that’s less codeto type.
By the way, aqspray object is an S4 object with twoslots:powers andcoeffs. Thepowers slot is a list of vector of exponents and thecoeffs slot is a character vector, whose each element iscoercable to abigq number by an application of thefunctiongmp::as.bigq. TheshowMonomialfunctions act only on thepowers slot.
When an arithmetic operation is performed between twoqspray objects, the show options of the first one arepassed to the result,if possible:
qspray+qlone(4)^99## -2*X1^4.X2^3.X3^4 - 4*X2^2.X3^2 + X4^99For example, this is not possible if you specify only three lettersfor the variables and you perform an operation with aqspray involving the fourth variable:
showQsprayOption(qspray,"showMonomial")<-showMonomialXYZ(c("a","b","c"))qspray## -2*a^4.b^3.c^4 - 4*b^2.c^2qspray+qlone(4)^99## -2*a1^4.a2^3.a3^4 - 4*a2^2.a3^2 + a4^99The package provides a function which returns the exact value of theintegral of a polynomial with rational coefficients over a simplex whosevertices have rational Cartesian coordinates:
# variablesx<-qlone(1); y<-qlone(2); z<-qlone(3)# polynomialP<- x^4+ y+2*x*y^2-3*z# simplex (tetrahedron) verticesv1<-c(1,1,1)v2<-c(2,2,3)v3<-c(3,4,5)v4<-c(3,2,1)# simplexS<-rbind(v1, v2, v3, v4)# integralintegratePolynomialOnSimplex(P, S)## Big Rational ('bigq') :## [1] 1387/42qsprayLet’s take aqspray polynomial:
f<-function(x, y, z) {4*x^2+ y/2-5*x*y*z/3}x<-qlone(1); y<-qlone(2); z<-qlone(3)P<-f(x, y, z)You can get a derivative of this polynomial:
derivQspray(P,i =2)# derivative w.r.t y## -5/3*x.z + 1/2You can permute the variables of this polynomial:
swapVariables(P,1,3)==f(z, y, x)## [1] TRUEYou can perform a change of variables on this polynomial:
changeVariables(P,list(x+1,2*x, y^2))==f(x+1,2*x, y^2)## [1] TRUEFinally, let us mention thegroebner function, whichcomputes a Gröbner basis of the ideal generated by a list ofqspray polynomials:
f<-qsprayMaker(string ="x^(3) - 2 x^(1,1)")g<-qsprayMaker(string ="x^(2,1) - 2 x^(0,2) + x^(1)")groebner(list(f, g))## [[1]]## x - 2*y^2#### [[2]]## y^3As an application of Gröbner bases, there is the functionisPolynomialOf. This function checks whether a polynomialcan be obtained by substituting the variables of a polynomial with somegiven polynomials: given aqspray polynomialQand someqspray polynomialsP1, …,Pn, does there exist a polynomial functionfsuch thatQ = f(P1, ..., Pn)? If this is true, theisPolynomialOf function also returnsf.
There are packages depending on theqspray package(some of them are not on CRAN yet):
polyhedralCubature:this package uses theintegratePolynomialOnSimplex functionto get the exact integral of a multivariate polynomial over apolytope.
jack: Jackpolynomials.
resultant:resultant, subresultants, and greatest common divisor of twoqspray polynomials.
ratioOfQsprays:fractions ofqspray polynomials.
symbolicQspray:multivariate polynomials whose coefficients areratioOfQsprays fractions of polynomials; they representmultivariate polynomials with parameters.
The three packagesjack,ratioOfQsprays andsymbolicQspray usesome C++ code based on the header file of the C++ code ofqspray. If you want to use it in your package too,include the following instruction in the DESCRIPTION file:
LinkingTo: Rcpp, RcppArmadillo, qsprayAnd include the following instruction in your C++ code:
#include"qspray.h"Then you can use theqspray header file in your C++code by using the namespaceQSPRAY.
The header file provides the templated classQspray. Anobject of typeQspray<T> represents a multivariatepolynomials whose coefficients are represented by the objects of thetypeT. For example, multivariate polynomials with numericcoefficients can be represented by the objects of typeQspray<double> (so I should have chosen another namesinceQ is here to indicate the field of rational numbers).The classQspray provides the operators==,!=,+,-,* and thepower for the objects of typeQspray<T> as long asthe operators==,!=,+,- and* are available for the typeT. So you don’t have to implement the comparison operatorsnor the arithmetic operations for theQspray<double>polynomials if you instantiate this type. The classQsprayalso provides a function to calculate derivatives. This class isincluded in the namespaceQSPRAY which also includes afunction performing the division of two multivariate polynomials but itis restricted to polynomials with rational coefficients. Anyway it wouldnot be a good idea to use the algorithm performed by this function forpolynomials whose coefficients type is not an “exact type”, such asdouble. If you want to useQspray<T>with an exact typeT and if you need the division, send mea few words about your use case and I will see whether I can help. Iwill probably remove the division from the namespaceQSPRAY. I originally included it to use it in theratioOfQsprays package, but I finally used the divisionprovided by theCGAL library instead, which isfaster.
A few words about the implementation. The classQspray<T> has only one member object: an object oftypePolynomial<T>, which is an alias of the typestd::unordered_map<std::vector<int>, T> (plus atemplate argument for the hasher). So aPolynomial<T>object is a map whose keys arestd::vector<int>objects and whose values areT objects. An element of thismap represents a term of the polynomial: a key represents a monomial,e.g. the vector{2,1,3} represents the monomialx^2*y*z^3, and the value attached to this key representsthe coefficient of this monomial. This way to represent a multivariatepolynomial has been copied from Robin Hankin’sspraypackage, without which theqspray package would havenever existed.