Fast expansion of a polynomial with R – part 2

[This article was first published on Saturn Elephant, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

In the previous post, I showed how to expand a polynomial with symbolic parameters with the help of the spray package. As I said, it has one problem: it doesn’t preserve the rational numbers in the polynomial expression.

I’m going to provide a solution here which overcomes this problem, with the help of the Ryacas package. In fact I provide a solution with the packages Ryacas and partitions, and then I improve it with the help of the spray package.

Here is the polynomial expression I use for the illustration:

expr <- "((x*x+y*y+1)*(a*x*x+b*y*y)+z*z*(b*x*x+a*y*y)-2*(a-b)*x*y*z-a*b*(x*x+y*y))^2-4*(x*x+y*y)*(a*x*x+b*y*y-x*y*z*(a-b))^2"

In fact, the equation expr(x,y,z) = 0 defines a solid Möbius strip. That is why I was interested in this expression, because I wanted to draw the solid Möbius strip with POV-Ray. It is nice, in spite of a sad artifact (please leave a comment if you know how to get rid of it):

Let’s assign this expression to Yacas and let’s have a look at the degrees of the three variables x, y and z:

yac_assign(expr, "POLY")
yac_str("Degree(POLY, x)")
## [1] "8"
yac_str("Degree(POLY, y)")
## [1] "8"
yac_str("Degree(POLY, z)")
## [1] "4"

The degrees are 8, 8 and 4 respectively. So we can get all possible combinations \((i,j,k)\) of \(x^iy^jz^k\) with the blockparts function of the partitions package:

combins <- blockparts(c(8L, 8L, 4L))
## [1]   3 405

There are \(405\) such combinations. Of course they don’t all appear in the polynomial, and that is the point we will improve later. For the moment we will iterate over all these combinations. Here is the function which takes one combination and returns the corresponding coefficient of \(x^iy^jz^k\) in terms of \(a\) and \(b\), written in POV-Ray syntax:

f <- function(combin) {
  xyz <- sprintf("xyz(%s): ", toString(combin))
  coef <- yac_str(sprintf(
    "ExpandBrackets(Coef(Coef(Coef(POLY, x, %d), y, %d), z, %d))",
    combin[1L], combin[2L], combin[3L]
  if(coef == "0") return(NULL)
  coef <- gsub("([ab])\\^(\\d+)", "pow(\\1,\\2)", x = coef)
  paste0(xyz, coef, ",")
# for example:
f(c(2L, 6L, 0L))
## [1] "xyz(2, 6, 0): 2*pow(b,2)+2*b*a,"

Then we get the POV-Ray code as follows:

povray <- apply(combins, 2L, f)
cat(povray, sep = "\n", file = "SolidMobiusStrip.pov")

The file SolidMobiusStrip.pov:

xyz(4, 0, 0): pow(a,2)*pow(b,2)-2*pow(a,2)*b+pow(a,2),
xyz(6, 0, 0): (-2)*pow(a,2)*b-2*pow(a,2),
xyz(8, 0, 0): pow(a,2),
xyz(2, 2, 0): 2*pow(b,2)*pow(a,2)-2*pow(b,2)*a+(-2)*b*pow(a,2)+2*b*a,
xyz(4, 2, 0): (-2)*pow(b,2)*a+(-4)*b*pow(a,2)-4*b*a-2*pow(a,2),
xyz(6, 2, 0): 2*pow(a,2)+2*a*b,
xyz(0, 4, 0): pow(b,2)*pow(a,2)-2*pow(b,2)*a+pow(b,2),
xyz(2, 4, 0): (-4)*pow(b,2)*a-2*pow(b,2)+(-2)*b*pow(a,2)-4*b*a,
xyz(4, 4, 0): pow(b,2)+4*b*a+pow(a,2),
xyz(0, 6, 0): (-2)*pow(b,2)*a-2*pow(b,2),
xyz(2, 6, 0): 2*pow(b,2)+2*b*a,
xyz(0, 8, 0): pow(b,2),
xyz(3, 1, 1): 4*pow(a,2)*b-4*pow(a,2)+(-4)*a*pow(b,2)+4*a*b,
xyz(5, 1, 1): 4*pow(a,2)-4*a*b,
xyz(1, 3, 1): 4*pow(a,2)*b+(-4)*a*pow(b,2)-4*a*b+4*pow(b,2),
xyz(3, 3, 1): 4*pow(a,2)-4*pow(b,2),
xyz(1, 5, 1): 4*a*b-4*pow(b,2),
xyz(4, 0, 2): (-2)*a*pow(b,2)+2*a*b,
xyz(6, 0, 2): 2*b*a,
xyz(2, 2, 2): (-2)*pow(b,2)*a+6*pow(b,2)+(-2)*b*pow(a,2)-8*b*a+6*pow(a,2),
xyz(4, 2, 2): (-2)*pow(a,2)+10*a*b-2*pow(b,2),
xyz(0, 4, 2): (-2)*b*pow(a,2)+2*b*a,
xyz(2, 4, 2): (-2)*pow(a,2)+10*a*b-2*pow(b,2),
xyz(0, 6, 2): 2*a*b,
xyz(3, 1, 3): 4*pow(b,2)-4*b*a,
xyz(1, 3, 3): (-4)*pow(a,2)+4*a*b,
xyz(4, 0, 4): pow(b,2),
xyz(2, 2, 4): 2*a*b,
xyz(0, 4, 4): pow(a,2)

Now we will restrict the \(405\) combinations. There are only \(29\) combinations of exponents in the polynomial expansion. How to get them? With spray. We don’t care if there are rational numbers in the polynomial because we will take the exponents only.

x <- lone(1L, 5L)
y <- lone(2L, 5L)
z <- lone(3L, 5L)
a <- lone(4L, 5L)
b <- lone(5L, 5L)
P <- ((x*x+y*y+1)*(a*x*x+b*y*y) + z*z*(b*x*x+a*y*y) - 
        2*(a-b)*x*y*z - a*b*(x*x+y*y))^2 - 
exponents <- index(P)
exponents_xyz <- unique(exponents[, c(1L, 2L, 3L)])
## [1] 29  3

Indeed, there are \(29\) combinations. Now we can proceed as before and get the POV-Ray within a couple of seconds:

povray <- apply(exponents_xyz, 1L, f)
cat(povray, sep = "\n", file = "SolidMobiusStrip.pov")
To leave a comment for the author, please follow the link and comment on their blog: Saturn Elephant. offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)