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