文档建议反馈控制台
首页
学习
活动
专区
圈层
工具
MCP广场
文章/答案/技术大牛
发布
社区首页 >问答首页 >用加拉加斯(Sympy)代替多项式的变量

问用加拉加斯(Sympy)代替多项式的变量
EN

Stack Overflow用户
提问于 2022-07-15 16:06:10
回答 4查看 115关注 0票数 0

在四变量x，y，z，w中有一个长多项式：

代码语言：javascript
复制
((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

我正在使用R.，我想使用加拉加斯包(Sympy的包装器)在进行变量更改后作为多项式获得这个表达式。即，我想替换x，y，z和w

代码语言：javascript
复制
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

分别使用。我没有运气就试了subs。这是我发现的唯一的工作方式：

代码语言：javascript
复制
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]"
)

但30分钟后，poly的计算仍未完成。有没有更有效的方法？

sympy
python
r
EN

回答 4

Stack Overflow用户

回答已采纳

发布于 2022-07-15 17:41:48

而不是生成完整的表达式，然后要求系数的展开，当你完成，你可以把它逐项。我已经用完整的表达式做过两次了，但是您可以使用toy表达式，然后是完整的表达式，作为注释：

代码语言：javascript
复制
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])))
票数 1
EN

Stack Overflow用户

发布于 2022-07-15 19:01:11

基于对注释中问题的澄清，我们假定所要将结果扩展为x、y、z和w中的多项式，使得系数是a，b，c，d中的多项式，输出是以a、b、c和d表示的每个单项多项式的指数和相应系数多项式的指数。

这使用字符替换，创建ch2，然后从mpoly调用mp来创建mpoly。在内部，它们是由一组单调词组成的。我们定义并使用f迭代列表，提取a，b，c和d中的系数和x，y，z，w中的单项式，用tapply来收集和求和相同x，y，z，w的所有a，b，c，d系数。result是a，b，c，d系数的特征向量，其名称为x，y，z，w单调式。

在一台相对较慢的笔记本电脑上，这需要286秒。

代码语言：javascript
复制
library(magrittr) 
library(mpoly)

pt <- proc.time()[[3]]

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_", .) %>%
  gsub("x_", "(a*x - b*y - c*z - d*w)", .) %>%
  gsub("y_", "(a*y + b*x + c*w - d*z)", .) %>%
  gsub("z_", "(a*z - b*w + c*x + d*y)", .) %>%
  gsub("w_", "(a*w + b*z - c*y + d*x)", .)

p <- mp(ch2)

# x is a component of an mpoly.  names are variables to extract 
# plus "coef" or if "coef" not among names then coef is set to 1.
# output is character vector.
f <- function(x, names) {
   names2 <- unique(c(names, "coef"))
   x0 <- x[names(x) %in% names2]
   if (!"coef" %in% names) x0["coef"] <- 1
   p <- structure(list(x0), class = "mpoly")
   print(p, stars = TRUE, silent = TRUE)
}

xx <- sapply(exponents(p), function(x) toString(x[c("x","y","z","w")]))
aa <- sapply(p, f, c("a", "b", "c", "d", "coef"))
r <- tapply(aa, xx, paste, collapse = "+")
r <- paste(names(r), unname(r), sep = ", ")

proc.time()[[3]] - pt
## [1] 285.94
票数 1
EN

Stack Overflow用户

发布于 2022-07-15 22:43:43

这也使用了Python中的SymPy，但可能更容易将其转换为R：

代码语言：javascript
复制
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)))

请注意，我使用小写poly而不是Poly，因为它更快。我还认为，在扩展sqrt(3)的过程中，用符号t代替poly会大大加快速度。最后，将gmpy2与SymPy一起安装将使这类事情变得更快。

这样，您就可以将其转换为字符串(听起来就像您想要做的那样)。结果有超过一百万个字符，但我将只显示前1000个字符：

代码语言：javascript
复制
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'
票数 0
EN
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：

https://stackoverflow.com/questions/72996803

复制
相关文章

相似问题

领券

腾讯云开发者

扫码关注腾讯云开发者

扫码关注腾讯云开发者

领取腾讯云代金券

热门产品

热门推荐

更多推荐

Copyright © 2013 - 2026 Tencent Cloud. All Rights Reserved. 腾讯云 版权所有 

深圳市腾讯计算机系统有限公司 ICP备案/许可证号：粤B2-20090059 粤公网安备44030502008569号

腾讯云计算（北京）有限责任公司 京ICP证150476号 |  京ICP备11018762号

问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档