Expanding a symbolic multivariate polynomial with R for usage in POV-Ray

Question

Consider the following trivariate polynomial with two parameters a and b:

P(x,y,z) = ((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 POV-Ray, I want to plot the algebraic isosurface of equation P(x,y,z)=0 for some values of a and b. In POV-Ray, one has to define the polynomial by listing its monomials where each monomial is given as follows:

xyz(i,j,k): coef

where i,j,k are the exponents and coef is the coefficient of x^i y^j z^k.

For example one has the monomials b^2 y^8 = b^2 x^0 y^8 z^0 and 2*b*a x^6 z^2 = 2*b*a x^6 y^0 z^2 and then they must be given as follows:

xyz(0, 8, 0): pow(b,2),
xyz(6, 0, 2): 2*b*a,
......

It is not funny to manually expand the given polynomial. I want to generate this POV-Ray code with R.

Answer 1

I found a solution with the packages Ryacas and partitions.

First, let's enter the polynomial in Ryacas:

library(Ryacas)
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"
yac_assign(expr, "POLY")

Let's see the degree of each variable:

yac_str("Degree(POLY, x)") # 8
yac_str("Degree(POLY, y)") # 8
yac_str("Degree(POLY, z)") # 4

So we need the exponents i,j,k with possible values for i and j between 0 and 8 and possible values for k between 0 and 4. We use the partitions package to generate the compositions of these three integers. For example, the compositions summing to 4 are obtained by:

library(partitions)
compositions(4, 3)
                                  
[1,] 4 3 2 1 0 3 2 1 0 2 1 0 1 0 0
[2,] 0 1 2 3 4 0 1 2 3 0 1 2 0 1 0
[3,] 0 0 0 0 0 1 1 1 1 2 2 2 3 3 4

Here is a function which takes a composition (i,j,k) and returns the coefficient of x^i y^j z^k prepended by xyz(i,j,k)::

f <- function(comp) {
  if(comp[3L] > 4L) return(NULL)
  xyz <- sprintf("xyz(%s): ", toString(comp))
  coef <- yac_str(sprintf(
    "ExpandBrackets(Coef(Coef(Coef(POLY, x, %d), y, %d), z, %d))",
    comp[1L], comp[2L], comp[3L]
  ))
  if(coef == "0") return(NULL)
  # replace a^p with pow(a,p) for POV-Ray, same for b^p:
  coef <- gsub("([ab])\\^(\\d+)", "pow(\\1,\\2)", x = coef)
  paste0(xyz, coef, ",")
}

Then we loop over the degree:

for(deg in 0L:8L){
  comps <- compositions(deg, 3L)
  povray <- apply(comps, 2L, f, simplify = FALSE)
  cat(
    unlist(povray), sep = "\n", file = "polynomial.pov", append = deg > 0L
  )
}

We obtain the desired code in the file polynomial.pov:

xyz(4, 0, 0): pow(a,2)*pow(b,2)-2*pow(a,2)*b+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(0, 4, 0): pow(b,2)*pow(a,2)-2*pow(b,2)*a+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(1, 3, 1): 4*pow(a,2)*b+(-4)*a*pow(b,2)-4*a*b+4*pow(b,2),
xyz(6, 0, 0): (-2)*pow(a,2)*b-2*pow(a,2),
xyz(4, 2, 0): (-2)*pow(b,2)*a+(-4)*b*pow(a,2)-4*b*a-2*pow(a,2),
xyz(2, 4, 0): (-4)*pow(b,2)*a-2*pow(b,2)+(-2)*b*pow(a,2)-4*b*a,
xyz(0, 6, 0): (-2)*pow(b,2)*a-2*pow(b,2),
xyz(4, 0, 2): (-2)*a*pow(b,2)+2*a*b,
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(0, 4, 2): (-2)*b*pow(a,2)+2*b*a,
xyz(5, 1, 1): 4*pow(a,2)-4*a*b,
xyz(3, 3, 1): 4*pow(a,2)-4*pow(b,2),
xyz(1, 5, 1): 4*a*b-4*pow(b,2),
xyz(3, 1, 3): 4*pow(b,2)-4*b*a,
xyz(1, 3, 3): (-4)*pow(a,2)+4*a*b,
xyz(8, 0, 0): pow(a,2),
xyz(6, 2, 0): 2*pow(a,2)+2*a*b,
xyz(4, 4, 0): pow(b,2)+4*b*a+pow(a,2),
xyz(2, 6, 0): 2*pow(b,2)+2*b*a,
xyz(0, 8, 0): pow(b,2),
xyz(6, 0, 2): 2*b*a,
xyz(4, 2, 2): (-2)*pow(a,2)+10*a*b-2*pow(b,2),
xyz(2, 4, 2): (-2)*pow(a,2)+10*a*b-2*pow(b,2),
xyz(0, 6, 2): 2*a*b,
xyz(4, 0, 4): pow(b,2),
xyz(2, 2, 4): 2*a*b,
xyz(0, 4, 4): pow(a,2)

Answer 2

Here is a solution with the spray package.

library(spray)

x <- lone(1, 5)
y <- lone(2, 5)
z <- lone(3, 5)
a <- lone(4, 5)
b <- lone(5, 5)

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


nterms <- nrow(P[["index"]])

coeffs <- P[["value"]]

XYZpowers <- P[["index"]][, c(1L, 2L, 3L)]
XYZ <- apply(XYZpowers, 1L, function(comp) {
  sprintf("xyz(%s): ", toString(comp))
})

ABpowers  <- P[["index"]][, c(4L, 5L)]

a <- lone(1, 2)
b <- lone(2, 2)
abPolys <- lapply(1L:nterms, function(i) {
  pows <- ABpowers[i, ]
  coef <- coeffs[i]
  coef * a^pows[1L] * b^pows[2L]
})

AB <- split(abPolys, XYZ)

ABgroups <- sapply(AB, function(poly) {
  Reduce(`+`, poly)
}, simplify = FALSE)

ABgroups <- Filter(Negate(is.empty), ABgroups)


asCharacter <- function(poly) {
  op <- options(sprayvars = letters)
  x <- capture.output(print_spray_polyform(poly))
  options(op)
  x
}

ABchar <- sapply(ABgroups, asCharacter, simplify = FALSE)

asPOVRay <- function(poly) {
  gsub("([ab])\\^(\\d+)", "pow(\\1,\\2)", x = poly)
}

ABpov <- sapply(ABchar, asPOVRay)

paste0(names(ABpov), unlist(ABpov))

#  [1] "xyz(0, 4, 0): +pow(b,2) +pow(a,2)*pow(b,2) -2*a*pow(b,2)"                
#  [2] "xyz(0, 4, 2): -2*pow(a,2)*b +2*a*b"                                      
#  [3] "xyz(0, 4, 4): +pow(a,2)"                                                 
#  [4] "xyz(0, 6, 0): -2*a*pow(b,2) -2*pow(b,2)"                                 
#  [5] "xyz(0, 6, 2): +2*a*b"                                                    
#  [6] "xyz(0, 8, 0): +pow(b,2)"                                                 
#  [7] "xyz(1, 3, 1): -4*a*b +4*pow(a,2)*b +4*pow(b,2) -4*a*pow(b,2)"            
#  [8] "xyz(1, 3, 3): +4*a*b -4*pow(a,2)"                                        
#  [9] "xyz(1, 5, 1): -4*pow(b,2) +4*a*b"                                        
# [10] "xyz(2, 2, 0): +2*a*b -2*pow(a,2)*b +2*pow(a,2)*pow(b,2) -2*a*pow(b,2)"   
# [11] "xyz(2, 2, 2): -2*a*pow(b,2) +6*pow(b,2) -2*pow(a,2)*b -8*a*b +6*pow(a,2)"
# [12] "xyz(2, 2, 4): +2*a*b"                                                    
# [13] "xyz(2, 4, 0): -4*a*b -2*pow(a,2)*b -4*a*pow(b,2) -2*pow(b,2)"            
# [14] "xyz(2, 4, 2): -2*pow(a,2) -2*pow(b,2) +10*a*b"                           
# [15] "xyz(2, 6, 0): +2*pow(b,2) +2*a*b"                                        
# [16] "xyz(3, 1, 1): +4*pow(a,2)*b +4*a*b -4*pow(a,2) -4*a*pow(b,2)"            
# [17] "xyz(3, 1, 3): +4*pow(b,2) -4*a*b"                                        
# [18] "xyz(3, 3, 1): +4*pow(a,2) -4*pow(b,2)"                                   
# [19] "xyz(4, 0, 0): -2*pow(a,2)*b +pow(a,2) +pow(a,2)*pow(b,2)"                
# [20] "xyz(4, 0, 2): -2*a*pow(b,2) +2*a*b"                                      
# [21] "xyz(4, 0, 4): +pow(b,2)"                                                 
# [22] "xyz(4, 2, 0): -4*a*b -4*pow(a,2)*b -2*pow(a,2) -2*a*pow(b,2)"            
# [23] "xyz(4, 2, 2): +10*a*b -2*pow(b,2) -2*pow(a,2)"                           
# [24] "xyz(4, 4, 0): +pow(b,2) +4*a*b +pow(a,2)"                                
# [25] "xyz(5, 1, 1): +4*pow(a,2) -4*a*b"                                        
# [26] "xyz(6, 0, 0): -2*pow(a,2) -2*pow(a,2)*b"                                 
# [27] "xyz(6, 0, 2): +2*a*b"                                                    
# [28] "xyz(6, 2, 0): +2*pow(a,2) +2*a*b"                                        
# [29] "xyz(8, 0, 0): +pow(a,2)" 

Answer 3

Those are indeed instances of the spray and partitions packages being used as intended. I do very much like the method of looping over the degree from 0-8 and then using compositions(), although it is possible to streamline this very slightly:

suppressMessages(library("partitions"))
summary(blockparts(rep(8,3),8,TRUE))
#>                                                 
#> [1,] 0 1 2 3 4 5 6 7 8 0 ... 0 1 2 0 1 0 0 1 0 0
#> [2,] 0 0 0 0 0 0 0 0 0 1 ... 0 0 0 1 1 2 0 0 1 0
#> [3,] 0 0 0 0 0 0 0 0 0 0 ... 6 6 6 6 6 6 7 7 7 8

^{Created on 2023-07-04 with reprex v2.0.2}

It is somewhat mortifying to me that I can't easily reproduce Stéphane's wonderful solution with mvp. Here I would observe that the applications above are a rare case where the symbolic nature of the mvp package is counterproductive, and the somewhat more rigid structure of spray is actually helpful. We see that x,y,z and a,b are confined to their own columns, and Stéphane exploited this immutability of the spray package.