Substitute the variables of a polynomial with caracas (Sympy)

I have a long polynomial in four variables x, y, z, w:

((x^2+y^2+z^2+w^2+145/3)^2-4*(9*z^2+16*w^2))^2*((x^2+y^2+z^2+w^2+145/3)^2+296*(x^2+y^2)-4*(9*z^2+16*w^2)) -16*(x^2+y^2)*(x^2+y^2+z^2+w^2+145/3)^2*(37*(x^2+y^2+z^2+w^2+145/3)^2-1369*(x^2+y^2)-7*(225*z^2+448*w^2)) -16*sqrt(3)/9*(x^3-3*x*y^2)*(110*(x^2+y^2+z^2+w^2+145/3)^3 -148*(x^2+y^2+z^2+w^2+145/3)*(110*x^2+110*y^2-297*z^2+480*w^2)) -64*(x^2+y^2)*(3*(729*z^4+4096*w^4)+168*(x^2+y^2)*(15*z^2-22*w^2)) +64*(12100/27*(x^3-3*x*y^2)^2 -7056*(3*x^2*y-y^3)^2) -592240896*z^2*w^2

I’m working with R. I want to use the caracas package (a wrapper of Sympy) to get this expression as a polynomial after doing a change of variables. Namely, I want to substitue x, y, z and w by

a*x - b*y - c*z - d*w,
a*y + b*x + c*w - d*z,
a*z - b*w + c*x + d*y,
a*w + b*z - c*y + d*x

respectively. I tried subs with no luck. Here is the only working way I found:

library(caracas)
def_sym(x, y, z, w, a, b, c, d)
X <- a*x - b*y - c*z - d*w
Y <- a*y + b*x + c*w - d*z
Z <- a*z - b*w + c*x + d*y
W <- a*w + b*z - c*y + d*x

expr <- ((X^2+Y^2+Z^2+W^2+145/3)^2-4*(9*Z^2+16*W^2))^2*((X^2+Y^2+Z^2+W^2+145/3)^2+296*(X^2+Y^2)-4*(9*Z^2+16*W^2)) -16*(X^2+Y^2)*(X^2+Y^2+Z^2+W^2+145/3)^2*(37*(X^2+Y^2+Z^2+W^2+145/3)^2-1369*(X^2+Y^2)-7*(225*Z^2+448*W^2)) -16*sqrt(3)/9*(X^3-3*X*Y^2)*(110*(X^2+Y^2+Z^2+W^2+145/3)^3 -148*(X^2+Y^2+Z^2+W^2+145/3)*(110*X^2+110*Y^2-297*Z^2+480*W^2)) -64*(X^2+Y^2)*(3*(729*Z^4+4096*W^4)+168*(X^2+Y^2)*(15*Z^2-22*W^2)) +64*(12100/27*(X^3-3*X*Y^2)^2 -7056*(3*X^2*Y-Y^3)^2) -592240896*Z^2*W^2

poly <- sympy_func(
  expr, "Poly", domain = "QQ[a,b,c,d]"
)

But after 30 minutes the computation of poly is not finished. Is there a more efficient way?

Answer

Rather than generating the full expression and then requesting coefficients of the expansion when you are done, you can take it term-by-term. I have done so twice with the full expression, but you can do so with the toy expression that is followed by your full expression as a comment:

from sympy import *
from sympy.parsing.sympy_parser import *
from sympy.abc import x,y,z,w,a,b,c,d
eq = x*y + 1#parse_expr('((x^2+y^2+z^2+w^2+145/3)^2-4*(9*z^2+16*w^2))^2*((x^2+y^2+z^2+w^2+145/3)^2+296*(x^2+y^2)-4*(9*z^2+16*w^2)) -16*(x^2+y^2)*(x^2+y^2+z^2+w^2+145/3)^2*(37*(x^2+y^2+z^2+w^2+145/3)^2-1369*(x^2+y^2)-7*(225*z^2+448*w^2)) -16*sqrt(3)/9*(x^3-3*x*y^2)*(110*(x^2+y^2+z^2+w^2+145/3)^3 -148*(x^2+y^2+z^2+w^2+145/3)*(110*x^2+110*y^2-297*z^2+480*w^2)) -64*(x^2+y^2)*(3*(729*z^4+4096*w^4)+168*(x^2+y^2)*(15*z^2-22*w^2)) +64*(12100/27*(x^3-3*x*y^2)^2 -7056*(3*x^2*y-y^3)^2) -592240896*z^2*w^2', transformations=T[:])
reps = {x: a*x - b*y - c*z - d*w,y:a*y + b*x + c*w - d*z,z:a*z - b*w + c*x + d*y,w:a*w + b*z - c*y + d*x}
eq = eq.xreplace(reps)
c = {}
for i in Add.make_args(eq):
    f = i.xreplace(reps).expand()
    for s in Add.make_args(f):
        co, mo = s.as_coeff_mul(x,y,z,w)
        c.setdefault(Mul._from_args(mo), []).append(co)

for k in c:
    print(k,Add(*c[k])))