辛普森不能把普罗迪变蓝

Question

我使用的是Sympy1.0和Python2.7。我要计算前100个整数的和：

此代码成功运行

import sympy as sy
from sympy.tensor import IndexedBase, Idx
import numpy as np

x = sy.IndexedBase('x')
i = sy.symbols('i', cls=Idx)
s = sy.Sum(x[i], (i, 0, 100))
s_lambda = sy.lambdify(sy.DeferredVector('x'), s, 'numpy')
s_lambda(np.arange(101))

并按预期给出5050。但是当我试图用Product而不是Sum来做同样的事情时：

^{pr2}$

我得到了NameError: global name 'Product' is not defined 我做错什么了？有没有一个变通办法来得到我想要的？在

编辑1： 如果我事先不知道Product的极限怎么办。比如说

import sympy as sy
from sympy.tensor import IndexedBase, Idx
import numpy as np

x = sy.IndexedBase('x')
i = sy.symbols('i', cls=Idx)
n = sy.symbols('n', integer=True, positive=True)
s = sy.Product(x[i], (i, 0, n))
s_lambda = sy.lambdify((sy.DeferredVector('x'), n) s.doit(), 'numpy')
s_lambda(np.arange(101), 5)

编辑2： 我想找个解决办法。NameError: global name 'Product' is not defined由于以下原因而出现错误： 在

lambdastr((sy.DeferredVector('x'), n), p)

这就产生了：

lambda x,n: (Product(x[i], (i, 0, n)))

而对于Sum，我们得到了一个有效的lambda函数：

lambda x,n: ((builtins.sum(x[i] for i in range(0, n+1))))

此时，问题围绕Product函数的定义展开。根据手册，我可以通过dict注入我对函数的定义 在

def my_prod(a, b):
    # my implementation
    pass

my_fun = {"Product" : my_prod}
f = sy.lambdify((sy.DeferredVector('x'), n), p, modules=['numpy', my_fun])
f([1,2,3,4,5], 2)

问题是，当我尝试使用lambdified函数f时，出现了list indices must be integers, not Symbol错误。我想这是由于i是一个符号，而它应该是一个整数。我不明白为什么它在尝试调用my_prod之前没有传递实际的integer，就像对Sum的情况一样。在

Answer 1

当Product的极限是预先知道的

您可以通过调用.doit()将Product扩展到其组件部分来解决此问题：

In [104]: s = sy.Product(x[i], (i, 1, 10)); s
Out[104]: Product(x[i], (i, 1, 10))

In [105]: s.doit()
Out[105]: x[1]*x[2]*x[3]*x[4]*x[5]*x[6]*x[7]*x[8]*x[9]*x[10]

---

例如

^{pr2}$

印刷品

3628800

---

但是，如果您将.doit()与sy.Product(x[i], (i, 1, 100))一起使用，那么您将得到一个算术溢出，因为np.arange(101)具有{}或{}（取决于您的操作系统）和产品100!

In [109]: math.factorial(100)
Out[109]: 93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000

太大，不能存储在int32或int64数组值中。在

In [118]: np.iinfo('int64').max
Out[118]: 9223372036854775807

In [119]: np.iinfo('int64').max < math.factorial(100)
Out[119]: True

因此

s = sy.Product(x[i], (i, 1, 100))
s_lambda = sy.lambdify(sy.DeferredVector('x'), s.doit(), 'numpy')
print(s_lambda(np.arange(101)))

引发

RuntimeWarning: overflow encountered in long_scalars

并错误地打印0。在

---

如果将输入从dtype int64数组更改为Pythonint的列表，则 可以正确计算产品：

import sympy as sy
from sympy.tensor import IndexedBase, Idx
import numpy as np

x = sy.IndexedBase('x')
i = sy.symbols('i', cls=Idx)
s = sy.Product(x[i], (i, 1, 100))
s_lambda = sy.lambdify(sy.DeferredVector('x'), s.doit(), 'numpy')
print(s_lambda(np.arange(101).tolist()))

印刷品

93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000

---

当Product的极限是不是预先知道的

周围的工作 （非洲）变得更加复杂。如果使用调试器跟踪代码路径 当使用Sum时，您会发现 ^{} 调用以将Sum(x[i], (i, 0, n))转换为表达式builtins.sum(x[i] for i in range(0, n+1))。在

如果我们将_print_Product方法添加到NumPyPrinter（一个LambdaPrinter的子类）， 然后我们可以得到lambdify来成功地将Product转换成NumPy可以计算的表达式：

import sympy as sy
from sympy.tensor import IndexedBase, Idx
import numpy as np
import sympy.printing.lambdarepr as SPL

def _print_Product(self, expr):
    loops = (
        'for {i} in range({a}, {b}+1)'.format(
            i=self._print(i),
            a=self._print(a),
            b=self._print(b))
        for i, a, b in expr.limits)
    return '(prod([{function} {loops}]))'.format(
        function=self._print(expr.function),
        loops=' '.join(loops))
SPL.NumPyPrinter._print_Product = _print_Product

x = sy.IndexedBase('x')
i = sy.symbols('i', cls=Idx)
n = sy.symbols('n', integer=True, positive=True)
s = sy.Product(x[i], (i, 1, n))
s_lambda = sy.lambdify((sy.DeferredVector('x'), n), s, 'numpy')
print(s_lambda(np.arange(101), 5))

印刷品

120