Substitute the variables of a polynomial with caracas (Sympy)-CodePudding

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'

CodePudding user response：

With the help of smichr's answer, I finally managed to get what I want in R. Here is an example shorter than my real case.

library(caracas)
sympy <- get_sympy()

# define the variables x,y,z and the constants a,b
# as well as auxiliary variables X,Y,Z
def_sym(x, y, z, a, b, X, Y, Z)

# define expression in terms of X,Y,Z
expr <- sympy$parse_expr("X**2   X*Z/3   Y   1/2") 

# define the substitutions in new variables
Xs <- a*x   b*y
Ys <- x   9*z
Zs <- a*y   z^2

# define the list of replacements
substitutions <- list(
  X = Xs$pyobj,
  Y = Ys$pyobj,
  Z = Zs$pyobj
)

# perform the substitutions
expr <- expr$subs(substitutions)

# 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){
  sprintf("xyz(%s): %s,", degree, povray[[degree]])
}, character(1L))
cat(polynomial, sep = "\n")
# xyz(0, 0, 0): 1/2,
# xyz(1, 0, 0): 1,
# xyz(2, 0, 0): a**2,
# xyz(0, 2, 0): a*b/3   b**2,
# xyz(1, 1, 0): a**2/3   2*a*b,
# xyz(0, 0, 1): 9,
# xyz(1, 0, 2): a/3,
# xyz(0, 1, 2): b/3,