Expanding a polynomial with 'caracas', part 2
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.