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