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?
CodePudding user response：
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])))
CodePudding user response：
This uses character substitution and mpoly package rather than caracas.
ch <- "((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"
ch2 <- ch |>
gsub("([xyzw])", "\\1_", x = _) |>
gsub("x_", "(a*x - b*y - c*z - d*w)", x = _) |>
gsub("y_", "(a*y b*x c*w - d*z)", x = _) |>
gsub("z_", "(a*z - b*w c*x d*y)", x = _) |>
gsub("w_", "(a*w b*z - c*y d*x)", x = _)
library(mpoly)
p <- mp(ch2)
CodePudding user response：
This is also using SymPy from Python but might be easier to translate to R:
from sympy import *
from sympy.abc import x,y,z,w,a,b,c,d,t
eq = sympify('((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')
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}
result = poly(eq.xreplace(reps).subs(sqrt(3), t), domain=QQ)
result = poly(result.as_expr().subs(t, sqrt(3)))
Note that I used lower case poly rather than Poly because it is faster. I think also that replacing sqrt(3) with a symbol t speeds things up a lot during the expansion in poly. Lastly installing gmpy2 alongside SymPy will make this sort of thing a lot faster.
With that you can then convert it to a string (it sounds like that's what you want to do). The result has over a million characters but I'll just show the first 1000:
In [22]: s = str(result)
In [23]: len(s)
Out[23]: 1603555
In [24]: s[:1000]
Out[24]: 'Poly(x**12*a**12 6*x**12*a**10*b**2 6*x**12*a**10*c**2 6*x**12*a**10*d**2 15*x**12*a**8*b**4 30*x**12*a**8*b**2*c**2 30*x**12*a**8*b**2*d**2 15*x**12*a**8*c**4 30*x**12*a**8*c**2*d**2 15*x**12*a**8*d**4 20*x**12*a**6*b**6 60*x**12*a**6*b**4*c**2 60*x**12*a**6*b**4*d**2 60*x**12*a**6*b**2*c**4 120*x**12*a**6*b**2*c**2*d**2 60*x**12*a**6*b**2*d**4 20*x**12*a**6*c**6 60*x**12*a**6*c**4*d**2 60*x**12*a**6*c**2*d**4 20*x**12*a**6*d**6 15*x**12*a**4*b**8 60*x**12*a**4*b**6*c**2 60*x**12*a**4*b**6*d**2 90*x**12*a**4*b**4*c**4 180*x**12*a**4*b**4*c**2*d**2 90*x**12*a**4*b**4*d**4 60*x**12*a**4*b**2*c**6 180*x**12*a**4*b**2*c**4*d**2 180*x**12*a**4*b**2*c**2*d**4 60*x**12*a**4*b**2*d**6 15*x**12*a**4*c**8 60*x**12*a**4*c**6*d**2 90*x**12*a**4*c**4*d**4 60*x**12*a**4*c**2*d**6 15*x**12*a**4*d**8 6*x**12*a**2*b**10 30*x**12*a**2*b**8*c**2 30*x**12*a**2*b**8*d**2 60*x**12*a**2*b**6*c**4 120*x**12*a**2*b**6*c**2*d**2 60'