Expanding a polynomial with ‘caracas’

[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.

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:

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

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(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)
    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.

To leave a comment for the author, please follow the link and comment on their blog: Saturn Elephant.

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.

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)