I wanted to plot an algebraic isosurface with POV-Ray but the expression of the polynomial defining the isosurface was very long (the polynomial had degree 12). Moreover there was a square root in the coefficients (\(\sqrt{3}\)) as well as \(\cos t\) and \(\sin t\), where \(t\) is a parameter I wanted to vary in order to make an animation. So I needed a tool able to expand a polynomial with some literal values in the coefficients. This is not possible with the Ryacas package.

I finally found this tool: the caracas package. It allows to use the Python library SymPy in R. I didn’t carefully read its documentation yet, I don’t know whether it has other features. But this feature is a great one.

Here is a small example:

library(caracas)
def_sym(x, y, z, a, b) # symbolic values
poly <- sympy_func(
  x^2 + a*x^2 + 2/3*y + b*y + x*z + a*x*z, "Poly", domain = "QQ[a,b]"
)
as.character(poly)

This gives:

"Poly((a + 1)*x^2 + (a + 1)*x*z + (b + 2/3)*y, x, y, z, domain='QQ[a,b]')"

That is great. Here QQ[a,b] is the field \(\mathbb{Q}[a,b]\). I lost a significant part of my knowledge in mathematics but I think this is a field. It doesn’t matter. Roughly speaking, this is the set of rational numbers to which we add the two elements \(a\) and \(b\). So there are treated as constants, as if they were some numbers.

To get a coefficient, for example the one of \(xz = x^1y^0z^1\):

sympy <- get_sympy()
sympy$Poly$nth(poly$pyobj, 1L, 0L, 1L)

This gives:

a + 1

Everything needed for writing the POV-Ray code was there. I wrote a small script in addition to generate this code. I show it below with the above small example:

library(caracas)
library(partitions) # to get the compositions of an integer, 
                    # representing the degrees with a given total
def_sym(x, y, z, a, b) 
poly <- sympy_func(
  x^2 + a*x^2 + 2/3*y + b*y + x*z + a*x*z, "Poly", domain = "QQ[a,b]"
)
sympy <- get_sympy()
f <- function(comp){
  xyz <- sprintf("xyz(%s): ", toString(comp))
  coef <- sympy$Poly$nth(poly$pyobj, comp[1L], comp[2L], comp[3L])
  if(coef == 0) return(NULL)
  paste0(xyz, coef, ",")
}
for(deg in 0L:2L){
  comps <- compositions(deg, 3L)
  povray <- apply(comps, 2L, f, simplify = FALSE)
  cat(
    unlist(povray), sep = "\n", file = "povray.txt", append = deg > 0L
  )
}

And here is the povray.txt file generated by this script:

xyz(0, 1, 0): b + 2/3,
xyz(2, 0, 0): a + 1,
xyz(1, 0, 1): a + 1,

One just has to remove the trailing comma, and this the desired POV-Ray code.

I won’t leave you without showing the animation:

Credit to ‘ICN5D’ for the isosurface.