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.

score 1 · Answer 1 · answered Jul 01 '23 at 04:26

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)

Heh, looks fantastic :). If you don't mind to add the rendered image, please? I was using POVray about 25 years ago, and now I feel great nostalgia for it :). — Grzegorz Sapijaszko, Jul 01 '23 at 19:14
@GrzegorzSapijaszko the image is here: https://stackoverflow.com/q/72316222/1100107 — Stéphane Laurent, Jul 02 '23 at 16:21
And for those of us stumbling to render the above in `povray`, either here or at Saturn, if you could point to a properly formed `that_annoying_line_artifact.pov`, that uses the output of this fascinating `Ryacas`, `partition`, `spray`, `povray` pipeline. — Chris, Jul 10 '23 at 18:43
@Chris It's probably on my [youtube channel](https://www.youtube.com/@stla3716/videos) with a link to the code. — Stéphane Laurent, Jul 10 '23 at 20:00
Unfortunately not, but perhaps, transitioning from version, camera, lights to action, are the above ex. `xyz(0,0,2): 2*b*a` evaluated to `xyz(0, 6, 2): 0.045000000000000005,` by povray when wrapped in `polynomial {8 ...`when presented in their polynomial form? The link above seems to suggest that evaluated values are in the `.pov` rather than the polynomial (29) strings. Strictly, none of my question here addresses finding which polynomial is leaving y and(x,z -though to lesser extent) unbounded to the mobius extent/skin, but being able to 'look' at it might help., — Chris, Jul 10 '23 at 22:31
@Chris It seems this is an antialisaing issue. When I change the antialiasing method (there are three methods), this is worse. — Stéphane Laurent, Jul 11 '23 at 16:57

score 0 · Answer 2 · answered Jul 02 '23 at 21:28

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)"

score 0 · Answer 3 · answered Jul 03 '23 at 20:44

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.

I think you mean `blockparts(c(8, 8, 4))`, since the degree of `z` is 4. — Stéphane Laurent, Jul 04 '23 at 08:37

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

3 Answers3