# Expanding a polynomial with ‘caracas’

**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:

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.

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