# Expanding a polynomial with ‘caracas’, part 2

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

Last month, I posted an
article
showing a way to expand a polynomial in R when the coefficients of the
polynomial contain some literal values, with the help of the
**caracas** package. Today I wanted to apply it with a
polynomial expression having about 500 characters, and highly
factorized. The method took more than 30 minutes, so I looked for a more
efficient one.

Thanks to some help on StackOverflow, I came to the following method which is more efficient. It consists in splitting the expression according to its additive terms and to work on each term, instead of expanding the whole polynomial. In the example below I take the polynomial expression defining the isosurface equation of a Solid Möbius strip.

library(caracas) sympy <- get_sympy() # define the variables x,y,z and the constants a,b def_sym(x, y, z, a, b) # define expression expr <- sympy$parse_expr( "((x*x+y*y+1)*(a*x*x+b*y*y)+z*z*(b*x*x+a*y*y)-2*(a-b)*x*y*z-a*b*(x*x+y*y))**2-4*(x*x+y*y)*(a*x*x+b*y*y-x*y*z*(a-b))**2" ) # extraction of monomials in the 'povray' list povray <- list() terms <- sympy$Add$make_args(expr) for(term in terms){ f <- term$expand() fterms <- sympy$Add$make_args(f) for(fterm in fterms){ decomp <- fterm$as_coeff_mul(x$pyobj, y$pyobj, z$pyobj) coef <- decomp[[1]] mono <- decomp[[2]] polexpr <- sympy$Mul$fromiter(mono) poly <- polexpr$as_poly(x$pyobj, y$pyobj, z$pyobj) degree <- toString(poly$monoms()[[1]]) if(degree %in% names(povray)){ povray[[degree]] <- sympy$Add(povray[[degree]], coef) }else{ povray[[degree]] <- coef } } } polynomial <- vapply(names(povray), function(degree){ coeff <- povray[[degree]] |> gsub("([ab])\\*\\*(\\d+)", "pow(\\1,\\2)", x = _) sprintf("xyz(%s): %s,", degree, coeff) }, character(1L)) cat(polynomial, sep = "\n", file = "SolidMobiusStrip.txt")

At the last step I use `gsub`

to replace the powers like
`a**2`

to their POV-Ray syntax `pow(a,2)`

. The
above code writes this POV-Ray code:

xyz(4, 0, 0): pow(a,2)*pow(b,2) - 2*pow(a,2)*b + pow(a,2), xyz(8, 0, 0): pow(a,2), xyz(0, 4, 0): pow(a,2)*pow(b,2) - 2*a*pow(b,2) + pow(b,2), xyz(0, 8, 0): pow(b,2), xyz(6, 0, 0): -2*pow(a,2)*b - 2*pow(a,2), xyz(0, 6, 0): -2*a*pow(b,2) - 2*pow(b,2), xyz(4, 4, 0): pow(a,2) + 4*a*b + pow(b,2), xyz(0, 4, 4): pow(a,2), xyz(4, 0, 4): pow(b,2), xyz(4, 2, 0): -4*pow(a,2)*b - 2*pow(a,2) - 2*a*pow(b,2) - 4*a*b, xyz(6, 2, 0): 2*pow(a,2) + 2*a*b, xyz(2, 4, 0): -2*pow(a,2)*b - 4*a*pow(b,2) - 4*a*b - 2*pow(b,2), xyz(2, 6, 0): 2*a*b + 2*pow(b,2), xyz(1, 3, 3): -4*pow(a,2) + 4*a*b, xyz(3, 1, 1): 4*pow(a,2)*b - 4*pow(a,2) - 4*a*pow(b,2) + 4*a*b, xyz(5, 1, 1): 4*pow(a,2) - 4*a*b, xyz(3, 3, 1): 4*pow(a,2) - 4*pow(b,2), xyz(2, 2, 0): 2*pow(a,2)*pow(b,2) - 2*pow(a,2)*b - 2*a*pow(b,2) + 2*a*b, xyz(4, 0, 2): -2*a*pow(b,2) + 2*a*b, xyz(0, 4, 2): -2*pow(a,2)*b + 2*a*b, xyz(6, 0, 2): 2*a*b, xyz(0, 6, 2): 2*a*b, xyz(2, 4, 2): -2*pow(a,2) + 10*a*b - 2*pow(b,2), xyz(4, 2, 2): -2*pow(a,2) + 10*a*b - 2*pow(b,2), xyz(1, 3, 1): 4*pow(a,2)*b - 4*a*pow(b,2) - 4*a*b + 4*pow(b,2), xyz(1, 5, 1): 4*a*b - 4*pow(b,2), xyz(3, 1, 3): -4*a*b + 4*pow(b,2), xyz(2, 2, 2): -2*pow(a,2)*b + 6*pow(a,2) - 2*a*pow(b,2) - 8*a*b + 6*pow(b,2), xyz(2, 2, 4): 2*a*b,

This is very fast for this example, but it still took 20 minutes with my case, which is a slight modification of an animation by ‘ICN5D’; here it is:

The difference with the original animation is that this one uses an isoclinic rotation for the animation.

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