There is something new in the qspray package: Gröbner bases. Gröbner bases have many applications, and we will see one of them now. My package qspray deals with multivariate polynomials with rational coefficients.

Consider the following system of three polynomial equations:

\[\begin{align} x^2 + y + z & = 1 \\ x + y^2 + z & = 1 \\ x + y + z^2 & = 1 \end{align}\]

The goal is to find all solutions of this system, that is to say, the simultaneous roots of the three polynomials \(x^2 + y + z – 1\), \(x + y^2 + z – 1\), \(x + y + z^2 – 1\). Not easy, is it?

Fact: a Gröbner basis of (the ideal generated by) a list of polynomials is a list of polynomials with the same simultaneous roots.

library(qspray)
x <- qlone(1); y <- qlone(2); z <- qlone(3)
f1 <- x^2 + y + z - 1
f2 <- x + y^2 + z - 1
f3 <- x + y + z^2 - 1
gb <- groebner(list(f1, f2, f3)) # Gröbner basis of f1, f2, f3
( gbstrings <- lapply(gb, prettyQspray, vars = c("x", "y", "z")) )
## [[1]]
## [1] "y + x + z^2 - 1"
## 
## [[2]]
## [1] "((-9)*z^3)/2 - 2*z^5 - z^6 + 6*z^4 - y + z + z^7/2 + y^2"
## 
## [[3]]
## [1] "y*z^2 - z^2/2 + z^4/2"
## 
## [[4]]
## [1] "z^6 - z^2 + 4*z^3 - 4*z^4"

Do you see something nice? The fourth polynomial, \(z^6 - z^2 + 4z^3 - 4z^4\), is univariate: it depends only on the variable \(z\). We could numerically solve it with the R function polyroot, but we want exact roots. It is easy to see that \(0\) and \(1\) are two roots of this polynomial. So we can divide this polynomial by \(z(z-1)\), and it will remain a factor of degree \(4\). But I know, you are lazy. Me too. That’s why I prefer to use the Ryacas package to perform this task:

library(Ryacas)
yac_str("Factor(z^6 - z^2 + 4*z^3 - 4*z^4)")
## [1] "(2*z+z^2-1)*(z-1)^2*z^2"

So it remains to find the roots of \(2z + z^2 - 1\). We use Ryacas again, as we’re lazy:

yac_str("PSolve(2*z + z^2 - 1, z)")
## [1] "{(Sqrt(8)-2)/2,(-(Sqrt(8)+2))/2}"

There are two roots: \(\sqrt{2} - 1\) and \(-\sqrt{2} - 1\).

But let’s consider the two friendly roots first: \(0\) and \(1\). If you plug \(z=0\) and \(z=1\) in the three other polynomials of the Gröbner basis, you easily find three simultaneous roots: \((1, 0, 0)\), \((0, 1, 0)\) and \((0, 0, 1)\).

You should do this exercise, because now we will plug \(z=\sqrt{2}-1\) and \(z=-\sqrt{2}-1\), and this will be less easy. In fact I will only treat \(z=\sqrt{2}-1\), and leave the other root as another exercise for you.

Let’s go. We plug \(z=\sqrt{2}-1\) in the third polynomial. Clearly, we will obtain an univariate polynomial in \(y\), and we will solve it. Or rather, Ryacas will solve it:

yac_str(sprintf("p3 := %s", gbstrings[[3L]]))
## [1] "y*z^2-z^2/2+z^4/2"
yac_str("PSolve(Eliminate(z, Sqrt(2)-1, p3), y)")
## [1] "-((Sqrt(2)-1)^4/2-(Sqrt(2)-1)^2/2)/(Sqrt(2)-1)^2"

Yacas is not always nice. It does not simplify this expression. So let’s look at its numerical approximation:

yac_str("N(PSolve(Eliminate(z, Sqrt(2)-1, p3), y))")
## [1] "0.4142135623"

We recognize \(\sqrt{2}-1\). Now we plug \(z=\sqrt{2}-1\) into the second polynomial and we solve the resulting polynomial of degree \(2\) in \(y\):

yac_str(sprintf("p2 := %s", gbstrings[[2L]]))
## [1] "((-9)*z^3)/2-2*z^5-z^6+6*z^4-y+z+z^7/2+y^2"
yac_str("N(PSolve(Eliminate(z, Sqrt(2)-1, p2), y))")
## [1] "{0.5857864376,0.4142135623}"

There are two roots: \(y=\sqrt{2}-1\) and the other one has no importance, since we didn’t find it before.

So finally we just have to plug \(y=\sqrt{2}-1\) and \(z=\sqrt{2}-1\) into the first polynomial and we solve the resulting univariate polynomial in \(x\):

yac_str(sprintf("p1 := %s", gbstrings[[1L]]))
## [1] "y+x+z^2-1"
yac_str(
  "N(PSolve(Eliminate(y, Sqrt(2)-1, Eliminate(z, Sqrt(2)-1, p1)), x))"
)
## [1] "0.4142135623"

Again? Well, we find \(x = \sqrt{2}-1\). Finally we found one simultaneous root: \[(x, y, z) = (\sqrt{2}-1, \sqrt{2}-1, \sqrt{2}-1).\] Now I leave you with the promised exercise: treat the case \(z=-\sqrt{2}-1\). I give you the solution, you will find: \[(x, y, z) = (-\sqrt{2}-1, -\sqrt{2}-1, -\sqrt{2}-1).\]

Isn’t it impressive? I started to learn about Gröbner bases only a few days ago, and I still have to learn. I hope some applications will be implemented in qspray in the future, and/or exposed on this blog.

The new version of qspray is not on CRAN yet. If you are in a hurry, you can install the development version hosted on Github.