--- title: "Expanding a polynomial with 'caracas', part 2" author: "Stéphane Laurent" date: '2022-07-16' 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 --- Last month, I posted an [article](https://laustep.github.io/stlahblog/posts/caracas01.html) 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. ``` {.r .numberLines} 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: ``` povray 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: ![](./figures/ICN5D_02.gif) The difference with the original animation is that this one uses an [isoclinic rotation](https://laustep.github.io/stlahblog/posts/slicedImplicitHypersurface.html) for the animation.