Fast expansion of a polynomial with R – part 2
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:
library(Ryacas)
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:
library(partitions) combins <- blockparts(c(8L, 8L, 4L)) dim(combins) ## [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.
library(spray)
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 -
4*(x*x+y*y)*(a*x*x+b*y*y-x*y*z*(a-b))^2
exponents <- index(P)
exponents_xyz <- unique(exponents[, c(1L, 2L, 3L)])
dim(exponents_xyz)
## [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")
R-bloggers.com 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.