--- title: "Expanding a polynomial with 'caracas'" author: "Stéphane Laurent" date: '2022-06-07' tags: R, povray, graphics, maths, python rbloggers: yes output: md_document: variant: markdown preserve_yaml: true html_document: highlight: kate keep_md: no highlighter: pandoc-solarized --- 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:  r 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:  r "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$:  r sympy <- get_sympy() sympy$Poly$nth(poly$pyobj, 1L, 0L, 1L)  This gives:  r 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:  r 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: ![](./figures/ICN5D_01.gif) Credit to '**ICN5D**' for the isosurface.